next up previous contents
Next: Producing Plots: Modifying the Up: Walks through XSPEC Previous: Simultaneous Fitting: Examples from

Using XSPEC to Simulate Data: an Example from ASCA

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

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 tex2html_wrap_inline10379 erg cm tex2html_wrap_inline12821 s tex2html_wrap_inline10317 , 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 tex2html_wrap_inline10373 cm tex2html_wrap_inline12821 . 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) 
  mo = phabs[1]( raymond[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.21
 Mod parameter  2 of component  2 raymond  kT       keV
    1.000      1.0000E-02  8.0000E-03  8.0000E-03   64.00       64.00
0.86
 Mod parameter  3 of component  2 raymond  Abundanc
    1.000     -1.0000E-03  0.0000E+00  0.0000E+00   5.000       5.000
0.27
 Mod parameter  4 of component  2 raymond  Redshift
   0.0000E+00 -1.0000E-03  0.0000E+00  0.0000E+00   2.000       2.000

 Mod parameter  5 of component  2 raymond  norm
    1.000      1.0000E-02  0.0000E+00  0.0000E+00  1.0000E+24  1.0000E+24

  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  mo = phabs[1]( raymond[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.2100     +/-  0.0000E+00
    2    2    2   raymond    kT       keV       0.8600     +/-  0.0000E+00
    3    3    2   raymond    Abundanc           0.2700     frozen
    4    4    2   raymond    Redshift           0.0000E+00 frozen
    5    5    2   raymond    norm                1.000     +/-  0.0000E+00
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
    3 variable fit parameters

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.2673     photons ( 4.7803E-10 ergs)cm**-2 s**-1 (  0.500-  4.500)
XSPEC> newpar 5 0.014
    3 variable fit parameters
XSPEC> flux
 Model flux  3.7424E-03 photons ( 6.6924E-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 tex2html_wrap_inline10377 to give the flux observed by BBXRT ( tex2html_wrap_inline10379 erg cm tex2html_wrap_inline12821 s tex2html_wrap_inline10317 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.rsp
              ... 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--the other quantities A, Bkg, cornorm can be left at their default values).

Use counting statistics in creating fake data? (y)
 Input optional fake file prefix (max 12 chars):
 Override default values for file #    1
 Fake data filename (sis_small.fak): ngc4472_sis.fak
 T, A, Bkg, cornorm ( 1.0000 , 1.0000 , 1.0000 , 0. ): 40000
 Net count rate (cts/s) for file   1  0.3503    +/-  2.9973E-03
   using response (RMF) file...       s0c1g0234p40e1_512_1av0_8i.rsp
 Chi-Squared =      168.1253     using   512 PHA bins.
 Reduced chi-squared =     0.3303052     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> ign 1-15 76-** Chi-Squared =      53.13371     using    60 PHA bins.
 Reduced chi-squared =     0.9321703     for     57 degrees of freedom
 Null hypothesis probability = 0.621

As expected, tex2html_wrap_inline10235 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

   47.3996     -3     0.2011      0.8696      0.2705      0.0000E+00
              1.3923E-02
   47.3954     -4     0.2013      0.8693      0.2707      0.0000E+00
              1.3928E-02
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               1     2     3     5
 1.45E-08 | -0.03 -0.01  0.03  1.00
 9.81E-06 |  0.41  0.91 -0.10  0.03
 8.11E-05 | -0.91  0.41  0.07 -0.03
 2.14E-04 | -0.10 -0.06 -0.99  0.03
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  mo = phabs[1]( raymond[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.2013     +/-  0.8399E-02
    2    2    2   raymond    kT       keV       0.8693     +/-  0.4775E-02
    3    3    2   raymond    Abundanc           0.2707     +/-  0.1453E-01
    4    4    2   raymond    Redshift           0.0000E+00 frozen
    5    5    2   raymond    norm               1.3928E-02 +/-  0.5187E-03
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      47.39543     using    60 PHA bins.
 Reduced chi-squared =     0.8463470     for     56 degrees of freedom
 Null hypothesis probability = 0.787
XSPEC> err 1 2 3
 Parameter   Confidence Range (     2.706)
     1    0.187992        0.215481
     2    0.861300        0.876959
     3    0.248081        0.296034

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.


next up previous contents
Next: Producing Plots: Modifying the Up: Walks through XSPEC Previous: Simultaneous Fitting: Examples from

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