Fitting Grating Data

Overview

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

Synopsis:

This thread provides a general introduction to fitting grating data in Sherpa. Loading and filtering data are covered, as well as defining instrument and source models.

Users working with HRC-S/LETG grating data will also find the Fitting Multiple Orders of HRC-S/LETG Data thread helpful for their analysis.

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 and Background Models
Examining Method & Statistic Settings
Fitting
Examining Fit Results
Saving and Quitting the Session
History
Images

Getting Started

Sample ObsID used: 459 (HETG/ACIS-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:
459_heg_m1_bin10.pha
459_heg_p1_bin10.pha
459_meg_m1_bin10.pha
459_meg_p1_bin10.pha

gARFs:
459_heg_m1.arf
459_heg_p1.arf
459_meg_m1.arf
459_meg_p1.arf

The spectrum that will be used in this session has been binned by a factor of 10.

Users may also choose to run the Create Grating RMFs for ACIS Observations thread. Creating observation-specific gRMFs is optional, and is discussed further in the Building the Instrument Responses section.

The data files are available in sherpa.tar.gz, as explained in the Sherpa Getting Started thread.

Reading the Spectrum Files

The data are input to Sherpa with the data command (a shorthand version of "read data"):

sherpa> data 1 459_heg_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 459_heg_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> data 3 459_meg_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 4 459_meg_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 spectra as follows:

HEG, -1 order = dataset 1
HEG, +1 order = dataset 2
MEG, -1 order = dataset 3
MEG, +1 order = dataset 4

Building the Instrument Responses

First, the instrument models are established by the rsp command. The arf parameter value is then set to the corresponding file for each order and arm:

sherpa> paramprompt off

sherpa> rsp[hegm1]
sherpa> rsp[hegp1]
sherpa> rsp[megm1]
sherpa> rsp[megp1]

sherpa> hegm1.arf = 459_heg_m1.arf
sherpa> hegp1.arf = 459_heg_p1.arf
sherpa> megm1.arf = 459_meg_m1.arf
sherpa> megp1.arf = 459_meg_p1.arf

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 gARF and spectrum for each order, via the instrument command:

sherpa> instrument 1 = hegm1
sherpa> instrument 2 = hegp1
sherpa> instrument 3 = megm1
sherpa> instrument 4 = megp1

The current definition of the instrument model may be examined using show instrument:

sherpa> show instrument 1
Instrument 1: rsp1d[hegm1]
    Param   Type  Value
    -----   ----  -----
 1    rmf string: "none" (N_E=8192,N_PHA=8192)
 2    arf string: "459_heg_m1.arf" (N_E=8192)

Notice that Sherpa has defined properties for the rmf parameter, even though we did not enter a file. Sherpa has support for "dummy" instruments: if data have been input and the instrument stack contains only an ARF, a dummy RMF will be created that maps the ARF bins to the data bins, if possible. ahelp instrument contains more information on "dummy" instruments.

The datasets may now be plotted:

sherpa> lplot 4 data 1 data 2 data 3 data 4

Figure 1 shows the resulting plot.

Filtering the Data

We choose to filter the data to focus on an area of interest:

sherpa> ignore allsets all
sherpa> notice allsets wave 1:15

The ignore command is used to ignore all the data in every dataset, then notice is used to select the desired regions. You may wish to adjust the limits to exclude more or less of your data.

Each filtered dataset may then be plotted:

sherpa> lplot 4 data 3 data 4 data 3 data 4

Notice that the plot now includes only the data in the specified wavelength regions. Figure 2 shows the resulting plot.

Defining the Source and Background Models

We plan on simultaneously fitting the background data (rather than subtracting it), so we need to create a model expression for the source and the background. We model this source with a broken power law (bpl1d) absorbed by the interstellar medium (atten). The background will be modeled by a one-dimensional power law (powlaw1d), also absorbed by the ISM (the same atten model).

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

sherpa> paramprompt on
Model parameter prompting is on

sherpa> atten[abs]
abs.hcol parameter value [1e+20] 
abs.heiRatio parameter value [0.1] 
abs.heiiRatio parameter value [0.01] 

sherpa> bpl1d[bpow]
bpow.gamma1 parameter value [0] 
bpow.gamma2 parameter value [0] 
bpow.eb parameter value [7.99625] 
bpow.ref parameter value [7] 
bpow.ampl parameter value [0.0238299] 

sherpa> powlaw1d[pow1d]
pow1d.gamma parameter value [1] 
pow1d.ref parameter value [7] 
pow1d.ampl parameter value [0.0238299]

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:
-------------------------------------------

Atten[abs]  (integrate: off)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1   hcol thawed      1e+20      1e+17      1e+24                      
 2heiRatio thawed      1e-01          0          1                      
 3heiiRatio thawed      1e-02          0          1                      

bpl1d[bpow]  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1 gamma1 thawed          0        -10         10                      
 2 gamma2 thawed          0        -10         10                      
 3     eb thawed     7.9963     1.0075     14.985                      
 4    ref frozen          7          1     14.985                      
 5   ampl thawed  2.383e-02  2.383e-04      2.383                      

powlaw[pow1d]  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1  gamma thawed          1        -10         10                      
 2    ref frozen          7          1     14.985                      
 3   ampl thawed  2.383e-02  2.383e-04      2.383

Next we modify the initial parameter value for abs.hcol:

sherpa> abs.hcol=1.81e20 
sherpa> freeze abs

The hydrogen column density (hcol) is set to the Galactic value. All the abs parameters are then frozen, which means they will not be allowed to vary during the fit.

Now that the model components have been established, the product of abs and bpow is assigned as the source model for all datasets:

sherpa> source 1:4 = abs*bpow

while the background model is set as the product of abs and pow1d:

sherpa> background 1:4 = abs*pow1d

Both model definitions can be listed with the show command:

sherpa> show source 1
Source 1: (abs * bpow)
Atten[abs]  (integrate: off)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1   hcol frozen   1.81e+20      1e+17      1e+24                      
 2heiRatio frozen      1e-01          0          1                      
 3heiiRatio frozen      1e-02          0          1                      
bpl1d[bpow]  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1 gamma1 thawed          0        -10         10                      
 2 gamma2 thawed          0        -10         10                      
 3     eb thawed     7.9963     1.0075     14.985                      
 4    ref frozen          7          1     14.985                      
 5   ampl thawed  2.383e-02  2.383e-04      2.383

sherpa> show background 1
Background 1: (abs * pow1d)
Atten[abs]  (integrate: off)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1   hcol frozen   1.81e+20      1e+17      1e+24                      
 2heiRatio frozen      1e-01          0          1                      
 3heiiRatio frozen      1e-02          0          1                      
powlaw[pow1d]  (integrate: on)
    Param   Type      Value        Min        Max                 Units
    -----   ----      -----        ---        ---                 -----
 1  gamma thawed          1        -10         10                      
 2    ref frozen          7          1     14.985                      
 3   ampl thawed  2.383e-02  2.383e-04      2.383

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.

Fitting

The datasets are now fit:

sherpa> fit
 LVMQT: V2.0
 LVMQT: initial statistic value = 1.06211e+10
 LVMQT: final statistic value = 789026 at iteration 114
          bpow.gamma1  0.420507     
          bpow.gamma2  -0.0786515     
          bpow.eb  11.7909     
          bpow.ampl  0.00225955     
         pow1d.gamma  0.276206     
         pow1d.ampl  0.000238299     

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

As the warning says, we reset the minimum boundary of powbkgmh.ampl and refit the data:

sherpa> pow1d.ampl.min=2.383e-10

sherpa> fit                     
 LVMQT: V2.0
 LVMQT: initial statistic value = 789026
 LVMQT: final statistic value = 1954.98 at iteration 6
          bpow.gamma1  0.414261     
          bpow.gamma2  -0.0607346     
          bpow.eb  11.7837     
          bpow.ampl  0.00237397     
         pow1d.gamma  0.21045     
         pow1d.ampl  9.79897e-06

To plot the fits:

sherpa> lplot 4 fit 3 fit 4 fit 3 fit 4

sherpa> d 1,3,4 ylabel ""
sherpa> title "3C 273 (ObsID 459)"

sherpa> d 1 label 12 0.075 "HEG -1"
sherpa> d 2 label 12 0.075 "HEG +1"
sherpa> d 3 label 12 0.125 "MEG -1"
sherpa> d 4 label 12 0.125 "MEG +1"
sherpa> d all l all green

sherpa> redraw

The ChIPS commands are used to add labels to the drawing areas. The plot is shown in Figure 3 .

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

sherpa> lplot 2 fit 1 delchi

This plot is shown in Figure 4 . After creating a plot, it may be saved as a PostScript file:

sherpa> print postfile 459_1_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: 561 data points -- 555 degrees of freedom.
 Statistic value       = 473.223
 Probability [Q-value] = 0.994862
 Reduced statistic     = 0.852654

Background for DataSet 1: 561 data points -- 559 degrees of freedom.
 Statistic value       = 114.649
 Probability [Q-value] = 1
 Reduced statistic     = 0.205097

DataSet 2: 561 data points -- 555 degrees of freedom.
 Statistic value       = 499.407
 Probability [Q-value] = 0.956151
 Reduced statistic     = 0.899832

Background for DataSet 2: 561 data points -- 559 degrees of freedom.
 Statistic value       = 98.6341
 Probability [Q-value] = 1
 Reduced statistic     = 0.176447

DataSet 3: 281 data points -- 275 degrees of freedom.
 Statistic value       = 279.114
 Probability [Q-value] = 0.419589
 Reduced statistic     = 1.01496

Background for DataSet 3: 281 data points -- 279 degrees of freedom.
 Statistic value       = 119.737
 Probability [Q-value] = 1
 Reduced statistic     = 0.429164

DataSet 4: 281 data points -- 275 degrees of freedom.
 Statistic value       = 260.716
 Probability [Q-value] = 0.722871
 Reduced statistic     = 0.948058

Background for DataSet 4: 281 data points -- 279 degrees of freedom.
 Statistic value       = 109.496
 Probability [Q-value] = 1
 Reduced statistic     = 0.392459

Total   : 3368 data points -- 3362 degrees of freedom.
 Statistic =         1954.98
 Probability =       1
 Reduced statistic = 0.581492

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.gamma1       0.414261  -0.00746802     +0.00746802   
          bpow.gamma2     -0.0607528  -0.119213       +0.119213     
          bpow.eb            11.7837  -0.243908       +0.243908     
          bpow.ampl       0.00237397  -8.36986e-06    +8.36986e-06  
         pow1d.gamma         0.21045  -0.0654054      +0.0654054    
         pow1d.ampl      9.79894e-06  -2.59661e-07    +2.59661e-07

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 459_fitting_session.shp

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

Finally, quit the session:

sherpa> quit

History

14 Jan 2005	updated for CIAO 3.2: minor changes to screen output
11 Jul 2005	overall revision to thread, changes to screen output
21 Dec 2005	reviewed for CIAO 3.3: no changes
01 Dec 2006	reviewed for CIAO 3.4: no changes