Introduction to Fitting PHA Spectra

Sherpa Threads (CIAO 4.0 Beta)

[Python Syntax]

Overview

Last Update: 29 Apr 2008 - show_all() command is available in CIAO 4.0

Synopsis:

The basic steps used in fitting spectral data are illustrated in this thread. The data used herein were created by running the Step-by-Step Guide to Creating ACIS Spectra for Pointlike Sources thread (or, equivalently, the psextract script).

Proceed to the HTML or hardcopy (PDF: A4 | letter) version of the thread.

Load the Spectrum & Instrument Responses
Filter the Data & Subtract the Background
Defining the Source Model
Fitting
Examining Fit Results
History
Images

Load the Spectrum & Instrument Responses

First, load the spectrum file:

sherpa> load_pha("3c273.pi")
statistical errors were found in file '3c273.pi'
but not used; to use them, re-read with use_errors=True
read ARF file 3c273.arf
read RMF file 3c273.rmf
statistical errors were found in file '3c273_bg.pi'
but not used; to use them, re-read with use_errors=True
read background file 3c273_bg.pi

Since the RESPFILE, ANCRFILE, and BACKFILE header keywords were updated in the spectrum file, the response files (RMF and ARF) and background file are automatically read in as well. If the default dataset id of "1" is to be used, it does not need to be explicitly included in the load function; only the data filenames are required in this case.

If Sherpa does not automatically read in the background and response files, read them manually:

sherpa> load_arf("3c273.arf")
sherpa> load_arf("3c273.rmf")
sherpa> load_bkg("3c273_bg.pi")

Use show_all() to get the status of the Sherpa session. Some additional commands are used to get the total number of counts and counts rates calculated from the data.

sherpa> show_all()
Optimization Method: LevMar
Statistic:           Chi2Gehrels

Data Set: 1
name           = 3c273.pi
channel        = Int32[1024]
counts         = Int32[1024]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = Int16[1024]
quality        = Int16[1024]
exposure       = 38564.6089269
backscal       = 2.52643646989e-06
areascal       = 1.0
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = [1]

RMF Data Set: 1:1
name     = 3c273.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = Int16[2002]
n_chan   = Int16[2002]
matrix   = Float32[61834]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

ARF Data Set: 1:1
name     = 3c273.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float32[1090]
bin_lo   = None
bin_hi   = None
exposure = 38564.1414549

Background Data Set: 1:1
name           = 3c273_bg.pi
channel        = Int32[1024]
counts         = Int32[1024]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 38564.6089269
backscal       = 1.87253514146e-05
areascal       = 1.0
grouped        = False
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Background RMF Data Set: 1:1
name     = 3c273.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = Int16[2002]
n_chan   = Int16[2002]
matrix   = Float32[61834]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

Background ARF Data Set: 1:1
name     = 3c273.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float32[1090]
bin_lo   = None
bin_hi   = None
exposure = 38564.1414549

# total counts (or values) in the data
sherpa> data_sum=calc_data_sum()
sherpa> print(data_sum)
736.0

# calculating count rate in cts/sec
sherpa> data_exp=get_data().exposure
sherpa> data_crate=data_sum/data_exp
sherpa> print(data_crate)
0.0190848557908

# total counts (or values) in the background data
sherpa> bkg_sum=sum(get_bkg().counts)
sherpa> print(bkg_sum)
216

# calculating background count rate in cts/sec
sherpa> bkg_exp=get_bkg().exposure
sherpa> bkg_crate=bkg_sum/bkg_exp
sherpa> print(bkg_crate)
0.00560099028644

Plot the data:

sherpa> plot_data()

The data are plotted in energy space - as seen in Figure 1 - since the instrument model provides the information necessary to compute the predicted counts for each bin. In general, the units of the x-axis are determined by the value in the units field of the data, which may be accessed by get_data().units.

Filter the Data & Subtract the Background

Looking at the plot, we decide to evaluate the fit on the energy range between 0.1 and 6.0 keV. The CIAO why topic on Choosing an Energy Filter has more information on this process. To apply this filter to the dataset we use the notice_id() function, which requires explicit input of the data set id. The first argument in notice_id() is the data set id, the second defines the lower energy of the range, and the third the higher energy of the range. The data between 0.1 and 6.0 keV will be noticed with this command:

sherpa> notice_id(1, 0.1, 6.0)

At this point, we also opt to subtract the background data:

sherpa> subtract()

Figure 2 shows the resulting plot.

Defining the Source Model

Before fitting the data, it is necessary to define a model that characterizes the source. We use a source model composed of two model components:

powlaw1d - a one-dimensional power law.
xsphabs - an XSpec photoelectric absorption model.

We define an expression that is the product of these two components. The hydrogen column density (nH) is set to the known Galactic value for the source and the parameter is frozen so that it will not be allowed to vary in the fit.

sherpa> set_source(xsphabs.abs1 * powlaw1d.p1)
sherpa> abs1.nH = 0.07
sherpa> freeze(abs1.nH)

The current source model definition may be displayed:

sherpa> print(get_model())
(abs1 * p1)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      frozen         0.07            0       100000      10^22
   p1.gamma     thawed            1          -10           10
   p1.ref       frozen            1      -1e+120       1e+120
   p1.ampl      thawed            1            0       1e+120

Note that Sherpa and XSPEC absorption models have to be multiplied by a model which has normalization and amplitude parameters, such as powlaw1d. They cannot be used as single model in the source expression. It may be necessary to modify the parameter values, since there is no guess function in the Beta version of Sherpa.

Fitting

Now we are ready to run the fit. The default statistics (chi2gehrels) and optimization method (levmar) are used:

sherpa> fit()
Statistic value = 37.9079 at function evaluation 4
Data points = 44
Degrees of freedom = 42
Probability [Q-value] = 0.651155
Reduced statistic  = 0.902569
   p1.gamma       2.15852
   p1.ampl        0.00022484

The fit information includes the statistic value for chi2gehrels, goodness-of-fit and reduced chi-square along with the best-fit parameter values of the photon index and amplitude.

The best fit model with the data and residuals may be plotted in the same window:

sherpa> plot_fit_delchi()

which creates Figure 3 . The errors are plotted as "delchi", the sigma residuals of the fit [(data - model)/error].

Examining Fit Results

Goodness of fit

In CIAO 3.4, the goodness command was used to get the chi-squared goodness-of-fit. This information is now reported with the best-fit values after a fit; it is no longer necessary to run a separate command. The get_fit_results() command allows access to this information after the fit has been performed:

sherpa> print(get_fit_results())
succeeded = True
parnames  = ('p1.gamma', 'p1.ampl')
parvals   = (2.1585155109, 0.000224840147863)
covarerr  = None
statval   = 37.9079137903
numpoints = 44
dof       = 42
qval      = 0.651155274054
rstat     = 0.90256937596
message   = both actual and predicted relative reductions in the sum of squares are at most ftol=1.19209e-07
nfev      = 4

# retrieve a single value
sherpa> print(get_fit_results().qval)
0.651155274054
sherpa> print(get_fit_results().rstat)
0.90256937596

The number of bins in the fit (numpoints), the number of degrees of freedom (dof, i.e. the number of bins minus the number of free parameters), and the statistic value (statval) are reported. If the chosen statistic is one of the chi-square statistics, as in this example, the reduced statistic (rstat, i.e. the statistic value divided by the number of degrees of freedom) and the probability (qval, Q-value) are included as well.

The calc_chisqr() command calculates the statistic contribution per bin:

sherpa> calc_chisqr()

array([  8.18672857e+00,   5.97520920e+00,   1.12656176e+00,
         2.70996433e-01,   1.52222564e-01,   6.41849654e-03,
         7.45752883e-01,   3.65878490e-01,   7.64892124e-01,
         7.54582962e-02,   6.15857733e-01,   2.33666307e+00,
         2.01512592e-01,   3.57967769e-02,   9.88219657e-01,
         3.98382525e-01,   1.35736379e-02,   1.26937320e+00,
         1.71508420e+00,   2.12311365e-01,   8.44121611e-02,
         1.43547662e-01,   5.85917785e-01,   4.80676462e-02,
         1.97969717e+00,   1.00172042e-01,   1.05894103e+00,
         4.39172443e-01,   3.28729744e-01,   6.84238565e-02,
         1.50504771e-01,   1.24953231e+00,   6.00567767e-01,
         5.35535306e-04,   8.38929631e-01,   1.14986971e+00,
         1.77134436e-01,   6.99634175e-02,   2.85814372e-01,
         2.05568965e+00,   2.35952177e-01,   2.13674053e-01,
         1.51009822e-01,   4.34761022e-01])

Confidence intervals

The covariance() command - which may be shortened to covar() - computes covariance matrices and provides an estimate of confidence intervals for the thawed parameters (see the related command projection() also):

sherpa> covar()
covariance 1-sigma bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   p1.gamma          2.15852   -0.0827856    0.0827856
   p1.ampl        0.00022484 -1.48256e-05  1.48256e-05

The output is the best-fit parameter value with positive and negative error estimates.

Flux and Counts

To calculate the flux of a PHA dataset, use the calc_photon_flux() and calc_energy_flux() commands. The flux may be calculated over the entire dataset or over a specific range. If a range is requested, the dataset id ("1" in this case) must be included as the first argument.

sherpa> calc_photon_flux()
0.00046995137298

sherpa> calc_photon_flux(1, 2., 10.) 
7.26548742347e-05

sherpa> calc_energy_flux()
9.61260829653e-13

sherpa> calc_energy_flux(1, 2., 10.)
4.55133191273e-13

To calculate the counts of a PHA dataset, use the calc_data_sum(), calc_model_sum(), or calc_source_sum() commands. As with the flux calculations, these commands may be given a range:

sherpa> calc_data_sum()
706.85714092

sherpa> calc_data_sum(1, 2., 10.)
306.230157017


sherpa> calc_model_sum()
638.456686064

sherpa> calc_model_sum(1, 2., 10.)
270.73162768


sherpa> calc_source_sum()
0.0469951376773

sherpa> calc_source_sum(1, 2., 10.)
0.0469951376773

History

14 Nov 2007	rewritten for CIAO 4.0 Beta 3
29 Apr 2008	`show_all()` command is available in CIAO 4.0