Introduction to Fitting PHA Spectra

Getting Started

Please follow the "Sherpa Threads: Getting Started" thread.

Load the Spectrum & Instrument Responses

First, load the spectrum file:

sherpa> DATA 3c273.pi
The inferred file type is PHA.  If this is not what you want, please 
specify the type explicitly in the data command.
WARNING: statistical errors specified in the PHA file.
         These are currently IGNORED.  To use them, type:
         READ ERRORS "<filename>[cols CHANNEL,STAT_ERR]" fitsbin
Background data are being input from:
  /data/sherpa/pha_intro/3c273_bg.pi
WARNING: statistical errors specified in the PHA file.
         These are currently IGNORED.  To use them, type:
         READ BERRORS "<filename>[cols CHANNEL,STAT_ERR]" fitsbin
RMF is being input from:
  /data/sherpa/pha_intro/3c273.rmf
ARF is being input from:
  /data/sherpa/pha_intro/3c273.arf
sherpa>

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.

In this case, Sherpa automatically defines the instrument model (RSP) for source and background data and adds it to the instrument model stack. If Sherpa does not automatically read in the response files, define the instrument model manually:

sherpa> INSTRUMENT SOURCE 1 = RSP[mydetector]("3c273.rmf", "3c273.arf")
The inferred file type is ARF.  If this is not what you want, please 
specify the type explicitly in the data command.

sherpa> INSTRUMENT BACK 1 = mydetector

Use SHOW to get the status of the Sherpa session:

sherpa> SHOW
...
------------------------------
Defined analysis model stacks:
------------------------------

instrument source 1 = AutoReadResponse
instrument back 1 = AutoReadResponse

------------------------------------
Defined instrument model components:
------------------------------------

rsp1d[AutoReadResponse]
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1    rmf string: "/data/sherpa/pha_intro/3c273.rmf" (N_E=1090,N_PHA=1024)
 2    arf string: "/data/sherpa/pha_intro/3c273.arf" (N_E=1090)

Plot the data:

sherpa> LPLOT 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.

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:

sherpa> NOTICE ENERGY 0.1:6.0

At this point, we also opt to subtract the background data (see the Fitting Spectral Data: With/Without Independent Background Responses thread for other options):

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.

After turning off parameter prompting, we define an expression that is the product of these two components:

sherpa> PARAMPROMPT OFF
Model parameter prompting is off
sherpa> SOURCE = xsphabs[abs]*powlaw1d[p1]
sherpa> abs.nh=0.07
sherpa> FREEZE abs

Based on the data that has been loaded, Sherpa estimates the initial value for each parameter. We manually set the value of the hydrogen column density (nH) to the known Galactic value for our source and freeze the parameter; this means that it will not be included in the fit.

The current source model definition may be displayed:

sherpa> SHOW SOURCE
Source 1: (abs * p1)
xsphabs[abs]  (XSPEC model name: phabs)  (integrate: off)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1     nH frozen      7e-02      1e-07         10    10**22 atoms/cm**2
powlaw[p1]  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1  gamma thawed          1        -10         10                      
 2    ref frozen          1     0.1248     5.9057                      
 3   ampl thawed 7.9257e-06 7.9257e-08 7.9257e-04                      

By default, integration of model values over each energy bin (delta E) is turned off for xsphabs and is turned on for powlaw1d. Since we are using a multiplicative source model expression, integration will be performed during fitting; see "ahelp integrate" for more information.

Note that Sherpa and XSPEC absorption models have to be multiplied by models which have the normalization and amplitude parameters, such as powlaw1d. These models cannot be used as single model in the source expression, while the models with the normalization can.

Fitting

Now we are ready to run the fit. We use the default statistics (CHI GEHRELS) and optimization method (Levenberg-Marquardt):

sherpa> FIT
 LVMQT: V2.0
 LVMQT: initial statistic value = 355.854
 LVMQT: final statistic value = 37.9079 at iteration 10
            p1.gamma  2.1585     
            p1.ampl  0.000224838     

The screen output shows the final value of the minimum statistics with the total number of iterations and the final values of the thawed parameters of the power law, namely photon index and amplitude at ~1 keV.

Plotting the fit and residuals:

sherpa> LPLOT 2 FIT DELCHI

produces this plot . The errors are plotted as DELCHI, the sigma residuals of the fit [(data - model)/error].

Examining Fit Results

There are several ways to calculate statistics on fit results in Sherpa. Here we show two of them: GOODNESS and COVARIANCE.

sherpa> GOODNESS
Goodness: computed with Chi-Squared Gehrels

DataSet 1: 44 data points -- 42 degrees of freedom.
 Statistic value       = 37.9079
 Probability [Q-value] = 0.651155
 Reduced statistic     = 0.902569

sherpa> COVAR

Computed for sherpa.cov.sigma = 1
        --------------------------------------------------------
        Parameter Name      Best-Fit Lower Bound     Upper Bound
        --------------------------------------------------------
            p1.gamma          2.1585  -0.0827851      +0.0827851    
            p1.ampl      0.000224838  -1.48256e-05    +1.48256e-05  

Script the Session

The commands used in this thread may also be executed via a Sherpa script, which we have named pha_intro.shp:

DATA 3c273.pi
NOTICE ENERGY 0.1:6.0
SUBTRACT
PARAMPROMPT OFF
SOURCE = xsphabs[abs]*powlaw1d[p1]
abs.nh=0.07
FREEZE abs
FIT
GOODNESS
COVAR

There are two ways to load this script:

• With the USE command in an existing session:

  sherpa> USE pha_intro.shp 
  

• When starting a new Sherpa session:

  unix% sherpa pha_intro.shp
  

Note that when the script ends, you are at the Sherpa prompt. This allows you to continue the analysis - plot, adjust parameters, fit the data again - if desired.