Replicating a Chandra Observation / Using ChaRT or SAOTrace

The example presented here was originally designed to see how well the marx PSF compares to that of an actual Chandra observation of a far off-axis source. We will also be comparing the marx intrinsic PSF to that of SAOTrace (SAOTrace is a high-fidelity ray-tracer for the Chandra mirrors that simulates many details that are treated statistically in marx to improve efficiency). Hence, this example will serve several purposes:

  • How to simulate an existing observation.

  • How to use SAOTrace dithered rays with marx.

  • To see how the SAOTrace and marx PSFs compare to the observed PSF for a far off-axis source.

In this example, we will use the 1999 Chandra calibration observation of LMC X-1 observed about 25 arc-minutes off-axis with the ACIS detector (ObsID 1068). We download the data with download_chandra_obsid and reprocess it with chandra_repro.

The figure shows the source as seen in ds9:

The Chandra PSF at this position is an ellipse with spokes.

ds9 image of the PSF

The crosshairs in ds9 image of the PSF mark the position where the support strut shadows meet as determined by eye: ra=84.91093 degrees, dec=-69.74348 degrees, (x,y)=(5256,6890). This coordinate will be used as the position of the source. For this observation, this coordinate falls on ACIS-4 (S0).

Creating the source spectral file

For an accurate simulation of an observed source, it is important to make the detector configuration match that of the observation as closely as possible. And since the PSF is energy dependent, a rough idea of the energy dependence of the incident source flux must be obtained. There are several approaches to get this input spectrum. One way would be to extract the spectrum, fit a model using XSPEC, ISIS, or Sherpa and write an input file from the model as we did in several previous examples.

Just to illustrate another approach, we will do something different here. (But do not worry if that looks like black magic to you; you can download the input spectrum file below or select another method to get the input spectrum. We just think it is instructive to show something slightly different in each example.)

Here, we will simply extract the counts in a region containing the observed PSF and divide it by the effective area. This is known as “flux-correction”, and involves no spectral fitting. Strictly speaking, the validity of this technique assumes that the spectrum does not vary much over the scale of the RMF.

The first step is to create the ARF to be used for the flux-correction. The creation of the ARF is straightforward via this Bourne shell script:




asphist infile="$asolfile" outfile="$aspfile" evtfile="$evtfile" clobber=yes

# It would be more accuarte to use mkwarf.
# However, that's slower and for our purposes an approximate solution is good enough.
mkarf detsubsys="ACIS-S0" grating="NONE" outfile="$arffile" \
  obsfile="$evtfile" asphistfile="$aspfile" \
  sourcepixelx=$X sourcepixely=$Y engrid="0.3:8.0:0.1" \
  maskfile=NONE pbkfile=NONE dafile=NONE \
  verbose=1 mode=h clobber=yes 

The next step is to extract the counts in the PSF and divide by the ARF to get the flux corrected spectrum. Here we show a way using only standard CIAO tools; there are some comments on what’s happening in the scripts.

First, create a region file containing the source events with dmellipse, and use dmextract to create the PI histogram.

dmcopy "1068/primary/acisf01068N003_evt2.fits.gz[EVENTS][bin x=4950:5530:1,y=6700:7100:1]" image.fits option=image clobber=yes
dmellipse image.fits psf90.reg 0.9 clobber=yes

# We have to use the same energy binning in energy space here that we used for the arf!
# So, first convert the PI to energy (to a precision that's good enough for this example.)
dmtcalc "1068/repro/acisf01068_repro_evt2.fits[EVENTS][sky=region(psf90.reg)]" obs1068_evt2_with_energy.fits  expr="energy=(float)pi*0.0149" clobber=yes
dmextract "obs1068_evt2_with_energy.fits[bin energy=.3:7.999:0.1]" obs1068.spec clobber=yes opt=generic

The energy binning of the ARF and extracted spectrum are the same. Thus, we can now calculate the corrected spectrum by dividing the flux by the effective area and exposure time. There is an extra factor of 0.9, because the ARF assumes that 100% of the PSF are included in the data extraction, while the ellipse above only contains 90% of the counts. Also, we divide by the exposure time (about 1760 s). We make intermediary files to copy columns back and forth and write the spectrum out as an ASCII table in the end.

dmcopy "obs1068.arf[cols energ_lo,energ_hi,specresp]" input_spec.fits clobber=yes
dmpaste input_spec.fits "obs1068.spec[cols counts]" input_spec_1.fits clobber=yes
dmtcalc input_spec_1.fits input_spec_2.fits  expr="flux=counts/((float)specresp * 0.9*1759.8)" clobber=yes
# devide by binwidth to turn the flux int oa flux DENSITY for marx
dmtcalc input_spec_2.fits input_spec_3.fits  expr="fluxdens=flux/0.1" clobber=yes
dmcopy "input_spec_3.fits[cols energ_hi,fluxdens]" "input_spec_marx.tbl[opt kernel=text/simple]" clobber=yes
dmcopy "input_spec_3.fits[cols energ_lo,energ_hi,flux]" "input_spec_saotrace.tbl[opt kernel=text/simple]" clobber=yes

# We now use a combination of some of the most obscure UNIX tools to
# bring the SAOTRACE input spectrum into the right format
# see
# Remove the leading "#" from the line with the column names
awk '{ sub(/\# ENERG_LO/,"ENERG_LO"); print }' < input_spec_saotrace.tbl > input_spec_saotrace.temp
# Then, replace spaces with tabs
tr ' ' \\t < input_spec_saotrace.temp > input_spec_saotrace.rdb

The resulting ASCII tables for marx and for SAOTrace/Chart with the spectrum will be input into both marx and SAOTrace. The spectrum is the same, but the format of the tables is two columns (energy, flux density) for marx and three columns (lower energy, upper energy, flux) for SAOTrace / ChaRT.

Running marx (without SAOTrace)

Running marx for this example does not differ much from any of the previous examples. We use the SpectrumFile parameter to input the source spectrum we estimated above and set the remaining parameters to match the setting of the observation for the pointing direction, exposure time, etc. To get the numbers we display the header (e.g. in ds9) and manually look for the required fits header keywords (e.g. EXPOSURE for ExposureTime, RA_NOM for RA_Nom, etc.). Note the TStart is given in seconds on the spacecraft clock as in the TSTART and TSTOP header keywords.

marx ExposureTime=1737 TStart=51709004.8724941462 \
  GratingType='NONE' DetectorType='ACIS-S' \
  RA_Nom=85.368613 Dec_Nom=-70.12585 Roll_Nom=116.86955 \
  SourceRA=84.91093 SourceDEC=-69.74348 \
  DitherModel='FILE' DitherFile='1068/repro/pcadf051708271N003_asol1.fits' \
  SourceFlux=-1 SpectrumType='FILE' SpectrumFile='input_spec_marx.tbl' \

marx2fits marx_only marx_only.fits

Running SAOTrace / Chart and marx

SAOTrace is a high fidelity ray-trace of the Chandra mirrors, which simulates details of the mirror roughness and alignment that are only approximated in marx because these details are not important for almost all Chandra simulations and they require a much longer run time.

SAOTrace is a separate, stand-alone program that you need to download and install following the instructions in the SAOTrace documentation at Alternatively, you can use ChaRT, which is a web interface to SAOTrace. Using ChaRT saves you the installation, but it is less flexible (e.g. it can only simulate single sources).

Below, we use the installed version of SAOTrace, but to follow this example, you can equally well upload input_spec_saotrace.rdb to ChaRT. SAOTrace itself is a very complex system of many different parts. Here, we will explain the parameters used in this example, but we need to refer to the SAOTrace documentation for details.

We define the source in a file called saotrace_source.lua:

ra_pnt = 85.368613
dec_pnt = -70.12585
roll_pnt = 116.86955

dither_asol_chandra{ file = "1068/repro/pcadf051708271N003_asol1.fits", 
                     ra = ra_pnt, dec = dec_pnt, roll = roll_pnt }

point{ position = { ra = 84.91093,
          dec = -69.74348,
          ra_aimpt = ra_pnt,
          dec_aimpt = dec_pnt,
       spectrum = { { file = "input_spec_saotrace.rdb", 
                      units = "photons/s/cm2", 
                      scale = 1, 
                      format = "rdb",
                      emin = "ENERG_LO",
                      emax = "ENERG_HI",
                      flux = "FLUX"} }

and run SAOTrace with the following command:

trace-nest                    \
   tag=saotrace                  \
   srcpars=saotrace_source.lua   \
   tstart=51709004.8724941462   \
   limit=1759.8607333824  \

This gives us a ray file called saotrace.fits. We take this file as input for marx using SourceType="SAOSAC" and SAOSACFile="saotrace.fits". All the remaining parameters for the exposure time, pointing direction etc. are the same as above:

marx ExposureTime=1737 TStart=51709004.8724941462 \
  GratingType='NONE' DetectorType='ACIS-S' \
  RA_Nom=85.368613 Dec_Nom=-70.12585 Roll_Nom=116.86955 \
  SourceRA=84.91093 SourceDEC=-69.74348 \
  DitherModel='FILE' DitherFile='1068/repro/pcadf051708271N003_asol1.fits' \
  SourceFlux=-1 SourceType='SAOSAC' SAOSACFile='saotrace.fits' \

marx2fits marx_saotrace marx_saotrace.fits

As a last step, we will apply a GTI filter to the output, because the asol file, which contains the “aspect solution” (describing where exactly Chandra points at any one time), can have a short period of time in the beginning, when the telescope is still slewing.

dmcopy "marx_saotrace.fits[@1068/repro/acisf01068_repro_evt2.fits[GTI4]]" marx_saotrace-filtered.fits

Comparing the results

We can now use ds9 to compare the observation with the two simulated event lists:

Three very similar PSFs.

ds9 image of the PSF in the observation (top left), the simulation using only marx (top right), and the simulation using SAOTrace to trace the mirror and marx as the instrument model (bottom left).

The structure of the PSFs is very similar, emphasizing how good both mirror models are. On closer inspection, there is a small shadow just above and to the right of the point where the support strut shadows meet. This feature is a little smaller in marx than in SAOTrace or the real data due to the simplification that the marx mirror model makes. However, very few point sources are observed long enough this far off-axis that these tiny differences actually matter.

So, we can see from this comparison that marx is the tool of choice for almost all Chandra simulations.