In [1]:
%matplotlib inline
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
#np.seterr(all='ignore')
from threeML import *
In [2]:
triggerName = 'bn090217206'
ra = 204.9
dec = -8.4
#Data are in the current directory
datadir = os.path.abspath('.')
#The .pha, .bak and .rsp files have been prepared with the Fermi
#official software. In the future it will be possible to create
#them directly from the plugin
#Create an instance of the GBM plugin for each detector
#Data files
obsSpectrum = os.path.join( datadir, "bn090217206_n6_srcspectra.pha{1}" )
bakSpectrum = os.path.join( datadir, "bn090217206_n6_bkgspectra.bak{1}" )
rspFile = os.path.join( datadir, "bn090217206_n6_weightedrsp.rsp{1}" )
#Plugin instance
NaI6 = OGIPLike( "NaI6", obsSpectrum, bakSpectrum, rspFile )
#Choose energies to use (in this case, I exclude the energy
#range from 30 to 40 keV to avoid the k-edge, as well as anything above
#950 keV, where the calibration is uncertain)
NaI6.set_active_measurements( "10.0-30.0", "40.0-950.0" )
#Now repeat for the other GBM detectors
obsSpectrum = os.path.join( datadir, "bn090217206_n9_srcspectra.pha{1}" )
bakSpectrum = os.path.join( datadir, "bn090217206_n9_bkgspectra.bak{1}" )
rspFile = os.path.join( datadir, "bn090217206_n9_weightedrsp.rsp{1}" )
#Plugin instance
NaI9 = OGIPLike( "NaI9", obsSpectrum, bakSpectrum, rspFile )
#Choose chanels to use
NaI9.set_active_measurements( "10.0-30.0", "40.0-950.0" )
obsSpectrum = os.path.join( datadir, "bn090217206_b1_srcspectra.pha{1}" )
bakSpectrum = os.path.join( datadir, "bn090217206_b1_bkgspectra.bak{1}" )
rspFile = os.path.join( datadir, "bn090217206_b1_weightedrsp.rsp{1}" )
#Plugin instance
BGO1 = OGIPLike( "BGO1", obsSpectrum, bakSpectrum, rspFile )
#Choose chanels to use (in this case, from 200 keV to 10 MeV)
BGO1.set_active_measurements( "200-10000" )
data_list = DataList( NaI6, NaI9, BGO1 )
In [3]:
#Let's use a Band model, a phenomenological model typically used for GRB
bb= Blackbody()
pl = Powerlaw()
comp_model = bb+pl
band = Band()
In [4]:
GRB = PointSource( triggerName, ra, dec, spectral_shape=comp_model )
model = Model( GRB )
jl = JointLikelihood( model, data_list )
res = jl.fit()
GRB2 = PointSource( triggerName+'_band', ra, dec, spectral_shape=band )
model2 = Model( GRB2 )
jl2 = JointLikelihood( model2, data_list )
res2 = jl2.fit()
In [5]:
#spec_plot_mle = SpectralPlotter(jl)
plot_point_source_spectra(jl.results,
flux_unit='erg2/(cm2 s keV)')
In [6]:
plot_point_source_spectra(jl.results,
use_components=True,
equal_tailed=False,
best_fit='median',
flux_unit='erg2/(cm2 s keV)')
In [7]:
plot_point_source_spectra(jl.results,
energy_unit='MeV',
ene_min=1e-1,
ene_max=1e4,
fit_cmap='jet',
contour_cmap='rainbow_r',
confidence_level=.99,
flux_unit='erg/(cm2 s MeV)',
legend_kwargs={'loc':'upper right'},
plot_style_kwargs={'linestyle':'--'},
contour_style_kwargs={'lw':0}
)
In [8]:
plot_point_source_spectra(jl.results,
jl2.results,
confidence_level=.90,
equal_tailed=False,
flux_unit='erg2/(cm2 s keV)')
In [9]:
plot_point_source_spectra(jl.results,
jl2.results,
confidence_level=0.68,
sum_sources=True,
flux_unit='erg2/(cm2 s keV)')
In [5]:
comp_model.K_1.prior = Log_uniform_prior(lower_bound = 1E-7, upper_bound = 1E-5)
comp_model.K_2.prior = Log_uniform_prior(lower_bound =1E-1,upper_bound = 1E2)
comp_model.index_2.set_uninformative_prior(Uniform_prior)
comp_model.kT_1.prior = Log_uniform_prior(lower_bound =1E1,upper_bound = 1E4)
bayes = BayesianAnalysis(model, data_list)
res= bayes.sample(20,100,500)
band.K.prior = Log_uniform_prior(lower_bound = 1E-3, upper_bound = 1)
band.xp.prior =Log_uniform_prior(lower_bound = 10, upper_bound = 700)
band.alpha.set_uninformative_prior(Uniform_prior)
band.beta.set_uninformative_prior(Uniform_prior)
bayes2 = BayesianAnalysis(model2, data_list)
res2= bayes2.sample(20,100,500)
In [12]:
plot_point_source_spectra(bayes.results,
fit_cmap='plasma',
contour_cmap='winter',
confidence_level=0.9,
flux_unit='erg/(cm2 s MeV)',
legend_kwargs={'loc':'lower left'},
plot_style_kwargs={'linestyle':'--'},
contour_style_kwargs={'lw':0}
)
In [19]:
plot_point_source_spectra(bayes.results,bayes2.results,
fit_cmap='plasma',
equal_tailed=False,
contour_cmap='winter',
confidence_level=.95,
flux_unit='erg2/(cm2 s MeV)',
legend_kwargs={'loc':'lower left'},
plot_style_kwargs={'linestyle':'--'},
contour_style_kwargs={'lw':0}
)
In [14]:
plot_point_source_spectra(bayes.results,bayes2.results,
fit_cmap='plasma',
contour_cmap='winter',
flux_unit='erg2/(cm2 s MeV)',
legend_kwargs={'loc':'lower left'},
plot_style_kwargs={'linestyle':'--'},
contour_style_kwargs={'lw':0},
sum_sources=True
)
In [15]:
plot_point_source_spectra(bayes.results,
fit_cmap='plasma_r',
contour_cmap='jet_r',
use_components=True,
flux_unit='1/(cm2 s MeV)',
legend_kwargs={'loc':'lower left'},
plot_style_kwargs={'linestyle':'--'},
contour_style_kwargs={'lw':0},
)
In [18]:
plot_point_source_spectra(bayes.results,jl.results,
confidence_level=.95,
equal_tailed=False,
flux_unit='erg2/(cm2 s MeV)',
plot_style_kwargs={'linestyle':'--'},
contour_style_kwargs={'lw':1,'linestyle':':','alpha':0.4},
fraction_of_samples=.01
)
In [6]:
_=bayes.corner_plot_cc(color_params='index2')
In [ ]: