Fitting Multiple Orders of HRC-S/LETG Data

Overview

Last Update: 1 Dec 2006 - reviewed for CIAO 3.4: no changes

Synopsis:

Because of the low energy resolution in the HRC-S, the PHA2 file contains two rows (negative and postive) containing all the spectral orders. While it is not possible to separate the overlapping orders, they can be modeled in Sherpa by defining the instrument response as a composite of the orders in which you are interested.

This thread uses response files (gRMFs and gARFs) built in CIAO to model and fit the first three positive and negative orders of the spectra.

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

Getting Started
Reading the Spectrum Files
Building the Instrument Responses
Filtering the Data
Defining the Source Model
Examining Method & Statistic Settings
Subtract the Background
Fitting
Examining Fit Results
Saving and Quitting the Session
History
Images

Getting Started

Sample ObsID used: 460 (LETG/HRC-S, 3C 273)

The files used in this example were created by following several of the CIAO Grating threads:

Here is a list of all the necessary files:

spectra:
460_leg_m1_bin10.pha
460_leg_p1_bin10.pha

gRMFs:
460_leg_-1.grmf
460_leg_-2.grmf
460_leg_-3.grmf
460_leg_1.grmf
460_leg_2.grmf
460_leg_3.grmf

gARFs:
460_LEG_-1.garf
460_LEG_-2.garf
460_LEG_-3.garf
460_LEG_1.garf
460_LEG_2.garf
460_LEG_3.garf

Reading the Spectrum Files

The spectra that will be used in this session have already been binned by a factor of 10. The data are input to Sherpa with the data command (a shorthand version of "read data"):

sherpa> data 1 460_leg_m1_bin10.pha 
The inferred file type is PHA.  If this is not what you want, please 
specify the type explicitly in the data command.
Warning: could not find SYS_ERR column
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: backgrounds UP and DOWN are being read from this file,
         and are being combined into a single background dataset.
Warning: could not find SYS_ERR column

sherpa> data 2 460_leg_p1_bin10.pha 
The inferred file type is PHA.  If this is not what you want, please 
specify the type explicitly in the data command.
Warning: could not find SYS_ERR column
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: backgrounds UP and DOWN are being read from this file,
         and are being combined into a single background dataset.
Warning: could not find SYS_ERR column

Sherpa now refers to the negative-order spectrum as dataset 1 and the positive-order spectrum as dataset 2.

Building the Instrument Responses

The individual instrument response files (gRMFs and gARFs) need to be read into Sherpa as file-based model components: frmf1d for the gRMFs and farf1d for the gARFs. We choose to name the models to match the order of the response, e.g. arfm1 for the -1 gARF.

sherpa> frmf1d[rmfm1](460_leg_-1.grmf)
sherpa> frmf1d[rmfm2](460_leg_-2.grmf)
sherpa> frmf1d[rmfm3](460_leg_-3.grmf)
sherpa> frmf1d[rmfp1](460_leg_1.grmf)
sherpa> frmf1d[rmfp2](460_leg_2.grmf)
sherpa> frmf1d[rmfp3](460_leg_3.grmf)

sherpa> farf1d[arfm1](460_LEG_-1.garf)
sherpa> farf1d[arfm2](460_LEG_-2.garf)
sherpa> farf1d[arfm3](460_LEG_-3.garf)
sherpa> farf1d[arfp1](460_LEG_1.garf)
sherpa> farf1d[arfp2](460_LEG_2.garf)
sherpa> farf1d[arfp3](460_LEG_3.garf)

This message will be printed after each gARF is entered:

The inferred file type is ARF.  If this is not what you want, please 
specify the type explicitly in the data command.

In order to convolve the input datasets with the response model components that have been established, they must be defined as the instrument models. This involves pairing up the gARFs and gRMFs for each order and summing them together, keeping the negative (dataset 1) and positive (dataset 2) separated:

sherpa> instrument 1 = arfm1*rmfm1 + arfm2*rmfm2 + arfm3*rmfm3
sherpa> instrument 2 = arfp1*rmfp1 + arfp2*rmfp2 + arfp3*rmfp3

For each dataset, the photon spectrum will be folded through the arf*rmf combination to produce a counts spectrum, e.g. arfm1*rmfm1 creates spectrum c1. The overall spectrum is then a sum of the components, i.e. c1+c2+c3.

The datasets may now be plotted:

sherpa> lplot 2 data 1 data 2

Figure 1 shows the resulting plot.

Filtering the Data

There are two plate gaps in the HRC-S detector: one at ~50 Å and the other at ~60 Å; see Figure 7.1 in the POG for the HRC-S detector layout. The effect of dithering into these gaps appear in negative-order and positive-order data, respectively, as a flat area of zero counts:

sherpa> d all limits x 45 70

The regions where the data are in a plate gap are evident in Figure 2 . These regions are ignored in the fitting so that the statistic calculations are not skewed:

sherpa> ignore 1 wave 49:56
sherpa> ignore 2 wave 59:66

You may wish to adjust the limits to exclude more or less data around this region. Any other desired filters may be applied to the data at this point as well.

Defining the Source Model

We plan on background-subtracting the data (rather than fitting it simultaneously), so it is only necessary to create a source model expression. We model this source with a broken power law (xsbknpower) absorbed by the interstellar medium (xswabs).

In this example, we choose to use the XSpec version of the models. These models expect that the x-values will always be energy bins. When the analysis setting is using non-energy bins (e.g., WAVE in this case) and an XSpec model is defined, Sherpa converts the bins to energy before sending them to the XSpec model. After the XSpec model finishes, Sherpa converts back to the original units. Sherpa also scales the model values appropriately (e.g., if counts/keV came out of the XSpec model, and Sherpa is working with wavelength bins, then Sherpa scales the output of the XSpec model to counts/Angstrom).

First, we set up each model component. The absorption model will be referred to as "abs", and the broken power law will be "bpow".

sherpa> xswabs[abs]
abs.nH parameter value [0.1] 

sherpa> xsbknpower[bpow]
bpow.PhoInd1 parameter value [1] 
bpow.BreakE parameter value [5] 
bpow.PhoInd2 parameter value [2] 
bpow.norm parameter value [0.0434012]

Note that since a dataset has already been input, Sherpa estimates the initial parameter values for the models based on the data. These values can also be listed with the show command:

sherpa> show models
-------------------------------------------
Defined source/background model components:
-------------------------------------------

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

xsbknpower[bpow]  (XSPEC model name: bknpower)  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1PhoInd1 thawed          1         -3         10                      
 2 BreakE thawed          5          0      1e+06                   keV
 3PhoInd2 thawed          2         -3         10                      
 4   norm thawed 4.3401e-02 4.3401e-04     4.3401 ph/cm^2/s/keV @ 1 keV

We choose to modify a few of the initial parameter values:

sherpa> abs.nH=1.81E-02
sherpa> freeze abs.nh

sherpa> bpow.2=1

bpow.2=1 sets the second parameter of xsbknpower[bpow], which is the break energy [keV], to a lower starting point. The hydrogen column density (nH) is set to the Galactic value and then frozen, which means it will not be allowed to vary during the fit. The rest of the parameters remain thawed.

Now that the model components have been established, the product of the two is assigned as the source model for both datasets:

sherpa> source 1,2 = (abs * bpow)
sherpa> show source 1  
Source 1: (abs * bpow)
xswabs[abs]  (XSPEC model name: wabs)  (integrate: off)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1     nH frozen   1.81e-02      1e-07         10            10^22/cm^2
xsbknpower[bpow]  (XSPEC model name: bknpower)  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1PhoInd1 thawed          1         -3         10                      
 2 BreakE thawed          5          0      1e+06                   keV
 3PhoInd2 thawed          2         -3         10                      
 4   norm thawed 4.3401e-02 4.3401e-04     4.3401 ph/cm^2/s/keV @ 1 keV

Examining Method & Statistic Settings

Next we check 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

For this fit, the default fitting and statistic settings will be used. More information is available from ahelp lev-mar and ahelp chigehrels. For a list of all the available methods and statistic settings, see ahelp method and ahelp statistic, respectively.

Subtract the Background

The final thing to do before fitting is perform background subtraction on the data:

sherpa> subtract 1,2

Fitting

The datasets are now fit:

sherpa> fit 1,2
 LVMQT: V2.0
 LVMQT: initial statistic value = 13186
 LVMQT: final statistic value = 1855.63 at iteration 11
          bpow.PhoInd1  2.22993     
          bpow.BreakE  0.749234  keV  
          bpow.PhoInd2  1.616     
          bpow.norm  0.0196918  ph/cm^2/s/keV @ 1 keV

To plot the fits:

sherpa> lplot 2 fit 1 fit 2

sherpa> d 1 label 125 0.04 "LEG, -1 order"
sherpa> l 1 size 1.5
sherpa> d 2 label 125 0.04 "LEG, +1 order"
sherpa> l 1 size 1.5
sherpa> redraw

The ChIPS commands are used to add labels to the drawing areas. The plot is shown in Figure 3 . Notice that Sherpa does not attempt to fit the regions that we chose to omit.

It is also useful to plot the fit with the residuals:

sherpa> lplot 2 fit 1 delchi

This plot is shown in Figure 4 . By omitting the regions of data over a plate gap, the residuals are contained within three sigma. This will improve the statistical calculations shown in the Examining Fit Results section.

After creating a plot, it may be saved as a PostScript file:

sherpa> print postfile 460_fit_delchi.ps

Examining Fit Results

There are several methods available in Sherpa for examining fit results. The goodness command reports information on the chi-square goodness-of-fit:

sherpa> goodness
Goodness: computed with Chi-Squared Gehrels

DataSet 1: 1582 data points -- 1578 degrees of freedom.
 Statistic value       = 919.714
 Probability [Q-value] = 1
 Reduced statistic     = 0.582835

DataSet 2: 1582 data points -- 1578 degrees of freedom.
 Statistic value       = 935.915
 Probability [Q-value] = 1
 Reduced statistic     = 0.593102

Total   : 3164 data points -- 3160 degrees of freedom.
 Statistic =         1855.63
 Probability =       1
 Reduced statistic = 0.587224

The uncertainty, covariance, and projection commands can be used to estimate confidence intervals for the thawed parameters:

sherpa> covariance

Computed for sherpa.cov.sigma = 1
        --------------------------------------------------------
        Parameter Name      Best-Fit Lower Bound     Upper Bound
        --------------------------------------------------------
          bpow.PhoInd1       2.22993  -0.0129966      +0.0129966    
          bpow.BreakE       0.749234  -0.0164652      +0.0164652    
          bpow.PhoInd2         1.616  -0.0151005      +0.0151005    
          bpow.norm        0.0196918  -0.000245127    +0.000245127

Saving and Quitting the Session

Before exiting Sherpa, you may wish to save the session in order to return to the analysis at a later point:

sherpa> save all 460_fitting_session.shp

All the information about the current session is written to 460_fitting_session.shp, an ASCII file. It may be loaded into Sherpa again with the use command.

Finally, quit the session:

sherpa> quit

History

03 May 2005	original version, new for CIAO 3.2
21 Dec 2005	reviewed for CIAO 3.3: changed filenames to make the CIAO threads that generated them
01 Dec 2006	reviewed for CIAO 3.4: no changes