In [1]:
from threeML import *
import matplotlib.pyplot as plt
%matplotlib inline
%matplotlib notebook
Here we can see a simple example of basic spectral analysis for FITS PHA files from the Fermi Gamma-ray Burst Monitor. FITS PHA files are read in with the OGIPLike plugin that supports reading TYPE I & II PHA files with properly formatted OGIP keywords.
Here, we examine a TYPE II PHA file. Since there are multiple spectra embedded in the file, we must either use the XSPEC style {spectrum_number} syntax to access the appropriate spectum or use the keyword spectrum_number=1.
In [2]:
triggerName = 'bn090217206'
ra = 204.9
dec = -8.4
#Data are in the current directory
datadir = os.path.abspath('.')
#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", observation=obsSpectrum, background=bakSpectrum, response=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" )
As we can see, the plugin probes the data to choose the appropriate likelihood for the given obseration and background data distribution.
In GBM, the background is estimated from a polynomial fit and hence has Gaussian distributed errors via error propagation.
We can also select the energies that we would like to use in a spectral fit. To understand how energy selections work, there is a detailed docstring:
In [5]:
NaI6.set_active_measurements?
Signature: NaI6.set_active_measurements(*args, **kwargs) Docstring: Set the measurements to be used during the analysis. Use as many ranges as you need, and you can specify either energies or channels to be used.
NOTE to Xspec users: while XSpec uses integers and floats to distinguish between energies and channels specifications, 3ML does not, as it would be error-prone when writing scripts. Read the following documentation to know how to achieve the same functionality.
They are specified as 'emin-emax'. Energies are in keV. Example:
set_active_measurements('10-12.5','56.0-100.0')
which will set the energy range 10-12.5 keV and 56-100 keV to be used in the analysis. Note that there is no difference in saying 10 or 10.0.
They are specified as 'c[channel min]-c[channel max]'. Example:
set_active_measurements('c10-c12','c56-c100')
This will set channels 10-12 and 56-100 as active channels to be used in the analysis
You can also specify mixed energy/channel selections, for example to go from 0.2 keV to channel 20 and from channel 50 to 10 keV:
set_active_measurements('0.2-c10','c50-10')
Use 'all' to select all measurements, as in:
set_active_measurements('all')
Use 'reset' to return to native PHA quality from file, as in:
set_active_measurements('reset')
Excluding measurements work as selecting measurements, but with the "exclude" keyword set to the energies and/or channels to be excluded. To exclude between channel 10 and 20 keV and 50 keV to channel 120 do:
set_active_measurements(exclude=["c10-20", "50-c120"])
Call this method more than once if you need to select and exclude. For example, to select between 0.2 keV and channel 10, but exclude channel 30-50 and energy , do:
set_active_measurements("0.2-c10",exclude=["c30-c50"])
To simply add or exclude channels from the native PHA, one can use the use_quailty option:
set_active_measurements("0.2-c10",exclude=["c30-c50"], use_quality=True)
This translates to including the channels from 0.2 keV - channel 10, exluding channels 30-50 and any channels flagged BAD in the PHA file will also be excluded.
:param args: :param exclude: (list) exclude the provided channel/energy ranges :param use_quality: (bool) use the native quality on the PHA file (default=False) :return: File: ~/coding/3ML/threeML/plugins/SpectrumLike.py Type: instancemethod
In [6]:
NaI6
Out[6]:
In [7]:
NaI6.display()
To examine our energy selections, we can view the count spectrum:
In [8]:
NaI6.view_count_spectrum()
Deselected regions are marked shaded in grey.
We have also view which channels fall below a given significance level:
In [9]:
NaI6.view_count_spectrum(significance_level=5)
In [3]:
#This declares which data we want to use. In our case, all that we have already created.
data_list = DataList( NaI6 )
In [4]:
powerlaw = Powerlaw()
In [5]:
GRB = PointSource( triggerName, ra, dec, spectral_shape=powerlaw )
In [6]:
model = Model( GRB )
In [8]:
jl = JointLikelihood( model, data_list, verbose=False )
jl.set_minimizer('ROOT')
res = jl.fit()
In [9]:
jl.minimizer
Out[9]:
In [30]:
res = jl.get_errors()
In [31]:
res = jl.get_contours(powerlaw.index,-1.3,-1.1,20)
In [32]:
res = jl.get_contours(powerlaw.index,-1.25,-1.1,60,powerlaw.K,1.8,3.4,60)
And to plot the fit in the data space we call the data spectrum plotter:
In [33]:
jl.restore_best_fit()
_=display_spectrum_model_counts(jl)
Or we can examine the fit in model space. Note that we must use the analysis_results of the joint likelihood for the model plotting:
In [35]:
plot_point_source_spectra(jl.results,flux_unit='erg/(s cm2 keV)')
We can go Bayesian too!
In [36]:
powerlaw.index.prior = Uniform_prior(lower_bound=-5.0, upper_bound=5.0)
powerlaw.K.prior = Log_uniform_prior(lower_bound=1.0, upper_bound=10)
bayes = BayesianAnalysis(model, data_list)
In [37]:
samples = bayes.sample(n_walkers=50,burn_in=100, n_samples=1000)
In [38]:
fig = bayes.corner_plot()
In [39]:
fig = bayes.corner_plot_cc()
In [40]:
plot_point_source_spectra(bayes.results, flux_unit='erg/(cm2 s keV)',equal_tailed=False)
In [ ]: