Marx 4.3.0

Simulating Pileup with marx

The purpose of this example is to show how to to use the marxpileup program to simulate the effects of pileup. marxpileup program is a post-processor that performs a frame-by-frame analysis of an existing marx simulation. Here the simulated data from the powerlaw simulation will be used. The spectrum used for that simulation had a somewhat large normalization to ensure that pileup would occur. If you have not already done so, then you should review that example first.

Assuming that the simulation to be processed is in the plaw/ subdirectory, marxpileup is run using

    marxpileup MarxOutputDir=plaw

This will create a directory called plaw/pileup/ and place the pileup results there. It also writes a brief summary to the display resembling

  Total Number Input: 8230
  Total Number Detected: 5047
  Efficiency: 6.132442e-01

This shows that because of pileup, nearly 40 percent of the events were lost.

The next step is to run marx2fits to create a corresponding level-2 file called plaw_pileup_evt2.fits.

    marx2fits --pileup plaw/pileup plaw_pileup_evt2.fits

Since this is a pileup simulation, the --pileup flag was passed to marx2fits.

As in the powerlaw example, CIAO tools may be used to produce a PHA file, and ARF, and an RMF for subsequent analysis. The same ARF and RMF that was used for the powerlaw simulation can be used here. However, it is necessary to create a new PHA file from the plaw_pileup_evt2.fits event file. For analyzing a piled observation of a near on-axis point source, it is recommended that the extraction region have a radius of 4 ACIS tangent plane pixels. A Bourne shell script that runs dmextract to produce the PHA file plaw_pileup_pha.fits may be found here.

Analyzing the simulated data

As before, isis will be used to analyze the piled spectrum.

    load_dataset ("plaw_pileup_pha.fits", "plaw_rmf.fits", "plaw_arf.fits");
    rebin_data (1, 30);
    fit_fun("phabs(1)*powerlaw(1)");
    set_par (1, 1, 1);
    fit_counts;
    list_par;
    set_kernel (1, "pileup");
    set_par (7, 0.9);
    fit_counts;
    list_par;
    group_data (1,4);
    eval_counts;
    plot_open("plawpileupfit.ps/CPS");
    rplot_counts (1);
    plot_close ();
    exit (0);

The fit produces a very reasonable chi-square per dof of about 1.1. The best fit parameter values may be displayed using the list_par command:

phabs(1)*powerlaw(1)
 idx  param             tie-to  freeze         value         min         max
  1  phabs(1).nH            0     1                1           0      100000  10^22
  2  powerlaw(1).norm       0     0     0.0004284556           0       1e+10
  3  powerlaw(1).PhoIndex   0     0         1.531699          -2           9

Note that while the chi-square is good, the parameters are different from the power-law parameters that went into the simulation. For example, the normalization is about half that expected, and the powerlaw index is somewhat low compared to the expected value of 1.8. In fact, the conf function shows that the upper 90 percent confidence limit on the index is less than 1.6.

Suspecting that this observation suffers from pileup, we enable the isis pileup kernel:

isis> set_kernel(1, "pileup");

Now list_par produces

isis> list_par;
phabs(1)*powerlaw(1)
 idx  param             tie-to  freeze         value         min         max
  1  phabs(1).nH            0     1                1           0      100000  10^22
  2  powerlaw(1).norm       0     0     0.0004284586           0       1e+10
  3  powerlaw(1).PhoIndex   0     0         1.531706          -2           9
  4  pileup(1).nregions     0     1                1           1          10
  5  pileup(1).g0           0     1                1           0           1
  6  pileup(1).alpha        0     0              0.5           0           1
  7  pileup(1).psffrac      0     1             0.95           0           1

The value of the psfrac parameter is a bit high is set to 0.90 using:

isis> set_par (7, 0.9);

Now the fit_counts command may be used to compute the fit using the pileup kernel. Keep in mind that the parameter space in the context of the pileup kernel is very complex with many local minima. For this reason it is strongly recommended that pileup fits be performed many times with different initial values for the parameters.

isis> fit_counts;
 Parameters[Variable] = 7[3]
            Data bins = 123
           Chi-square = 92.95992
   Reduced chi-square = 0.774666
0
isis> list_par;
phabs(1)*powerlaw(1)
 idx  param             tie-to  freeze         value         min         max
  1  phabs(1).nH            0     1                1           0      100000  10^22
  2  powerlaw(1).norm       0     0     0.0009427394           0       1e+10
  3  powerlaw(1).PhoIndex   0     0         1.703912          -2           9
  4  pileup(1).nregions     0     1                1           1          10
  5  pileup(1).g0           0     1                1           0           1
  6  pileup(1).alpha        0     0        0.2173629           0           1
  7  pileup(1).psffrac      0     1              0.9           0           1

Under the pileup kernel the powerlaw normalization is very close to what was expected. In fact, for an on-axis point source, the extraction region used contains about 95 percent of the flux; hence the expect powerlaw normalization is about 0.95*0.001. The powerlaw index is a bit lower than expected, but much closer than that predicted by the standard linear kernel. The upper 90 percent confidence interval on the index is 1.78.

Here is a plot of the resulting fit:

This page was last updated Feb 21, 2008 by John E. Davis.
Technical questions should be addressed to marx-help at space mit edu.