Simultaneously Fitting Two Datasets

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 dataset id. Most Sherpa commands have the option of specifying the id of the dataset 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 stored in the FITS files pi2278.fits and pi2286.fits. Dataset 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 dataset 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 dataset 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 powerlaw, so a powerlaw component (powlaw1d, called pl1) and an absorbing component (xswabs, called abs1) are defined. The same model expression is assigned to each dataset:

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

Data above 5.5 keV is ignored in both datasets 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 datasets 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
Statistic value = 5.73184 at function evaluation 29
Data points = 14
Degrees of freedom = 11
Probability [Q-value] = 0.890644
Reduced statistic  = 0.521076
   abs1.nh        1.21008
   pl1.gamma      2.23979
   pl1.ampl       4.32024e-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"]);

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 powerlaw slope isn't, we could define a new powerlaw component and link the slopes together. The result is another free parameter in the fit: the normalization of the second powerlaw component.

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

Now fit the data again:

sherpa> fit(1,2);
Statistic value = 5.7289 at function evaluation 31
Data points = 14
Degrees of freedom = 10
Probability [Q-value] = 0.837503
Reduced statistic  = 0.57289
   abs1.nh        1.20707
   pl1.gamma      2.23648
   pl1.ampl       4.273e-05
   pl2.ampl       4.33222e-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> print(pl1*abs1);
(pl1 * abs1)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   pl1.gamma    thawed      2.23648          -10           10
   pl1.ref      frozen            1      -1e+120       1e+120
   pl1.ampl     thawed    4.273e-05            0       1e+120
   abs1.nh      thawed      1.20707            0       100000      10^22

sherpa> print(pl2*abs1);
(pl2 * abs1)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   pl2.gamma    linked      2.23648          -10           10
   pl2.ref      frozen            1      -1e+120       1e+120
   pl2.ampl     thawed  4.33222e-05            0       1e+120
   abs1.nh      thawed      1.20707            0       100000      10^22

This thread is complete, so we can exit the Sherpa session:

sherpa> quit();