# Fitting Multiple Orders of HRC-S/LETG Data

Sherpa Threads (CIAO 4.12 Sherpa v1)

## Overview

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

**Last Update:** 13 Dec 2019 -
replaced ChIPS commands and figures with Matplotlib equivalent.

## Contents

**Getting Started****Reading the Spectrum Files****Loading the Instrument Responses****Plotting Orders****Filtering the Data****Defining the Source Model****Examining Method & Statistic Settings****Subtracting the Background and Grouping the Data****Fitting****Examining Fit Results****Calculate Flux****Saving and Quitting the Session****Scripting It****History**-
**Images**- Figure 1: Plotting the two orders
- Figure 2: Modifying the plot of orders
- Figure 3: Regions where the data dithered into a plate gap
- Figure 4: Fit to the LEG +/- 1,2,3 orders of ObsID 460
- Figure 5: Fit and residuals for the negative-order spectrum
- Figure 6: Plotting one order
- Figure 7: Plotting a combination of orders

## Getting Started

**Download the sample data:** 460 (LETG/HRC-S, 3C 273)

unix% download_chandra_obsid 460

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.pharmfs:460_leg_-1.grmf 460_leg_-2.grmf 460_leg_-3.grmf 460_leg_1.grmf 460_leg_2.grmf 460_leg_3.grmfarfs: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") WARNING: systematic errors were not found in file '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 dataset from file 460_leg_m1_bin10.pha read background_down into a dataset from file 460_leg_m1_bin10.pha sherpa> load_pha(2, "460_leg_p1_bin10.pha") WARNING: systematic errors were not found in file '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 dataset from file 460_leg_p1_bin10.pha read background_down into a dataset 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

To establish the instrument response, the individual instrument response files (gRMFs
and gARFs) need to be read into *Sherpa* for each of the six
spectral orders (+/- 1,2,3).
If the ARF and RMF filenames are recorded 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.

**This warning does not apply to users of CIAO 4.12.1.**

In Sherpa 4.12, a bug with load_multi_arfs and load_multi_rmfs prevents fit from running. These functions can be patched using patch.py before loading the responses into the Sherpa session:

sherpa> %run -i patch.py 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])

where the contents of the patch is the Python function:

unix% cat patch.py #!/usr/bin/env python # patch the startup method of the MultiResponseSumModel *before* any # instances have been created # import numpy from sherpa.models.model import CompositeModel from sherpa.astro.instrument import MultiResponseSumModel def startup_monkey(self, cache): pha = self.pha if numpy.iterable(pha.mask): pha.notice_response(True) self.channel = pha.get_noticed_channels() self.mask = pha.get_mask() self._get_noticed_energy_list() CompositeModel.startup(self, cache) MultiResponseSumModel.startup = startup_monkey

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

OR

sherpa> for num in range(1,4): load_arf(1, "460_LEG_-{}.garf".format(num), num) load_rmf(1, "460_leg_-{}.grmf".format(num), num) sherpa> for num in range(1,4): load_arf(2, "460_LEG_{}.garf".format(num), num) load_rmf(2, "460_leg_{}.grmf".format(num), num)

OR

sherpa> arf_ids = range(1, 4) sherpa> rmf_ids = arf_ids[:] sherpa> rmfs_1 = ["460_leg_-{}.grmf".format(id) for id in rmf_ids] sherpa> arfs_1 = ["460_LEG_-{}.garf".format(id) for id in arf_ids] sherpa> rmfs_2 = ["460_leg_{}.grmf".format(id) for id in rmf_ids] sherpa> arfs_2 = ["460_LEG_{}.garf".format(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 established in the
*Sherpa* session with the
list_data_ids and list_response_ids commands. For viewing detailed
information about all loaded data sets and associated ARF and
RMF responses, the show_all,
show_data, and show_bkg
commands are available (which provide the option
to save the output data to a file with the 'outfile' argument)
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> show_data() Data Set: 1 Filter: 0.0603-12.1032 Energy (keV) Bkg Scale: 0.2 Noticed Channels: 1-16384 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 = energy rate = True plot_fac = 0 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 = UInt64[16384] f_chan = UInt32[16384] n_chan = UInt32[16384] matrix = Float64[1640230] offset = 1 e_min = Float64[16384] e_max = Float64[16384] ethresh = 1e-10 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 = 39910.3861311 ethresh = 1e-10 RMF Data Set: 1:2 name = 460_leg_-2.grmf detchans = 16384 energ_lo = Float64[16384] energ_hi = Float64[16384] n_grp = UInt64[16384] f_chan = UInt32[16384] n_chan = UInt32[16384] matrix = Float64[1655030] offset = 1 e_min = Float64[16384] e_max = Float64[16384] ethresh = 1e-10 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 = 39910.3861311 ethresh = 1e-10 RMF Data Set: 1:3 name = 460_leg_-3.grmf detchans = 16384 energ_lo = Float64[16384] energ_hi = Float64[16384] n_grp = UInt64[16384] f_chan = UInt32[16384] n_chan = UInt32[16384] matrix = Float64[1655236] offset = 1 e_min = Float64[16384] e_max = Float64[16384] ethresh = 1e-10 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 = 39910.3861311 ethresh = 1e-10 Background Data Set: 1:1 Filter: 0.0603-12.1032 Energy (keV) Noticed Channels: 1-16384 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 = 5.0 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [] background_ids = [] Background Data Set: 1:2 Filter: 0.0603-12.1032 Energy (keV) Noticed Channels: 1-16384 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 = 5.0 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [] background_ids = [] Data Set: 2 Filter: 0.0603-12.1032 Energy (keV) Bkg Scale: 0.2 Noticed Channels: 1-16384 name = 460_leg_p1_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 = energy rate = True plot_fac = 0 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 = UInt64[16384] f_chan = UInt32[16384] n_chan = UInt32[16384] matrix = Float64[1643640] offset = 1 e_min = Float64[16384] e_max = Float64[16384] ethresh = 1e-10 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 = 39910.3861311 ethresh = 1e-10 RMF Data Set: 2:2 name = 460_leg_2.grmf detchans = 16384 energ_lo = Float64[16384] energ_hi = Float64[16384] n_grp = UInt64[16384] f_chan = UInt32[16384] n_chan = UInt32[16384] matrix = Float64[1655295] offset = 1 e_min = Float64[16384] e_max = Float64[16384] ethresh = 1e-10 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 = 39910.3861311 ethresh = 1e-10 RMF Data Set: 2:3 name = 460_leg_3.grmf detchans = 16384 energ_lo = Float64[16384] energ_hi = Float64[16384] n_grp = UInt64[16384] f_chan = UInt32[16384] n_chan = UInt32[16384] matrix = Float64[1655558] offset = 1 e_min = Float64[16384] e_max = Float64[16384] ethresh = 1e-10 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 = 39910.3861311 ethresh = 1e-10 Background Data Set: 2:1 Filter: 0.0603-12.1032 Energy (keV) Noticed Channels: 1-16384 name = 460_leg_p1_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 = 5.0 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [] background_ids = [] Background Data Set: 2:2 Filter: 0.0603-12.1032 Energy (keV) Noticed Channels: 1-16384 name = 460_leg_p1_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 = 5.0 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [] background_ids = []

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.

## Plotting Orders

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> get_analysis() 'wavelength'

The data may now be plotted:

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

Figure 1 shows the resulting plot.

### Figure 1: Plotting the two orders

We can use Matplotlib 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"] = "solid" sherpa> prefs["marker"] = "None" sherpa> prefs["yerrorbars"] = False sherpa> plot('data', 1, 'data', 2)

### Figure 2: Modifying the plot of orders

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

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

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 subtracting the background from 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 are defined such that the x-values are always
interpreted as 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_source(xswabs.abs1*xsbknpower.bpow1)
sherpa> set_source(2, abs1*bpow1)
sherpa> show_model()
Model: 1
(39939.274500594 * 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.274500594 * 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.norm = 0.0434012
sherpa> show_model()
Model: 1
(39939.274500594 * 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.274500594 * 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 had not 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. To have
*Sherpa* automatically query for initial parameter values
when a model is established, set
'paramprompt(True)' (it is 'False' by
default).

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 Chi Squared with Gehrels variance. The variance is estimated from the number of counts in each bin, but unlike `Chi2DataVar`, the Gaussian approximation is not used. This makes it more-suitable for use with low-count data. The standard deviation for each bin is calculated using the approximation from [1]_: sigma(i,S) = 1 + sqrt(N(i,s) + 0.75) where the higher-order terms have been dropped. This is accurate to approximately one percent. For data where the background has not been subtracted then the error term is: sigma(i) = sigma(i,S) whereas with background subtraction, sigma(i)^2 = sigma(i,S)^2 + [A(S)/A(B)]^2 sigma(i,B)^2 Notes ----- The accuracy of the error term when the background has been subtracted has not been determined. A preferable approach to background subtraction is to model the background as well as the source signal. References ---------- .. [1] "Confidence limits for small numbers of events in astrophysical data", Gehrels, N. 1986, ApJ, vol 303, p. 336-346. http://adsabs.harvard.edu/abs/1986ApJ...303..336G 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.

## Subtracting the Background and Grouping the Data

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

sherpa> subtract() sherpa> subtract(2)

It is
important to note that when fitting data with
χ^{2} 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.
Since the original data was grouped, regrouping the data
clears the filter and will be required to be set again;
however, if the original data had been ungrouped, then
the ignore filter would be retained.

sherpa> group_counts(1, 1) sherpa> group_counts(2, 1) sherpa> ignore_id(1,'49.:56.') #re-establish the wavelength filter after re-grouping sherpa> ignore_id(2, '59.:66.')

## Fitting

The data sets are now simultaneously fit:

sherpa> fit() Datasets = 1, 2 Method = neldermead Statistic = chi2xspecvar Initial fit statistic = 38106.3 Final fit statistic = 22570.5 at function evaluation 368 Data points = 23212 Degrees of freedom = 23208 Probability [Q-value] = 0.998575 Reduced statistic = 0.972532 Change in statistic = 15535.8 bpow1.PhoIndx1 1.89984 bpow1.BreakE 0.746243 bpow1.PhoIndx2 1.66696 bpow1.norm 0.0205259

To plot the fits:

sherpa> plot("fit", 1, "fit", 2) sherpa> plt.suptitle("3C 273 (ObsID 460)") sherpa> plt.sca(plt.gcf().axes[0]) sherpa> plt.title("") sherpa> plt.text(x=125, y=0.04, s="LEG, -1 order",color="green") sherpa> plt.sca(plt.gcf().axes[1]) sherpa> plt.title("") sherpa> plt.text(x=125, y=0.04, s="LEG, +1 order",color="green")

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

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

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

```
sherpa> show_filter()
Data Set Filter: 1
1.0062-48.9812,56.0062-183.5551 Wavelength (Angstrom)
Data Set Filter: 2
1.0123-58.9812,66.0062-189.6353 Wavelength (Angstrom)
```

By omitting the regions of data over a plate gap, the residuals are contained within approximately 3σ. 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.

### 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> plt.savefig("460_fit_delchi.ps")

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 *Matplotlib* pyplot commands to assign a different color to
each order to separate them visually.

To plot individual orders:

sherpa> plot_data() sherpa> plot_model(overplot=True) # overplot hi-resolution model [orders 1-3] as histogram sherpa> plot_order(1,3,overplot=True) # overplot 3rd order contribution sherpa> plt.ycale("log") # set y-axis to log scale sherpa> plt.yscale("linear") # set y-axis to linear scale

Figure 6 displays the resulting plot.

### Figure 6: Plotting one order

To plot a combination of orders with colors:

sherpa> plot_data() sherpa> plot_data() sherpa> plot_model(overplot=True) sherpa> plt.gca().lines[-1].set_color("red") sherpa> plot_order(1,1,overplot=True) sherpa> plt.gca().lines[-1].set_color("darkred") sherpa> plt.gca().lines[-1].set_label(r"$1^{st}$") sherpa> plot_order(1,2,overplot=True) sherpa> plt.gca().lines[-1].set_color("pink") sherpa> plt.gca().lines[-1].set_label(r"$2^{nd}$") sherpa> plot_order(1,3,overplot=True) sherpa> plt.gca().lines[-1].set_color("orange") sherpa> plt.gca().lines[-1].set_label(r"$3^{rd}$") sherpa> plt.yscale("log") sherpa> set_plot_title("Model Orders") sherpa> plt.legend(loc="upper right") sherpa> plt.ylabel(r"$\mathrm{normalized\ counts}\ \mathrm{sec^{-1}}\ \mathrm{\AA^{-1}}$") sherpa> plt.xlabel(r"$m\lambda [\AA]$") sherpa> plt.savefig("460_leg_m1_bin10.png")

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

### Figure 7: Plotting a combination of orders

## Examining Fit Results

The χ^{2} goodness-of-fit information is reported with the
best-fit values after a fit is performed; and the
get_fit_results and show_fit
commands allow subsequent access to this information post-fitting:

sherpa> show_fit() Optimization Method: NelderMead name = simplex ftol = 1.19209289551e-07 maxfev = None initsimplex = 0 finalsimplex = 0 step = None iquad = 1 verbose = 0 Statistic: Chi2XspecVar Chi Squared with data variance (XSPEC style). The variance in each bin is estimated from the data value in that bin. See also `Chi2DataVar`, `Chi2Gehrels`, and `Chi2ModVar`. The calculation of the variance is the same as `Chi2DataVar` except that if the number of counts in a bin is less than 1 then the variance for that bin is set to 1. See Also -------- Chi2DataVar, Chi2Gehrels, Chi2ModVar Fit:Datasets = 1, 2 Method = neldermead Statistic = chi2xspecvar Initial fit statistic = 38106.3 Final fit statistic = 22570.5 at function evaluation 368 Data points = 23212 Degrees of freedom = 23208 Probability [Q-value] = 0.998575 Reduced statistic = 0.972532 Change in statistic = 15535.8 bpow1.PhoIndx1 1.89984 bpow1.BreakE 0.746243 bpow1.PhoIndx2 1.66696 bpow1.norm 0.0205259

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 χ^{2} 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() array([ 0.1509434 , 0.62745098, 11.13245033, ..., 0.3170667 , 3.19345547, 2.04647078])

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 confidence method (conf):

sherpa> covar() Datasets = 1, 2 Confidence Method = covariance Iterative Fit Method = None Fitting Method = neldermead Statistic = chi2xspecvar covariance 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- bpow1.PhoIndx1 1.89984 -0.0161525 0.0161525 bpow1.BreakE 0.746243 -0.0221191 0.0221191 bpow1.PhoIndx2 1.66696 -0.0126385 0.0126385 bpow1.norm 0.0205259 -0.000273059 0.000273059 sherpa> set_conf_opt("fast", False) # to force conf() to use current method, by default set to "False" sherpa> conf() bpow1.PhoIndx1 lower bound: -0.021153 bpow1.PhoIndx2 lower bound: -0.0151387 bpow1.norm lower bound: -0.000348323 bpow1.BreakE lower bound: -0.0377974 bpow1.PhoIndx1 upper bound: 0.0190064 bpow1.BreakE upper bound: 0.0749116 bpow1.norm upper bound: 0.000418121 bpow1.PhoIndx2 upper bound: 0.0138105 Datasets = 1, 2 Confidence Method = confidence Iterative Fit Method = None Fitting Method = neldermead Statistic = chi2xspecvar confidence 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- bpow1.PhoIndx1 1.89984 -0.021153 0.0190064 bpow1.BreakE 0.746243 -0.0377974 0.0749116 bpow1.PhoIndx2 1.66696 -0.0151387 0.0138105 bpow1.norm 0.0205259 -0.000348323 0.000418121

## 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.077839687728646093 sherpa> calc_energy_flux(id=1) 1.755224843601643e-10

## 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 %run -i 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)

(Note that restoring a *Sherpa* session from such a file
could be problematic since it may include syntax errors,
unwanted fitting trials, et cetera.)

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

10 Jan 2010 | updated for CIAO 4.2 |

13 Jul 2010 | updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread. |

22 Jan 2012 | reviewed for CIAO 4.4 (no changes) |

13 Dec 2012 | reviewed for CIAO 4.5: group commands no longer clear the existing data filter |

03 Mar 2014 | reviewed for CIAO 4.6: no changes |

26 Feb 2015 | updated for CIAO 4.7: no content change |

14 Dec 2015 | updated for CIAO 4.8, removed references to CIAO 3.4. |

14 Nov 2016 | reviewed for CIAO 4.9; no content change, updated fit results. |

30 May 2018 | reviewed for CIAO 4.10, no content change |

11 Dec 2018 | reviewed for CIAO 4.11, no content change |

13 Dec 2019 | replaced ChIPS commands and figures with Matplotlib equivalent. |