Fitting Multiple Orders of HRC-S/LETG Data

Sherpa Threads (CIAO 4.1)

[Python Syntax]

Overview

Last Update: 29 Apr 2009 - new script command is available with CIAO 4.1.2

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) created with 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

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

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

arfs:
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 load_pha command (a subset of the load_data command):

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

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

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

Loading the Instrument Responses

The individual instrument response files (gRMFs and gARFs) need to be read into Sherpa for each of the six spectral orders (+/- 1,2,3). (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 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/load_rmf or load_multi_arfs/load_multi_rmfs commands.

Several options are available for loading multiple spectral responses:

sherpa> load_multi_arfs(1, ["460_LEG_-1.garf","460_LEG_-2.garf","460_LEG_-3.garf"], [1,2,3])
sherpa> load_multi_rmfs(1, ["460_leg_-1.grmf","460_leg_-2.grmf","460_leg_-3.grmf"], [1,2,3])

sherpa> load_multi_arfs(2, ["460_LEG_1.garf","460_LEG_2.garf","460_LEG_3.garf"], [1,2,3])
sherpa> load_multi_rmfs(2, ["460_leg_1.grmf","460_leg_2.grmf","460_leg_3.grmf"], [1,2,3])

sherpa> for num in range(1,4):
            load_arf(1, "460_LEG_-%i.garf" % num, num)
            load_rmf(1, "460_leg_-%i.grmf" % num, num)

sherpa> for num in range(1,4):
            load_arf(2, "460_LEG_%i.garf" % num, num)
            load_rmf(2, "460_leg_%i.grmf" % num, num)

sherpa> arf_ids = range(1, 4)
sherpa> rmf_ids = arf_ids[:]

sherpa> rmfs_1 = ["460_leg_-%i.grmf" % id for id in rmf_ids]
sherpa> arfs_1 = ["460_LEG_-%i.garf" % id for id in arf_ids] 

sherpa> rmfs_2 = ["460_leg_%i.grmf" % id for id in rmf_ids]
sherpa> arfs_2 = ["460_LEG_%i.garf" % id for id in arf_ids]

sherpa> load_multi_arfs(1, arfs_1, arf_ids) 
sherpa> load_multi_rmfs(1, rmfs_1, rmf_ids)

sherpa> load_multi_arfs(2, arfs_2, arf_ids) 
sherpa> load_multi_rmfs(2, rmfs_2, rmf_ids)

We may list the data set and response IDs with the following list commands or with get_data, and return information about each ARF and RMF with the get_arf and get_rmf commands.

sherpa> list_data_ids()
        [1, 2]

sherpa> list_response_ids()
        [1, 2, 3]

sherpa> print(get_data())
name           = 460_leg_m1_bin10.pha
channel        = Float64[16384]
counts         = Float64[16384]
staterror      = None
syserror       = None
bin_lo         = Float64[16384]
bin_hi         = Float64[16384]
grouping       = Int16[16384]
quality        = Int16[16384]
exposure       = 39939.2745006
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = wavelength
response_ids   = [1, 2, 3]
background_ids = [1, 2]

sherpa> print(get_arf(1,3))
name     = 460_LEG_-3.garf
energ_lo = Float64[16384]
energ_hi = Float64[16384]
specresp = Float64[16384]
bin_lo   = Float64[16384]
bin_hi   = Float64[16384]
exposure = 39911.9152023

The response orders are stored in the order in which they are loaded, noted by the 'resp_id' value. The first response in the list of response IDs is treated as the primary (usually resp_id=1), i.e. it is used as the mapping from wavelength or energy units to spectral channel. The primary response is used as the independent variable grid for plotting, as well as for translating the y-axis to the proper units (keV or Angstroms). During fitting, the calculation of the source model uses the respective ARF grid from each of the responses. Similarly, the respective response is used when folding through the response matrix.

The current definition of the instrument response for all data sets may be examined at once with the show_all or show_data commands (which also provide the option of saving the data to a file):

sherpa> show_all()   
Data Set: 1
Filter: 205.7375-56.1124,48.8624-1.0244 Wavelength (Angstrom)
Noticed Channels: 1-11980,12551-16384
name           = 460_leg_m1_bin10.pha
channel        = Int16[16384]
counts         = Int16[16384]
staterror      = None
syserror       = None
bin_lo         = Float64[16384]
bin_hi         = Float64[16384]
grouping       = Int16[16384]
quality        = Int16[16384]
exposure       = 39939.2745006
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = wavelength
response_ids   = [1, 2, 3]
background_ids = [1, 2]

RMF Data Set: 1:1
name     = 460_leg_-1.grmf
detchans = 16384
energ_lo = Float64[16384]
energ_hi = Float64[16384]
n_grp    = Int16[16384]
f_chan   = UInt32[16384]
n_chan   = UInt32[16384]
matrix   = Float64[1640230]
offset   = 1
e_min    = Float64[16384]
e_max    = Float64[16384]

ARF Data Set: 1:1
name     = 460_LEG_-1.garf
energ_lo = Float64[16384]
energ_hi = Float64[16384]
specresp = Float64[16384]
bin_lo   = Float64[16384]
bin_hi   = Float64[16384]
exposure = 39911.9152023

RMF Data Set: 1:2
name     = 460_leg_-2.grmf
detchans = 16384
energ_lo = Float64[16384]
energ_hi = Float64[16384]
n_grp    = Int16[16384]
f_chan   = UInt32[16384]
n_chan   = UInt32[16384]
matrix   = Float64[1655030]
offset   = 1
e_min    = Float64[16384]
e_max    = Float64[16384]
 
ARF Data Set: 1:2
name     = 460_LEG_-2.garf
energ_lo = Float64[16384]
energ_hi = Float64[16384]
specresp = Float64[16384]
bin_lo   = Float64[16384]
bin_hi   = Float64[16384]
exposure = 39911.9152023

RMF Data Set: 1:3
name     = 460_leg_-3.grmf
detchans = 16384
energ_lo = Float64[16384]
energ_hi = Float64[16384]
n_grp    = Int16[16384]
f_chan   = UInt32[16384]
n_chan   = UInt32[16384]
matrix   = Float64[1655236]
offset   = 1
e_min    = Float64[16384]
e_max    = Float64[16384]

ARF Data Set: 1:3
name     = 460_LEG_-3.garf
energ_lo = Float64[16384]
energ_hi = Float64[16384]
specresp = Float64[16384]
bin_lo   = Float64[16384]
bin_hi   = Float64[16384]
exposure = 39911.9152023

Background Data Set: 1:1
name           = 460_leg_m1_bin10.pha
channel        = Int16[16384]
counts         = Int16[16384]
staterror      = None
syserror       = None
bin_lo         = Float64[16384]
bin_hi         = Float64[16384]
grouping       = Int16[16384]
quality        = Int16[16384]
exposure       = 39939.2745006
backscal       = 5.0
areascal       = None
grouped        = True
subtracted     = False
units          = wavelength
response_ids   = []
background_ids = []

Background Data Set: 1:2
name           = 460_leg_m1_bin10.pha
channel        = Int16[16384]
counts         = Int16[16384]
staterror      = None
syserror       = None
bin_lo         = Float64[16384]
bin_hi         = Float64[16384]
grouping       = Int16[16384]
quality        = Int16[16384]
exposure       = 39939.2745006
backscal       = 5.0
areascal       = None
grouped        = True
subtracted     = False
units          = wavelength
response_ids   = []
background_ids = []

Data Set: 2
Filter: 205.7375-56.1124,48.8624-1.0244 Wavelength (Angstrom)
Noticed Channels: 1-11980,12551-16384
name           = 460_leg_p1_bin10.pha
channel        = Int16[16384]
counts         = Int16[16384]
staterror      = None
syserror       = None
bin_lo         = Float64[16384]
bin_hi         = Float64[16384]
grouping       = Int16[16384]
quality        = Int16[16384]
exposure       = 39939.2745006
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = wavelength
response_ids   = [1, 2, 3]
background_ids = [1, 2]

RMF Data Set: 2:1
name     = 460_leg_1.grmf
detchans = 16384
energ_lo = Float64[16384]
energ_hi = Float64[16384]
n_grp    = Int16[16384]
f_chan   = UInt32[16384]
n_chan   = UInt32[16384]
matrix   = Float64[1643640]
offset   = 1
e_min    = Float64[16384]
e_max    = Float64[16384]

ARF Data Set: 2:1
name     = 460_LEG_1.garf
energ_lo = Float64[16384]
energ_hi = Float64[16384]
specresp = Float64[16384]
bin_lo   = Float64[16384]
bin_hi   = Float64[16384]
exposure = 39911.9152023

RMF Data Set: 2:2
name     = 460_leg_2.grmf
detchans = 16384
energ_lo = Float64[16384]
energ_hi = Float64[16384]
n_grp    = Int16[16384]
f_chan   = UInt32[16384]
n_chan   = UInt32[16384]
matrix   = Float64[1655295]
offset   = 1
e_min    = Float64[16384]
e_max    = Float64[16384]

ARF Data Set: 2:2
name     = 460_LEG_2.garf
energ_lo = Float64[16384]
energ_hi = Float64[16384]
specresp = Float64[16384]
bin_lo   = Float64[16384]
bin_hi   = Float64[16384]
exposure = 39911.9152023

RMF Data Set: 2:3
name     = 460_leg_3.grmf
detchans = 16384
energ_lo = Float64[16384]
energ_hi = Float64[16384]
n_grp    = Int16[16384]
f_chan   = UInt32[16384]
n_chan   = UInt32[16384]
matrix   = Float64[1655558]
offset   = 1
e_min    = Float64[16384]
e_max    = Float64[16384]

ARF Data Set: 2:3
name     = 460_LEG_3.garf
energ_lo = Float64[16384]
energ_hi = Float64[16384]
specresp = Float64[16384]
bin_lo   = Float64[16384]
bin_hi   = Float64[16384]
exposure = 39911.9152023

Background Data Set: 2:1
name           = 460_leg_p1_bin10.pha
channel        = Int16[16384]
counts         = Int16[16384]
staterror      = None
syserror       = None
bin_lo         = Float64[16384]
bin_hi         = Float64[16384]
grouping       = Int16[16384]
quality        = Int16[16384]
exposure       = 39939.2745006
backscal       = 5.0
areascal       = None
grouped        = True
subtracted     = False
units          = wavelength
response_ids   = []
background_ids = []

Background Data Set: 2:2
name           = 460_leg_p1_bin10.pha
channel        = Int16[16384]
counts         = Int16[16384]
staterror      = None
syserror       = None
bin_lo         = Float64[16384]
bin_hi         = Float64[16384]
grouping       = Int16[16384]
quality        = Int16[16384]
exposure       = 39939.2745006
backscal       = 5.0
areascal       = None
grouped        = True
subtracted     = False
units          = wavelength
response_ids   = []
background_ids = []

Plotting Orders with ChIPS

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 for data set 1:

sherpa> set_analysis("wave")

sherpa> print(get_data().units)
wavelength

The data may now be plotted:

sherpa> plot("data", 1, "data", 2)

Figure 1 shows the resulting plot.

[Plot of HRC LEG +/- 1 orders for 3C 273]

Figure 1: Plotting the two orders

Plot of HRC LEG +/- 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.)

We can use ChIPS to adjust the features of the plot, such as setting colors, changing font styles, and repositioning objects. Below, get_data_plot_prefs is used to to connect the data points with a solid line, and to remove data point symbols and y-axis error bars from the plot. The resulting plot is displayed in Figure 2.

sherpa> prefs = get_data_plot_prefs()

sherpa> prefs["linestyle"] = chips_solid

sherpa> prefs["symbolstyle"] = chips_none

sherpa> prefs["yerrorbars"] = 0

sherpa> plot_data()

Figure 2: Modifying the plot of orders with ChIPS

Plot of HRC LEG +/- 1 orders, modified with ChIPS.

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. The regions where the data are in a plate gap are evident in Figure 3.

[Plot displaying HRC-S plate gaps in the positive and negative orders]

[Print media version: Plot displaying HRC-S plate gaps in the positive and negative orders]

Figure 3: Regions where the data dithered into a plate gap

There are two plate gaps in the HRC-S detector: one at ~50 angstroms and the other at ~60 angstroms. The effect of dithering into these gaps can be seen in the negative-order (top frame) and positive-order data (bottom frame).

These regions are ignored in the fitting so that the statistic calculations are not skewed:

sherpa> ignore_id(1, 49., 56.)
sherpa> ignore_id(2, 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 with the source data), 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).

The absorption model will be referred to as "abs1", and the broken power law will be "bpow1"; the product of the two is assigned as the source model for the data sets:

sherpa> set_model(xswabs.abs1*xsbknpower.bpow1)
sherpa> set_model(2, abs1*bpow1)

sherpa> show_model()
Model: 1
(39939.2745006 * MultiResponseSumModel(apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1)),apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1)),apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed            1            0       100000 10^22 atoms / cm^2
   bpow1.phoindx1 thawed            1           -2            9           
   bpow1.breake thawed            5         0.01        1e+06        keV
   bpow1.phoindx2 thawed            2           -2            9           
   bpow1.norm   thawed            1            0        1e+24           

Model: 2
(39939.2745006 * MultiResponseSumModel(apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1)),apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1)),apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed            1            0       100000 10^22 atoms / cm^2
   bpow1.phoindx1 thawed            1           -2            9           
   bpow1.breake thawed            5         0.01        1e+06        keV
   bpow1.phoindx2 thawed            2           -2            9           
   bpow1.norm   thawed            1            0        1e+24           


sherpa> abs1.nh = 0.1 


sherpa> bpow1.phoindx1 = 1 
sherpa> bpow1.breake = 5
sherpa> bpow1.phoindx2 = 2 
sherpa> bpow1.norm = 0.0434012
sherpa> 

sherpa> show_model()
Model: 1
(39939.2745006 * MultiResponseSumModel(apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1)),apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1)),apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed          0.1            0       100000 10^22 atoms / cm^2
   bpow1.phoindx1 thawed            1           -2            9           
   bpow1.breake thawed            5         0.01        1e+06        keV
   bpow1.phoindx2 thawed            2           -2            9           
   bpow1.norm   thawed    0.0434012            0        1e+24           

Model: 2
(39939.2745006 * MultiResponseSumModel(apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1)),apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1)),apply_rmf(apply_arf)((xswabs.abs1 * xsbknpower.bpow1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed          0.1            0       100000 10^22 atoms / cm^2
   bpow1.phoindx1 thawed            1           -2            9           
   bpow1.breake thawed            5         0.01        1e+06        keV
   bpow1.phoindx2 thawed            2           -2            9           
   bpow1.norm   thawed    0.0434012            0        1e+24

Note that if we hadn't set initial parameter values for the model, we could have used the guess() Sherpa function to estimate the initial parameter values for each model component separately, based on the data input to the session.

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

sherpa> abs1.nh = 1.81e-02
sherpa> freeze(abs1.nh)

sherpa> bpow1.breake = 1

We choose to set the break energy [keV] to a lower starting point, and 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.

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

sherpa> set_method("neldermead")
sherpa> set_method_opt("finalsimplex", 0)

sherpa> set_stat("chi2xspecvar")

For this fit, we decide to change the statistic setting from the default (Chi2Gehrels) to Chi2XSpecvar, and the fitting optimization method from LevMar to Nelder-Mead. For a list of all the available methods and statistic settings with explanations, see the Sherpa pages Optimization Methods and Statistcs. To change the current method and statistic, use set_method and set_stat.

It is important to note that when fitting data with Χ² statistics, the data must be binned so that no channels/bins have zero counts - see the Sherpa group functions for the available binning options. If the number of counts in a bin is less than 1, the variance is set to 1 with Chi2XSpecvar. We can use the group_counts function after subtracting the background to specify that at least one count should be contained in each data bin.

Subtract the Background

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

sherpa> subtract()
sherpa> subtract(2)

sherpa> group_counts(1, 1)
sherpa> group_counts(2, 1)

Fitting

The data sets are now simultaneously fit:

sherpa> fit()
 Solar Abundance Vector set to angr:  Anders E. & Grevesse N. Geochimica et Cosmochimica Acta 53, 197 (1989)
 Cross Section Table set to bcmc:  Balucinska-Church and McCammon, 1998
Datasets              = 1, 2
Method                = neldermead
Statistic             = chi2xspecvar
Initial fit statistic = 27612.1
Final fit statistic   = 23722.8 at function evaluation 777
Data points           = 23853
Degrees of freedom    = 23849
Probability [Q-value] = 0.717606
Reduced statistic     = 0.994709
Change in statistic   = 3889.31
   bpow1.phoindx1   1.8749      
   bpow1.breake   0.795073    
   bpow1.phoindx2   1.66592     
   bpow1.norm     0.0209281

To plot the fits:

sherpa> plot("fit", 1, "fit", 2)

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

sherpa> current_plot("plot1")
sherpa> add_label(125, 0.04, "LEG, -1 order")
sherpa> set_label(["color","green"])

sherpa> current_plot("plot2")
sherpa> add_label(125, 0.04, "LEG, +1 order")
sherpa> set_label(["color","green"])

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

[Plot of the model fit to the spectral orders, with plot title and labels added.]

Figure 4: Fit to the LEG +/- 1,2,3 orders of ObsID 460

Plot of the model fit to the spectral orders, with plot title and labels added.

Notice that Sherpa does not attempt to fit the regions that we chose to omit.

sherpa>  show_filter()
Data Set Filter: 1
1.0244-48.8624,56.1124-205.7375 Wavelength (Angstrom)

Data Set Filter: 2
1.0244-58.8624,66.1124-205.7375 Wavelength (Angstrom)

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

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

sherpa> plot_fit_delchi()

The plot of the negative order spectrum is shown in Figure 5.

[Print media version: Fit and residuals for the negative-order spectrum]

Figure 5: Fit and residuals for the negative-order spectrum

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

sherpa> print_window("460_fit_delchi")

Finally, we may use the plot_order and plot_model commands to plot spectral orders individually or simultaneously. In the case of the latter, we may use ChIPS commands to assign a different color to each order to separate them visually.

To plot individual orders:

sherpa> plot_data()

sherpa> plot_model(overplot=1)  # overplot hi-resolution model [orders 1-3] as histogram

sherpa> plot_order(1,3,overplot=1)    # overplot 3rd order contribution

sherpa> set_histogram("line.color=blue")

sherpa> log_scale(Y_AXIS)    # set y-axis to log scale 

sherpa> linear_scale(Y_AXIS)  # set y-axis to linear scale

Figure 6 displays the resulting plot.

[Plot of negative orders with model contribution from orders 1 and 3 shown in blue]

[Print media version: Plot of negative orders with model contribution from orders 1 and 3 shown in blue]

Figure 6: Plotting one order

Plot of negative spectral orders on a linear scale, with the convolved model shown in red and the contribution for order 3 shown in blue.

To plot a combination of orders with colors:

sherpa> plot_data()

sherpa> plot_model(overplot=1)

sherpa> plot_order(1,1,overplot=1)
sherpa> set_histogram("line.color=darkred")

sherpa> plot_order(1,2,overplot=1)
sherpa> set_histogram("line.color=pink")

sherpa> plot_order(1,3,overplot=1)
sherpa> set_histogram("line.color=orange")

sherpa> log_scale(Y_AXIS)


sherpa> title=("Model Orders +
              "{\\color{darkred}1, }" +
              "{\\color{pink}2, }" +
              "{\\color{orange}3, }")

sherpa> set_plot_title(title)

sherpa> set_plot_ylabel("normalized counts sec^{-1} \\AA^{-1}")

sherpa> set_plot_xlabel("m\\lambda [\\AA]")

sherpa> print_window("460_leg_m1_bin10", ["format","png","clobber","true"])

sherpa> plot_data()

sherpa> plot_model(overplot=1)  # overplot hi-resolution model [orders 1-3] as histogram

sherpa> get_data().response_ids = [1,2]    # force session to use only
orders 1,2

sherpa> plot_model(overplot=1)             # overplot sum of orders 1 and 2

sherpa> set_histogram("line.color=blue")

sherpa> log_scale(Y_AXIS)  

sherpa> linear_scale(Y_AXIS)

sherpa> get_data().response_ids = list_response_ids() # restore all orders

Figure 7 displays the plot which results from the first plotting example.

[Plot of negative orders with convolved model contribution from orders 1, 2, 3 shown in different colors]

[Print media version: Plot of negative orders with convolved model contribution from orders 1, 2, 3 shown in different colors]

Figure 7: Plotting a combination of orders

Plot of negative spectral orders, with the convolved model contributions from orders 1, 2, and 3 shown in dark red, pink, and orange, respectively.

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 get_fit_results and show_fit commands 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
Method                = neldermead
Statistic             = chi2xspecvar
Initial fit statistic = 27612.1
Final fit statistic   = 23722.8 at function evaluation 777
Data points           = 23853
Degrees of freedom    = 23849
Probability [Q-value] = 0.717606
Reduced statistic     = 0.994709
Change in statistic   = 3889.31
   bpow1.phoindx1   1.8749      
   bpow1.breake   0.795073    
   bpow1.phoindx2   1.66592     
   bpow1.norm     0.0209281

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> calc_chisqr()

 [ 0.,  0.,  0., 0., ... , 2.67718857e-03   3.55541447e-01
 1.58184491e-01   1.25163916e-02   2.65451065e+00   1.62082315e+00   3.17077833e-01   3.19350908e+00   2.04649545e+00]

The covar (covariance) command can be used to estimate confidence intervals for the thawed parameters - though this method may not constrain parameters where the parameter space is not smooth. In this case, we can try the projection methodproj (projection):

sherpa> covar()
WARNING: hard minimum hit for parameter bpow1.phoindx1
WARNING: hard maximum hit for parameter bpow1.phoindx1
Datasets              = 1, 2
Confidence Method     = covariance
Fitting Method        = neldermead
Statistic             = chi2xspecvar
covariance 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   bpow1.phoindx1       1.8749        -----        -----
   bpow1.breake     0.795073   -0.0222528    0.0222528
   bpow1.phoindx2      1.66592  -0.00963342   0.00963342
   bpow1.norm      0.0209281 -1.44041e-05  1.44041e-05

sherpa> set_proj_opt("fast", False)   # to force proj() to use current method 
sherpa> proj()
Datasets              = 1, 2
Confidence Method     = projection
Fitting Method        = neldermead
Statistic             = chi2xspecvar
projection 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   bpow1.phoindx1       1.8749   -0.0162188    0.0188155
   bpow1.breake     0.795073   -0.0738628    0.0633429
   bpow1.phoindx2      1.66592   -0.0147633    0.0149894
   bpow1.norm      0.0209281 -0.000403257  0.000320916

Calculate Flux

The functions calc_photon_flux and calc_energy_flux can be used to return the total integrated photon or energy flux over the full range of orders, computed over the combined high resolution ARF grid. The photon flux returned is in photons/cm^2/sec, and the energy flux in ergs/cm^2/sec.

sherpa> calc_photon_flux(id=1)
        0.080791838

sherpa> calc_energy_flux(id=1)
        3.0232320043116234e-09

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

All the information about the current session is written to 460_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.py is a Python script which performs the primary commands used above; it can be executed by typing execfile("fit.py") 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

06 Aug 2008	updated for CIAO 4.1
04 Dec 2008	`plot_order()` is available in Sherpa 4.1
13 Dec 2008	`neldermead` method used in place of `levmar`
29 Apr 2009	new `script` command is available with CIAO 4.1.2