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
Note, that the parameter
SourceFlux sets the normalization of the flux; if the
normalization of the model file should be used, set
The next step is to run marx in the desired configuration. Some
prefer to use tools such as
pset to update the
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
This tells marx that the simulation is to continue until the absolute
value of that number of rays have been detected.
SourceFlux parameter may be used to indicate the integrated flux of
the spectrum. The value of
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
specified by the
OutputDir parameter. After the simulation has
completed, a standard Chandra event file may be created using the
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
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.
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
#!/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.