Fitting Grating Data
Sherpa Threads (CIAO 4.14 Sherpa)
Overview
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.
Last Update: 29 Mar 2029 - reviewed for CIAO 4.14, updated fit results.
Contents
- 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
Download the sample data: 459 (HETG/ACIS-S, 3C 273)
unix% download_chandra_obsid 459
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 source data are input to Sherpa with the load_pha command, where any associated background or response files will automatically be read in if the corresponding filenames are recorded in the header of the source data files:
sherpa> load_pha(1, "459_heg_m1_bin10.pha") WARNING: systematic errors were not found in file '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") WARNING: systematic errors were not found in file '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") WARNING: systematic errors were not found in file '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") WARNING: systematic errors were not found in file '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
Note that there are two background data sets associated with each Chandra source grating data set, 'background up' and 'background down', which Sherpa assigns to background IDs 1 and 2. These contain data extracted from background regions adjacent to the source region. The two background data sets are added and scaled to the source data when being subtracted with the Sherpa subtract command, or modeled with set_bkg_model.
In order to work with only one of the two background data sets associated with a source, the get_data command may be used to remove the unneeded background, as shown below:
sherpa> print(get_data(1).background_ids) [1,2] sherpa> get_data(1).background_ids = [1] # set background 1 as the only background to fit.
Loading the Instrument Responses
The instrument response is established when the appropriate response files (ARF, RMF) are input to the Sherpa session. If the names of the ARF and RMF response files are recorded in the header of the PHA file, Sherpa will load them automatically when the PHA file is read with load_pha; 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") sherpa> load_arf(1, "459_heg_m1.arf", bkg_id=1) sherpa> load_arf(1, "459_heg_m1.arf", bkg_id=2) sherpa> load_arf(2, "459_heg_p1.arf", bkg_id=1) sherpa> load_arf(2, "459_heg_p1.arf", bkg_id=2) sherpa> load_arf(3, "459_meg_m1.arf", bkg_id=1) sherpa> load_arf(3, "459_meg_m1.arf", bkg_id=2) sherpa> load_arf(4, "459_meg_p1.arf", bkg_id=1) sherpa> load_arf(4, "459_meg_p1.arf", bkg_id=2)
The current definition of the instrument response, along with information about the loaded data sets, may be examined using the commands show_data and show_bkg:
sherpa> show_data() Data Set: 1 Filter: 0.5772-12.3984 Energy (keV) Bkg Scale 1: 0.124414 Bkg Scale 2: 0.124414 Noticed Channels: 1-8192 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.608926889 backscal = 1.0 areascal = 1.0 grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [1, 2] 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.284798602 ethresh = 1e-10 Background Data Set: 1:1 Filter: 0.5772-12.3984 Energy (keV) Noticed Channels: 1-8192 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.608926889 backscal = 4.0188284 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [] 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.284798602 ethresh = 1e-10 Background Data Set: 1:2 Filter: 0.5772-12.3984 Energy (keV) Noticed Channels: 1-8192 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.608926889 backscal = 4.0188284 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [] Background ARF Data Set: 1:2 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.284798602 ethresh = 1e-10 Data Set: 2 Filter: 0.5772-12.3984 Energy (keV) Bkg Scale 1: 0.124414 Bkg Scale 2: 0.124414 Noticed Channels: 1-8192 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.608926889 backscal = 1.0 areascal = 1.0 grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [1, 2] 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.013342202 ethresh = 1e-10 Background Data Set: 2:1 Filter: 0.5772-12.3984 Energy (keV) Noticed Channels: 1-8192 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.608926889 backscal = 4.0188284 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [] 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.013342202 ethresh = 1e-10 Background Data Set: 2:2 Filter: 0.5772-12.3984 Energy (keV) Noticed Channels: 1-8192 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.608926889 backscal = 4.0188284 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [] Background ARF Data Set: 2:2 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.013342202 ethresh = 1e-10 Data Set: 3 Filter: 0.2955-12.3984 Energy (keV) Bkg Scale 1: 0.124414 Bkg Scale 2: 0.124414 Noticed Channels: 1-8192 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.608926889 backscal = 1.0 areascal = 1.0 grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [1, 2] 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.285573929 ethresh = 1e-10 Background Data Set: 3:1 Filter: 0.2955-12.3984 Energy (keV) Noticed Channels: 1-8192 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.608926889 backscal = 4.0188284 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [] 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.285573929 ethresh = 1e-10 Background Data Set: 3:2 Filter: 0.2955-12.3984 Energy (keV) Noticed Channels: 1-8192 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.608926889 backscal = 4.0188284 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [] Background ARF Data Set: 3:2 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.285573929 ethresh = 1e-10 Data Set: 4 Filter: 0.2955-12.3984 Energy (keV) Bkg Scale 1: 0.124414 Bkg Scale 2: 0.124414 Noticed Channels: 1-8192 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.608926889 backscal = 1.0 areascal = 1.0 grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [1, 2] 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.014116893 ethresh = 1e-10 Background Data Set: 4:1 Filter: 0.2955-12.3984 Energy (keV) Noticed Channels: 1-8192 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.608926889 backscal = 4.0188284 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [] 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.014116893 ethresh = 1e-10 Background Data Set: 4:2 Filter: 0.2955-12.3984 Energy (keV) Noticed Channels: 1-8192 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.608926889 backscal = 4.0188284 areascal = None grouped = True subtracted = False units = energy rate = True plot_fac = 0 response_ids = [1] background_ids = [] Background ARF Data Set: 4:2 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.014116893 ethresh = 1e-10
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.5772-12.3984 Energy (keV) Data Set Filter: 2 0.5772-12.3984 Energy (keV) Data Set Filter: 3 0.2955-12.3984 Energy (keV) Data Set Filter: 4 0.2955-12.3984 Energy (keV) sherpa> set_analysis("wave") sherpa> show_filter() Data Set Filter: 1 1.0000-21.4800 Wavelength (Angstrom) Data Set Filter: 2 1.0000-21.4800 Wavelength (Angstrom) Data Set Filter: 3 1.0000-41.9600 Wavelength (Angstrom) Data Set Filter: 4 1.0000-41.9600 Wavelength (Angstrom)
The data may now be plotted:
sherpa> plot("data", 1, "data", 2, "data", 3, "data", 4)
Figure 1 shows the resulting plot.
Figure 1: Plotting the four orders
Filtering the Data
We choose to filter the data to focus on an area of interest:
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.
Figure 2: Filtering the data sets
Defining the Source and Background Models
We plan on simultaneously fitting the background data with the source data (rather than subtracting it), so we need to create a model expression for both the source and background. When modeling (or subtracting) the grating background, Sherpa adds the two background up/down counts arrays and scales the result to the associated source data.
We model this source with a broken power-law (xsbknpower) absorbed by the interstellar medium (ISM) (atten). The background will be modeled by a one-dimensional power-law (xspowerlaw), also absorbed by the ISM (the same atten model).
First, we set up 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 is "bpow1", and the 1-D power-law is "pow1d":
sherpa> create_model_component("atten", "abs1") <Atten model instance 'atten.abs1'> sherpa> abs1.hcol = 1e+20 sherpa> abs1.heiRatio = 0.1 sherpa> abs1.heiiRatio = 0.01 sherpa> create_model_component("xsbknpower", "bpow1") <XSbknpower model instance 'xsbknpower.bpow1'> sherpa> bpow1.PhoIndx1 = 0 sherpa> bpow1.PhoIndx2 = 0 sherpa> bpow1.BreakE = 1.55 # ~8 Angstroms sherpa> bpow1.norm = 0.001 sherpa> create_model_component("xspowerlaw","pow1d") <XSpowerlaw model instance 'xspowerlaw.pow1d'> sherpa> pow1d.PhoIndex = 1 sherpa> pow1d.norm = 1e-5
We freeze the normalization parameters for bpow1 and pow1d (bpow1.ref and pow1d.ref) without changing the default values. 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. To have Sherpa automatically query for initial parameter values when a model is established, set 'paramprompt(True)' (it is 'False' by default).
The model parameter values can be listed with print() command (note that show_model/show_source is appropriate once the full model expression has been assigned to the data with set_source:
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) xsbknpower.bpow1 Param Type Value Min Max Units ----- ---- ----- --- --- ----- bpow1.PhoIndx1 thawed 0 -3 10 bpow1.BreakE thawed 1.55 0.01 1e+06 keV bpow1.PhoIndx2 thawed 0 -3 10 bpow1.norm thawed 0.001 0 1e+24 sherpa> print(pow1d) xspowerlaw.pow1d Param Type Value Min Max Units ----- ---- ----- --- --- ----- pow1d.PhoIndex thawed 1 -3 10 pow1d.norm thawed 1e-05 0 1e+24
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> mdl = abs * bpow1 sherpa> set_source(1, mdl) sherpa> set_source(2, mdl) sherpa> set_source(3, mdl) sherpa> set_source(4, mdl)
while the background model is set as the product of abs1 and pow1d:
sherpa> bmdl = abs1 * pow1d sherpa> set_bkg_model(1, bmdl, 1) sherpa> set_bkg_model(1, bmdl, 2) sherpa> set_bkg_model(2, bmdl, 1) sherpa> set_bkg_model(2, bmdl, 2) sherpa> set_bkg_model(3, bmdl, 1) sherpa> set_bkg_model(3, bmdl, 2) sherpa> set_bkg_model(4, bmdl, 1) sherpa> set_bkg_model(4, bmdl, 2)
The source and background model definitions can be listed with show_model and show_bkg_model:
sherpa> show_model() Model: 1 apply_arf((38564.608926889 * ((atten.abs1 * xsbknpower.bpow1) + (0.24882873824620128 * (atten.abs1 * xspowerlaw.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.PhoIndx1 thawed 0 -3 10 bpow1.BreakE thawed 1.55 0 1e+06 keV bpow1.PhoIndx2 thawed 0 -3 10 bpow1.norm thawed 0.001 0 1e+24 pow1d.PhoIndex thawed 1 -3 10 pow1d.norm thawed 1e-05 0 1e+24 Model: 2 {same as above} Model: 3 {same as above} Model: 4 {same as above} sherpa> show_bkg_model() Background Model: 1:1 apply_arf((38564.608926889 * (atten.abs1 * xspowerlaw.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.PhoIndex thawed 1 -3 10 pow1d.norm thawed 1e-05 0 1e+24 Background Model: 1:2 {same as above} Background Model: 2:2 {same as above} Background Model: 3:2 {same as above} Background Model: 4:2 {same as above}
Examining Method & Statistic Settings
Next we check the current method and statistics settings:
sherpa> show_method() Optimization Method: LevMar name = levmar ftol = 1.1920928955078125e-07 xtol = 1.1920928955078125e-07 gtol = 1.1920928955078125e-07 maxfev = None epsfcn = 1.1920928955078125e-07 factor = 100.0 numcores = 1 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 A(B) is the off-source "area", which could be the size of the region from which the background is extracted, or the length of a background time segment, or a product of the two, etc.; and A(S) is the on-source "area". These terms may be defined for a particular type of data: for example, PHA data sets A(B) to `BACKSCAL * EXPOSURE` from the background data set and A(S) to `BACKSCAL * EXPOSURE` from the source data set. See Also -------- Chi2DataVar, Chi2ModVar, Chi2XspecVar 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
The Sherpa default fitting statistic and optimization method are χ^{2}-Gehrels and Levenberg-Marquardt, respectively. For this fit, we will use the Nelder-Mead Simplex method and the CStat statistic; this is because χ^{2}-Gehrels 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 (we first fit the background datasets and then the combined source and background data although we could have just called fit to do both):
sherpa> fit_bkg() Datasets = 1, 2, 3, 4 Method = neldermead Statistic = cstat Initial fit statistic = 22665.8 Final fit statistic = 2913.89 at function evaluation 278 Data points = 2528 Degrees of freedom = 2526 Probability [Q-value] = 9.45287e-08 Reduced statistic = 1.15356 Change in statistic = 19751.9 pow1d.PhoIndex 1.66743 pow1d.norm 0.000467667 sherpa> fit() Datasets = 1, 2, 3, 4 Method = neldermead Statistic = cstat Initial fit statistic = 273796 Final fit statistic = 4442.47 at function evaluation 1566 Data points = 3792 Degrees of freedom = 3786 Probability [Q-value] = 4.18275e-13 Reduced statistic = 1.17339 Change in statistic = 269354 bpow1.PhoIndx1 2.14733 bpow1.BreakE 1.11393 bpow1.PhoIndx2 1.56252 bpow1.norm 0.0247635 pow1d.PhoIndex 1.66745 pow1d.norm 0.000467698
To plot the fits:
sherpa> plot("fit", 1, "fit", 2, "fit", 3, "fit", 4) WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with cstat WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with cstat WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with cstat WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with cstat sherpa> plt.suptitle("3C 273 (ObsID 459)") sherpa> for ax in plt.gcf().axes: plt.sca(ax) ; plt.title("") sherpa> plt.gcf().axes[0].text(x=10,y=0.1,s="HEG -1",color="lime") sherpa> plt.gcf().axes[1].text(x=10,y=0.1,s="HEG +1",color="lime") sherpa> plt.gcf().axes[2].text(x=10,y=0.15,s="MEG -1",color="lime") sherpa> plt.gcf().axes[3].text(x=10,y=0.15,s="MEG +1",color="lime")
The matplotlib functionality used to add a title and labels to the drawing area are also shown. The plot is shown in Figure 3.
Figure 3: Results of Simultaneous Fit
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('chi2datavar')
sherpa> plot_fit_delchi()
This plot is shown in Figure 4.
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> plt.savefig("all.ps")
Examining Fit Results
The χ^{2} goodness-of-fit is reported with the best-fit values after a fit and 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.1920928955078125e-07 maxfev = None initsimplex = 0 finalsimplex = 9 step = None iquad = 1 verbose = 0 reflect = True Statistic: Chi2DataVar Chi Squared with data variance. The variance in each bin is estimated from the data value in that bin. If the number of counts in each bin is large, then the shape of the Poisson distribution from which the counts are sampled tends asymptotically towards that of a Gaussian distribution, with variance sigma(i)^2 = N(i,S) + [A(S)/A(B)]^2 N(i,B) where N is the number of on-source (and off-source) bins included in the fit. The background term appears only if an estimate of the background has been subtracted from the data. A(B) is the off-source "area", which could be the size of the region from which the background is extracted, or the length of a background time segment, or a product of the two, etc.; and A(S) is the on-source "area". These terms may be defined for a particular type of data: for example, PHA data sets A(B) to `BACKSCAL * EXPOSURE` from the background data set and A(S) to `BACKSCAL * EXPOSURE` from the source data set. See Also -------- Chi2Gehrels, Chi2ModVar, Chi2XspecVar Fit:Datasets = 1, 2, 3, 4 Method = neldermead Statistic = cstat Initial fit statistic = 273796 Final fit statistic = 4442.47 at function evaluation 1566 Data points = 3792 Degrees of freedom = 3786 Probability [Q-value] = 4.18275e-13 Reduced statistic = 1.17339 Change in statistic = 269354 bpow1.PhoIndx1 2.14733 bpow1.BreakE 1.11393 bpow1.PhoIndx2 1.56252 bpow1.norm 0.0247635 pow1d.PhoIndex 1.66745 pow1d.norm 0.000467698
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> print(calc_chisqr()) [2.92161203e-01 2.64440162e-01 1.79178565e+00 8.42223735e-02 3.62919502e-01 4.81765977e-01 2.60517362e-03 1.18987414e+00 ... 7.84338821e-01 1.73202191e+00 2.79090040e-02 4.75389283e-01 8.01371384e-01 4.97691406e-02 1.12280469e-01 3.12539503e-01]
The confidence (conf), covariance (covar) and projection (proj) commands can be used to estimate confidence intervals for the thawed parameters:
sherpa> set_stat("cstat") sherpa> conf() bpow1.PhoIndx2 lower bound: -0.0218748 pow1d.PhoIndex lower bound: -0.0806091 bpow1.PhoIndx2 upper bound: 0.0120026 bpow1.PhoIndx1 lower bound: -0.51134 bpow1.norm -: WARNING: The confidence level lies within (2.321579e-02, 2.398965e-02) bpow1.norm lower bound: -0.00116079 bpow1.PhoIndx1 +: WARNING: The confidence level lies within (2.409754e+00, 2.413448e+00) bpow1.PhoIndx1 upper bound: 0.264274 bpow1.BreakE lower bound: -0.0274067 pow1d.PhoIndex upper bound: 0.081217 pow1d.norm lower bound: -2.81036e-05 pow1d.norm upper bound: 2.95137e-05 bpow1.BreakE upper bound: 0.317126 bpow1.norm +: WARNING: The confidence level lies within (2.556533e-02, 2.556458e-02) bpow1.norm upper bound: 0.000801444 Datasets = 1, 2, 3, 4 Confidence Method = confidence Iterative Fit Method = None Fitting Method = neldermead Statistic = cstat confidence 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- bpow1.PhoIndx1 2.14733 -0.51134 0.264274 bpow1.BreakE 1.11393 -0.0274067 0.317126 bpow1.PhoIndx2 1.56252 -0.0218748 0.0120026 bpow1.norm 0.0247635 -0.00116079 0.000801444 pow1d.PhoIndex 1.66745 -0.0806091 0.081217 pow1d.norm 0.000467698 -2.81036e-05 2.95137e-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") sherpa> save("459_fitting_session.ascii")
The save function records all the information about the current session to the binary file 459_fitting_session.save, and the save_all function records the session settings to an editable ASCII file.
To restore the session that was saved to the binary file 459_fitting_session.save or ASCII file 459_fitting_session.ascii:
sherpa> restore("session1.save")
sherpa> %run -i session1.ascii
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
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 |
08 Jan 2010 | updated for CIAO 4.2 |
13 Jul 2010 | updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread. |
19 Aug 2011 | updated the Reading the Spectrum Files section with information on the two background up/down data sets included in Chandra grating PHA files |
20 Jan 2012 | reviewed for CIAO 4.4 (no changes) |
10 Dec 2013 | reviewed for CIAO 4.6: no changes |
25 Feb 2015 | updated for CIAO 4.7, no content changes |
14 Dec 2015 | reviewed for CIAO 4.8, no content changes |
15 Nov 2016 | reviewed for CIAO 4.9, no content changes, updated fit results. |
29 May 2018 | reviewed for CIAO 4.10, no content changes |
11 Dec 2018 | reviewed for CIAO 4.11, no content changes |
13 Dec 2019 | reviewed for CIAO 4.12, ChIPS figures and commands replaced with Matplotlib equivalent. Use xsbknpower for the broken-power law model instead of bpl1d and xspowerlaw instead of powlaw1d, since it returns better results if the model incorporates an RMF. |
21 Dec 2020 | Updated for CIAO 4.13: re-generated the plots as the display of PHA plots has changed, and changed the noticed data range to 2-12.5 Angstroms (it was 1-15). |
29 Mar 2029 | reviewed for CIAO 4.14, updated fit results. |