Basic PHA FITS Analysis with 3ML

Giacomo Vianello (Stanford University) giacomov@stanford.edu

IPython Notebook setup.

This is needed only if you are using the IPython Notebook on your own computer, it is NOT needed if you are on threeml.stanford.edu. This line will activate the support for inline display of matplotlib images:

In [2]:

from threeML import *

import matplotlib.pyplot as plt

%matplotlib inline
%matplotlib notebook

Configuration read from /Users/jburgess/.threeML/threeML_config.yml

Simple standard analysis¶

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 [4]:

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

Auto-probed noise models:
- observation: poisson
- background: gaussian
Range 10.0-30.0 translates to channels 4-20
Range 40.0-950.0 translates to channels 26-125
Now using 117 channels out of 128

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


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

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.

Energy selections:

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.

Channel selections:

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

Mixed channel and energy selections:

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 measurements (i.e., reset to initial state):

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

Exclude measurements:

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

Select and exclude:

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

Using native PHA quality:

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

Investigating the contents of the data¶

We can examine some quicklook properties of the plugin by executing it or calling its display function:

In [6]:

NaI6

Out[6]:

bak file                 /Users/jburgess/coding/3ML/examples/bn09021720...
pha file                 /Users/jburgess/coding/3ML/examples/bn09021720...
n. channels                                                            128
total rate                                                         1549.26
total bkg. rate                                                    993.421
total bkg. rate error                                              6.15657
exposure                                                           19.9127
bkg. exposure                                                           20
significance                                                       54.4521
is poisson                                                            True
bkg. is poisson                                                      False
response                 /Users/jburgess/coding/3ML/examples/bn09021720...

In [7]:

NaI6.display()

	0
bak file	/Users/jburgess/coding/3ML/examples/bn09021720...
pha file	/Users/jburgess/coding/3ML/examples/bn09021720...
n. channels	128
total rate	1549.26
total bkg. rate	993.421
total bkg. rate error	6.15657
exposure	19.9127
bkg. exposure	20
significance	54.4521
is poisson	True
bkg. is poisson	False
response	/Users/jburgess/coding/3ML/examples/bn09021720...

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)

channels below the significance threshold shown in red

Setup for spectral fitting¶

Now we will prepare the plugin for fitting by:

Creating a DataList
Selecting a spectral shape
Creating a likelihood model
Building a joint analysis object

In [10]:

#This declares which data we want to use. In our case, all that we have already created.

data_list = DataList( NaI6 )

In [11]:

powerlaw = Powerlaw()

In [26]:

GRB = PointSource( triggerName, ra, dec, spectral_shape=powerlaw )

In [27]:

model = Model( GRB )

In [28]:

jl = JointLikelihood( model, data_list, verbose=False )
res = jl.fit()

Best fit values:

	Value	Unit
bn090217206.spectrum.main.Powerlaw.K	2.53 +/- 0.21	1 / (cm2 keV s)
bn090217206.spectrum.main.Powerlaw.index	-1.183 +/- 0.015

Correlation matrix:

1.00	-0.98
-0.98	1.00

Values of -log(likelihood) at the minimum:

	-log(likelihood)
NaI6	869.098854
total	869.098854

Examining the fitted model¶

We can now look at the asymmetric errors, countours, etc.

In [30]:

res = jl.get_errors()

Name	Value	Unit
bn090217206.spectrum.main.Powerlaw.K	2.53 -0.20 +0.21	1 / (cm2 keV s)
bn090217206.spectrum.main.Powerlaw.index	-1.183 +/- 0.015

In [31]:

res = jl.get_contours(powerlaw.index,-1.3,-1.1,20)

Widget Javascript not detected.  It may not be installed properly. Did you enable the widgetsnbextension? If not, then run "jupyter nbextension enable --py --sys-prefix widgetsnbextension"

In [32]:

res = jl.get_contours(powerlaw.index,-1.25,-1.1,60,powerlaw.K,1.8,3.4,60)

Widget Javascript not detected.  It may not be installed properly. Did you enable the widgetsnbextension? If not, then run "jupyter nbextension enable --py --sys-prefix widgetsnbextension"

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

Widget Javascript not detected.  It may not be installed properly. Did you enable the widgetsnbextension? If not, then run "jupyter nbextension enable --py --sys-prefix widgetsnbextension"

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)

Widget Javascript not detected.  It may not be installed properly. Did you enable the widgetsnbextension? If not, then run "jupyter nbextension enable --py --sys-prefix widgetsnbextension"
Widget Javascript not detected.  It may not be installed properly. Did you enable the widgetsnbextension? If not, then run "jupyter nbextension enable --py --sys-prefix widgetsnbextension"

Mean acceptance fraction: 0.71602

Maximum a posteriori probability (MAP) point:

	Value	Unit
bn090217206.spectrum.main.Powerlaw.K	2.52 -0.17 +0.22	1 / (cm2 keV s)
bn090217206.spectrum.main.Powerlaw.index	-1.183 -0.016 +0.013

Values of -log(posterior) at the minimum:

	-log(posterior)
NaI6	-869.501547
total	-869.501547

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)

Widget Javascript not detected.  It may not be installed properly. Did you enable the widgetsnbextension? If not, then run "jupyter nbextension enable --py --sys-prefix widgetsnbextension"

In [ ]: