Simulating an LETG/ACIS-S observation

In this example, we simulate an LETG/ACIS-S observation with marx. The primary goal of this exercise is to see how well marx can simulate this instrument combination by performing a spectral analysis of the result. For simplicity, we assume an on-axis point source and an absorbed powerlaw spectrum with a column density of 10^19 atoms/cm^2, a spectral index of 1.8. Rather than specifying an exposure time for the observation as in the other examples, here we specify the desired number of detected events. In this case, the simulation will be run until 1,000,000 (1e6) events have been detected.

Creating an input spectrum from a Sherpa model

The first step is to create a 2-column text file that tabulates the absorbed powerlaw 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. For a change, we will use the Sherpa program in this example.

We chose the specific physical model (a positive powerlaw with an unrealistically high flux) because we want to construct a really detailed picture of the features seen in the PSF and we need a large number of photons over a wide range of energies:

import numpy as np

# set source properties
set_source(xsphabs.a * xspowerlaw.p)
a.nH = 0.001
p.PhoIndex = 1.8
p.norm = 0.1

# get source
my_src = get_source()

# set energy grid
bin_width = 0.01
energies = np.arange(0.03, 12., bin_width)

# evaluate source on energy grid
flux = my_src(energies)

# Sherpa uses the convention that an energy array holds the LOWER end of the bins;
# Marx that it holds the UPPER end of a bin.
# Thus, we need to shift energies and flux by one bin.
# Also, we need to divide the flux by the bin width to obtain the flux density.
save_arrays("source_flux.tbl", [energies[1:], flux[:-1] / bin_width], ["keV","photons/s/cm**2/keV"], ascii=True)

More details about the format of the marx input spectrum can be found at SpectrumFile. Note, that the parameter SourceFlux sets the normalization of the flux; if the normalization of the model file should be used, set SourceFlux=-1.

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="source_flux.tbl" \
       NumRays=-1000000 ExposureTime=0 TStart=2012.5 \
       OutputDir=letgplaw GratingType="LETG" DetectorType="ACIS-S" \
       DitherModel="INTERNAL" RA_Nom=30 Dec_Nom=40 Roll_Nom=50 \
       SourceRA=30 SourceDEC=40 \
       Verbose=yes mode=h

Note the use of a negative value of the NumRays parameter. This tells marx that the simulation is to continue until the absolute value of that number of rays have been detected. 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. The results of the simulation will be written to a subdirectory called letgplaw, as specified by the OutputDir parameter. After the simulation has completed, a standard Chandra event file may be created using the marx2fits program:

marx2fits letgplaw letgplaw_evt2.fits

The fits file letgplaw_evt1.fits can be further processed with standard CIAO tools. As some of these tools require the aspect history, the marxasp program is used to create an aspect solution file that matches the simulation:

marxasp MarxDir="letgplaw" OutputFile="letgplaw_asol1.fits"

It is interesting to look at the event file with a viewer such as ds9. Here we use dmimg2jpg, which displays an RGB image of the events projected to the sky plane.

RGB image of the event file

The energy bands were color-coded as follows: red: 0.0 to 0.5 keV, green: 0.5 to 2.0 keV, blue: 2.0 to 12 keV.

This image RGB image of the event file shows what the pattern of events looks like near zeroth order. There are a number of features in this image that are worth pointing out. The thick blue line in the middle of the fan-like structure that extends from the top left of the figure to the bottom right corresponds to diffracted photons from the primary grating bars of the LEG. The events in this line are the subject of the spatial extraction step when creating the event histograms for spectral analysis, as described below. The other blue lines in the fan-like structure correspond to photons that have also been diffracted by the fine support structure of the LEG. This fine support structure forms a grating that diffracts in a direction orthogonal to the primary diffraction direction. The events in the multicolored line that runs from the bottom left to the upper right consists of two contributions. The first is from photons that arrived at the detector while the image was being read out between frames. The second is from the undiffracted photons from the primary grating that have been diffracted from the LEG fine support structure. The red star-like pattern in the center of the figure was produced by photons that have diffracted from the LEG’s coarse support structure, whose 2 mm period diffracts only the very lowest energy photons through an appreciable angle. Finally, if one looks closely enough, the shadowing of the mirror support struts may be seen near zeroth order.

Analyzing the simulated data

The (type-II) PHA file, grating ARFs and RMFs can be produced using the standard CIAO tools using the simulated event file letgplaw_evt2.fits and the aspect solution file letgplaw_asol1.fits as inputs. Here is a Bourne shell script that does this:

#!/bin/sh

X=4096.5
Y=4096.5

Verbose=1

# marx2fits generated file
evt2file="letgplaw_evt2.fits"

# These files are generated by this script
evt1afile="letgplaw_evt1a.fits"
reg1afile="letgplaw_reg1a.fits"
asolfile="letgplaw_asol1.fits"
pha2file="letgplaw_pha2.fits"

tg_create_mask infile="$evt2file" outfile="$reg1afile" \
  use_user_pars=yes last_source_toread=A \
  sA_id=1 sA_zero_x=$X sA_zero_y=$Y grating_obs=header_value \
  clobber=yes verbose=$Verbose

tg_resolve_events infile="$evt2file" outfile="$evt1afile" \
     regionfile="$reg1afile" osipfile="CALDB" acaofffile="$asolfile" \
     verbose=$Verbose clobber=yes

tgextract infile="$evt1afile" outfile="$pha2file" tg_order_list="-1,+1" \
  ancrfile=none respfile=NONE inregion_file=none clobber=yes \
  tg_srcid_list=1 outfile_type=pha_typeII tg_part_list='header_value'

mkgrmf grating_arm="LEG" order=-1 outfile="leg-1_rmf.fits" srcid=1 detsubsys="ACIS-S3" \
  threshold=1e-06 obsfile="$pha2file" regionfile="$pha2file" \
  wvgrid_arf="compute" wvgrid_chan="compute" clobber=yes

fullgarf phafile="$pha2file"  pharow=1 evtfile="$evt2file" \
  asol="$asolfile" engrid="grid(leg-1_rmf.fits[cols ENERG_LO,ENERG_HI])" \
  maskfile=NONE dafile=NONE dtffile=NONE badpix=NONE maskfile=NONE rootname="leg" ardlibqual=";UNIFORM;bpmask=0"

mkgrmf grating_arm="LEG" order=+1 outfile="leg+1_rmf.fits" srcid=1 detsubsys="ACIS-S3" \
   threshold=1e-06 obsfile="$pha2file" regionfile="$pha2file" \
   wvgrid_arf="compute" wvgrid_chan="compute" clobber=yes

fullgarf phafile="$pha2file"  pharow=2 evtfile="$evt2file" \
   asol="$asolfile" engrid="grid(leg+1_rmf.fits[cols ENERG_LO,ENERG_HI])" \
   maskfile=NONE dafile=NONE dtffile=NONE badpix=NONE maskfile=NONE rootname="leg" ardlibqual=";UNIFORM;bpmask=0"

In particular, this script constructs an order-sorted level 1.5 file from which plus and minus first order events are extracted, and creates the corresponding ARFs and RMFs for the extracted spectra.

For this example, we wish to verify that the marx simulation is consistent with the input spectrum. To this end, we use Sherpa to fit an absorbed powerlaw to the spectrum. The spectral fits to the minus first order histograms can be carried out using the following commands (positive order works similarly):

# Sherpa code

load_pha('letgplaw_pha2.fits')
load_arf('legLEG_-1_garf.fits')
load_rmf('leg-1_rmf.fits')
set_source(xsphabs.a * xspowerlaw.xp)

group_counts(25)
fit()

# Do not display a window
# Delete this line for interactive use!
# set_preferences(['window.display', False])

plot_fit_delchi()
xscale('log')
savefig('m1_fit.png', bbox_inches='tight')
savefig('m1_fit.eps', bbox_inches='tight')

This command sequence produced the following parameter values and a reduced chi-square around 1:

Reduced statistic     = 1.08292
Change in statistic   = 2.84957e+07
   a.nH           0
   xp.PhoIndex    1.75727
   xp.norm        0.10434

The reduced chi-squared value indicates that we found an acceptable fit and all parameter values are close to the original input values (the nH is so small in the input that the fit may give a zero value as fit result). Note that your results may vary slightly if you run this example, because marx is a Monte-Carlo simulation based on random numbers.

Simulated Spectrum in order -1 and a best-fit model.