Fitting PHA Data with Multi-Component Source Models

Getting Started

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

Reading Data & Instrument Responses Into Sherpa

In this thread, we wish to fit spectral data from the FITS PHA datafile source_pi.fits. This dataset is input into Sherpa with the READ command:

sherpa> DATA source_pi.fits
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

Note that the READ DATA command may be shortened to just DATA.

Now the dataset may be plotted:

sherpa> LPLOT DATA

Note that these data are plotted in bin (i.e. channel) space, because there is no information about RMF/ARF attached to the header of this file. A standard PHA file contains number of counts in each PHA channel, which Sherpa reads. Therefore only counts per channel will be plotted with the above command. For details on Sherpa plotting see Data Visualization thread.

In order to assign the instrument properties (such as for example RMF and ARF files) an instrument model has to be defined. The Sherpa model name for an instrument response model, defined by the standard instrument response files, is RSP. The RSP model component is established for use and is named acis:

sherpa> RSP[acis]
acis.rmf parameter value ["none"] rmf.fits
acis.arf parameter value ["none"] arf.fits
The inferred file type is ARF.  If this is not what you want, please 
specify the type explicitly in the data command.

Notice that in the square bracket we named the response model as acis. Sherpa prompts the user for input of the model's default parameter values, e.g. two filenames for rmf and arf. Because the parameters belong to the acis model they are named acis.rmf and acis.arf. We set the rmf and arf parameter values to the proper RMF and ARF filenames for this dataset. It is possible, however, to establish a response model using only an RMF or an ARF file.

In order to assign acis model to the data and convolve the PHA data with the response, the instrument model has to be defined. Note that this allows for multiple models to be established during the Sherpa session, while only the one assigned to instrument used during fitting.

sherpa> INSTRUMENT = acis

The current definition of Sherpa's instrument model may be examined using SHOW INSTRUMENT:

sherpa> SHOW INSTRUMENT
Instrument 1: rsp1d[acis]
    Param   Type  Value
    -----   ----  -----
 1    rmf string: "rmf.fits" (N_E=198,N_PHA=1024)
 2    arf string: "arf.fits" (N_E=495)

This output shows that rmf.fits and arf.fits currently define the instrument model. Check RSP for requirement on the binning in RMF and ARF files.

The input dataset may be plotted again:

sherpa> LPLOT DATA

The data are now plotted in energy space since the instrument model is providing the information necessary for the computation of the number of predicted counts in each energy bin.

Filtering Spectral Data

Sherpa has many filtering options with IGNORE and NOTICE commands. During this Sherpa session, we would like to fit only the data between 0.5 and 8.0 keV. To apply this filter to the dataset:

sherpa> IGNORE ENERGY :0.5, 8:

Note: this command may be entered only after an instrument model has been defined.

To remove this filter:

sherpa> NOTICE ALL

The filter may then be re-applied (or a different filter may be applied):

sherpa> IGNORE ENERGY :0.5, 8:

The filtered dataset may then be plotted:

sherpa> LPLOT DATA

Figure 1  shows the resulting plot. Notice that the plot now includes only the data in the specified energy region.

Establishing Model Components

We will fit the spectral data using a source model expression that involves multiple model components.

The first of these model components is the one-dimensional power-law model, POWLAW1D. The POWLAW1D model component is established and is named p1 defined in the square bracket:

sherpa> PARAMPROMPT OFF 
Model parameter prompting is off 
sherpa> POWLAW1D[p1]

Since a dataset has been input, Sherpa estimates the initial parameter values (and the minimum and maximum for their ranges) for this model based on the data. The command PARAMPROMPT OFF cancels prompting for immediate changes to these model parameter value estimates.

The second of the model components is the XSpec photoelectric absorption model, called XSWABS in Sherpa. Please see the XSpec User's Guide, or ahelp xs, or ahelp xswabs for more information about this model. Note that all XSpec models are available from within Sherpa and their XSpec names must be preceded by XS.

The XSWABS model component is established and is named abs1:

sherpa> XSWABS[abs1]

Again, since a dataset has been previously input, Sherpa estimates the initial parameter values based on the data.

Defining a Multi-Component Source Model Expression

We will now fit the spectral data using an expression that is the product of the model components established above. The SOURCE defines the final model expression used for fitting the data:

sherpa> SOURCE = abs1*p1

The current definition of Sherpa's source model expression including initial parameter values for each model component may be examined using SHOW SOURCE:

sherpa> SHOW SOURCE
Source 1: (abs1 * p1)
xswabs[abs1]  (XSPEC model name: absw)  (integrate: off)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1     nH thawed      1e-01      1e-07         10            10^22/cm^2
powlaw[p1]  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1  gamma thawed          1        -10         10                      
 2    ref frozen          1     0.5183     7.9789                      
 3   ampl thawed 7.1163e-06 7.1163e-08 7.1163e-04                      

This output shows that abs * p1 is currently defined as the source model expression. By default, integration over an energy bin is turned off for the abs1 model component and is turned on for the p1 model component. However, since we are using a multiplicative source model expression to fit binned PHA data, integration of source expression over each energy bin will be performed during fitting.

Modifying Method & Statistic Settings

The SHOW command may be used to view the current method and statistics settings:

sherpa> SHOW METHOD
Optimization Method: Levenberg-Marquardt

      Name       Value         Min         Max                     Description
      ----       -----         ---         ---                     -----------
 1   iters        2000           1       10000    Maximum number of iterations
 2     eps       1e-03       1e-09           1               Absolute accuracy
 3   smplx           0           0           1  Refine fit with simplex (0=no)
 4 smplxep           1       1e-04        1000    Switch-to-simplex eps factor
 5 smplxit           3           1          20  Switch-to-simplex iters factor

sherpa> SHOW STATISTIC
Statistic:           Chi-Squared Gehrels

We will change the statistic to Cash for fitting this data:

sherpa> STATISTIC CASH

Note that truncation is turned on when the Cash statistic is used. This setting prevents negative model-predicted data values from affecting the convergence process.

Further details about the Levenberg-Marquardt optimization method which is the default method in Sherpa are available by typing:

sherpa> ahelp lev-mar

Fitting

The dataset is then fit:

sherpa> FIT
WARNING: the Levenberg-Marquardt optimization method works
         less robustly when the Cash or cstat statistic is used.
         Consider using the Powell or Simplex method instead.
 LVMQT: V2.0
 LVMQT: initial statistic value = 15306.3
 LVMQT: final statistic value = -6708.78 at iteration 31
          abs1.nH  1e-07  10^22/cm^2  
            p1.gamma  -0.373229     
            p1.ampl  0.000113554     

WARNING:
   The value of abs1.nH is equal to the abs1.nH.min limit boundary.
   You may wish to consider changing min/max values and refitting.

Notice that the first warning recommends that we use a different optimization method. Before refitting the data, we reset the initial parameter values and switch to the Powell method:

sherpa> RESET
sherpa> METHOD POWELL

And now run the fit again:

sherpa> fit
 powll: v1.2
 powll:  initial statistic value =     1.53063E+04
 powll:     converged to minimum =    -7.62300E+03 at iteration =     12
 powll:    final statistic value =    -7.62300E+03
          abs1.nH  1.48605  10^22/cm^2  
            p1.gamma  0.815361     
            p1.ampl  0.000711629     

WARNING:
   The value of p1.ampl within 0.01% of the p1.ampl.max limit boundary.
   You may wish to consider changing min/max values and refitting.

The new warning refers to the final value of the amplitude parameter being to close to the boundary. We will expand allow min/max range for this parameter and refit. First, we examine the current parameter values for the model p1:

sherpa> SHOW p1
powlaw[p1]  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1  gamma thawed     0.8154        -10         10                      
 2    ref frozen          1     0.5183     7.9789                      
 3   ampl thawed 7.1163e-04 7.1163e-08 7.1163e-04                      

As noted in the warning, we change the Max setting for the amplitude (p1.ampl):

sherpa> p1.ampl.max = 0.1

And fit the data a final time:

sherpa> FIT
 powll: v1.2
 powll:  initial statistic value =    -7.62300E+03
 powll:     converged to minimum =    -7.70549E+03 at iteration =      7
 powll:    final statistic value =    -7.70549E+03
          abs1.nH  2.37765  10^22/cm^2  
            p1.gamma  1.50086     
            p1.ampl  0.00200489     

Use LPLOT to plot the fit and residuals. Also see Customizing Sherpa plots thread to learn how to use Slang to make nice plots.

sherpa> LPLOT 2 FIT RESIDUALS
  ==> Error bars computed using Chi Gehrels.
  ==> Error bars computed using Chi Gehrels.
sherpa> PRINT POSTFILE sherpa.spectra.2.ps

Figure 2  shows the resulting plot. The features seen in the residuals could be due to the real source emission being different than the assumed power law model. This is in fact true for our data set. The spectrum is from a supernova remnant usually characterized by thermal emission accompanying with several emission lines. Also some features could be due to calibration uncertainties. We do not discuss further analysis here, which in normal data analysis process would consist of changing the source expression to include preferable plasma emission model and refiting.

Examining Fit Results

Several algorithms are available in Sherpa for examining fit results such as covariance, projection, uncertainty, interval-projection, region-projection A detailed description of how to use different options is available in several threads: Step-by-Step Guide to Estimating Errors and Confidence Levels and Estimating Errors and Confidence Levels. Also Sherpa S-lang module functions allows for easy access to the results in both scripts and on the command line. They are described in Accessing fit results using S-Lang thread.

Checking Sherpa Session Status

The final overall status of this Sherpa session may be viewed as follows:

sherpa> SHOW

Optimization Method: Powell
Statistic:           Cash 

-----------------
Input data files:
-----------------

Data 1: source_pi.fits pha.
Total Size: 1024 bins (or pixels)
Dimensions: 1
Total counts (or values): 3328
Exposure: 7854.47 sec
Count rate: 0.424 cts/sec
  Backscal: 2.366406e-06

Current filters for dataset 1:
notice source 1 all
IGNORE source 1 ENERGY : 0.5 , 8 : 
Noticed filter size: 512 bins 

Sum of data within filter: 3283

------------------------------
Defined analysis model stacks:
------------------------------

source 1 = (abs1 * p1)
instrument source 1 = acis

-------------------------------------------
Defined source/background model components:
-------------------------------------------

powlaw[p1]  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1  gamma thawed     1.5009        -10         10                      
 2    ref frozen          1     0.5183     7.9789                      
 3   ampl thawed 2.0049e-03 7.1163e-08      1e-01                      

xswabs[abs1]  (XSPEC model name: absw)  (integrate: off)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1     nH thawed     2.3776      1e-07         10            10^22/cm^2

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

rsp1d[acis]
    Param   Type  Value
    -----   ----  -----
 1    rmf string: "rmf.fits" (N_E=198,N_PHA=1024)
 2    arf string: "arf.fits" (N_E=495)

Saving a Sherpa Session

To save a Sherpa session so that we may return again later, use the SAVE command:

sherpa> SAVE ALL session1.shp
sherpa> EXIT

where session1.shp is the filename for an ASCII file to which a record the current session will be written, in the form of a Sherpa script.

Using Sherpa Scripts & Restoring a Session

Sherpa commands may either be issued on the Sherpa command line or from a file by the USE command.

To restore the session that was saved to the file session1.shp:

sherpa> USE session1.shp
The inferred file type is ARF.  If this is not what you want, please 
specify the type explicitly in the data command.
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
WARNING: model truncation turned on.  Override with TRUNCATE OFF.
sherpa>

One may verify that the session has been properly restored by comparing the SHOW output from the Checking Sherpa Session Status section.