Using XSPEC to Simulate Data: an Example from ASCA

In several cases, analyzing simulated data is a powerful tool to demonstrate feasibility. For example:

• To support an observing proposal. That is, to demonstrate what constraints a proposed observation would yield.

• To support a hardware proposal. If a response matrix is generated, it can be used to demonstrate what kind of science could be done with a new instrument.

• To support a theoretical paper. A theorist could write a paper describing a model, and then show how these model spectra would appear when observed. This, of course, is very like the first case.

Here, we'll use XSPEC to see how an ASCA observation of the elliptical galaxy NGC 4472 can constrain the condition of the hot gas. The first step is to define a model on which to base the simulation. The way XSPEC creates simulated data is to take the current model, convolve it with the current response matrix, while adding noise appropriate to the integration time specified. Once created, the simulated data can be analyzed in the same way as real data to derive confidence limits.

We begin by looking in the literature for the best estimate of the NGC 4472 spectrum. BBXRT observed the galaxy in 1990 and the results were published in Serlemitsos et al., (1993). They found a flux in the 0.5-4.5 keV range of $6.7\times 10^{-12}$ erg cm^-2 s^-1, a temperature range of 0.74 < kT < 0.98, an abundance range (as a fraction of solar) of 0.09 < A < 0.46 and a column range of $5.0\times 10^{20} < N_H < 3.7\times 10^{21}$ cm^-2. A Raymond-Smith spectral model was found to give a good fit. We specify this model at first with the median parameter values, except for the normalization of the Raymond-Smith, which we leave at its default value of unity at first (but adjust later):

XSPEC>mo pha(ray)
  Model:  phabs[1]( raymond[2] )
Input parameter value, delta, min, bot, top, and max values for ...
Current:          1     0.001         0         0     1E+05     1E+06
phabs:nH>0.21
Current:          1      0.01     0.008     0.008        64        64
raymond:kT>0.86
Current:          1    -0.001         0         0         5         5
raymond:Abundanc>0.27
Current:          0    -0.001         0         0         2         2
raymond:Redshift>/*
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( raymond[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.2100     +/-      0.
    2    2    2   raymond    kT       keV       0.8600     +/-      0.
    3    3    2   raymond    Abundanc           0.2700     frozen
    4    4    2   raymond    Redshift               0.     frozen
    5    5    2   raymond    norm                1.000     +/-      0.
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------

We now can derive the correct normalization by using the commands dummyrsp, flux and newpar. That is, we'll determine the flux of the model with the normalization of unity (this requires a response matrix to cover the BBXRT band--we use a dummy response here). We then work out the new normalization and reset it:

XSPEC> dummy 0.5 4.5
XSPEC>flux 0.5 4.5
 Model flux  0.2802     photons ( 4.9626E-10 ergs)cm**-2 s**-1 (  0.500-  4.500)
XSPEC> newpar 5 0.014
    3 variable fit parameters
XSPEC>flux
 Model flux  3.9235E-03 photons ( 6.9476E-12 ergs)cm**-2 s**-1 (  0.500-  4.500)

Here, we have changed the value of the normalization (the fifth parameter) from 1 to $6.7\times 10^{-12}/4.78\times 10^{-10} = 0.014$ to give the flux observed by BBXRT ( $6.7\times 10^{-12}$ erg cm^-2 s^-1 in the energy range 0.5-4.5).

The simulation is initiated with the command fakeit. If the argument none is given, the user will be prompted for the name of the response matrix. If no argument is given, the current response will be used:

XSPEC>fakeit none
For fake data, file #   1 needs response file: s0c1g0234p40e1_512_1av0_8i
              ... and ancillary response file: none

There then follows a series of prompts asking the user to specify whether he or she wants counting statistics (yes!), the name of the fake data file (ngc4472_sis.fak in our example), and the integration time T (40,000 seconds - cornorm can be left at its default value).

Use counting statistics in creating fake data? (y) /
Input optional fake file prefix (max 4 chars): /
 Fake data filename (s0c1g0234p40e1_512_1av0_8i.fak) [/ to use default]: ngc4472_sis.fak
 T, cornorm  (1, 0): 40000
 Net count rate (cts/s) for file   1  0.3563    +/-  3.0221E-03
   using response (RMF) file...       s0c1g0234p40e1_512_1av0_8i.rsp
 Chi-Squared =      188.6545     using   512 PHA bins.
 Reduced chi-squared =     0.3706375     for    509 degrees of freedom
 Null hypothesis probability =  1.00

We now have created a file containing a simulated spectrum of NGC 4472. As is usual before fitting, we need to check which channels to ignore. This time, we'll examine the actual numbers of counts in each channel and reject those that have fewer than 20 per channel. We use iplot counts and see that our criterion requires us to ignore channels 1-15 and 76-512:

XSPEC>ignore 1-15 76-** 
 Chi-Squared =      63.30437     using    60 PHA bins.
 Reduced chi-squared =      1.110603     for     57 degrees of freedom
 Null hypothesis probability = 0.264

As expected, $\chi^2$ is reasonable even before fitting because the model and the data have the same shape. But the point of this simulation is to determine confidence ranges. First, we thaw the value of the abundance (fixed by default), fit and then use the error command:

XSPEC> thaw 3
 Number of variable fit parameters =    4
XSPEC>fit
 Chi-Squared    Lvl  Fit param # 1     2           3           4
                 5
   55.3176     -3     0.2309      0.8569      0.2772          0.
              1.4322E-02
   55.2946     -4     0.2320      0.8565      0.2784          0.
              1.4322E-02
   55.2945     -5     0.2321      0.8565      0.2784          0.
              1.4322E-02
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               1     2     3     5
 1.51E-08 | -0.03 -0.01  0.03  1.00
 1.02E-05 |  0.39  0.91 -0.11  0.03
 9.98E-05 | -0.91  0.40  0.11 -0.03
 2.32E-04 | -0.15 -0.06 -0.99  0.03
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( raymond[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.2321     +/-  0.9426E-02
    2    2    2   raymond    kT       keV       0.8565     +/-  0.5048E-02
    3    3    2   raymond    Abundanc           0.2784     +/-  0.1510E-01
    4    4    2   raymond    Redshift               0.     frozen
    5    5    2   raymond    norm               1.4322E-02 +/-  0.5423E-03
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      55.29454     using    60 PHA bins.
 Reduced chi-squared =     0.9874024     for     56 degrees of freedom
 Null hypothesis probability = 0.502
XSPEC>err 1 2 3
 Parameter   Confidence Range (     2.706)
     1    0.217009        0.248102
     2    0.847909        0.864666
     3    0.254807        0.304992

These confidence ranges show that an ASCA observation would definitely constrain the parameters, especially the column and abundance, more tightly than the original BBXRT observation. Of course, whether these constraints are sufficient depends on the theories being tested. When producing and analyzing simulated data, it is crucial to keep in mind the purpose of the proposed observation, for the potential parameter space that can be covered with simulations is almost limitless.