Simulating an ACIS Spectrum for an On-axis Point Source, Including Pileup for Proposal Purposes

Run Sherpa and Read Instrument Response Files

Run sherpa:

unix% sherpa

-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Program
-----------------------------------------------------
Version: CIAO 3.4

Type AHELP SHERPA for overview.
Type EXIT, QUIT, or BYE to leave the program.

Notes:
    Temporary files for visualization will be written to the directory: 
    /tmp
    To change this so that these files are not deleted when you exit Sherpa,
    edit $ASCDS_WORK_PATH in your 'ciao' setup script.

    Abundances set to Anders and Grevesse

The CXC provides a set of instrument responses for use in proposal planning. For this thread we use responses appropriate for ACIS-S on-axis. These are the Area Response File ``aciss_aimpt_cy10.arf'' and the Redistribution Matrix ``aciss_aimpt_cy10.rmf''. The Area Response File contains the effective area information and the Redistribution Matrix gives the PHA to energy mapping.

The following command defines a name ``aciss'' for the responses (instrument files).

sherpa> rsp[aciss](aciss_aimpt_cy10.rmf,aciss_aimpt_cy10.arf)
The inferred file type is ARF.  If this is not what you want, please 
specify the type explicitly in the data command.

These two responses can be referenced by their name ``aciss'', as shown below.

sherpa> show aciss
rsp1d[aciss]
    Param   Type  Value
    -----   ----  -----
 1    rmf string: "aciss_aimpt_cy10.rmf" (N_E=1024,N_PHA=1024)
 2    arf string: "aciss_aimpt_cy10.arf" (N_E=1024)

The following command assigns the instruments files ``aciss'' to dataset number 1.

sherpa> instrument 1 = aciss

Define the Model

In Sherpa, model components are defined and labeled, allowing different model components to be applied differently to each dataset. The following command creates a model named ``pl''. This is the XSPEC powerlaw function, with parameters photon index and normalization. Note that the normalization is currently set to 1.0 and will be adjusted later so that the fake spectrum has the correct flux. Sherpa will prompt for the parameters:

sherpa> xspowerlaw[pl]
pl.PhoIndx parameter value [1] 1.6
pl.norm parameter value [1] 

Similarly, the following creates a model of the galactic absorption, g_abs, using the result obtained from colden:

sherpa> xswabs[g_abs] 
g_abs.nH parameter value [0.1] 0.026

Then tell Sherpa that the model source for dataset 1 is the galactic absorption modifying the power law

sherpa> source 1 = g_abs*pl

Run FAKEIT

Having defined a model for dataset 1, we can now create a fake datafile using the command FAKEIT. FAKEIT creates a simulated 1-D dataset from the source model and and instrument model, and then adds Poisson noise to the modeled data. FAKEIT has two parameters:

sherpa> show fakeit
Fakeit exposure time: 1 seconds.
Fakeit backscale: 1

The exposure time should be set to the required value (exposure time in seconds.) The backscale parameter is taken from the header of the background PHA file (if one is in use), but can also be set manually.

Set the exposure time to be 80 ks. This is an initial guess and might need to be changed later:

sherpa> fakeit time=80000

Now run FAKEIT:

sherpa> fakeit

The simulated dataset can be inspected using ``show data''. Note that the data is NOT written to disk until the ``write'' command is given.

sherpa> show data
Y Column: Counts/sec/keV
  Dimensions: 1
  Total Size: 1024 bins (or pixels)
  Axis: 0; Name: Energy (keV)
  Length: 1024 bins (or pixels)
File Name: fake
SubSection (if any): 
File Type: PHA
Exposure:       80000 sec
Background Scale:           1
[0.010950] =          0
[0.021900] =          0
[0.036500] =          0
[0.051100] =          0
[0.065700] =          0
[0.080300] =          0
[0.094900] =          0
[0.109500] =      4.393
[0.124100] =    14.7817
[0.138700] =    15.8579
[0.153300] =    18.1618
[0.167900] =    25.9375
[0.182500] =    43.8142

---output truncated---

Normalize the Simulated Data

The simulated data file currently does not have the correct normalization--the flux of the simulated data is incorrect because the power law normalization was arbitrarily set to 1.0. To correct the flux we need to adjust the normalization. First compute the flux predicted by the model in units of photons/cm**2/s and ergs/cm**2/s, respectively, using the FLUX and EFLUX commands:

sherpa> flux(2:10)
Flux for source dataset 1: 0.679321 photons/cm**2/s
sherpa> eflux(2:10)
Flux for source dataset 1: 4.7691e-09 ergs/cm**2/s

The 2-10 keV flux of this source has been measured at 2.70E-14 ergs s-1 cm-2 in the low state. Therefore the correct normalization is 2.7e-14/4.7691e-09=5.66e-6. The normalization of the model can now be adjusted and FAKEIT run again:

sherpa> pl.norm=5.66e-6
sherpa> fakeit

Now check that the flux is correct using FLUX and EFLUX as above:

sherpa> flux(2:10)
Flux for source dataset 1: 3.84496e-06 photons/cm**2/s
sherpa> eflux(2:10)
Flux for source dataset 1: 2.69931e-14 ergs/cm**2/s

The flux is now correct.

The total number of counts in the simulated dataset created by fakeit can be derived using the DCOUNTS command (dcounts). The following gives the number of counts between 0.3 and 8keV:

sherpa> dcounts(0.3:8)
Observed counts for source dataset 1: 388 counts

sherpa> mcounts(0.3:8)
Model counts for source dataset 1: 391.111 counts

Fit the Simulated Spectrum and Calculate Uncertainties

This is a small number of counts for a spectral fit, but it should be possible to constrain the spectral slope. To estimate the uncertainties in the slope when the source is in the low state, freeze the galactic absorption and fit the fake spectrum.

sherpa> freeze  g_abs.nh
sherpa>  fit
 LVMQT: V2.0
 LVMQT: initial statistic value = 63.7788
 LVMQT: final statistic value = 56.2718 at iteration 4
            pl.PhoIndx  1.52874     
            pl.norm  4.2475e-06  photons/keV/cm**2/s at 1 keV  

 

Now use the covariance command to get 1 sigma uncertainties. Setting the set_verbose parameter allows you to view covariance matrix and eigenvector information:

sherpa> set_verbose(2)
sherpa> covariance
Information Matrix (Second Derivatives of Fit Statistic): 

     pl.PhoIndx      pl.norm    
         20.1981    -1.92991e+06
    -1.92991e+06     3.91217e+12

Eigenvectors (Principal Axes) of the Covariance Matrix: 

     pl.PhoIndx      pl.norm    
               1     -4.9331e-07
      4.9331e-07               1

Eigenvalues of the Covariance Matrix: 

       0.0519588     2.55613e-13

Covariance Matrix (Inverse of Information Matrix): 

     pl.PhoIndx      pl.norm    
       0.0519588     2.56318e-08
     2.56318e-08     2.68257e-13

Covariance Matrix Determinant (Product of Eigenvalues): 1.32813e-14


Computed for sherpa.cov.sigma = 1
        --------------------------------------------------------
        Parameter Name      Best-Fit Lower Bound     Upper Bound
        --------------------------------------------------------
            pl.PhoIndx       1.52874  -0.227945       +0.227945     
            pl.norm       4.2475e-06  -5.17936e-07    +5.17936e-07  
   

With this many counts, the photon index can be measured to about 30%, as required. This means that even if the source is in the lowest observed state the slope can be meaningfully constrained. If the source is brighter, then the uncertainties will be reduced.

Alternate Strategy: Normalize Spectrum Using Predicted Count Rate

The spectrum can also be normalized to the Chandra predicted count rate. In this case the steps are:

• Obtain predicted Chandra count rate (e.g. from PIMMS ).
• Follow steps outlined above to get a fake spectrum with normalization of 1.0.
• Use the Sherpa SHOW command to get the count rate in the fake spectrum with normalization=1.0.
• Adjust the normalization so that the count rate in the new spectrum is correct.

Getting Started

In the high state the flux is 100 times that of the low state and pileup must be considered. If you haven't followed the first part of the thread, repeat the first four steps as in the low state case above (Read Responses through Normalize Simulated Data) using a power law normalization of 100 times the low state normalization (5.66e-4). If you have followed the steps above, reset the slope to 1.6, set the normalization to 5.66e-4, and rerun FAKEIT.

sherpa> pl.phoindx=1.6
sherpa> pl.norm=5.66e-4
sherpa> fakeit

Check that normalization is as expected:

sherpa> flux(2:10)
Flux for source dataset 1: 0.000384496 photons/cm**2/s
sherpa> eflux(2:10)
Flux for source dataset 1: 2.69931e-12 ergs/cm**2/s

Define Pileup Model

Pileup can be included in the simulation by defining a pileup model, which is treated in the same way as any other spectral model. The model, jdpileup, was developed by John Davis at MIT. The parameter `alpha' is the probability that two piled photons are assigned a ``bad'' grade and rejected. It should be left at the default value of 0.5 and frozen. One means of pileup mitigation is to change the frame time parameter, ftime. It is strongly recommended that the other four parameters be frozen at the default values unless the proposer is thoroughly familiar with the pileup model. A detailed discussion is provided by Davis (2001, ApJ, 562, 575).

Here we define a pileup model (jdp).

sherpa> jdpileup[jdp]
jdp.alpha parameter value [0.5] 
jdp.g0 parameter value [1] 
jdp.f parameter value [0.95] 
jdp.n parameter value [1] 
jdp.ftime parameter value [3.241] 
jdp.fracexp parameter value [0.987] 

Here we use the command ``pileup'' to tell sherpa which pileup model to use (jdp, defined above) and to freeze alpha.

sherpa> pileup = jdp
sherpa> freeze jdp.alpha

Simulate Spectrum with Pileup

Now simulate the piled spectrum and plot the results (ignoring energies higher than 8 keV):

sherpa> fakeit
sherpa> ignore energy 8:
sherpa> lplot  data 1

The plot is shown here . To illustrate the effect of pileup we have also plotted the spectrum without pileup. The unpiled spectrum has a peak count rate of slightly more than 0.3 counts/sec/keV. In contrast, the piled spectrum has a peak count rate of 0.1 counts/sec/keV because many photons are rejected and/or many low energy photons are mistakenly counted as high energies. Additional information about pileup is available in the ACIS chapter of the POG.

Pilup information for the simulated spectrum can be displayed with the SHOW command.

sherpa> show pileup
   1: 0.327826  0.660215
   2: 0.255685  0.257465
   3: 0.132946  0.0669357
   4: 0.0518451  0.0130515
   5: 0.0161744  0.00203588
   6: 0.00420503  0.000264644
   7: 0.000937051  2.94866e-05
   8: 0.000182711  2.87473e-06
   *** pileup fraction: 0.339785

The last line gives the pileup fraction. Lines 1-7 give further details of the pileup statistics. Line 1 indicates that 33% of the frames contained 1 photon in the pileup region, and 66% of the events were single photon events. Line 2 shows that 26% of the frames contained 2 photons, and 26% of the events were 2 photon events. For further details please see the ISIS manual.

Pileup Mitigation

There are a number of ways to configure ACIS to minimize pileup. These are discussed in the POG and include offset pointing, using a transmission grating, using CC mode, alternating exposures, and reducing the frame time with a subarray. For the case of JXXXX.X-XXX, we do not know in advance whether the spectrum will be piled. The source might actually be quite faint. In this case, we wouldn't use the first four options because they are not optimized for faint sources. We are therefore left with only two choices -- either use a subarray or do not take any steps to mitigate pileup. Most of the ACIS field-of-view is lost if a subarray is chosen. This would be a big disadvantage if there were other sources in the field which formed part of the proposal. Since JXXXX.X-XXX is a point source, a subarray is a good choice.

Sometimes, regardless of whether or not a subarray has been used, the pileup fraction is larger than anticipated. In this case, it is possible to estimate the pileup and recover the source spectrum using the pileup model described earlier in this section. Please see the CIAO "Fitting Spectral Data: Using A Pileup Model" thread on fitting a piled spectrum for details.

Simulating Pileup with Subarray

A subarray reduces pileup by shortening the frame time. With a shorter frame time, the probability that two photons are detected as a single event decreases. Here we simulate a 1/8 subarray, which correspons to a frame time of 0.4 sec. Set the frame time to be 0.4 sec and then run FAKEIT again:

sherpa> jdp.ftime=0.4
sherpa> fakeit
sherpa> show pileup
   1: 0.158805  0.952677
   2: 0.0152864  0.045852
   3: 0.000980973  0.00147123
   *** pileup fraction: 0.0473232

The pileup fraction is now reduced significantly. This can also be seen in the simulated spectrum .

Finally, check the number of counts between 0.3 - 8 keV for the simulated source spectrum (High state):

sherpa> dcounts(0.3:8)
Observed counts for source dataset 1: 34252 counts