next up previous contents
Next: PIMMS, VIEWING and ROLL Up: ASCA AO7 Appendix E Previous: Feasibility: SIS+XRT

Subsections

Simulations

Introduction

Proposers are required, as a minimum, to estimate the SIS count rates of their targets. Full simulations are optional. But whether they perform simulations or not, the claims proposers make will be subjected to a thorough feasibility assessment by the GOF staff. Many previous proposals, unsupported by simulations, were rejected outright on feasibility grounds, regardless of the how interesting the science appeared. Clearly, including simulations in proposals is a good idea.

This chapter is arranged as follows. The next section contains feasibility recommendations based on our experience of the the first three ASCA AOs. Using PIMMS to estimate GIS and SIS count rates is described next, followed by a description of how to use the spectral fitting program XSPEC to produce and analyze simulated ASCA data. Both these last two sections contain several, fully worked examples.

Helpful hints for good feasibilities

Common mistakes

Most of the infeasible proposals from previous AO periods had one or more of the following failings:

Do I really have do a full-blown simulation?

Not necessarily. In some cases you may only require an estimate of the SIS count rate in the proposal. If you do not do simulations you must demonstrate your feasibility in some other fashion, or else run the risk of proposing an infeasible observation. Infeasible observations will be caught by either the technical reviewers or by memebers of the peer-review panel, ultimately resulting in rejection. Proposers are warned especially to take care with SIS feasibility since the energy resolution has degraded to the extent that it is only a factor of $\sim 2$ better than the GIS in the Fe K band. SIS matrices with degraded energy resolution are now available with the files you will download in §10.4 (please consult the 'README' that comes with the files for details).

In some cases, depending on the proposed science, you may use one of the two alternatives below, to a full-blown simulation (certainly in most cases for the GIS).

Estimating source count rate using PIMMS

As a first step in investigating the feasibility of a proposed observation, the count rate in the ASCA instruments should be estimated. We have provided PIMMS (Portable, Interactive, Multi-Mission Simulator) for this purpose. The following examples show actual PIMMS sessions (some details are different in newer versions).

PIMMS terminology

PIMMS uses the command GO for actual execution, while other commands are mostly used for setting up various parameters. This approach allows users to repeat similar calculations using a slightly different set of parameters. PIMMS uses the following terms:

Example I: estimating ASCA count rates

In this example, the default spectral model is used to estimate the ASCA SIS count rate (the default product). The only place where the default set-up is not used is to specify the conversion from Einstein IPC count rate.

*** PIMMS Version 1.0b (first release) ***
    Reading mission directory, please wait
* PIMMS simulation product is COUNT
        count rates for various instruments or intrinsic fluxes can be 
estimated
   <--- Use 'PRODUCT' to simulate images
* Current spectral model is BREMSSTRAHLUNG, kT=  10.0000 keV; NH 
=  1.000E+21
   <--- Use 'MODEL' command to change
* By default, input rate is taken to be
 Flux (  2.000- 10.000 keV) in ergs/cm/cm/s
   <--- Use 'FROM' command to change the default
* Simulation product will be
 Count rate in ASCA SIS
   <--- Use 'INSTRUMENT' command to switch to another instrument
PIMMS > go 1.2 einstein ipc
* For thermal Bremsstrahlung model with kT= 10.0000 keV; NH =  
1.000E+21
  and  1.200E+00 cps in EINSTEIN IPC
  (Model normalization =  7.059E-03)
* PIMMS predicts  1.956E+00 cps with ASCA SIS
* An exposure of      12.78s is required for a 5-sigma detection

              Percent Telemetry Full
                 Faint Mode        Bright Mode   Fast Mode
               4CCD  2CCD  1CCD  4CCD  2CCD  1CCD  1CCD
Bit Rate High   12     6     3     3     2     1     0
Bit Rate Med.   98    49    24    24    12     6     3
Bit Rate Low    **    **    98    98    49    24    12
* Recommended mode (assuming a point source)
  is Bright (4-CCD) mode
  Medium bit rate is sufficient for your observation
PIMMS > quit

Example II: estimating ASCA count rates

In this example, a 1.0-keV Raymond-Smith model is used, with two different values of interstellar absorption. After estimating ASCA SIS count rate, we switch to ASCA GIS count rate. The conversion is from EXOSAT ME count rate. Here, we use the FROM command to change the default rather than explicitly specifying within the GO command.

*** PIMMS Version 1.0b (first release) ***
    Reading mission directory, please wait
* PIMMS simulation product is COUNT
        count rates for various instruments or intrinsic fluxes can be 
estimated
   <--- Use 'PRODUCT' to simulate images
* Current spectral model is BREMSSTRAHLUNG, kT=  10.0000 keV; NH 
=  1.000E+21
   <--- Use 'MODEL' command to change
* By default, input rate is taken to be
 Flux (  2.000- 10.000 keV) in ergs/cm/cm/s
   <--- Use 'FROM' command to change the default
* Simulation product will be
 Count rate in ASCA SIS
   <--- Use 'INSTRUMENT' command to switch to another instrument
PIMMS > from exosat me
PIMMS > model rs 1.0 5e19
PIMMS > go 1.0
* For Raymond Smith model with kT=  1.0000 keV; NH =  5.000E+19
  and  1.000E+00 cps in EXOSAT ME
  (Model normalization =  3.535E-02)
* PIMMS predicts  3.072E+00 cps with ASCA SIS
* An exposure of       8.14s is required for a 5-sigma detection

              Percent Telemetry Full
                 Faint Mode        Bright Mode   Fast Mode
               4CCD  2CCD  1CCD  4CCD  2CCD  1CCD  1CCD
Bit Rate High   19    10     5     5     2     1     1
Bit Rate Med.   **    77    38    38    19    10     5
Bit Rate Low    **    **    **    **    77    38    19
* Recommended mode (assuming a point source)
  is Bright (4-CCD) mode
  Medium bit rate is sufficient for your observation
PIMMS > model rs 1.0 1e20
PIMMS > go 1.0
* For Raymond Smith model with kT=  1.0000 keV; NH =  1.000E+20
  and  1.000E+00 cps in EXOSAT ME
  (Model normalization =  3.543E-02)
* PIMMS predicts  3.042E+00 cps with ASCA SIS
* An exposure of       8.22s is required for a 5-sigma detection

              Percent Telemetry Full
                 Faint Mode        Bright Mode   Fast Mode
               4CCD  2CCD  1CCD  4CCD  2CCD  1CCD  1CCD
Bit Rate High   19    10     5     5     2     1     1
Bit Rate Med.   **    76    38    38    19    10     5
Bit Rate Low    **    **    **    **    76    38    19
* Recommended mode (assuming a point source)
  is Bright (4-CCD) mode
  Medium bit rate is sufficient for your observation
PIMMS > inst asca gis
PIMMS > go 1.0
* For Raymond Smith model with kT=  1.0000 keV; NH =  1.000E+20
  and  1.000E+00 cps in EXOSAT ME
  (Model normalization =  3.543E-02)
* PIMMS predicts  1.686E+00 cps with ASCA GIS
PH mode can be used at low bit rate
  (56% of telemetry will be used with this)
PIMMS > quit

Using XSPEC for Spectral Simulations

Getting started

XSPEC is a large, sophisticated and powerful program equipped with many built-in models and extensive on-line help. But novice users can produce and analyze simulated data without being experts. In fact, the GOF has made the process a bit easier and more accurate: ``dummy'' PHA files are now available which include real background. To obtain the requisite files, FTP anonymously to heasarc.gsfc.nasa.gov and get all the files in the directory:

asca/nra_info/simulations
The README file explains what the files are.

For simulations, the basic XSPEC plan of action is to:

More details and explanation follow. Here, what appears on the screen during the sample XSPEC session is rendered in typewriter font. Some of the screen output might be slightly different, depending on what sort of computer you're using. For more information about XSPEC, please consult the on-line help or the XSPEC manual available at the URL

http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/manual.html.

Note that there were many significant changes between XSPEC v10.0 and previous versions. The examples that follow refer to XSPEC v10.0. They will work with earlier versions with only occassional changes in command format. Again consult the XSPEC on-line help.

Reading in the dummy data

This is simple. At the XSPEC> prompt, just use the data command followed by the name of the dummy data file. In this example, we'll read in both the GIS and SIS data files to simulate simultaneous fitting:

XSPEC> data g2_dummy.pha s0c1_dummy.pi
   .
   .
   .
 Net count rate (cts/cm^2/s) for file   1  2.3788E-02+/-  
    1.2142E-03( 51.0% total)
 Net count rate (cts/cm^2/s) for file   2  3.9099E-02+/-  
    9.0386E-04( 65.1% total)
(You can ignore the error messages about mismatching header keywords-- they're harmless). The headers of the dummy data files contain the names of the appropriate GIS and SIS response and background files g2.rsp & g2_back.pha; s0c1.rsp & s0c1_back.pi, respectively.10.2 This set of response and background files will be used for all subsequent fitting and simulating, unless you read in a different set.

Defining a model

In XSPEC, models are built up from the suite of model components. To see a list of the models available in XSPEC, type:

 XSPEC> model ?
And for more information about the model command, type:
 XSPEC> help model
The example we'll be working through is to determine what sort of constraints an ASCA observation of a high-redshift quasar would yield. Simulated data should be treated the same way as real data, e.g., rebinned when necessary and fitted with more than one model. If this example seems rather ``involved'', it's partly because real analysis is also involved. At first, we'll assume that the spectrum is a power-law (redshifted, of course) with two absorption components: one due to the Galaxy, the other intrinsic to the quasar. In XSPEC, this model is specified like this:
 XSPEC> model wabs zwabs zpowerl
For more information about these models, type, for example:
 XSPEC> help model zwabs
Once the model is specified, the user will be prompted for the parameters. In this example, we're interested in a quasar at a redshift of 2.5 and with absorption components of $7\times10^{19}$ cm-2 and $1\times10^{22}$ cm-2 due to the Galaxy and the quasar, respectively. We don't know the slope of the quasar spectrum, but will assume it's the same as nearby Seyfert I AGN, i.e., a ``canonical'' photon index of 1.7. Note that to enter the default values, hit return:
 Input parameter value, delta, min, bot, top, and max values for ...
 Mod parameter  1 of component  1  1 wabs     nH 10^22
    1.000      1.0000E-03      0.          0.      1.0000E+05  1.0000E+06
 0.007
 Mod parameter  2 of component  2  1 zwabs    nH 10^22
    1.000      1.0000E-03      0.          0.      1.0000E+05  1.0000E+06

 Mod parameter  3 of component  2  2 zwabs    redshift
       0.     -1.0000E-02      0.          0.       10.00       10.00
 2.5
 Mod parameter  4 of component  3  1 zpowerlw PhoIndex
    1.000      1.0000E-02  -10.00      -9.000       9.000       10.00
 1.7
 Mod parameter  5 of component  3  2 zpowerlw Redshift
       0.     -1.0000E-02      0.          0.       10.00       10.00
 2.5
 Mod parameter  6 of component  3  0 zpowerlw norm
    1.000      1.0000E-03      0.          0.      1.0000E+05  1.0000E+06

  ----------------------------------------------------------------------
  mo = wabs[1] zwabs[2] (zpowerlw[3])
  Model Fit Model   Component     Parameter   Value
  par   par comp  index  type
    1    1    1    1  wabs        nH 10^22  7.000000E-03 +/-      0.
    2    2    2    1  zwabs       nH 10^22   1.00000     +/-      0.
    3    3    2    2  zwabs       redshift   2.50000     frozen
    4    4    3    1  zpowerlw    PhoIndex   1.70000     +/-      0.
    5    5    3    2  zpowerlw    Redshift   3.50000     frozen
    6    6    3    0  zpowerlw    norm       1.00000     +/-      0.
  ----------------------------------------------------------------------
    4 variable fit parameters
  Chi-Squared =     6.9425E+08 using   768 PHA bins.
  Reduced chi-squared =     9.0871E+05
At this stage, the user need only specify the initial parameter values: the other quantities (delta, min, bot, top, and max) are required by the fitting process but can be left at their default values for now. Notice that the normalization has not been changed from its default value of unity, which accounts for the very large $\chi^2$. To get the right value for the quasar, we need the expected count rate (from PIMMS) or the flux (from elsewhere). In this example, we'll assume that we expect the 2-10 keV flux to be $1.2\times 10^{-12}$ erg cm-2 s-1. First, we derive the current model flux with the flux command which integrates the model over the energy bounds specified:
XSPEC> flux 2 10
 Model flux 4.5878E-02 photons (4.8259E-10 ergs)cm**-2 s**-1
   (  2.000- 10.000)
Given that the model currently has a normalization of unity, flux implies that normalization should be $1.2\times 10^{-12}/4.8259\times 10^{-12}$, i.e., 0.0025. To change the value of any parameter, use the newpar command, the arguments of which are the parameter number (6 in our example) and the new value (0.0025). Like this:
XSPEC> newpar 6 0.0025
We're now ready to make some simulated ASCA data.

Faking data

The fakeit command initiates the production of simulated data. In this process, the current model, including background, is convolved with the instrument response to produce a fake data file. Counting statistics are also included and the current background file subtracted. Here, since we've read in two dummy files, fakeit will produce two fake files. We'll name these files g2_quasar.fak and s0c1_quasar.fak and integrate for 100 ks:

 XSPEC> fakeit
 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 (g2_dummy.fak): g2_quasar.fak
  T, A, Bkg, cornorm ( 1.00000E+05, 1.0000 , 2.87320E-02, 1.0000 ): 100000
  *Prefixed fake file name is too long:s0c1_dummy will be truncated.
 Override default values for file #    2
  Fake data filename (s0c1_dummy.fak): s0c1_quasar.fak
  T, A, Bkg, cornorm ( 1.00000E+05, 1.0000 , 4.53320E-02, 1.0000 ): 100000
  Net count rate (cts/cm^2/s) for file   1  3.9570E-02+/-  9.3291E-04( 77.7% tot
al)
  Net count rate (cts/cm^2/s) for file   2  4.8399E-02+/-  9.5381E-04( 69.8% tot
al)
  Chi-Squared =      902.5     using   768 PHA bins.
  Reduced chi-squared =      1.181
To see what these fake spectra look like, first set the plot device appropriately using the cpd (change plot device) command. (In my session, I'm using the Tektronics emulator, so the device is /te. To see the list of available plotting devices, type cpd ?):
 XSPEC> cpd /te
To plot the data, type:
 XSPEC> plot data
Other things can be plotted. To see the options, type plot ? Notice that when a fake spectrum is plotted, the originating model is also plotted, along with the residuals of the ``fit''. If you want a plot in energy rather than channels, use the command set energy.

At last, we can answer some questions of scientific feasibility by fitting the fake data.

Fitting the fake data

Before fitting the data, it is prudent to check which channels have data in them. This is because empty or sparsely filled channels can give misleading behavior of the fit statistic and model parameters. For a discussion of this issue, those interested should consult Nousek & Shue (1989, ApJ, 342, 1207, and references therein).

The contents of the fake data file can be viewed by running the program fdump which is part of the ftools package. This is done most conveniently by spawning out a process within XSPEC. (To spawn for versions of XSPEC earlier than 10.0, use a dollar sign).

 XSPEC> fdump g2_quasar.fak
 Spawning...
Name of optional output file[STDOUT]
Names of columns[-]
Lists of rows[-]
SIMPLE  =                    T / file does conform to FITS standard
BITPIX  =                    8 / number of bits per data pixel     
   .
   .
   .
After scrolling by the header information, we reach the ASCII dump of the extension which contains the PHA information. The first column is the row number; the second is the channel number; the third is the number of counts in the channel.
        CHANNEL COUNTS
                COUNTS
      1      1    0.0000000E+00
      2      2    3.0000000E+00
      3      3    3.0000000E+00
      4      4    4.0000000E+00
           .
           .
           .
    254    254    0.0000000E+00
    255    255    0.0000000E+00
    256    256    0.0000000E+00
For our file, g2_quasar.fak, about half the channels are inadequately filled, which is hardly surprising given the weakness of the target. If these were real data, the next step would be to rebin: simulated data should be treated the same way. The ftool grppha can do this:
XSPEC> $ grppha
Please enter PHA filename[dmp] g2_quasar.fak
Please enter output filename[] g2_quasar_b.fak
   .
   .
   .
GRPPHA[] group min 20
 Spare channels  238: 256
GRPPHA[] exit
 Exiting, Changes written to file : g2_quasar_b.fak
Here, we have rebinned--i.e., ``grouped'' together--channels such that each contains a minimum of 20 counts. The rebinned file is named g2_quasar_b.fak. The ``spare channels'' are those at the end of the channel range which, even when grouped together, still contain fewer than our minimum of 20. These spare channels have been flagged as bad by grppha and can be conveniently ignored in fitting (see below). For more information about grppha, try its on-line help. After doing the same thing to s0c1_quasar.fak, we read in the new, rebinned files:
XSPEC> data g2_quasar_b.fak s0c1_quasar_b.fak
To ignore the spare, ``bad'' channels, we use the ignore bad command:
XSPEC> ignore bad

We are now almost ready to fit the data. Before doing so, we must decide on what model to fit to the data. This depends on the scientific purpose of the proposed observation. Here, we have assumed that our quasar is really like a redshifted Seyfert I AGN with intrinsic absorption: but can we really identify the two absorption components?: can we constrain the slope of the continuum? In general, when fitting models to spectral data, it's best to proceed incrementally, i.e., first to keep some parameters fixed (``frozen'' in XSPEC terms) and then, once convergence has been achieved, to free the remaining parameters (``thaw'' in XSPEC terms). The model we used to create the simulated data is still the current model. We'll begin by first freezing the Galactic absorption:

 XSPEC> freeze 1
 Number of variable fit parameters =    3
In the model, the three variable parameters are the intrinsic absorption and the slope and normalization of the power law. The redshifts of the absorption and power law are frozen, by default, at their initial values. Fitting is set in motion using the command fit:
 XSPEC> fit
  Chi-Squared  Lvl  Fit param # 1     2           3           4
                  5           6
    159.41     -3     7.0000E-02   1.240       2.500       1.808
                2.500      2.9586E-03
    158.16     -4     7.0000E-02   1.301       2.500       1.820
                2.500      3.0896E-03
    158.16     -5     7.0000E-02   1.300       2.500       1.820
                2.500      3.0927E-03
   ---------------------------------------------------------------------------
   mo = wabs[1] zwabs[2] (zpowerlw[3])
   Model Fit Model   Component     Parameter   Value
   par   par comp
     1    1    1       wabs        nH 10^22  7.000000E-02 frozen
     2    2    2       zwabs       nH 10^22   1.30035     +/- 0.35751
     3    3    2       zwabs       redshift   2.50000     frozen
     4    4    3       zpowerlw    PhoIndex   1.82042     +/- 0.46748E-01
     5    5    3       zpowerlw    Redshift   2.50000     frozen
     6    6    3       zpowerlw    norm      3.092726E-03 +/- 0.30527E-03
   ---------------------------------------------------------------------------
  Chi-Squared =      158.2     using   277 PHA bins.
  Reduced chi-squared =     0.5772

Not surprisingly, the fit is good ($\chi^2$ is 158.2 for 274 degrees of freedom) and the best-fitting parameters are consistent with the original values used to generate the fake data. But one of the points of a simulation is to see how accurately the parameters can be constrained. Although the result of fit includes errors, these errors are approximate. To get statistically meaningful errors, the error command should be used. The default is 90 per cent confidence for one interesting parameter. In the following the default is used for illustration but for realistic estimates you should use the $\Delta \chi^{2}$appropriate for the number of interesting parameters in your model. For parameters 2 & 4: (intrinsic absorption column and slope):

 XSPEC> error 2 4
 Fit Index   Confidence Range ( 2.706)
      2    0.766446         1.92867
      4     1.74461         1.89992
This shows that with a 100-ks ASCA observation using one GIS and one SIS, the intrinsic absorption can be constrained such that $7.7\times 10^{21} < N_{\rm H} < 1.9\times 10^{22}$ and that $1.7 < \Gamma < 1.9$, provided the Galactic column is known a priori. If we allow the Galactic column to be a free parameter, we get a slightly different result, but the parameters are still constrained:
 XSPEC> thaw 1
  Number of variable fit parameters =    4
 XSPEC> fit
  Chi-Squared  Lvl  Fit param # 1     2           3           4
                  5           6
    157.04     -2     2.7383E-02   1.956       2.500       1.803
                2.500      2.9457E-03
    156.65     -2     4.6359E-03   2.285       2.500       1.787
                2.500      2.8365E-03
    156.62     -1     2.4715E-03   2.316       2.500       1.785
                2.500      2.8270E-03
    156.60     -1     4.9028E-04   2.344       2.500       1.784
                2.500      2.8170E-03
    156.56     -2     0.0000E+00   2.320       2.500       1.778
                2.500      2.7799E-03
    156.56     -3     0.0000E+00   2.301       2.500       1.775
                2.500      2.7626E-03
    156.56      0     0.0000E+00   2.301       2.500       1.775
                2.500      2.7626E-03
   ---------------------------------------------------------------------------
   mo = wabs[1] zwabs[2] (zpowerlw[3])
   Model Fit Model   Component     Parameter   Value
   par   par comp
     1    1    1       wabs        nH 10^22  0.000000E+00 +/- -1.0000
     2    2    2       zwabs       nH 10^22   2.30107     +/- 0.37728
     3    3    2       zwabs       redshift   2.50000     frozen
     4    4    3       zpowerlw    PhoIndex   1.77465     +/- 0.47094E-01
     5    5    3       zpowerlw    Redshift   2.50000     frozen
     6    6    3       zpowerlw    norm      2.762617E-03 +/- 0.27726E-03
   ---------------------------------------------------------------------------
  Chi-Squared =      156.6     using   277 PHA bins.
  Reduced chi-squared =     0.5735
 XSPEC> err 1 2 4
  *WARNING*: Parameter sigma indicates possible error:     -1.0000
  Parameter pegged at hard limit    0.0000E+00 with delta ftstat=  0.
      1    0.000000E+00    9.777232E-02 
  Fit Index   Confidence Range (     2.706)
      2    0.830053         2.97709
 Number of trials exceeded - last iteration delta =   2.4414E-04
  Continue fitting? (Y)
      4     1.69781         1.87123


next up previous contents
Next: PIMMS, VIEWING and ROLL Up: ASCA AO7 Appendix E Previous: Feasibility: SIS+XRT
Michael Arida
1998-06-02