In [3]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
%matplotlib notebook
from threeML import *
In [2]:
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','b0'],destination_directory=data_dir_gbm,compress_tte=True)
src_selection = "0.-5."
nai3 = FermiGBMTTELike('NAI3',
os.path.join(data_dir, "glg_tte_n3_bn080916009_v01.fit.gz"),
"-10-0,100-200",
src_selection,
rsp_file=os.path.join(data_dir, "glg_cspec_n3_bn080916009_v00.rsp2"))
bgo0 = FermiGBMTTELike('BGO0',
os.path.join(data_dir, "glg_tte_b0_bn080916009_v01.fit.gz"),
"-10-0,100-200",
src_selection,
rsp_file=os.path.join(data_dir, "glg_cspec_b0_bn080916009_v00.rsp2"))
nai3.set_active_measurements("8.0-30.0", "40.0-950.0")
bgo0.set_active_measurements("250-43000")
In [4]:
nai3.create_time_bins(start=0,stop=10,dt=2,method='constant')
bgo0.read_bins(nai3)
nai3_intervals = nai3.get_ogip_from_binner()
bgo0_intervals = bgo0.get_ogip_from_binner()
all_data = nai3_intervals + bgo0_intervals
In [5]:
triggerName = 'bn080916009'
ra = 121.8
dec = -61.3
data_list = DataList( *all_data )
band = Band()
GRB = PointSource( triggerName, ra, dec, spectral_shape=band )
model = Model( GRB )
time = IndependentVariable("time",13.0, unit='s')
model
law = Powerlaw(K=500,index=-1)
#law = Line(a=-0.02,b=-2.0)
model.add_independent_variable(time)
model.link(model.bn080916009.spectrum.main.Band.xp,time,law)
mean_time = [np.mean([x,y]) for x,y in zip(nai3.bins[0],nai3.bins[1])]
ff = step_generator(mean_time,band.alpha)
model.link(model.bn080916009.spectrum.main.Band.alpha,time,ff)
ff = step_generator(mean_time,band.K)
model.link(model.bn080916009.spectrum.main.Band.K,time,ff)
ff = step_generator(mean_time,band.beta)
model.link(model.bn080916009.spectrum.main.Band.beta,time,ff)
for n,b,t in zip(nai3_intervals,bgo0_intervals,mean_time):
n.external_property(time,t)
b.external_property(time,t)
In [6]:
model
Out[6]:
In [7]:
jl = JointLikelihood( model, data_list, verbose=False )
#jl.set_minimizer("PYOPT","COBYLA")
res = jl.fit()
In [ ]:
cleanup_downloaded_GBM_data(gbm_data)