next up previous contents
Next: Using XSPEC to Simulate Up: Walks through XSPEC Previous: Fitting Models to Data:

Simultaneous Fitting: Examples from Einstein and Ginga

XSPEC has the very useful facility of allowing models to be fitted simultaneously to more than one data file. It is even possible to group files together and to fit different models simultaneously. Reasons for fitting in this manner include:

Other scenarios are possible--the important thing is to recognize the flexibility of XSPEC in this regard.

As an example of the first case, we'll fit two spectra derived from two separate Einstein Solid State Spectrometer (SSS) observations of the cooling-flow cluster Abell 496. Although the two observations were carried out on consecutive days (in August 1979), the response is different, due to the variable build-up of ice on the detector. This problem bedeviled analysis during the mission; however, it has now been calibrated successfully and is incorporated into the response matrices associated with the spectral files in the HEASARC archive. The SSS also provides an example of how background correction files are used in XSPEC.

To fit the same model with the same set of parameters to more than one data file, simply enter the names of the data files after the data command:

XSPEC> data sa496b.pha sa496c.pha
 Net count rate (cts/s) for file   1  0.7806    +/-  9.3808E+05( 86.9% total)
   using response (RMF) file...       sa496b.rsp
   using background file...           sa496b.bck
   using correction file...           sa496b.cor
 Net count rate (cts/s) for file   2  0.8002    +/-  9.3808E+05( 86.7% total)
   using response (RMF) file...       sa496c.rsp
   using background file...           sa496c.bck
   using correction file...           sa496c.cor
 Net correction flux for file  1=    8.4469E-04
 Net correction flux for file  2=    8.7577E-04
   2 data sets are in use

As the messages indicate, XSPEC also has read in the associated:

These files are all listed in the headers of the data files (sa496b.pha & sa496c.pha).

To ignore channels, the file number (1 & 2 in this example) precedes the range of channels to be ignored. Here, we wish to ignore, for both files, channels 1-15 and channels 100-128. This can be done by specifying the files one after the other with the range of ignored channels:

XSPEC> ignore 1:1-15 1:100-128 2:1-15 2:100-128

or by specifying the range of file number with the channel range:

XSPEC> ignore 1-2:1-15 100-128

In this example, we'll fit a cooling-flow model under the reasonable assumption that the small SSS field of view sees mostly just the cool gas in the middle of the cluster. We'll freeze the values of the maximum temperature (the temperature from which the gas cools) and of the abundance to the values found by instruments such as the Ginga LAC and the EXOSAT ME, which observed the entire cluster. The minimum gas temperature is frozen at 0.1 keV; the ``slope'' is frozen at zero (isobaric cooling) and the normalization is given an initial value of 100 solar masses per year:

XSPEC> mo pha(cflow)
  mo = phabs[1]( cflow[2] )
 Input parameter value, delta, min, bot, top, and max values for ...
 Mod parameter  1 of component  1 phabs    nH       10^22
    1.000      1.0000E-03  0.0000E+00  0.0000E+00  1.0000E+05  1.0000E+06
0.045
 Mod parameter  2 of component  2 cflow    slope
   0.0000E+00  1.0000E-02  -5.000      -5.000       5.000       5.000
0,0
 Mod parameter  3 of component  2 cflow    lowT     keV
   0.1000      1.0000E-03  8.0800E-02  8.0800E-02   79.90       79.90
0.1,0
 Mod parameter  4 of component  2 cflow    highT    keV
    4.000      1.0000E-03  8.0800E-02  8.0800E-02   79.90       79.90
4,0
 Mod parameter  5 of component  2 cflow    Abundanc
    1.000      1.0000E-02  0.0000E+00  0.0000E+00   5.000       5.000
0.5,0
 Mod parameter  6 of component  2 cflow    Redshift
   0.0000E+00 -0.1000      0.0000E+00  0.0000E+00   100.0       100.0
0.032
 Mod parameter  7 of component  2 cflow    norm
    1.000      1.0000E-02  0.0000E+00  0.0000E+00  1.0000E+24  1.0000E+24
100
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  mo = phabs[1]( cflow[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     4.5000E-02 +/-  0.0000E+00
    2    2    2   cflow      slope              0.0000E+00 frozen
    3    3    2   cflow      lowT     keV       0.1000     frozen
    4    4    2   cflow      highT    keV        4.000     frozen
    5    5    2   cflow      Abundanc           0.5000     frozen
    6    6    2   cflow      Redshift           3.2000E-02 frozen
    7    7    2   cflow      norm                100.0     +/-  0.0000E+00
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
    2 variable fit parameters
 Chi-Squared =      2705.736     using   168 PHA bins.
 Reduced chi-squared =      16.29961     for    166 degrees of freedom
 Null hypothesis probability = 0.000E+00
XSPEC> fit
 Chi-Squared    Lvl  Fit param # 1     2           3           4
                 5           6           7
   401.110     -3     0.1889      0.0000E+00  0.1000       4.000
              0.5000      3.2000E-02   285.1
   364.046     -4     0.2298      0.0000E+00  0.1000       4.000
              0.5000      3.2000E-02   316.1
   363.536     -5     0.2351      0.0000E+00  0.1000       4.000
              0.5000      3.2000E-02   319.8
   363.530     -6     0.2358      0.0000E+00  0.1000       4.000
              0.5000      3.2000E-02   320.2
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               1     7
 3.04E-05 | -1.00  0.00
 3.38E+01 |  0.00 -1.00
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  mo = phabs[1]( cflow[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.2358     +/-  0.8515E-02
    2    2    2   cflow      slope              0.0000E+00 frozen
    3    3    2   cflow      lowT     keV       0.1000     frozen
    4    4    2   cflow      highT    keV        4.000     frozen
    5    5    2   cflow      Abundanc           0.5000     frozen
    6    6    2   cflow      Redshift           3.2000E-02 frozen
    7    7    2   cflow      norm                320.2     +/-   5.816
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      363.5299     using   168 PHA bins.
 Reduced chi-squared =      2.189939     for    166 degrees of freedom
 Null hypothesis probability = 8.325E-17

As we can see, tex2html_wrap_inline10235 is not good, but the high statistic could be because we have yet to adjust the correction file. Correction files in XSPEC take into account detector features that cannot be completely prescribed ab initio and which must be fitted at the same time as the model. Einstein SSS spectra, for example, have a background feature the level of which varies unpredictably. Its spectral form is contained in the correction file, but its normalization is determined by fitting. This fitting is set in motion using the command recornrm ( reset the correction-file normalization):

XSPEC> reco 1
  File #     Correction
      1      0.4022 +/-  0.0673
 After correction norm adjustment  0.402 +/- 0.067
 Chi-Squared =      327.7646     using   168 PHA bins.
 Reduced chi-squared =      1.974485     for    166 degrees of freedom
 Null hypothesis probability = 1.091E-12
XSPEC> reco 2
  File #     Correction
      2      0.4775 +/-  0.0597
 After correction norm adjustment  0.477 +/- 0.060
 Chi-Squared =      263.8333     using   168 PHA bins.
 Reduced chi-squared =      1.589336     for    166 degrees of freedom
 Null hypothesis probability = 2.025E-06

This process is iterative, and, in order to work, must be used in tandem with fitting the model. Successive fits and recorrections are applied until the fit is stable, i.e., until further improvement in tex2html_wrap_inline10235 no longer results. Of course, this procedure is only worthwhile when the model gives a reasonably good account of the data. Eventually, we end up at:

XSPEC> fit
 Chi-Squared    Lvl  Fit param # 1     2           3           4
                 5           6           7
   189.931     -3     0.2791      0.0000E+00  0.1000       4.000
              0.5000      3.2000E-02   300.4
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               1     7
 4.85E-05 | -1.00  0.00
 3.84E+01 |  0.00 -1.00
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  mo = phabs[1]( cflow[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.2791     +/-  0.1072E-01
    2    2    2   cflow      slope              0.0000E+00 frozen
    3    3    2   cflow      lowT     keV       0.1000     frozen
    4    4    2   cflow      highT    keV        4.000     frozen
    5    5    2   cflow      Abundanc           0.5000     frozen
    6    6    2   cflow      Redshift           3.2000E-02 frozen
    7    7    2   cflow      norm                300.4     +/-   6.198
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      189.9309     using   168 PHA bins.
 Reduced chi-squared =      1.144162     for    166 degrees of freedom
 Null hypothesis probability = 9.831E-02

The final value of tex2html_wrap_inline10235 is much better than the original, but is not quite acceptable. However, the current model has only two free parameters: further explorations of parameter space would undoubtedly improve the fit.

We'll leave this example and move on to look at another kind of simultaneous fitting: one where the same model is fitted to two different data files. This time, not all the parameters will be identical. The massive X-ray binary Centaurus X-3 was observed with the LAC on Ginga in 1989. Its flux level before eclipse was much lower than the level after eclipse. Here, we'll use XSPEC to see whether spectra from these two phases can be fitted with the same model, which differs only in the amount of absorption. This kind of fitting relies on introducing an extra dimension, the group, to the indexing of the data files. The files in each group share the same model but not necessarily the same parameter values, which may be shared as common to all the groups or varied separately from group to group. Although each group may contain more than one file, there is only one file in each of the two groups in this example. Groups are specified with the data command, with the group number preceding the file number, like this:

XSPEC> da 1:1 losum 2:2 hisum
 Net count rate (cts/s) for file   1   140.1    +/-  0.3549
   using response (RMF) file...       ginga_lac.rsp
 Net count rate (cts/s) for file   2   1371.    +/-   3.123
   using response (RMF) file...       ginga_lac.rsp
   2 data sets are in use

Here, the first group makes up the file losum.pha, which contains the spectrum of all the low, pre-eclipse emission. The second group makes up the second file, hisum.pha, which contains all the high, post-eclipse emission. Note that file number is ``absolute'' in the sense that it is independent of group number. Thus, if there were three files in each of the two groups (lo1.pha, lo2.pha, lo3.pha, hi1.pha, hi2.pha & hi3.pha, say), rather than one, the six files would be specified as follows:

XSPEC> da 1:1 lo1 1:2 lo2 1:3 lo3 2:4 hi1 2:5 hi2 2:6 hi3

The ignore command works, as usual, on file number, and does not take group number into account. So, to ignore channels 1-3 and 37-48 of both files:

XSPEC> ignore 1-2:1-3 37-48

The model we'll use at first to fit the two files is an absorbed power law with a high-energy cut-off:

XSPEC> mo phabs * highecut (po)

After defining the model, the user is prompted for two sets of parameter values, one for the first group of data files (losum.pha), the other for the second group (hisum.pha). Here, we'll enter the absorption column of the first group as tex2html_wrap_inline10355 cm tex2html_wrap_inline12821 and enter the default values for all the other parameters in the first group. Now, when it comes to the second group of parameters, we enter a column of tex2html_wrap_inline12819 cm tex2html_wrap_inline12821 and then enter defaults for the other parameters. The rule being applied here is as follows: to tie parameters in the second group to their equivalents in the first group, take the default when entering the second-group parameters; to allow parameters in the second group to vary independently of their equivalents in the first group, enter different values explicitly. Like this:

XSPEC> mo phabs * highecut (po)
  mo = phabs[1]*highecut[2]( powerlaw[3] )
 Input parameter value, delta, min, bot, top, and max values for ...
 Mod parameter  1 of component  1 phabs    nH       10^22     for data group 1
    1.000      1.0000E-03  0.0000E+00  0.0000E+00  1.0000E+05  1.0000E+06
100 
 Mod parameter  2 of component  2 highecut cutoffE  keV       for data group 1
    10.00      1.0000E-02  1.0000E-04  1.0000E-02   100.0       200.0

 Mod parameter  3 of component  2 highecut foldE    keV       for data group 1
    15.00      1.0000E-02  1.0000E-04  1.0000E-02   100.0       200.0

 Mod parameter  4 of component  3 powerlaw PhoIndex           for data group 1
    1.000      1.0000E-02  -3.000      -2.000       9.000       10.00

 Mod parameter  5 of component  3 powerlaw norm               for data group 1
    1.000      1.0000E-02  0.0000E+00  0.0000E+00  1.0000E+24  1.0000E+24

 Mod parameter  6 of component  4 phabs    nH       10^22     for data group 2
    100.0      1.0000E-03  0.0000E+00  0.0000E+00  1.0000E+05  1.0000E+06
1
 Mod parameter  7 of component  5 highecut cutoffE  keV       for data group 2
    10.00      1.0000E-02  1.0000E-04  1.0000E-02   100.0       200.0

 Mod parameter  8 of component  5 highecut foldE    keV       for data group 2
    15.00      1.0000E-02  1.0000E-04  1.0000E-02   100.0       200.0

 Mod parameter  9 of component  6 powerlaw PhoIndex           for data group 2
    1.000      1.0000E-02  -3.000      -2.000       9.000       10.00

 Mod parameter 10 of component  6 powerlaw norm               for data group 2
    1.000      1.0000E-02  0.0000E+00  0.0000E+00  1.0000E+24  1.0000E+24

  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  mo = phabs[1]*highecut[2]( powerlaw[3] )
  Model Fit Model Component  Parameter  Unit     Value                    Data
  par   par comp                                                          group
    1    1    1   phabs      nH       10^22      100.0     +/-  0.0000E+00   1
    2    2    2   highecut   cutoffE  keV        10.00     +/-  0.0000E+00   1
    3    3    2   highecut   foldE    keV        15.00     +/-  0.0000E+00   1
    4    4    3   powerlaw   PhoIndex            1.000     +/-  0.0000E+00   1
    5    5    3   powerlaw   norm                1.000     +/-  0.0000E+00   1
    6    6    4   phabs      nH       10^22      1.000     +/-  0.0000E+00   2
    7    2    5   highecut   cutoffE  keV        10.00     = par   2         2
    8    3    5   highecut   foldE    keV        15.00     = par   3         2
    9    4    6   powerlaw   PhoIndex            1.000     = par   4         2
   10    5    6   powerlaw   norm                1.000     = par   5         2
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
    6 variable fit parameters
 Chi-Squared =     1.9666626E+07 using    66 PHA bins.
 Reduced chi-squared =      327777.1     for     60 degrees of freedom
 Null hypothesis probability = 0.000E+00

Notice how the summary of the model, displayed immediately above, is different now that we have two groups, as opposed to one (as in all the previous examples). We can see that of the 10 model parameters, 6 are free (i.e., 4 of the second group parameters are tied to their equivalents in the first group). Fitting this model results in a huge tex2html_wrap_inline10235 (not shown here), because our assumption that only a change in absorption can account for the spectral variation before and after eclipse is clearly wrong. Perhaps scattering also plays a role in reducing the flux before eclipse. This could be modeled (simply at first) by allowing the normalization of the power law to be smaller before eclipse than after eclipse. To decouple tied parameters, we change the parameter value in the second group to a value--any value--different from that in the first group (changing the value in the first group has the effect of changing both without decoupling). As usual, the newpar command is used:

XSPEC> newpar 10 1
    7 variable fit parameters
 Chi-Squared =     1.9666626E+07 using    66 PHA bins.
 Reduced chi-squared =      333332.7     for     59 degrees of freedom
 Null hypothesis probability = 0.000E+00

After fitting (not shown), this decoupling reduces tex2html_wrap_inline10235 by a factor of six to 15,478, but this is still too high. Indeed, this simple attempt to account for the spectral variability in terms of ``blanket'' cold absorption and scattering does not work. More sophisticated models, involving additional components and partial absorption, should be investigated.


next up previous contents
Next: Using XSPEC to Simulate Up: Walks through XSPEC Previous: Fitting Models to Data:

Keith Arnaud (kaa@genji.gsfc.nasa.gov)
Wed May 28 10:59:33 EDT 1997