Using A Pileup Model

Remove the acis_detect_afterglow Correction

An afterglow is the residual charge from the interaction of a cosmic ray in a CCD. Standard data processing (SDP, aka "the pipeline") uses the tool acis_detect_afterglow to flag possible cosmic ray events in the level 1 event file; these are then filtered out in the level 2 event file. It has been determined, however, that 3-5 % of the valid source photons may be rejected from diffracted spectra. These rejections, though a small fraction of the total events, are systematic and non-uniform.

In order to accurately model the pileup, it is necessary to remove this correction so that no source photons are eliminated from the event file. Read the Pileup Talk for more information, in particular Data Preparation and Caveats. The Remove the acis_detect_afterglow correction thread shows how to undo the afterglow filtering.

Important Note

A new, more precise method for identifying afterglow events was introduced to SDP at version DS 7.4.0. If your data was processed with DS 7.4.0 or later, acis_detect_afterglow was not run in the pipeline. The new hot pixel tools are more judicious with respect to throwing away piled source events. Users should be able to analyze the data without having to unfilter events first.

The Get Started section of the Remove the acis_detect_afterglow correction thread shows how to check the processing version used for your data.

Multi-Component Source Model

Since we have background-subtracted the data (rather than fitting it simultaneously), it is only necessary to create a source model expression. We model this source with a power law (xspowerlaw) absorbed by the interstellar medium (xswabs).

The absorption model will be referred to as "abs1", and the broken power law will be "power1"; the product of the two is assigned as the source model for the data set:

sherpa> set_model(xswabs.abs1*xspowerlaw.power1)

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((20651.1490345 * (xswabs.abs1 * xspowerlaw.power1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed            1            0       100000 10^22 atoms / cm^2
   power1.phoindex  thawed            1           -2            9           
   power1.norm      thawed            1            0        1e+24           

sherpa> guess(get_model())
WARNING: No guess found for apply_rmf(apply_arf((20651.1490345 * (xswabs.abs1 * xspowerlaw.power1))))

sherpa> guess(power1)

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((20651.1490345 * (xswabs.abs1 * xspowerlaw.power1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed            1            0       100000 10^22 atoms / cm^2
   power1.phoindex thawed            1           -2            9           
   power1.norm  thawed  0.000116019  1.16019e-07     0.116019           

Note that since the model contains multiple components, the Sherpa guess() command cannot be used to estimate the initial parameter values for the full model expression based on the data; each model component must be considered separately.

If you know the Galactic column density for your source, you may set abs1.nh and freeze it:

sherpa> abs1.nh = <value>
sherpa> freeze(abs1.nh) 

In this thread, however, we allow it to vary during the fit.

sherpa> thaw(abs1.nh)
 
sherpa> show_model()
Model: 1
apply_rmf(apply_arf((20651.1490345 * (xswabs.abs1 * xspowerlaw.power1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed            1            0       100000 10^22 atoms / cm^2
   power1.phoindex thawed            1           -2            9           
   power1.norm  thawed  0.000116019  1.16019e-07     0.116019           

Pileup Model

It is also necessary to define the pileup model (jdpileup) immediately after defining the source model:

sherpa> set_pileup_model(jdpileup.jdp)

sherpa> print(get_pileup_model())
jdpileup.jdp
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   jdp.alpha    thawed          0.5            0            1           
   jdp.g0       frozen            1  1.17549e-38            1           
   jdp.f        thawed         0.95          0.9            1           
   jdp.n        frozen            1  1.17549e-38          100           
   jdp.ftime    frozen        3.241  1.17549e-38            5        sec
   jdp.fracexp  frozen        0.987            0            1           
   jdp.nterms   frozen           30            1          100                         

Read carefully the following information on how to set the pileup model parameters:

• alpha

  alpha parameterizes "grade migration" in the detector, and represents the probability, per photon count greater than one, that the piled event is not rejected by the spacecraft software as a "bad event". Specifically, if n photons are piled together in a single frame, the probability of them being retained (as a single photon event with their summed energy) is given by alpha^(n-1).

  alpha is the parameter most likely to be allowed to vary in a fit.

• g0:

  Grade correction for single photon detection, i.e. a fraction g0 of single photon events will be retained as good grades.

  In practice, this should be frozen to unity (the default) in any fit.

• f:

  The fraction of events in the source extraction region to which pileup will be applied. A typical value is approximately 0.95, and it should always be > 0.85. f is the second most likely parameter, after alpha, to be allowed to vary in a fit.

  It is recommended that the lower limit of this parameter be set to 0.85 before fitting (the default is 0.9):

  sherpa> jdp.f.min=0.85
  

• n:

  Divide the model counts among n regions, to which the pileup model will be applied independently. This should be approximately the number of 3x3 pixel islands in the extracted spectra; for point sources, it is unity (the default).

  This should remain a frozen parameter in any fit.

• ftime

  The frame time parameter (ftime) should be set to the good exposure time per frame, which is recorded in the header keyword EXPTIME of the event file. (Note that EXPTIME is used instead of TIMEDEL because the latter includes the transfer time, which ftime should not.)

  The CIAO command may be executed from within Sherpa, as long as it is escaped with an exclamation point (!):

  sherpa> !dmkeypar acisf01618_evt2.fits EXPTIME echo+
  3.2
  
  sherpa> jdp.ftime=3.2
  

• fracexp

  This parameter is a fraction <= 1 that multiplies the frame exposure time to create a shorter, effective frame exposure time.

  Users should set and freeze fracexp to unity, and set the frame exposure time via the ftime parameter only:

  sherpa> jdp.fracexp=1   
  

  Information on using this parameter to model novel deadtime effects is available in the The Chandra ABC Guide to Pileup.

• nterms:

  The nterms parameter represents the maximum number of photons considered for pileup in a single frame.

  This should be left frozen at a value of 5. In other words, the expansion of the model will include terms corresponding to 0, 1, 2, 3, 4, and 5 photon events landing in the same extraction region during the same frame time.

  sherpa> jdp.nterms=5   
  

The pileup model parameters are now set as shown:

sherpa> print(get_pileup_model())
jdpileup.jdp
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   jdp.alpha    thawed          0.5            0            1           
   jdp.g0       frozen            1  1.17549e-38            1           
   jdp.f        thawed         0.95         0.85            1           
   jdp.n        frozen            1  1.17549e-38          100           
   jdp.ftime    frozen          3.2  1.17549e-38            5        sec
   jdp.fracexp  frozen            1            0            1           
   jdp.nterms   frozen            5            1          100   

Fit using the Neldermead method

We set the optimization method to 'neldermead', the fit statistic to 'chi2datavar', and start the fit:

sherpa> set_method("neldermead")
sherpa> set_stat("chi2datavar")

sherpa> fit()
 Solar Abundance Vector set to angr:  Anders E. & Grevesse N. Geochimica et Cosmochimica Acta 53, 197 (1989)
 Cross Section Table set to bcmc:  Balucinska-Church and McCammon, 1998
Dataset               = 1
Method                = neldermead
Statistic             = chi2datavar
Initial fit statistic = 1727.08
Final fit statistic   = 77.7473 at function evaluation 828
Data points           = 42
Degrees of freedom    = 37
Probability [Q-value] = 0.000101442
Reduced statistic     = 2.10128
Change in statistic   = 1649.34
   jdp.alpha      0.432791   
   jdp.f          0.928049     
   abs1.nh        5.94924    
   power1.phoindex   1.35462     
   power1.norm    0.00195233  

The parameter space of models which include a pileup model component is very complex and presents a challenge to the search and fitting algorithm. By using the Neldermead method, we are settling for finding a local minimum; now one should examine the solution, check the error bars, find a new minimum, refit, etc. Another option is to use a different optimization method, one which will start from several initial parameter settings and explore in a more consistent way the entire parameter space; the "moncar" method that uses random starting locations can be used in Sherpa.

Fit using the Monte Carlo method

Now we change the method to Monte Carlo and re-run the fit to see how the results change. Note that we start at the results obtained by Neldermead, but you could reset the model parameters before continuing with the reset() command, if desired.

sherpa> set_method("moncar")
sherpa> fit()
Dataset               = 1
Method                = moncar
Statistic             = chi2datavar
Initial fit statistic = 77.7473
Final fit statistic   = 77.7473 at function evaluation 9237
Data points           = 42
Degrees of freedom    = 37
Probability [Q-value] = 0.000101442
Reduced statistic     = 2.10128
Change in statistic   = 1.44529e-06
   jdp.alpha      0.432506    
   jdp.f          0.92814     
   abs1.nh        5.94892     
   power1.phoindex   1.35457     
   power1.norm    0.00195322  

This table compares the results from the two runs with different optimization methods:

Calculate the parameter uncertainties

At this point, we continue with the first set of fit results obtained by using the Neldermead optimization method. For the best-fit values to really mean something, we need to find the uncertainties on the parameters. The covariance and projection commands can be used to estimate confidence intervals for the thawed parameters:

sherpa> set_method("neldermead")
sherpa> set_proj_opt("fast", False)
sherpa> proj()
WARNING: hard maximum hit for parameter jdp.alpha
Dataset               = 1
Confidence Method     = projection
Fitting Method        = neldermead
Statistic             = chi2datavar
projection 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   jdp.alpha        0.432791    -0.289777        -----
   jdp.f            0.928049    -0.303487    0.0316072
   abs1.nh           5.94924     -0.56033     0.588849
   power1.phoindex      1.35462    -0.234803     0.239429
   power1.norm    0.00195233  -0.00113728    0.0159352

The table lists the best-fit values for the model parameters, as well as the upper and lower bounds on each of them; 'jdp.alpha' does not have an upper bound, so it is unconstrained towards the hard maximum. Notice that the proj() command must be preceded by 'set_proj_opt('fast', False)' to use the current optimization method (Neldermead) instead of LevMar, the default method used by projection.

Finally, plot the fit and the residuals of the fit:

sherpa> plot_fit_delchi()

The final plot is show in Figure 3.

Figure 3: Plot of the fit and residuals

Plot of the model fit to the source spectrum, with residuals.