Using Spectral Fits to Calculate Fluxes

3ML provides a module to calculate the integral flux from a spectral fit and additionally uses the covariance matrix or posteror to calculate the error in the flux value for the integration range selected


In [1]:
%pylab inline

from threeML import *


Populating the interactive namespace from numpy and matplotlib
Configuration read from /Users/jburgess/.threeML/threeML_config.yml

Data setup

Using GRB 080916C as an example, we will fit two models to the time-integrated spectrum to demostrate the flux calculations capabilites.


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")


WARNING UserWarning: No TLMIN keyword found. This DRM does not follow OGIP standards. Assuming TLMIN=1


WARNING RuntimeWarning: Maximum MC energy is smaller than maximum EBOUNDS energy


WARNING RuntimeWarning: Maximum MC energy is smaller than maximum EBOUNDS energy


WARNING RuntimeWarning: Maximum MC energy is smaller than maximum EBOUNDS energy


WARNING RuntimeWarning: Maximum MC energy is smaller than maximum EBOUNDS energy


WARNING RuntimeWarning: Maximum MC energy is smaller than maximum EBOUNDS energy


WARNING RuntimeWarning: Minimum MC energy is larger than minimum EBOUNDS energy


WARNING RuntimeWarning: Maximum MC energy is smaller than maximum EBOUNDS energy


WARNING RuntimeWarning: Minimum MC energy is larger than minimum EBOUNDS energy


WARNING RuntimeWarning: Maximum MC energy is smaller than maximum EBOUNDS energy


WARNING RuntimeWarning: Minimum MC energy is larger than minimum EBOUNDS energy


WARNING RuntimeWarning: Maximum MC energy is smaller than maximum EBOUNDS energy


WARNING RuntimeWarning: Minimum MC energy is larger than minimum EBOUNDS energy

Model setup

We will fit two models: a Band function and a CPL+Blackbody


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 )

Fitting

MLE

We fit both models using MLE


In [4]:
jl1 = JointLikelihood( model1, data_list, verbose=False )

res = jl1.fit()



jl2 = JointLikelihood( model2, data_list, verbose=False )

res = jl2.fit()


Best fit values:

Value Unit
bn080916009.spectrum.main.Band.K (1.49 +/- 0.10) x 10^-2 1 / (cm2 keV s)
bn080916009.spectrum.main.Band.alpha -1.07 +/- 0.05
bn080916009.spectrum.main.Band.xp (5.0 +/- 1.0) x 10^2 keV
bn080916009.spectrum.main.Band.beta -2.05 +/- 0.12
Correlation matrix:

1.000.93-0.970.56
0.931.00-0.900.49
-0.97-0.901.00-0.65
0.560.49-0.651.00
Values of -log(likelihood) at the minimum:

-log(likelihood)
BGO0 1197.745203
NAI3 1091.776951
total 2289.522155
Best fit values:

Value Unit
bn080916009.spectrum.main.composite.K_1 8.8 +/- 0.7 1 / (cm2 keV s)
bn080916009...index_1 -1.576 +/- 0.019
bn080916009.spectrum.main.composite.K_2 (4.5 +/- 0.7) x 10^-6 1 / (cm2 keV3 s)
bn080916009.spectrum.main.composite.kT_2 (5.04 +/- 0.24) x 10 keV
Correlation matrix:

1.00-0.91-0.330.41
-0.911.000.17-0.35
-0.330.171.00-0.94
0.41-0.35-0.941.00
Values of -log(likelihood) at the minimum:

-log(likelihood)
BGO0 1210.246464
NAI3 1095.307328
total 2305.553792

Flux caluclation

Total flux

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')


flux negative error positive error
bn080916009: total 1.91114162965e-06 erg / (cm2 s) 1.57453116206e-06 erg / (cm2 s) 2.26762602494e-06 erg / (cm2 s)
bn080916009_2: total 1.87661507183e-06 erg / (cm2 s) 1.63222306088e-06 erg / (cm2 s) 2.15504525488e-06 erg / (cm2 s)

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)


flux negative error positive error
bn080916009: total 176.087354288 1 / (cm2 s) 0.442687488649 1 / (cm2 s) 451.886189451 1 / (cm2 s)
bn080916009_2: total 39434350432.4 1 / (cm2 s) 8970680196.46 1 / (cm2 s) 58995295825.8 1 / (cm2 s)

Components

If we want to look at component fluxes, we examine our second fit.

We can first look at the total flux:


In [5]:
res = calculate_point_source_flux(10,40000,
                                  jl1.results,jl2.results,
                                  flux_unit='erg/(s cm2)',
                                  energy_unit='keV',use_components=True)


flux negative error positive error
bn080916009: total 1.9039848819e-06 erg / (cm2 s) 1.53646152236e-06 erg / (cm2 s) 2.29025979966e-06 erg / (cm2 s)
bn080916009_2: Blackbody 2.98162435614e-07 erg / (cm2 s) 2.32816767523e-07 erg / (cm2 s) 3.74874099063e-07 erg / (cm2 s)
bn080916009_2: Powerlaw 1.57446724856e-06 erg / (cm2 s) 1.3498807165e-06 erg / (cm2 s) 1.85325378616e-06 erg / (cm2 s)

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'])


flux negative error positive error
bn080916009: total 1.9098744467e-06 erg / (cm2 s) 1.50734491539e-06 erg / (cm2 s) 2.27706381595e-06 erg / (cm2 s)
bn080916009_2: Blackbody 2.9672438355e-07 erg / (cm2 s) 2.26109650025e-07 erg / (cm2 s) 3.69757209064e-07 erg / (cm2 s)
bn080916009_2: total 1.87978597056e-06 erg / (cm2 s) 1.61169204266e-06 erg / (cm2 s) 2.12809837574e-06 erg / (cm2 s)

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!

Bayesian fitting

Now we will look at the results when a Bayesian fit is performed.

We set our priors and then sample:


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)


Mean acceptance fraction: 0.575

Maximum a posteriori probability (MAP) point:

Value Unit
bn080916009.spectrum.main.composite.K_1 8.7 -0.4 +0.8 1 / (cm2 keV s)
bn080916009...index_1 -1.574 -0.022 +0.013
bn080916009.spectrum.main.composite.K_2 (4.5 -0.6 +0.7) x 10^-6 1 / (cm2 keV3 s)
bn080916009.spectrum.main.composite.kT_2 (5.03 -0.20 +0.25) x 10 keV
Values of -log(posterior) at the minimum:

-log(posterior)
BGO0 -1207.673512
NAI3 -1092.482153
total -2300.155665

Flux Calculation

Total Flux

Just as with MLE, we pass the BayesianAnalysis object to the SpectralFlux class.

Now the propagation of fluxes is done using the posterior of the analysis.


In [10]:
res = calculate_point_source_flux(10,40000,
                                  bayes.results,
                                  flux_unit='erg/(s cm2)',
                                  energy_unit='keV')


flux negative error positive error flux distribution
bn080916009: total 1.86758978808e-06 erg / (cm2 s) 1.64122284622e-06 erg / (cm2 s) 2.12871511049e-06 erg / (cm2 s) [1.80524710826e-06 erg / (cm2 s), 2.0143686888...

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)


Components

We can also look at components as before. A dictionary of sources is returned, each containing Dataframes of the components information and distributions.


In [12]:
res = calculate_point_source_flux(10,40000,
                                  bayes.results,
                                  flux_unit='erg/(s cm2)',
                                  energy_unit='keV',
                                  use_components=True)


flux negative error positive error flux distribution
bn080916009: Blackbody 3.02645553975e-07 erg / (cm2 s) 2.41284530882e-07 erg / (cm2 s) 3.84229537662e-07 erg / (cm2 s) [2.73586762211e-07 erg / (cm2 s), 3.2080879109...
bn080916009: Powerlaw 1.55673566178e-06 erg / (cm2 s) 1.35453225651e-06 erg / (cm2 s) 1.79666593764e-06 erg / (cm2 s) [1.48721358571e-06 erg / (cm2 s), 1.7530923838...

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')


flux negative error positive error
bn080916009: total 1.88373471868e-06 erg / (cm2 s) 1.53858737596e-06 erg / (cm2 s) 2.2728799556e-06 erg / (cm2 s)
bn080916009_2: total 1.86461731624e-06 erg / (cm2 s) 1.63196809404e-06 erg / (cm2 s) 2.12230564683e-06 erg / (cm2 s)
flux negative error positive error flux distribution
bn080916009: total 1.85951690412e-06 erg / (cm2 s) 1.64495301858e-06 erg / (cm2 s) 2.11183597466e-06 erg / (cm2 s) [1.86976760607e-06 erg / (cm2 s), 1.4308721697...