Fitting Grating Data

Sherpa Threads (CIAO 4.1)

[S-Lang Syntax]

Overview

Last Update: 28 Apr 2009 - replaced use of atten model with Sherpa user model "atten_wave"; new script command is available with CIAO 4.1.2

Synopsis:

This thread provides a general introduction to fitting grating data in Sherpa. Loading and filtering data are covered, as well as defining instrument responses 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
Loading 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
Scripting It
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.garf
459_heg_p1.garf
459_meg_m1.garf
459_meg_p1.garf

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

Users may also choose to run the ACIS Grating RMFs thread. Creating observation-specific gRMFs is optional, and is discussed further in the Loading 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 load_pha command:

sherpa> load_pha(1, "459_heg_m1_bin10.pha");
statistical errors were found in file '459_heg_m1_bin10.pha' 
but not used; to use them, re-read with use_errors=True
read background_up into a dataset from file 459_heg_m1_bin10.pha
read background_down into a dataset from file 459_heg_m1_bin10.pha

sherpa> load_pha(2,"459_heg_p1_bin10.pha");
statistical errors were found in file '459_heg_p1_bin10.pha' 
but not used; to use them, re-read with use_errors=True
read background_up into a dataset from file 459_heg_p1_bin10.pha
read background_down into a dataset from file 459_heg_p1_bin10.pha

sherpa> load_pha(3, "459_meg_m1_bin10.pha");
statistical errors were found in file '459_meg_m1_bin10.pha' 
but not used; to use them, re-read with use_errors=True
read background_up into a dataset from file 459_meg_m1_bin10.pha
read background_down into a dataset from file 459_meg_m1_bin10.pha

sherpa> load_pha(4, "459_meg_p1_bin10.pha");
statistical errors were found in file '459_meg_p1_bin10.pha' 
but not used; to use them, re-read with use_errors=True
read background_up into a dataset from file 459_meg_p1_bin10.pha
read background_down into a dataset from file 459_meg_p1_bin10.pha

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

Loading the Instrument Responses

Note that it is no longer necessary to build instrument models as in Sherpa3.4 - the instrument response is automatically established when the appropriate response files (ARF, RMF) are input to the Sherpa session. If the instrument response files are written in the header of the PHA file, Sherpa will load them automatically; if not, they need to be loaded manually with the load_arf and load_rmf commands. Since we are working with grating data in this example, we load only the ARF files corresponding to each of the four orders and the associated backgrounds:

sherpa> load_arf(1, "459_heg_m1.arf");
sherpa> load_arf(2, "459_heg_p1.arf");
sherpa> load_arf(3, "459_meg_m1.arf");
sherpa> load_arf(4, "459_meg_p1.arf");

load_arf(1, "459_heg_m1.arf", bkg_id=1);
load_arf(1, "459_heg_m1.arf", bkg_id=2);

load_arf(2, "459_heg_p1.arf", bkg_id=1);
load_arf(2, "459_heg_p1.arf", bkg_id=2);

load_arf(3, "459_meg_m1.arf", bkg_id=1);
load_arf(3, "459_meg_m1.arf", bkg_id=2);

load_arf(4, "459_meg_p1.arf", bkg_id=1);
load_arf(4, "459_meg_p1.arf", bkg_id=2);

The current definition of the instrument response may be examined using the command show_all:

sherpa> show_all();
Data Set: 1
Filter: 0.5775-12.3676 Energy (keV)
Noticed Channels: 1.0-8192.0
name           = 459_heg_m1_bin10.pha
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = [1, 2]

RMF Data Set: 1:1
None

ARF Data Set: 1:1
name     = 459_heg_m1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 38565.2847986

Background Data Set: 1:1
name           = 459_heg_m1_bin10.pha
channel        = Float64[8192]
counts         = Int16[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Background RMF Data Set: 1:1
None

Background ARF Data Set: 1:1
name     = 459_heg_m1.arf
energ_lo = Float64[8192]
exposure = 38565.2847986

Background Data Set: 1:1
name           = 459_heg_m1_bin10.pha
channel        = Float64[8192]
counts         = Int16[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Background RMF Data Set: 1:1
None

Background ARF Data Set: 1:1
name     = 459_heg_m1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 38565.2847986

Background Data Set: 1:2
name           = 459_heg_m1_bin10.pha
channel        = Float64[8192]
counts         = Int16[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Data Set: 2
Filter: 0.5775-12.3676 Energy (keV)
Noticed Channels: 1.0-8192.0
name           = 459_heg_p1_bin10.pha
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = [1, 2]

RMF Data Set: 2:1
None

ARF Data Set: 2:1
name     = 459_heg_p1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 38563.0133422

Background Data Set: 2:1
name           = 459_heg_p1_bin10.pha
channel        = Float64[8192]
counts         = Int16[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Background RMF Data Set: 2:1
None

Background ARF Data Set: 2:1
name     = 459_heg_p1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 38563.0133422

Background Data Set: 2:2
name           = 459_heg_p1_bin10.pha
channel        = Float64[8192]
counts         = Int16[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Data Set: 3
Filter: 0.2957-12.3370 Energy (keV)
Noticed Channels: 1.0-8192.0
name           = 459_meg_m1_bin10.pha
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = [1, 2]

RMF Data Set: 3:1
None

ARF Data Set: 3:1
name     = 459_meg_m1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 38565.2855739

Background Data Set: 3:1
name           = 459_meg_m1_bin10.pha
channel        = Float64[8192]
counts         = Int16[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Background RMF Data Set: 3:1
None

Background ARF Data Set: 3:1
name     = 459_meg_m1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 38565.2855739

Background Data Set: 3:2
name           = 459_meg_m1_bin10.pha
channel        = Float64[8192]
counts         = Int16[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Data Set: 4
Filter: 0.2957-12.3370 Energy (keV)
Noticed Channels: 1.0-8192.0
name           = 459_meg_p1_bin10.pha
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = [1, 2]

RMF Data Set: 4:1
None

ARF Data Set: 4:1
name     = 459_meg_p1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 38563.0141169

Background Data Set: 4:1
name           = 459_meg_p1_bin10.pha
channel        = Float64[8192]
counts         = Int16[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Background RMF Data Set: 4:1
None

Background ARF Data Set: 4:1
name     = 459_meg_p1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 38563.0141169

Background Data Set: 4:2
name           = 459_meg_p1_bin10.pha
channel        = Float64[8192]
counts         = Int16[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Int16[8192]
quality        = Int16[8192]
exposure       = 38564.6089269
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Before plotting the data with the plot command, ensure that the units field of each data set is set to "wavelength" with the set_analysis command, as in the following example:

sherpa> show_filter(); 

Data Set Filter: 1
0.5775-12.3676 Energy (keV)

Data Set Filter: 2
0.5775-12.3676 Energy (keV)

Data Set Filter: 3
0.2957-12.3370 Energy (keV)

Data Set Filter: 4
0.2957-12.3370 Energy (keV)

sherpa> set_analysis("wave");

sherpa> show_filter();    

Data Set Filter: 1
1.0025-21.4675 Wavelength (Angstrom)

Data Set Filter: 2
1.0025-21.4675 Wavelength (Angstrom)

Data Set Filter: 3
1.0050-41.9350 Wavelength (Angstrom)

Data Set Filter: 4
1.0050-41.9350 Wavelength (Angstrom)

The data may now be plotted:

sherpa>  plot("data", 1 , "data", 2, "data", 3, "data", 4);

Figure 1 shows the resulting plot.

[Plot of ACIS HEG and MEG +/- 1 orders for 3C 273]

Figure 1: Plotting the four orders

Plot of ACIS HEG and MEG +/- 1 orders for 3C 273

Users conducting wavelength analysis are cautioned that model evaluation for PHA data sets in Sherpa 4.1 produces model values based on an energy grid, regardless of whether the 'set_analysis("wavelength")' command has been issued. This is done to ensure that XSpec models (those with prefix "xs") receive the intended energy grid for evaluation. As a result, model parameters and plots associated with non-XSpec models will reflect values computed on the energy grid and not the wavelength grid. Refer to the bug page for set_analysis for more information, including a recommended work-around.

Filtering the Data

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

sherpa> ignore()
sherpa> notice(1., 15.)

The ignore command is used to ignore all the data in every data set, then notice is used to select the desired data range in all data sets. You may wish to adjust the limits to exclude more or less of your data. To apply ignore() or notice() data constraints to a specific data set, use the ignore_id and notice_id commands, respectively.

Each filtered data set may then be plotted:

sherpa> plot("data", 1, "data", 2, "data", 3, "data", 4);

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

[Plot of ACIS HEG and MEG +/- 1 orders of 3C 273, restricted to the wavelength range 1.0 - 15.0 angstroms]

Figure 2: Filtering the data sets

Plot of ACIS HEG and MEG +/- 1 orders of 3C 273, restricted to the wavelength range 1.0 - 15.0 angstroms.

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

Note: Currently, the Sherpa atten model erroneously returns model parameter values in energy units, as the internal grid of this model is defined in Angstroms; this will be fixed in a future release. This is the case even if the command set_analysis("wave") has been issued. In the meantime, the workaround involves defining a Sherpa user model which allows atten to be used in wavelength analysis. The user model "atten_wave" is defined as follows, and is assigned an ID of "abs1"; for more on Sherpa user models, see the thread "Sherpa User Models."

set_analysis("wave");

hc = 12.39841874;  #in [keV-Angstrom]
variable dummy = atten.dummy;
dummy.integrate=0;

define atten_wave(p, elo, ehi)
{
    variable hc = 12.39841874;
    variable lo = hc/ehi;
    variable hi = hc/elo;
    return py_call(dummy.calc,p, lo, hi);
}

load_user_model(&atten_wave, "abs1");
add_user_pars("abs1",
              ["hcol","heiRatio","heiiRatio"],
              [dummy.hcol.val,dummy.heiRatio.val,dummy.heiiRatio.val],
              [dummy.hcol.min,dummy.heiRatio.min,dummy.heiiRatio.min],
              [dummy.hcol.max,dummy.heiRatio.max,dummy.heiiRatio.max],
              [dummy.hcol.units,dummy.heiRatio.units,dummy.heiiRatio.units],
              [0,0,0]);

Next, we set up the rest of the model components with create_model_component, and set some initial parameter values. The absorption model is referred to as "abs1", the broken power law will be "bpow1", and the 1-D power law will be "pow1d":


sherpa> abs1.hcol = 1e+20;
sherpa> abs1.heiRatio = 0.1; 
sherpa> abs1.heiiRatio = 0.01;

sherpa> create_model_component("bpl1d", "bpow1");
sherpa> bpow1.gamma1 = 0;
sherpa> bpow1.gamma2 = 0; 
sherpa> bpow1.eb = 7.99625; 
freeze(bpow1.ref);
sherpa> bpow1.ampl = 0.001;

sherpa> create_model_component("powlaw1d","pow1d");
sherpa> pow1d.gamma = 1; 
sherpa> pow1d.ampl = 1e-5;

Note that the normalization parameters for each model component will be 1 keV by default, which we want, so we freeze 'bpow1.ref' and 'pow1d.ref' without changing the value (model parameter values are defined in energy space, even if the 'set_analysis("wave")' setting has been established). For the "bpow1" and "pow1d" parameters for which we did set initial values, we could have used the Sherpa guess() function to estimate reasonable starting values, based on the data input to the Sherpa session. These values can be listed with print() command:

 sherpa> print(abs1);
atten.abs1
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.hcol     thawed        1e+20        1e+17        1e+24           
   abs1.heiRatio thawed          0.1            0            1           
   abs1.heiiRatio thawed         0.01            0            1           

sherpa> print(bpow1);
bpl1d.bpow1
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   bpow1.gamma1  thawed            0          -10           10           
   bpow1.gamma2  thawed            0          -10           10           
   bpow1.eb      thawed      7.99625            0        1e+03           
   bpow1.ref     frozen            1     -3.4e+38      3.4e+38           
   bpow1.ampl    thawed        0.001        1e-20      3.4e+38           

sherpa> print(pow1d);
powlaw1d.pow1d
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   pow1d.gamma  thawed            1          -10           10           
   pow1d.ref    frozen            1     -3.4e+38      3.4e+38           
   pow1d.ampl   thawed        1e-05            0      3.4e+38

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

sherpa> abs1.hcol = 1.81e20 
sherpa> freeze(abs1)

The hydrogen column density (hcol) is set to the Galactic value. All the abs1 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 abs1 and bpow1 is assigned as the source model for all data sets :

sherpa> set_model(1, abs1*bpow1);
sherpa> set_model(2, abs1*bpow1);
sherpa> set_model(3, abs1*bpow1);
sherpa> set_model(4, abs1*bpow1);

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

sherpa> set_bkg_model(1, abs1*pow1d, 1);
sherpa> set_bkg_model(1, abs1*pow1d, 2);
sherpa> set_bkg_model(2, abs1*pow1d, 1);
sherpa> set_bkg_model(2, abs1*pow1d, 2);
sherpa> set_bkg_model(3, abs1*pow1d, 1);
sherpa> set_bkg_model(3, abs1*pow1d, 2);
sherpa> set_bkg_model(4, abs1*pow1d, 1);
sherpa> set_bkg_model(4, abs1*pow1d, 2);

The source and background model definitions can be listed with get_model() and get_bkg_model() :

sherpa> print(get_model());
apply_arf((38564.6089269 * ((atten.abs1 * bpl1d.bpow1) + scale_factor * ((atten.abs1 * powlaw1d.pow1d) + (atten.abs1 * powlaw1d.pow1d)))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.hcol    frozen     1.81e+20        1e+17        1e+24           
   abs1.heiRatio frozen          0.1            0            1           
   abs1.heiiRatio frozen         0.01            0            1           
   bpow1.gamma1 thawed            0          -10           10           
   bpow1.gamma2 thawed            0          -10           10           
   bpow1.eb     thawed      7.99625            0         1000           
   bpow1.ref    frozen            1 -3.40282e+38  3.40282e+38           
   bpow1.ampl   thawed        0.001        1e-20  3.40282e+38           
   pow1d.gamma  thawed            1          -10           10           
   pow1d.ref    frozen            1 -3.40282e+38  3.40282e+38           
   pow1d.ampl   thawed        1e-05            0  3.40282e+38          

sherpa> print(get_bkg_model());
apply_arf((38564.6089269 * (atten.abs1 * powlaw1d.pow1d)))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.hcol    frozen     1.81e+20        1e+17        1e+24           
   abs1.heiRatio frozen          0.1            0            1           
   abs1.heiiRatio frozen         0.01            0            1           
   pow1d.gamma  thawed            1          -10           10           
   pow1d.ref    frozen            1 -3.40282e+38  3.40282e+38           
   pow1d.ampl   thawed        1e-05            0  3.40282e+38

Examining Method & Statistic Settings

Next we check the current method and statistics settings:

sherpa> show_method();
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

sherpa> show_stat();
Statistic: Chi2Gehrels

The Sherpa default fitting statistic and optimization method are Chi2Gehrels and LevMar, respectively. For this fit, we will use the Neldermead method and the CStat statistic - Chi2Gehrels could bias the fit results and yield an artificially low reduced statistic for this data, and the CStat (and Cash) statistic is appropriate for simultaneously fitting source and background data. For a list of all the available methods and statistic settings, see the Sherpa Statistics and Optimization Methods pages.

To change the current method and statistic, we use set_method and set_stat.

sherpa> set_method('neldermead');
sherpa> set_stat("cstat");

Fitting

The data sets are now fit:

sherpa> fit();   
Datasets              = 1, 2, 3, 4
Method                = neldermead
Statistic             = cstat
Initial fit statistic = 336328
Final fit statistic   = 5870.13 at function evaluation 1423
Data points           = 5052
Degrees of freedom    = 5046
Probability [Q-value] = 3.0515e-15
Reduced statistic     = 1.16332
Change in statistic   = 330458
   bpow1.gamma1   1.588       
   bpow1.gamma2   9.99999     
   bpow1.eb       7.27685     
   bpow1.ampl     0.023751    
   pow1d.gamma    1.66834     
   pow1d.ampl     0.000477663

To plot the fits:

sherpa> plot("fit", 1, "fit", 2, "fit", 3, "fit", 4);
WARNING: unable to calculate errors using current statistic: cstat
WARNING: unable to calculate errors using current statistic: cstat
WARNING: unable to calculate errors using current statistic: cstat
WARNING: unable to calculate errors using current statistic: cstat

sherpa> current_plot("all");
sherpa> set_plot_title("3C 273 (ObsID 459)");

sherpa> current_plot("plot1");
sherpa> add_label(12, 0.075, "HEG -1");
sherpa> set_label({"color","green"});
sherpa> current_plot("plot2");
sherpa> add_label(12, 0.075, "HEG +1");
sherpa> set_label({"color","green"});
sherpa> current_plot("plot3");
sherpa> add_label(12, 0.125, "MEG -1");
sherpa> set_label({"color","green"});
sherpa> current_plot("plot4");
sherpa> add_label(12, 0.125, "MEG +1");
sherpa> set_label({"color","green"});

The ChIPS commands are used to add a title and labels to the drawing area. The plot is shown in Figure 3.

[Labeled plots of the simultaneous fit on ACIS HEG and MEG +/- 1 orders of 3C 273]

Figure 3: Results of Simultaneous Fit

Labeled plots of the simultaneous fit on ACIS HEG and MEG +/- 1 orders of 3C 273.

Note that the cstat statistic does not calculate errors for the data points. Since it is useful to do so, we change the fit statistic to something suitable for calculating errors, and view the residuals of the fit with the plot_fit_delchi command:

sherpa> set_stat('chi2xspecvar'); 
sherpa> plot_fit_delchi();

This plot is shown in Figure 4.

[Plot of the fit and residuals for the HEG -1 order spectrum]

[Print media version: Plot of the fit and residuals for the HEG -1 order spectrum]

Figure 4: Fit and residuals for the HEG -1 order spectrum

After creating a plot, it may be saved as a PostScript file; in this example, "all.ps" is returned:

sherpa> print_window("all");

Examining Fit Results

In CIAO 3.4, the GOODNESS command was used to get the chi-squared goodness-of-fit. This information is now reported with the best-fit values after a fit; it is no longer necessary to run a separate command. The commands get_fit_results and show_fit allow access to this information after the fit has been performed:

sherpa> show_fit();
Optimization Method: NelderMead
name         = simplex
ftol         = 1.19209289551e-07
maxfev       = None
scale        = 3.0
initsimplex  = 0
finalsimplex = 9
step         = None
seed         = 7151
iquad        = 1
verbose      = 0

Statistic: Chi2XspecVar

Fit:Datasets              = 1, 2, 3, 4
Method                = neldermead
Statistic             = cstat
Initial fit statistic = 336328
Final fit statistic   = 5870.13 at function evaluation 1423
Data points           = 5052
Degrees of freedom    = 5046
Probability [Q-value] = 3.0515e-15
Reduced statistic     = 1.16332
Change in statistic   = 330458
   bpow1.gamma1   1.588       
   bpow1.gamma2   9.99999     
   bpow1.eb       7.27685     
   bpow1.ampl     0.023751    
   pow1d.gamma    1.66834     
   pow1d.ampl     0.000477663

The number of bins in the fit (Data points), the number of degrees of freedom (i.e. the number of bins minus the number of free parameters), and the final fit statistic value are reported. If the chosen statistic is one of the chi-square statistics, as in this example, the reduced statistic, i.e. the statistic value divided by the number of degrees of freedom, and the probability (Q-value) are included as well.

The calc_chisqr command calculates the statistic contribution per bin; in this example, the results for data set 1 are returned:

sherpa> print(calc_chisqr());

 [ 0.,  0.,  0., 0., ... , 2.42532909e-08,   3.06920209e-09,   1.39313146e-09,
         1.18245634e-10,   4.34819228e-11,   1.37782796e-12]

The covariance (covar) and projection (proj) commands can be used to estimate confidence intervals for the thawed parameters:

sherpa> set_stat("cstat");
sherpa> covar();
Datasets              = 1, 2, 3, 4
Confidence Method     = covariance
Fitting Method        = neldermead
Statistic             = cstat
covariance 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   bpow1.gamma1        1.588  -0.00632353   0.00632353
   bpow1.gamma2      9.99999     -1.28663      1.28663
   bpow1.eb          7.27685   -0.0726007    0.0726007
   bpow1.ampl       0.023751 -0.000116491  0.000116491
   pow1d.gamma       1.66834   -0.0325026    0.0325026
   pow1d.ampl    0.000477663  -1.1824e-05   1.1824e-05

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("459_fitting_session.save");

All the information about the current session is written to 459_fitting_session.save, a binary file. It may be loaded into Sherpa again with the restore() command.

Finally, quit the session:

sherpa> quit

Scripting It

The file fit.sl is an S-Lang script which performs the primary commands used above; it can be executed by typing execfile("fit.sl") 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.

History

18 Jul 2008	updated for CIAO 4.1
04 Dec 2008	`set_analysis()`, `show_data()`, and `show_fit()` are available in Sherpa 4.1
12 Dec 2008	`create_model_component` is available in Sherpa 4.1
28 Apr 2009	replaced use of `atten` model with Sherpa user model "`atten_wave`"; new `script` command is available with CIAO 4.1.2