Fitting PHA Data with Multi-Component Source Models
OverviewLast Update: 1 Dec 2006 - reviewed for CIAO 3.4: no changes Synopsis: This thread shows how to do a basic spectral fit with the appropriate response files. If you are interested in including the background in the fit as well, see the Independent Background Responses thread. |
Contents
- Getting Started
- Reading Data & Instrument Responses Into Sherpa
- Filtering Spectral Data
- Establishing Model Components
- Defining a Multi-Component Source Model Expression
- Modifying Method & Statistic Settings
- Fitting
- Examining Fit Results
- Checking Sherpa Session Status
- Saving a Sherpa Session
- Using Sherpa Scripts & Restoring a Session
- Summary
- History
- Images
Getting Started
Please follow the "Sherpa Threads: Getting Started" thread.
Reading Data & Instrument Responses Into Sherpa
In this thread, we wish to fit spectral data from the FITS PHA datafile source_pi.fits. This dataset is input into Sherpa with the READ command:
sherpa> DATA source_pi.fits The inferred file type is PHA. If this is not what you want, please specify the type explicitly in the data command. WARNING: statistical errors specified in the PHA file. These are currently IGNORED. To use them, type: READ ERRORS "<filename>[cols CHANNEL,STAT_ERR]" fitsbin
Note that the READ DATA command may be shortened to just DATA.
Now the dataset may be plotted:
sherpa> LPLOT DATA
Note that these data are plotted in bin (i.e. channel) space, because there is no information about RMF/ARF attached to the header of this file. A standard PHA file contains number of counts in each PHA channel, which Sherpa reads. Therefore only counts per channel will be plotted with the above command. For details on Sherpa plotting see Data Visualization thread.
In order to assign the instrument properties (such as for example RMF and ARF files) an instrument model has to be defined. The Sherpa model name for an instrument response model, defined by the standard instrument response files, is RSP. The RSP model component is established for use and is named acis:
sherpa> RSP[acis] acis.rmf parameter value ["none"] rmf.fits acis.arf parameter value ["none"] arf.fits The inferred file type is ARF. If this is not what you want, please specify the type explicitly in the data command.
Notice that in the square bracket we named the response model as acis. Sherpa prompts the user for input of the model's default parameter values, e.g. two filenames for rmf and arf. Because the parameters belong to the acis model they are named acis.rmf and acis.arf. We set the rmf and arf parameter values to the proper RMF and ARF filenames for this dataset. It is possible, however, to establish a response model using only an RMF or an ARF file.
In order to assign acis model to the data and convolve the PHA data with the response, the instrument model has to be defined. Note that this allows for multiple models to be established during the Sherpa session, while only the one assigned to instrument used during fitting.
sherpa> INSTRUMENT = acis
The current definition of Sherpa's instrument model may be examined using SHOW INSTRUMENT:
sherpa> SHOW INSTRUMENT Instrument 1: rsp1d[acis] Param Type Value ----- ---- ----- 1 rmf string: "rmf.fits" (N_E=198,N_PHA=1024) 2 arf string: "arf.fits" (N_E=495)
This output shows that rmf.fits and arf.fits currently define the instrument model. Check RSP for requirement on the binning in RMF and ARF files.
The input dataset may be plotted again:
sherpa> LPLOT DATA
The data are now plotted in energy space since the instrument model is providing the information necessary for the computation of the number of predicted counts in each energy bin.
Filtering Spectral Data
Sherpa has many filtering options with IGNORE and NOTICE commands. During this Sherpa session, we would like to fit only the data between 0.5 and 8.0 keV. To apply this filter to the dataset:
sherpa> IGNORE ENERGY :0.5, 8:
Note: this command may be entered only after an instrument model has been defined.
To remove this filter:
sherpa> NOTICE ALL
The filter may then be re-applied (or a different filter may be applied):
sherpa> IGNORE ENERGY :0.5, 8:
The filtered dataset may then be plotted:
sherpa> LPLOT DATA
Figure 1 shows the resulting plot. Notice that the plot now includes only the data in the specified energy region.
Establishing Model Components
We will fit the spectral data using a source model expression that involves multiple model components.
The first of these model components is the one-dimensional power-law model, POWLAW1D. The POWLAW1D model component is established and is named p1 defined in the square bracket:
sherpa> PARAMPROMPT OFF Model parameter prompting is off sherpa> POWLAW1D[p1]
Since a dataset has been input, Sherpa estimates the initial parameter values (and the minimum and maximum for their ranges) for this model based on the data. The command PARAMPROMPT OFF cancels prompting for immediate changes to these model parameter value estimates.
The second of the model components is the XSpec photoelectric absorption model, called XSWABS in Sherpa. Please see the XSpec User's Guide, or ahelp xs, or ahelp xswabs for more information about this model. Note that all XSpec models are available from within Sherpa and their XSpec names must be preceded by XS.
The XSWABS model component is established and is named abs1:
sherpa> XSWABS[abs1]
Again, since a dataset has been previously input, Sherpa estimates the initial parameter values based on the data.
Defining a Multi-Component Source Model Expression
We will now fit the spectral data using an expression that is the product of the model components established above. The SOURCE defines the final model expression used for fitting the data:
sherpa> SOURCE = abs1*p1
The current definition of Sherpa's source model expression including initial parameter values for each model component may be examined using SHOW SOURCE:
sherpa> SHOW SOURCE Source 1: (abs1 * p1) xswabs[abs1] (XSPEC model name: absw) (integrate: off) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 nH thawed 1e-01 1e-07 10 10^22/cm^2 powlaw[p1] (integrate: on) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 gamma thawed 1 -10 10 2 ref frozen 1 0.5183 7.9789 3 ampl thawed 7.1163e-06 7.1163e-08 7.1163e-04
This output shows that abs * p1 is currently defined as the source model expression. By default, integration over an energy bin is turned off for the abs1 model component and is turned on for the p1 model component. However, since we are using a multiplicative source model expression to fit binned PHA data, integration of source expression over each energy bin will be performed during fitting.
Modifying Method & Statistic Settings
The SHOW command may be used to view the current method and statistics settings:
sherpa> SHOW METHOD Optimization Method: Levenberg-Marquardt Name Value Min Max Description ---- ----- --- --- ----------- 1 iters 2000 1 10000 Maximum number of iterations 2 eps 1e-03 1e-09 1 Absolute accuracy 3 smplx 0 0 1 Refine fit with simplex (0=no) 4 smplxep 1 1e-04 1000 Switch-to-simplex eps factor 5 smplxit 3 1 20 Switch-to-simplex iters factor sherpa> SHOW STATISTIC Statistic: Chi-Squared Gehrels
We will change the statistic to Cash for fitting this data:
sherpa> STATISTIC CASH
Note that truncation is turned on when the Cash statistic is used. This setting prevents negative model-predicted data values from affecting the convergence process.
Further details about the Levenberg-Marquardt optimization method which is the default method in Sherpa are available by typing:
sherpa> ahelp lev-mar
Fitting
The dataset is then fit:
sherpa> FIT WARNING: the Levenberg-Marquardt optimization method works less robustly when the Cash or cstat statistic is used. Consider using the Powell or Simplex method instead. LVMQT: V2.0 LVMQT: initial statistic value = 15306.3 LVMQT: final statistic value = -6708.78 at iteration 31 abs1.nH 1e-07 10^22/cm^2 p1.gamma -0.373229 p1.ampl 0.000113554 WARNING: The value of abs1.nH is equal to the abs1.nH.min limit boundary. You may wish to consider changing min/max values and refitting.
Notice that the first warning recommends that we use a different optimization method. Before refitting the data, we reset the initial parameter values and switch to the Powell method:
And now run the fit again:
sherpa> fit powll: v1.2 powll: initial statistic value = 1.53063E+04 powll: converged to minimum = -7.62300E+03 at iteration = 12 powll: final statistic value = -7.62300E+03 abs1.nH 1.48605 10^22/cm^2 p1.gamma 0.815361 p1.ampl 0.000711629 WARNING: The value of p1.ampl within 0.01% of the p1.ampl.max limit boundary. You may wish to consider changing min/max values and refitting.
The new warning refers to the final value of the amplitude parameter being to close to the boundary. We will expand allow min/max range for this parameter and refit. First, we examine the current parameter values for the model p1:
sherpa> SHOW p1 powlaw[p1] (integrate: on) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 gamma thawed 0.8154 -10 10 2 ref frozen 1 0.5183 7.9789 3 ampl thawed 7.1163e-04 7.1163e-08 7.1163e-04
As noted in the warning, we change the Max setting for the amplitude (p1.ampl):
sherpa> p1.ampl.max = 0.1
And fit the data a final time:
sherpa> FIT powll: v1.2 powll: initial statistic value = -7.62300E+03 powll: converged to minimum = -7.70549E+03 at iteration = 7 powll: final statistic value = -7.70549E+03 abs1.nH 2.37765 10^22/cm^2 p1.gamma 1.50086 p1.ampl 0.00200489
Use LPLOT to plot the fit and residuals. Also see Customizing Sherpa plots thread to learn how to use Slang to make nice plots.
sherpa> LPLOT 2 FIT RESIDUALS ==> Error bars computed using Chi Gehrels. ==> Error bars computed using Chi Gehrels. sherpa> PRINT POSTFILE sherpa.spectra.2.ps
Figure 2 shows the resulting plot. The features seen in the residuals could be due to the real source emission being different than the assumed power law model. This is in fact true for our data set. The spectrum is from a supernova remnant usually characterized by thermal emission accompanying with several emission lines. Also some features could be due to calibration uncertainties. We do not discuss further analysis here, which in normal data analysis process would consist of changing the source expression to include preferable plasma emission model and refiting.
Examining Fit Results
Several algorithms are available in Sherpa for examining fit results such as covariance, projection, uncertainty, interval-projection, region-projection A detailed description of how to use different options is available in several threads: Step-by-Step Guide to Estimating Errors and Confidence Levels and Estimating Errors and Confidence Levels. Also Sherpa S-lang module functions allows for easy access to the results in both scripts and on the command line. They are described in Accessing fit results using S-Lang thread.
Checking Sherpa Session Status
The final overall status of this Sherpa session may be viewed as follows:
sherpa> SHOW Optimization Method: Powell Statistic: Cash ----------------- Input data files: ----------------- Data 1: source_pi.fits pha. Total Size: 1024 bins (or pixels) Dimensions: 1 Total counts (or values): 3328 Exposure: 7854.47 sec Count rate: 0.424 cts/sec Backscal: 2.366406e-06 Current filters for dataset 1: notice source 1 all IGNORE source 1 ENERGY : 0.5 , 8 : Noticed filter size: 512 bins Sum of data within filter: 3283 ------------------------------ Defined analysis model stacks: ------------------------------ source 1 = (abs1 * p1) instrument source 1 = acis ------------------------------------------- Defined source/background model components: ------------------------------------------- powlaw[p1] (integrate: on) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 gamma thawed 1.5009 -10 10 2 ref frozen 1 0.5183 7.9789 3 ampl thawed 2.0049e-03 7.1163e-08 1e-01 xswabs[abs1] (XSPEC model name: absw) (integrate: off) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 nH thawed 2.3776 1e-07 10 10^22/cm^2 ------------------------------------ Defined instrument model components: ------------------------------------ rsp1d[acis] Param Type Value ----- ---- ----- 1 rmf string: "rmf.fits" (N_E=198,N_PHA=1024) 2 arf string: "arf.fits" (N_E=495)
Saving a Sherpa Session
To save a Sherpa session so that we may return again later, use the SAVE command:
where session1.shp is the filename for an ASCII file to which a record the current session will be written, in the form of a Sherpa script.
Using Sherpa Scripts & Restoring a Session
Sherpa commands may either be issued on the Sherpa command line or from a file by the USE command.
To restore the session that was saved to the file session1.shp:
sherpa> USE session1.shp The inferred file type is ARF. If this is not what you want, please specify the type explicitly in the data command. The inferred file type is PHA. If this is not what you want, please specify the type explicitly in the data command. WARNING: statistical errors specified in the PHA file. These are currently IGNORED. To use them, type: READ ERRORS "<filename>[cols CHANNEL,STAT_ERR]" fitsbin WARNING: model truncation turned on. Override with TRUNCATE OFF. sherpa>
One may verify that the session has been properly restored by comparing the SHOW output from the Checking Sherpa Session Status section.
Summary
This thread is complete, so we can exit the Sherpa session:
sherpa> EXIT Goodbye.
History
14 Jan 2005 | reviewed for CIAO 3.2: no changes |
21 Dec 2005 | reviewed for CIAO 3.3: no changes |
01 Dec 2006 | reviewed for CIAO 3.4: no changes |