In [1]:
%pylab inline
from threeML import *
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','b0'],destination_directory=data_dir_gbm,compress_tte=True)
src_selection = '0-71'
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,100-200",
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",
verbose=False)
nai3.set_active_measurements("8.0-30.0", "40.0-950.0")
bgo0.set_active_measurements("250-43000")
In [3]:
triggerName = 'bn080916009'
ra = 121.8
dec = -61.3
data_list = DataList(nai3,bgo0 )
band = Band()
GRB1 = PointSource( triggerName, ra, dec, spectral_shape=band )
model1 = Model( GRB1 )
pl_bb= Powerlaw() + Blackbody()
GRB2 = PointSource( triggerName, ra, dec, spectral_shape=pl_bb )
model2 = Model( GRB2 )
In [4]:
jl1 = JointLikelihood( model1, data_list, verbose=False )
res = jl1.fit()
jl2 = JointLikelihood( model2, data_list, verbose=False )
res = jl2.fit()
The JointLikelihood objects are passed to the SpectralFlux class. Then either model_flux or component_flux are called depending on the flux desired.
The astropy system of units is used to specfiy flux units and an error is raised if the user selects an improper unit. The integration range is specified and the unit for this range can be altered.
In [5]:
res = calculate_point_source_flux(10,40000,jl1.results,jl2.results,flux_unit='erg/(s cm2)',energy_unit='keV')
A panadas DataFrame is returned with the sources (a fitting object can have multiple sources) flux, and flux error.
We can also change to photon fluxes by specifying the proper flux unit (here we changed to m^2). Here, the integration unit is also changed.
In [6]:
res = calculate_point_source_flux(10,40000,jl1.results,jl2.results,flux_unit='1/(s cm2)',energy_unit='Hz',equal_tailed=False)
In [5]:
res = calculate_point_source_flux(10,40000,
jl1.results,jl2.results,
flux_unit='erg/(s cm2)',
energy_unit='keV',use_components=True)
Then we can look at our component fluxes. The class automatically solves the error propagation equations to properly propagate the parameter errors into the components
In [7]:
res = calculate_point_source_flux(10,40000,jl1.results,jl2.results,flux_unit='erg/(s cm2)',
energy_unit='keV',
equal_tailed=False,
use_components=True, components_to_use=['Blackbody','total'])
A dictionary of sources is return that contains pandas DataFrames listing the fluxes and errors of each componenet.
NOTE: With proper error propagation, the total error is not always the sqrt of the sum of component errors squared!
In [8]:
pl_bb.K_1.prior = Log_uniform_prior(lower_bound = 1E-1, upper_bound = 1E2)
pl_bb.index_1.set_uninformative_prior(Uniform_prior)
pl_bb.K_2.prior = Log_uniform_prior(lower_bound = 1E-6, upper_bound = 1E-3)
pl_bb.kT_2.prior = Log_uniform_prior(lower_bound = 1E0, upper_bound = 1E4)
In [9]:
bayes = BayesianAnalysis(model2,data_list)
_=bayes.sample(30,100,500)
In [10]:
res = calculate_point_source_flux(10,40000,
bayes.results,
flux_unit='erg/(s cm2)',
energy_unit='keV')
Once again, a DataFrame is returned. This time, it contains the mean flux from the distribution, the specfied level (default is 0.05) credible regions and the flux distribution itself.
One can plot the distribtuion:
In [11]:
from astropy.visualization import quantity_support
quantity_support()
_=hist(res[1]['flux distribution'][0],bins=20)
In [12]:
res = calculate_point_source_flux(10,40000,
bayes.results,
flux_unit='erg/(s cm2)',
energy_unit='keV',
use_components=True)
We can easily now visulaize the flux distribtuions from the individual components.
In [13]:
_=hist(log10(res[1]['flux distribution'][0].value),bins=20)
_=hist(log10(res[1]['flux distribution'][1].value),bins=20)
In [14]:
res = calculate_point_source_flux(10,40000,
bayes.results,jl1.results,jl2.results,
flux_unit='erg/(s cm2)',
energy_unit='keV')