Simultaneously Fitting Two Data Sets

Statistical Issues: Background Subtraction

A typical dataset may contain multiple spectra, one of which contains contributions from the source and the background, and one or more which contain background counts alone. (The background itself may contain contributions from the cosmic X-ray background, the particle background, an so on, but we ignore this complication.)

The proper way to treat background data is to model them. However, many X-ray astronomers subtract background data from the raw data.

Why should one not subtract background?

• It reduces the amount of statistical information in the analysis - the final fit parameter values will be a less accurate estimate of the true values.
• The background subtracted data are not Poisson-distributed; one cannot fit them with the Poisson likelihood or the Cash statistic, even in the low-count limit. For example, subtracting a background can give negative counts; this is definitely not Poissonian!
• Fluctuations, particularly in the vicinity of localized features, can adversely affect analysis.

Reading FITS Data

Simultaneous fitting is identical to the basic method of fitting data, except for the specification of the data set id. Most Sherpa commands have the option of specifying the id of the data set on which to invoke the command. By default, this number is assumed to be 1; the id is assigned when the data are loaded into Sherpa.

The spectral data from the two observations of the source are stored in the FITS files pi2278.fits and pi2286.fits. Data Set pi2278.fits is input with the load_pha() command and given id "1":

sherpa> load_pha(1, "pi2278.fits")
statistical errors were found in file 'pi2278.fits'
but not used; to use them, re-read with use_errors=True

The second data set is loaded similarly, but is given id "2":

sherpa> load_pha(2, "pi2286.fits")
statistical errors were found in file 'pi2286.fits'
but not used; to use them, re-read with use_errors=True

Defining Instrument Responses

Sherpa will automatically read in any instrument response files that are defined in the header of the spectrum (keywords "RESPFILE" and "ANCRFILE").

In this example, the response files were not defined automatically, so we load them with load_arf() and load_rmf(). The data set id of the corresponding PI file is given in each command:

sherpa> load_arf(1, "arf2278.fits")
sherpa> load_rmf(1, "rmf2278.fits")

sherpa> load_arf(2, "arf2286.fits")
sherpa> load_rmf(2, "rmf2286.fits")

Defining Source Models

We want to model each of the spectra with an absorbed power-law, so a power-law component (powlaw1d, called pl1) and an absorbing component (xswabs, called abs1) are defined. The same model expression is assigned to each data set:

sherpa> set_source(1, xswabs.abs1 * powlaw1d.pl1)
sherpa> set_source(2, abs1 * pl1)

We can use the Sherpa 4.1 guess command to guess the initial parameter values and ranges for the power-law model component (as model parameter values are not automatically guessed in Sherpa 4.1 as they were in previous versions of the software):

sherpa> guess(pl1)

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((11619.0814305 * (xswabs.abs1 * powlaw1d.pl1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed            1            0       100000 10^22 atoms / cm^2
   pl1.gamma    thawed            1          -10           10           
   pl1.ref      frozen            1 -3.40282e+38  3.40282e+38           
   pl1.ampl     thawed  6.42531e-06  6.42531e-08  0.000642531           

Model: 2
apply_rmf(apply_arf((11619.0814107 * (xswabs.abs1 * powlaw1d.pl1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed            1            0       100000 10^22 atoms / cm^2
   pl1.gamma    thawed            1          -10           10           
   pl1.ref      frozen            1 -3.40282e+38  3.40282e+38           
   pl1.ampl     thawed  6.42531e-06  6.42531e-08  0.000642531    

(The guess command makes an initial guess at parameter values to ensure convergence, but it is always a good idea to check that the initial range of values (soft limits) is sensible for the data being fit.)

Data above 5.5 keV is ignored in both data sets to cut out energies with no data. While this isn't strictly necessary for this data, it does produce a nicer plot. Ignoring the higher energies may be especially useful for ACIS-S data, where the high-energy background is stronger.

The comma after the value indicates an open range, i.e. everything above "5.5" will be ignored:

sherpa> ignore(5.5,)

Fit the Data

In order to fit both data sets simultaneously, fit() is called with both ids:

sherpa> fit(1,2)
 Solar Abundance Vector set to angr:  Anders E. & Grevesse N. Geochimica et Cosmochimica Acta 53, 197 (1989)
 Cross Section Table set to bcmc:  Balucinska-Church and McCammon, 1998
Data Set              = 1, 2
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 27.6983
Final fit statistic   = 5.73184 at function evaluation 29
Data points           = 14
Degrees of freedom    = 11
Probability [Q-value] = 0.890644
Reduced statistic     = 0.521076
Change in statistic   = 21.9665
   abs1.nh        1.2101   
   pl1.gamma      2.3981   
   pl1.ampl       4.32033e-05

Plot both fits, as shown in Figure 1, and create a postscript version:

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

sherpa> print_window("two_fits",["format","ps"])

Figure 1: Two spectra, fit with same model

Link Parameters and Fit Again

If the second spectrum is similar but not identical to the first spectrum, a different model component could be used in one of the source definitions. Additionally, parameter linking can be used to fix parameter values to each other.

For example, if the intensity of the source is expected to change and the power-law slope is not, we could define a new power-law component and link the slopes together. The result is another free parameter in the fit: the normalization of the second power-law component.

sherpa> set_source(2, abs1 * powlaw1d.pl2)
sherpa> guess(pl2)
sherpa> pl2.gamma = pl1.gamma

Now fit the data again:

sherpa> fit(1,2)
Data Set              = 1, 2
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 28.4437
Final fit statistic   = 5.7289 at function evaluation 31
Data points           = 14
Degrees of freedom    = 10
Probability [Q-value] = 0.837503
Reduced statistic     = 0.57289
Change in statistic   = 22.7148
   abs1.nh        1.20703     
   pl1.gamma      2.23639     
   pl1.ampl       4.27257e-05 
   pl2.ampl       4.33182e-05 

Only the amplitude of the new component (pl2.ampl) is fit. In this case, the amplitudes of the two power laws are nearly identical. The final model parameter values may be printed:

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((11619.0814305 * (xswabs.abs1 * powlaw1d.pl1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed      1.20703            0       100000 10^22 atoms / cm^2
   pl1.gamma    thawed      2.23639          -10           10           
   pl1.ref      frozen            1 -3.40282e+38  3.40282e+38           
   pl1.ampl     thawed  4.27257e-05  6.42531e-08  0.000642531           

Model: 2
apply_rmf(apply_arf((11619.0814107 * (xswabs.abs1 * powlaw1d.pl2))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed      1.20703            0       100000 10^22 atoms / cm^2
   pl2.gamma    linked      2.23639          expr: pl1.gamma           
   pl2.ref      frozen            1 -3.40282e+38  3.40282e+38           
   pl2.ampl     thawed  4.33182e-05  6.42531e-08  0.000642531   

This thread is complete, so we can exit the Sherpa 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 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.