Introduction to Fitting PHA Spectra

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_rmf("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()
Data Set: 1
Filter: 0.1248-12.4100 Energy (keV)
Noticed Channels: 1.0-1024.0
name           = 3c273.pi
channel        = Float64[1024]
counts         = Float64[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   = UInt32[2002]
n_chan   = UInt32[2002]
matrix   = Float64[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 = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 38564.1414549

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



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

 
sherpa> data_cnt_rate = calc_data_sum()/get_exposure(id=1)  # calculating count rate in cts/sec 
sherpa> print(data_cnt_rate)
0.0190848557908


sherpa> bkg_sum = calc_data_sum(bkg_id=1)   # total counts (or values) in the background data
sherpa> print(bkg_sum)
216


sherpa> bkg_cnt_rate = calc_data_sum(bkg_id=1)/get_exposure(bkg_id=1)  # calculating background count rate in cts/sec
 
sherpa> print(bkg_cnt_rate)
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 with 'print(get_data().units)' or show_filter, and modified with set_analysis.

Figure 1: Plot of source spectrum

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

sherpa> plot_data()

Figure 2 shows the resulting plot.

Figure 2: Source spectrum, filtered and background-subtracted

The plot axis scaling may be changed to log by calling log_scale (and changed back to linear with linear_scale):

sherpa> log_scale() #Set both axes to log scale
sherpa> log_scale(XY_AXIS) 

sherpa> log_scale(X_AXIS) #Set only x-axis to log scale

sherpa> log_scale(Y_AXIS) #Set only y-axis to log scale

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> show_model()
Model: 1
apply_rmf(apply_arf((38564.6089269 * (xsphabs.abs1 * powlaw1d.p1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      frozen         0.07            0       100000 10^22 atoms / cm^2
   p1.gamma     thawed            1          -10           10
   p1.ref       frozen            1 -3.40282e+38  3.40282e+38
   p1.ampl      thawed            1            0  3.40282e+38

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 models in the source expression. It may be necessary to modify the parameter values, since the Sherpa guess functionality does not apply to absorption models. However, we can use this command to guess the initial parameter values and ranges for the power-law model component (as model parameter values are not automatically guessed in Sherpa 4.1 as they were in previous versions of the software):

sherpa> guess(p1)

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((38564.6089269 * (xsphabs.abs1 * powlaw1d.p1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      frozen         0.07            0       100000 10^22 atoms / cm^2
   p1.gamma     thawed            1          -10           10           
   p1.ref       frozen            1 -3.40282e+38  3.40282e+38           
   p1.ampl      thawed  0.000148802  1.48802e-06    0.0148802           

The guess command makes an initial guess at parameter values to ensure convergence, but it is always a good idea to check that the initial range of values (soft limits) is sensible for the data being fit.

Fitting

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

sherpa> fit()
Solar Abundance Vector set to angr:  Anders E. & Grevesse N. Geochimica et Cosmochimica Acta 53, 197 (1989)
 Cross Section Table set to bcmc:  Balucinska-Church and McCammon, 1998
Dataset               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 6.83386e+10
Final fit statistic   = 37.9079 at function evaluation 22
Data points           = 44
Degrees of freedom    = 42
Probability [Q-value] = 0.651155
Reduced statistic     = 0.902569
Change in statistic   = 6.83386e+10
   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 "sigma", the sigma residuals of the fit [(data - model)/error].

Figure 3: Fit and sigma residuals

Examining Fit Results

sherpa> show_fit()

Optimization Method: LevMar
name    = levmar
ftol    = 1.19209289551e-07
xtol    = 1.19209289551e-07
gtol    = 1.19209289551e-07
maxfev  = None
epsfcn  = 1.19209289551e-07
factor  = 100.0
verbose = 0

Statistic: Chi2Gehrels

Fit:Dataset               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 6.83386e+10
Final fit statistic   = 37.9079 at function evaluation 22
Data points           = 44
Degrees of freedom    = 42
Probability [Q-value] = 0.651155
Reduced statistic     = 0.902569
Change in statistic   = 6.83386e+10
   p1.gamma       2.15852
   p1.ampl        0.00022484

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

sherpa> calc_chisqr()

array([  8.18672858e+00,   5.97520920e+00,   1.12656172e+00,
         2.70996397e-01,   1.52222572e-01,   6.41849665e-03,
         7.45752882e-01,   3.65878489e-01,   7.64892123e-01,
         7.54582959e-02,   6.15857733e-01,   2.33666307e+00,
         2.01512590e-01,   3.57967845e-02,   9.88219709e-01,
         3.98382524e-01,   1.35736379e-02,   1.26937320e+00,
         1.71508420e+00,   2.12311366e-01,   8.44121608e-02,
         1.43547663e-01,   5.85917784e-01,   4.80676459e-02,
         1.97969717e+00,   1.00172041e-01,   1.05894104e+00,
         4.39172444e-01,   3.28729745e-01,   6.84238569e-02,
         1.50504771e-01,   1.24953231e+00,   6.00567768e-01,
         5.35535361e-04,   8.38929628e-01,   1.14986971e+00,
         1.77134437e-01,   6.99634183e-02,   2.85814374e-01,
         2.05568964e+00,   2.35952178e-01,   2.13674054e-01,
         1.51009824e-01,   4.34761024e-01])

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

sherpa> calc_photon_flux()
0.00046995137298

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

sherpa> calc_energy_flux()
9.61260829653e-13

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

sherpa> calc_data_sum()
706.85714092

sherpa> calc_data_sum(2., 10.)
306.230157017


sherpa> calc_model_sum()
638.45693375313488

sherpa> calc_model_sum(2., 10.)
270.73187526053857


sherpa> calc_source_sum()
0.0469964088469405

sherpa> calc_source_sum(2., 10.)
0.0072654873922453405

Scripting It

The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing execfile("fit.py") on the Sherpa command line.

The Sherpa script command may be used to save everything typed on the command line in a Sherpa session:

sherpa> script(filename="sherpa.log", clobber=False)

The CXC is committed to helping Sherpa users transition to new syntax as smoothly as possible. If you have existing Sherpa scripts or save files, submit them to us via the CXC Helpdesk and we will provide the CIAO/Sherpa 4.1 syntax to you.