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.
Most of the infeasible proposals from previous AO periods had one or more of the following failings:
No: all that is required in the proposal is an estimate of the SIS count rate. But without a simulation, you run the risk of proposing an infeasible observation. For those reluctant to invest the time in doing a full-blown simulation, there are two alternative--but not always applicable--routes for establishing feasibility:
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 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:
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
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
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/simulationsThe README file explains what the files are.
For simulations, the basic XSPEC plan of action is to:
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.
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. This set of response and background files will be used for all subsequent fitting and simulating, unless you read in a different set.
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 modelThe 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 zpowerlFor more information about these models, type, for example:
XSPEC> help model zwabsOnce 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 cm and cm 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+05At 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 . 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 erg cm s. 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 , 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.0025We're now ready to make some simulated ASCA 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.181To 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 /teTo plot the data, type:
XSPEC> plot dataOther 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.
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+00For 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.fakHere, 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.fakTo 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 = 3In 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 ( 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 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.89992This shows that with a 100-ks ASCA observation using one GIS and one SIS, the intrinsic absorption can be constrained such that and that , 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