Fitting PHA Data with Multi-Component Source Models

Getting Started

Please follow the "Obtaining data used in Sherpa threads" thread.

Reading Data & Instrument Responses Into Sherpa

In this thread, we wish to fit spectral data from the FITS PHA data file source_pi.fits. This data set is input to Sherpa with the load_pha function:

sherpa> load_pha("source_pi.fits")
read ARF file arf.fits
read RMF file rmf.fits

Note that it is no longer necessary to build instrument models as in Sherpa 3.4 - the instrument response is automatically established when the appropriate response files (ARF, RMF) are input to the Sherpa session. If the instrument response files are written in the header of the PHA file, as in this example, Sherpa will load them automatically; if not, they need to be loaded manually with the load_arf and load_rmf commands.

Now the data set may be plotted:

sherpa> plot_data()

The data are plotted in energy space, because information about the RMF and ARF are attached to the header of this file (the instrument response provides the information necessary for the computation of the number of predicted counts in each energy bin). If the response files are not automatically loaded with a PHA file, Sherpa will plot the number of counts in each PHA channel.

The current definition of the instrument response may be examined using show_all:

sherpa> show_all()
Data Set: 1
Filter: 0.0110-14.9431 Energy (keV)
Noticed Channels: 1-1024
name           = source_pi.fits
channel        = Int32[1024]
counts         = Int32[1024]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 56494.4943719
backscal       = 9.36267570731e-05
areascal       = 1.0
grouped        = False
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

RMF Data Set: 1:1
name     = rmf.fits
detchans = 1024
energ_lo = Float64[1070]
energ_hi = Float64[1070]
n_grp    = Int16[1070]
f_chan   = UInt32[1344]
n_chan   = UInt32[1344]
matrix   = Float64[382884]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

ARF Data Set: 1:1
name     = arf.fits
energ_lo = Float64[1070]
energ_hi = Float64[1070]
specresp = Float64[1070]
bin_lo   = None
bin_hi   = None
exposure = 56494.4943719

This output shows that rmf.fits and arf.fits currently define the instrument response.

Filtering Spectral Data

Sherpa has many filtering options with the ignore and notice functions. 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 data set, we use ignore:

sherpa> ignore(":0.5, 8.0:")

Note that the units of the arguments supplied to the ignore / notice functions must match the units of the data set. The data units can be changed with the get_data or set_analysis functions:

sherpa> set_analysis("energy")  [OR 'get_data().units=energy']

sherpa> print(get_data().units)
energy

To remove this filter:

sherpa> notice()

The filter may then be re-applied (or a different filter may be applied):

sherpa> ignore(":0.5, 8.0:")

The filtered data set may then be plotted:

sherpa> plot_data()

Figure 1 shows the resulting plot. Notice that the plot now includes only the data in the specified energy region.

Figure 1: Source spectrum (0.5-8.0 keV)

Defining a Multi-Component Source Model Expression

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 second of the model components is the XSpec photoelectric absorption model, called XSWABS in Sherpa. Please see the XSpec User's Guide 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.

It is no longer necessary to establish each model component separately before defining a multi-component source model expression, as in Sherpa 3.4; both steps may now be done simultaneously with the set_model function. Below, the POWLAW1D and XSWABS model components are established as p1 and abs1, respectively, and the multi-component source model expression is defined:

sherpa> set_model(xswabs.abs1*powlaw1d.p1)

We define an expression that is the product of the two model components. The function set_model defines the final model expression used for fitting the data.

The guess function may be used to estimate the initial parameter values for this model, as well as the minima and maxima for their ranges; where guess() is not able to estimate initial parameter values, we set our own:

sherpa> guess(p1)
sherpa> guess(abs1)
WARNING: No guess found for xswabs.abs1
sherpa> abs1.nH = 1e-07

The current definition of the Sherpa source model expression, including initial parameter values for each model component, may be examined using show_model:

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((56494.4943719 * (xswabs.abs1 * powlaw1d.p1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed        1e-07            0        1e+05 10^22 atoms / cm^2
   p1.gamma     thawed            1          -10           10           
   p1.ref       frozen            1     -3.4e+38      3.4e+38           
   p1.ampl      thawed  1.02908e-05        1e-07        0.001     


sherpa> print(abs1.integrate)
True

sherpa> print(p1.integrate)
True          

This output shows that abs * p1 is currently defined as the source model expression. By default, integration over an energy bin is turned on for both the abs1 and p1 model components; a change to the integration flag will not change the fit. Since we are using a multiplicative source model expression to fit binned PHA data, integration of the source expression over each energy bin will be performed during fitting.

Modifying Method & Statistic Settings

The show_method and show_stat functions may be used to view 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

We will change the statistic to Chi2XSpecVar for fitting this data:

sherpa> set_stat("chi2xspecvar")

Chi2XspecVar is the XSpec version of the χ² fit statistic with data variance.

Fitting

The data set is then fit:

sherpa> fit()
Datasets: (1,)
LevMar:
Initial fit statistic = 1784.24
Final fit statistic   = 403.989 at function evaluation 21
Data points           = 512
Degrees of freedom    = 509
Probability [Q-value] = 0.999791
Reduced statistic     = 0.793692
Change in statistic   = 1380.25
   abs1.nh        1e-07       
   p1.gamma       1.64249     
   p1.ampl        2.01787e-05 

Next, we decide to try a different optimization method. Before refitting the data, we reset the initial parameter values and switch to the Nelder-Mead / Simplex method:

sherpa> reset(get_model())
sherpa> set_method("simplex")

Now we run the fit again:

sherpa> fit()
Datasets: (1,)
NelderMead:
Initial fit statistic = 1784.24
Final fit statistic   = 390.678 at function evaluation 236
Data points           = 512
Degrees of freedom    = 509
Probability [Q-value] = 0.999971
Reduced statistic     = 0.76754
Change in statistic   = 1393.56
   abs1.nh        0.10621     
   p1.gamma       2.16445     
   p1.ampl        2.94145e-05 

Note: If the output returned by the fit function includes a warning that a final parameter value is too close to a boundary, the min/max range for the parameter may be expanded:

sherpa> model.parameter.max = desired upper bound
sherpa> model.parameter.min = desired lower bound

For example:
sherpa> p1.ampl.max = 0.003; 

See the section Examining Fit Results for a list of Sherpa functions which provide detailed information on the quality of parameter values used in a fit.

Use plot_fit_resid to plot the fit and residuals, and the Chips function print_window to save the plot as a PostScript file:

sherpa> plot_fit_resid()

sherpa> print_window("sherpa.spectra.2")

Figure 2 shows the resulting plot.

Figure 2: Power-law + absorption fit with residuals

The features seen in the residuals could be due to the real source emission being different than the assumed power law model, or a result of calibration uncertainties. We do not discuss further scientific analysis here, which would involve changing the source expression to include a preferable plasma emission model and refitting.

Examining Fit Results

Several algorithms are available in Sherpa for examining fit results, such as covariance (covar()), projection (proj()), interval-projection (int_proj()), and region-projection (reg_proj()). Figure 3 displays a contour plot of confidence regions produced by the function reg_proj:

sherpa> reg_proj(abs1.nH, p1.gamma)

Figure 3: Confidence region as a contour plot

In CIAO 3.4, the goodness command was used to get the chi-squared goodness-of-fit. This information is now reported with the best-fit values after a fit; it is no longer necessary to run a separate command. The function get_fit_results allows access to this information after the fit has been performed.

Checking Sherpa Session Status

The final overall status of this Sherpa session may be viewed as follows:

sherpa> show_all()
Data Set: 1
Filter: 0.5183-7.9789 Energy (keV)
Noticed Channels: 36-547
name           = source_pi.fits
channel        = Int32[1024]
counts         = Int32[1024]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 56494.4943719
backscal       = 9.36267570731e-05
areascal       = 1.0
grouped        = False
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

RMF Data Set: 1:1
name     = rmf.fits
detchans = 1024
energ_lo = Float64[1070]
energ_hi = Float64[1070]
n_grp    = Int16[1070]
f_chan   = UInt32[1344]
n_chan   = UInt32[1344]
matrix   = Float64[382884]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

ARF Data Set: 1:1
name     = arf.fits
energ_lo = Float64[1070]
energ_hi = Float64[1070]
specresp = Float64[1070]
bin_lo   = None
bin_hi   = None
exposure = 56494.4943719

Model: 1
apply_rmf(apply_arf((56494.4943719 * (xswabs.abs1 * powlaw1d.p1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed      0.10621            0        1e+05 10^22 atoms / cm^2
   p1.gamma     thawed      2.16445          -10           10           
   p1.ref       frozen            1     -3.4e+38      3.4e+38           
   p1.ampl      thawed  2.94145e-05      3.1e-07       0.0031           

Optimization Method: NelderMead
name         = simplex
ftol         = 1.19209289551e-07
maxfev       = None
scale        = 3.0
initsimplex  = 0
finalsimplex = 1
step         = None
seed         = 7151
iquad        = 1
verbose      = 0

Statistic: Chi2XspecVar

Fit:Datasets: (1,)
NelderMead:
Initial fit statistic = 1784.24
Final fit statistic   = 390.678 at function evaluation 236
Data points           = 512
Degrees of freedom    = 509
Probability [Q-value] = 0.999971
Reduced statistic     = 0.76754
Change in statistic   = 1393.56
   abs1.nh        0.10621     
   p1.gamma       2.16445     
   p1.ampl        2.94145e-05 

Saving a Sherpa Session

Before exiting Sherpa, you may wish to save the session in order to return to the analysis at a later point:

sherpa> save("session1.save")

All the information about the current session is written to session1.save, a binary file.

To restore the session that was saved to the file session1.save:

sherpa> restore("session1.save")

One may verify that the session has been properly restored by comparing the show_all() output from the Checking Sherpa Session Status section.

Scripting It

The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing execfile("fit.py") on the Sherpa command line.

The Sherpa script command may be used to save everything typed on the command line in a Sherpa session:

sherpa> script(filename="sherpa.log", clobber=False)

The CXC is committed to helping Sherpa users transition to new syntax as smoothly as possible. If you have existing Sherpa scripts or save files, submit them to us via the CXC Helpdesk and we will provide the CIAO/Sherpa 4.1 syntax to you.