Fitting Multiple Orders of HRCS/LETG Data
Sherpa Threads (CIAO 4.9 Sherpa v1)
Overview
Synopsis:
Because of the low energy resolution in the HRCS, 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: 14 Nov 2016  reviewed for CIAO 4.9; no content change, updated fit results.
Contents
 Getting Started
 Reading the Spectrum Files
 Loading the Instrument Responses
 Plotting Orders with ChIPS
 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 with ChIPS
 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 negativeorder spectrum
 Figure 6: Plotting one order
 Figure 7: Plotting a combination of orders
Getting Started
Download the sample data: 460 (LETG/HRCS, 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.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") 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, reread 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, reread 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 negativeorder spectrum as data set 1 and the positiveorder 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.
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] dict_keys([1, 2, 3]) # Python3 sherpa> show_data() Data Set: 1 Filter: 0.060312.1032 Energy (keV) Bkg Scale: 0.2 Noticed Channels: 116384 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] 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 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] 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 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] 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 Background Data Set: 1:1 Filter: 0.060312.1032 Energy (keV) Noticed Channels: 116384 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.060312.1032 Energy (keV) Noticed Channels: 116384 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.060312.1032 Energy (keV) Bkg Scale: 0.2 Noticed Channels: 116384 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] 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 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] 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 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] 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 Background Data Set: 2:1 Filter: 0.060312.1032 Energy (keV) Noticed Channels: 116384 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.060312.1032 Energy (keV) Noticed Channels: 116384 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 yaxis 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 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> get_analysis() 'wavelength'
The data may now be plotted:
sherpa> plot("data", 1, "data", 2)
Figure 1 shows the resulting plot.
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 yaxis 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', 1, 'data', 2)
Filtering the Data
There are two plate gaps in the HRCS detector: one at ~50 Å and the other at ~60 Å; see Figure 7.1 in the POG for the HRCS detector layout. The effect of dithering into these gaps appear in negativeorder and positiveorder data, respectively, as a flat area of zero counts. The regions where the data are in a plate gap are evident in Figure 3.
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 powerlaw (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 xvalues are always interpreted as energy bins. When the analysis setting is using nonenergy 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 powerlaw 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.81e02 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.19209289551e07 xtol = 1.19209289551e07 gtol = 1.19209289551e07 maxfev = None epsfcn = 1.19209289551e07 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 moresuitable for use with lowcount 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 higherorder 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. 336346. 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 NelderMead. 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.') #reestablish the wavelength filter after regrouping 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 [Qvalue] = 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> 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.
Notice that Sherpa does not attempt to fit the regions that we chose to omit.
sherpa> show_filter() Data Set Filter: 1 1.006248.9812,56.0062183.5551 Wavelength (Angstrom) Data Set Filter: 2 1.012358.9812,66.0062189.6353 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.
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 hiresolution model [orders 13] 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 yaxis to log scale sherpa> linear_scale(Y_AXIS) # set yaxis to linear scale
Figure 6 displays the resulting plot.
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(r"normalized counts sec^{1} \\AA^{1}") sherpa> set_plot_xlabel(r"m\\lambda [\\AA]") sherpa> print_window("460_leg_m1_bin10", ["format","png","clobber","true"])
OR
sherpa> plot_data() sherpa> plot_model(overplot=1) # overplot hiresolution model [orders 13] 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.
Examining Fit Results
The χ^{2} goodnessoffit information is reported with the bestfit values after a fit is performed; and the get_fit_results and show_fit commands allow subsequent access to this information postfitting:
sherpa> show_fit() Optimization Method: NelderMead name = simplex ftol = 1.19209289551e07 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. 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 [Qvalue] = 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 (Qvalue), 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 methodconf:
sherpa> covar() Datasets = 1, 2 Confidence Method = covariance Iterative Fit Method = None Fitting Method = neldermead Statistic = chi2xspecvar covariance 1sigma (68.2689%) bounds: Param BestFit 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 1sigma (68.2689%) bounds: Param BestFit 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.755224843601643e10
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 exec(open("fit.py").read()) 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.)
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.9 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 
10 Jan 2010  updated for CIAO 4.2 
13 Jul 2010  updated for CIAO 4.2 Sherpa v2: removal of SLang 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. 