Fermi GBM data is in a format that does not lend itself to being used with standard software such as XSPEC. However, the FermiGBMLikeTTE in 3ML plugin allows the user to work directly with the TTE data in its native format
FermiGBMLikeTTE provides the following functionality
It creates a standard 3ML Model and therefore can be used like any other plugin without using specical tools to create PHA files.
Import 3ML as always to make sure you have the plugin
In [1]:
%matplotlib inline
%matplotlib notebook
from threeML import *
get_available_plugins()
We will look at GRB080916C as a test case
FermiGBM_TTE_Like takes as arguments:
FermiGBM_TTE_Like will attempt to find the best background polynomial order via a LRT. The background is fit with an Poisson likehood via method developed by Giacomo V.
In [2]:
# os.path.join is a way to generate system-independent
# paths (good for unix, windows, Mac...)
data_dir = os.path.join('gbm','bn080916009')
trigger_number = 'bn080916009'
# Download the data
data_dir_gbm = os.path.join('gbm',trigger_number)
gbm_data = download_GBM_trigger_data(trigger_number,detectors=['n3','n4','b0'],
destination_directory=data_dir_gbm,compress_tte=True)
In [4]:
plot_tte_lightcurve(os.path.join(data_dir, "glg_tte_n3_bn080916009_v01.fit.gz"),stop=100)
In [ ]:
src_selection = "0.-71."
# We start out with a bad background interval to demonstrate a few features
nai3 = FermiGBMTTELike('NAI3',
os.path.join(data_dir, "glg_tte_n3_bn080916009_v01.fit.gz"),
os.path.join(data_dir, "glg_cspec_n3_bn080916009_v00.rsp2"),
src_selection,
"-10-0,50-200", # setting this badly on purpose!
poly_order=2)
nai4 = FermiGBMTTELike('NAI4',
os.path.join(data_dir_gbm, "glg_tte_n4_bn080916009_v01.fit.gz"),
os.path.join(data_dir_gbm, "glg_cspec_n4_bn080916009_v00.rsp2"),
src_selection,
"-10-0,120-200",
poly_order=-1,
verbose=False)
bgo0 = FermiGBMTTELike('BGO0',
os.path.join(data_dir, "glg_tte_b0_bn080916009_v01.fit.gz"),
os.path.join(data_dir, "glg_cspec_b0_bn080916009_v00.rsp2"),
src_selection,
"-10-0,100-200")
You will get a warning about GBM DRMS because they lack the keyword to indicate the first channel of the DRM. It is ok because the normal choice is 1
The TTE class build upon the generic EventList class which can be used to get some information on or selections.
If you are connected to the internet, timing information for other instruments can be obtained:
In [3]:
nai3.display()
In [4]:
nai3
Out[4]:
Let's look at the lightcurve of NAI3 to check out background fit:
In [5]:
nai3.view_lightcurve(-10,200.,.5)
Oy! That is not so nice! Luckily, we can simply select another interval!
We can set the background polynomial order ourself, or use LRT to decide which order is best. If the order is set to -1, then it is automatically determined. Setting the order will initiate a refit of the background using the last imput background interval selections.
In [6]:
nai3.set_background_interval("-10.--1.","120-200") # You can select as many as required!
In [7]:
nai3.background_poly_order = -1 # Auto determine!
nai3.view_lightcurve(-10,200.,.5)
It is also possible to select multiple (disjoint only) src intervals
In [8]:
nai3.set_active_time_interval('0-10','15-50')
nai3.view_lightcurve(-10,200.,.5)
Selecting non-disjoint intervals will result in the intervals being merged into one:
In [9]:
nai3.set_active_time_interval('0-10','5-50')
nai3.view_lightcurve(-10,200.,.5)
In [10]:
# go back to our original selection
nai3.set_active_time_interval(src_selection)
If we want to examine the light curve in different energy ranges, we can supply an argument:
In [11]:
nai3.view_lightcurve(-10,200.,.5,energy_selection='50-300')
In [12]:
nai3.view_lightcurve(-10,200.,.5,energy_selection='50-300, 500-900')
We need to select the energies we would like to fit over. GBM has over/underflow channels which must be exlcuded from the fit. This is not always at the same energy, so we need to check. FermiGBM_TTE_Like (and FermiGBMLike ) allow you to plot the count spectra so you can see what you will be excluding in the fit.
In [13]:
nai3.set_active_measurements("10.0-30.0", "40.0-900.0")
nai3.view_count_spectrum()
We can even checkout which channels have significant counts above the background:
In [14]:
bgo0.set_active_measurements("250-42500")
bgo0.view_count_spectrum(significance_level=1)
In [15]:
triggerName = 'bn080916009'
ra = 121.8
dec = -61.3
data_list = DataList(nai3,nai4,bgo0 )
band = Band()
GRB = PointSource( triggerName, ra, dec, spectral_shape=band )
model = Model( GRB )
model.display()
In [16]:
jl = JointLikelihood( model, data_list, verbose=False )
res = jl.fit()
In [17]:
sp = SpectralPlotter(jl)
_=sp.plot_model(y_unit='erg2/(cm2 s keV)', num_ene=300)
We can examine our fit with the data:
In [18]:
_ = display_spectrum_model_counts(jl,min_rate=10,step=False)
In [19]:
res = jl.get_contours(band.xp,250,500,20)
In [20]:
res = jl.get_contours(band.xp,250,500,25,band.alpha,-1.11,-.8,25)
In [21]:
nai3.create_time_bins(start=0,stop=10,method='significance',sigma=30)
In [22]:
nai3.view_lightcurve(use_binner=True,stop=100)
We must define functions to get the data and the model for the JointLikelihoodSet. For the TTE data, this means simply seting the bins we want to use.
The easiest way to acheive this is with a binner object and then a iterating through the text_bins list that is generated:
In [44]:
def data_getter(id):
new_data = []
for data in data_list.values():
data.set_active_time_interval(nai3.text_bins[id])
new_data.append(data)
return DataList(*new_data)
def model_getter(id):
return clone_model(model)
Now we simply pass these functions to the JointLikelihoodSet and we can go
In [45]:
jl_set = JointLikelihoodSet(data_getter=data_getter,
model_getter=model_getter,
n_iterations=len(nai3.text_bins))
res, lh = jl_set.go(compute_covariance=True)
In [46]:
res
Out[46]:
In [47]:
lh
Out[47]:
In [48]:
# First define priors
# We can do it explicitly like this:
# (be careful not to choose the boundaries outside of the allowed value
# for the parameter, according to the min_value and max_value properties)
band.K.prior = Log_uniform_prior(lower_bound=1e-4, upper_bound=3)
band.xp.prior = Log_uniform_prior(lower_bound=10, upper_bound=1e5)
# or use the set_uninformative_prior method, which will use as lower_bound
# and upper_bound the current boundaries for the parameter. Such boundaries
# must exists and be finite
band.alpha.set_uninformative_prior(Uniform_prior)
band.beta.set_uninformative_prior(Uniform_prior)
bayes = BayesianAnalysis(model, data_list)
In [49]:
samples = bayes.sample(n_walkers=50,burn_in=200, n_samples=1000)
In [50]:
fig = bayes.corner_plot(plot_contours=True, plot_density=False)
In [18]:
samples = bayes.sample_multinest(n_live_points=400,resume=False)
In [19]:
fig = bayes.corner_plot_cc()
In [51]:
# equal-tailed credible regions
eq_tail = bayes.get_credible_intervals()
In [52]:
# highest denisty intervals
hdi = bayes.get_highest_density_interval()
In [53]:
bayes.get_effective_free_parameters()
Out[53]:
In [54]:
nai3.write_pha('nai3',overwrite=True)
In [55]:
# Create the time bins
nai3.create_time_bins(0,10,method='constant',dt=1)
# Write them to a PHA file
nai3.write_pha_from_binner('nai3_binner',overwrite=True)
In [56]:
saved_pha = OGIPLike('test',pha_file='nai3_binner.pha{1}')
saved_pha.view_count_spectrum()
In [57]:
cleanup_downloaded_GBM_data(gbm_data)
In [ ]: