Simulating a user-defined CCD spectrum with ACIS¶
The purpose of this example is to show how to use marx to simulate an ACIS observation of a point source with a user-specified spectrum. For simplicity, suppose that we wish to simulate a 3000 ksec observation of an on-axis point source whose spectrum is represented by an absorbed powerlaw, with a spectral index of 1.8, and a column density of 10^22 atoms/cm^2. The normalization of the powerlaw will be set to 0.001 photons/keV/cm^2/s at 1 keV. The large exposure time was chosen to illustrate the consistency of the marx ray trace with that of the underlying calibration data.
Creating the spectral file¶
The first step is to create a 2-column text file that tabulates the flux [photons/sec/keV/cm^2] (second column) as a function of energy [keV] (first column). The easiest way to create such a file is to make use of a spectral modeling program such as ISIS, Sherpa or XSpec. The rest of this tutorial is given in the context of ISIS.
In ISIS the absorbed powerlaw model is specified using:
fit_fun("phabs(1)*powerlaw(1)");
set_par(1, 1);
set_par(2, 0.001);
set_par(3, 1.8);
save_par ("plaw.p");
list_par;
The model parameters lists by ISIS should look like this:
phabs(1)*powerlaw (1)
idx param tie-to freeze value min max
1 phabs(1).nH 0 0 1 0 100000 10^22
2 powerlaw(1).norm 0 0 0.001 0 0.01
3 powerlaw(1).PhoIndex 0 0 1.8 1 3
The next step is to convert the parameter file plaw.p to the
spectrum file that marx expects. The marxflux script may be used to
create a file called plawflux.tbl in the appropriate format via:
marxflux -e '[0.3:12.0:0.003]' plaw.p plawflux.tbl
This script requires ISIS to be installed and linked to at least
version 2.1 of the S-Lang library. (For XSpec users, marx is distributed with a
script called xspec2marx that may be used to create such a file,
and Sherpa users find instructions here:
http://cxc.harvard.edu/sherpa/threads/marx/ ).
The plawflux.tbl file is input to marx using the following
marx.par parameters:
SpectrumType=FILE
SpectrumFile=plawflux.tbl
SourceFlux=-1
The SpectrumType parameter is set to FILE to indicate that
marx is to read the spectrum from the file specified by the
SpectrumFile parameter. The SourceFlux parameter may be
used to indicate the integrated flux of the spectrum. The value of -1
as given above means that the integrated flux is to be taken from the
file.
Running marx¶
The next step is to run marx in the desired configuration. Some
prefer to use tools such as pset to update the marx.par
file and then run marx. Here, the parameters will be explicitly
passed to marx via the command line:
marx SourceFlux=-1 SpectrumType="FILE" SpectrumFile="plawflux.tbl" \
ExposureTime=3000000 TStart=2012.5 \
OutputDir=plaw GratingType="NONE" DetectorType="ACIS-S" \
DitherModel="INTERNAL" RA_Nom=30 Dec_Nom=40 Roll_Nom=50 \
SourceRA=30 SourceDEC=40 \
Verbose=1 mode=h
This will run the simulation and place the results in a subdirectory
called plaw. The results may be converted to a standard Chandra
level-2 fits file by the marx2fits program:
marx2fits plaw plaw_evt2.fits
The resulting fits file plaw_evt2.fits may be further processed
with standard CIAO tools. As some of these tools require the aspect
history, the marxasp program will be used to create an aspect
solution file that matches the simulation:
marxasp MarxDir="plaw" OutputFile="plaw_asol1.fits"
Creating CIAO-based ARFs and RMFs for MARX Simulations¶
Armed with the simulated event file plaw_evt2.fits and the aspect
solution file plaw_asol1.fits, a PHA file, ARF and RMF may be
made using the standard CIAO tools. First, extract the PHA data. The asphist CIAO
tool may be used to create an aspect histogram from an aspect solution
file, as seen in this Bourne shell script:
# Extraction region parameters: circle of radius R at (X,Y)
X=4096.5
Y=4096.5
R=20.0;
phagrid="pi=1:1024:1"
evtfile="plaw_evt2.fits"
asolfile="plaw_asol1.fits"
phafile="plaw_pha.fits"
rmffile="plaw_rmf.fits"
arffile="plaw_arf.fits"
asphistfile="plaw_asp.fits"
asphist infile="$asolfile" outfile="$asphistfile" evtfile="$evtfile" \
verbose=1 mode=h clobber=yes
dmextract infile="$evtfile[sky=Circle($X,$Y,$R)][bin $phagrid]" \
outfile="$phafile" verbose=1 mode=h clobber=yes
While marx strives to accurately model of the Chandra Observatory, there are some differences that need to be taken into account when processing marx generated data with CIAO. As described in Current caveats for MARX page, marx does not incorporate the non-uniformity maps for the ACIS detector, nor does it employ the spatially varying CTI-corrected pha redistribution matrices (RMFs). As such there will be systematic differences between the marx ACIS effective area and that of the mkarf CIAO tool when run in its default mode. Similarly, the mapping from energy to pha by marx will be different from that predicted by mkacisrmf.
Creating an ARF to match a marx simulation¶
As mentioned above, marx does not implement the ACIS QE uniformity maps. The following commands show how to set the mkarf parameters to produce an ARF that is consistent with the marx effective area (this continues the Bourne shell script from the last paragraph):
# Select a CCD.
ccdid=7;
detname="ACIS-$ccdid;UNIFORM;bpmask=0"
grating="NONE"
# For ACIS-I, use engrid="0.3:11.0:0.003". This reflects a limitation of mkrmf.
engrid="0.3:12.0:0.003"
mkarf mirror="hrma" detsubsys="$detname" grating="$grating" \
outfile="$arffile" \
obsfile="$evtfile" engrid="$engrid" asphistfile="$asphistfile" \
sourcepixelx=$X sourcepixely=$Y \
maskfile=NONE pbkfile=NONE dafile=NONE \
verbose=1 mode=h clobber=yes
Notice the settings detsubsys="ACIS-7;uniform;bpmask=0" and maskfile=NONE
pbkfile=NONE dafile=NONE that are different from the way a normal Chandra
observation would be processed.
Creating an RMF to match a marx simulation¶
marx maps energies to phas using the FEF gaussian parametrization utilized by the mkrmf CIAO tool. The newer mkacisrmf tool uses a more complicated convolution model that does not appear to permit a fast, memory-efficient random number generation algorithm that marx would require. In contrast, a gaussian-distributed random number generator is all that is required to produce pha values that are consistent with mkrmf generated responses.
The Chandra CALDB includes several FEF files. The one that
marx currently employs is
acisD2000-01-29fef_phaN0005.fits, and is located in
the $CALDB/data/chandra/acis/cpf/fefs/ directory. This file must
be specified as the infile mkrmf parameter.
For a point source simulation, look at the marx2fits generated event
file and find the average chip coordinates. The CIAO dmstat tool
may be used for this purpose. The CCD_ID and the mean chip
coordinates are important for the creation of the filter that will be
passed to mkrmf to select the appropriate FEF tile. For simplicity
suppose that the mean point source detector location is at (308,494)
on ACIS-7. Then run mkrmf using (continuing the shell script form above):
# Mean chip position
chipx=221;
chipy=532;
fef="$CALDB/data/chandra/acis/fef_pha/acisD2000-01-29fef_phaN0005.fits"
cxfilter="chipx_hi>=$chipx,chipx_lo<=$chipx"
cyfilter="chipy_hi>=$chipy,chipy_lo<=$chipy"
mkrmf infile="$fef[ccd_id=$ccdid,$cxfilter,$cyfilter]" \
outfile="$rmffile" axis1="energy=$engrid" \
axis2="$phagrid" thresh=1e-8 mode=h clobber=yes verbose=1
Notice, how the FEF file is set explicitly to
fef="$CALDB/data/chandra/acis/cpf/fefs/acisD2000-01-29fef_phaN0005.fits"
and that the correct tile of that FEF file is chosen with
fef[ccd_id=7,chipx_hi>=221,chipx_lo<=221,chipy_hi>=532,chipy_lo<=532]
See the CIAO threads for other ways of running mkrmf. The important thing is to specify the correct FEF and tile.
Analyzing the simulated data¶
The PHA, ARF, and RMF files may be used in a spectral modelling program such as ISIS to see whether or not one can reach the desired science goal from the simulated observation. For this example, the goal is to verify that the marx simulation is consistent with the input spectrum. To this end, ISIS will be used to fit an absorbed powerlaw to the pha spectrum. The figure below showing the resulting fit was created via the following ISIS script:
load_dataset ("plaw_pha.fits", "plaw_rmf.fits", "plaw_arf.fits");
rebin_data (1, 30);
xnotice_en (1, 0.3, 11);
fit_fun("phabs(1)*powerlaw(1)");
set_par(1,1,1);
fit_counts;
list_par;
plot_open("plawfit.ps/CPS");
xrange(0.3,11);
rplot_counts (1);
plot_close ();
This script produced the following parameter values with a reduced chi-square of 1.2:
Parameters[Variable] = 3[2]
Data bins = 600
Chi-square = 722.847
Reduced chi-square = 1.208774
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.001003397 0 1e+10
3 powerlaw(1).PhoIndex 0 0 1.826904 -2 9
The rplot_counts may be used to produce a plot of the resulting
fit.
The residuals show that the model is systematically high for some energies. The reason for this can be traced back to energy-dependent scattering where photons are scattered outside the extraction region. The CIAO effective area does not include this loss factor, and as a result, this omission appears in the residuals. This effect is apparant because of the enormous number of counts in this simulation.
Simulating pile-up in an ACIS CCD spectrum¶
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.
We continue the example from above, where the spectrum had a
somewhat large normalization to ensure that pileup would occur.
Creating an event file with pile-up¶
marxpileup can be run on the eventfile generated by marx 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.
In this case we can use the same ARF and RMF as above. 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 is pileup_ciao.sh.
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;
which gives:
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:
Plot of the pileup spectrum¶
