Marx 5.0.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. It 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: 770961
  Total Number Detected: 482771
  Efficiency: 6.261938e-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);
    xnotice_en (1, 0.3, 11.0);
    fit_fun("phabs(1)*powerlaw(1)");
    set_par (1, 1, 1);		       %  freeze nH
    fit_counts;
    list_par;

The fit produces a rather large chi-square per dof of more than 7.5:

 Parameters[Variable] = 3[2]
            Data bins = 659
           Chi-square = 4946.838
   Reduced chi-square = 7.529434
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.0004742587           0       1e+10
  3  powerlaw(1).PhoIndex   0     0         1.593586          -2           9

Note also that the parameters are different from the power-law parameters that went into the simulation. For example, the normalization is less than half of what was expected, and the powerlaw index is somewhat low compared to the expected value of 1.8. In fact, the conf function shows that the 99 percent confidence interval on the powerlaw index is from 1.59 to 1.62.

Suspecting that this observation suffers from pileup, we enable the isis pileup kernel, which introduces a few additional parameters:

    set_kernel (1, "pileup");
    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.0004742587           0       1e+10
  3  powerlaw(1).PhoIndex   0     0         1.593586          -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

Given the fact that the PSF fraction varies with off-axis angle and spectral shape, we will allow it to vary during the fit:

    thaw (7);

As before, the fit_counts command may be used to compute the best fit parameters. However, the parameter space in the context of the pileup kernel is very complex with many local minima, even for a function as simple as an absorbed powerlaw. For this reason it is strongly recommended that pileup fits be performed many times with different initial values for the parameters to better explore the parameter space. The easiest way to do this in isis is to use the fit_search function, which randomly samples the parameter space by choosing parameter values uniformly distributed between their minimum and maximum parameter values. Before starting, let's set the powerlaw normalization's maximum value to something reasonable:

    set_par (2; max=0.01);
    list_par;

The above produces:

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.0004742587           0        0.01
  3  powerlaw(1).PhoIndex   0     0         1.593586          -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     0             0.95           0           1

The next step is to use fit_search to perform the minimization. Here we will explore the parameter space by having fit_search run the fit_counts function at 100 randomly sampled locations. It should be noted that isis will run a number of the fits in parallel by distributing the computations across the available CPU cores.

    s = fit_search (100, &fit_counts);
    eval_counts;
    list_par;

This gives:

 Parameters[Variable] = 7[4]
            Data bins = 659
           Chi-square = 751.8835
   Reduced chi-square = 1.147914
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.0007544917           0        0.01
  3  powerlaw(1).PhoIndex   0     0         1.790701          -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.4686956           0           1
  7  pileup(1).psffrac      0     0         0.775919           0           1

We see that the reduced chi-square is near what one would expect for a good fit and that the powerlaw index is very close to the expected value. We can use the conf function to obtain its confidence interval:

    conf(3);

This indicates that the 90 percent confidence interval on the powerlaw index to be between 1.78 and 1.80.

The powerlaw normalization may appears to be a bit lower than expected. For a near on-axis point source, the 2 arc-second radius extraction region used contains roughly 90 percent of the flux. Hence the expected powerlaw normalization should be somewhere in the neighborhood of 0.0009. However, it is much less constrained, as can be seen by computing its 90 percent confidence interval, which runs from from 0.006 to 0.003.

Here is a plot of the resulting fit, as produced by the rplot_counts function:

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