HETG simulation of an extended source

The purpose of this example, contributed by Dan Dewey, is to show how to use marx to simulate an HETG observation of a point and a simple extended source with a user-specified spectrum. This example, like many of the other examples, uses marx, ISIS and CIAO, which should be installed on your system to run this example.

Create the spectral file

Use MARX to simulate an HETG observation of a powerlaw with two added lines (For example, this could be Fe and Ne flourescence lines). We will use ISIS to create a file plaw_hetg.par which describes this spectral model. This is similar to Simulating a user-defined CCD spectrum with ACIS except that the source is brighter and has a less steep powerlaw.

fit_fun("phabs(1)*powerlaw(1)+gauss(1)+gauss(2)");
list_par;  % This will show the model parameters.

% Set the model parameters:

set_par(1, 0.80);
set_par(2, 0.020);
set_par(3, 1.20);

% Set the Gaussian line parameters.
% It's also possible to set the parameters by name instead of their number.
% Note that the "center" must be specified in Angstroms,
% so keV values are converted using Const_hc [ keV A ].

% Fe line and Ne line centers

set_par("gauss(1).center", Const_hc/6.4038);
set_par("gauss(2).center", Const_hc/0.8486);

% Set the areas and widths of the lines:

set_par("gauss(1).area", 3.0e-4);
set_par("gauss(2).area", 1.0e-4);

% Make the lines have small but non-zero widths
% of 0.020 A FWHM, or sigma = 0.020/2.35 = 0.0085 A.
% These correspond to velocity widths of:
% c * 0.020/1.936 =  ~ 3100 km/s FWHM for Fe line
% c * 0.020/14.61 =  ~  410 km/s FWHM for Ne line

set_par("gauss(1).sigma", 0.0085);
set_par("gauss(2).sigma", 0.0085);

% Save the model/parameters to a file and exit isis.

save_par ("plaw_hetg.par");

Now, use marxflux to evaluate the plaw_hetg.par file to create a marx spectrum file. Here a fine binning is used having bins with dE/E = 0.0003 (v_bin = 90 km/s). This gives a high resolution spectrum across the whole 1 to 40 A (0.31 to 12.4 keV) range suitable for use with the HETG as well as, e.g., future microcalorimeter instruments.

marxflux -e '_A(exp([log(1.0):log(40.0):0.0003]))' plaw_hetg.par plawflux_hetg.tbl

The file plawflux_hetg.tbl can now be used by marx to define the spectrum.

Setup and run marx

Roughly these steps are the same as for the previous examples, except that the HETG is inserted by using GratingType="HETG".

It can be convenient to use pset to set the parameters in a local par file, so that method is demonstrated and used here.

Set up the various marx parameters as desired. Note that you can edit and re-do the following lines for new or modified simulations. As an example, the first simulation here is for a POINT Source; it is followed by a few additional lines to do a second simulation with a DISK Source.

In the working directory paste these sets of lines to the unix prompt to run a first simulation using a point source:

# Make a par file for my marx simulation:
cp marx.par mysim.par

# set the spectrum file to use:
pset mysim.par SpectrumType="FILE"
pset mysim.par SpectrumFile="plawflux_hetg.tbl"
pset mysim.par SourceFlux=-1

# Set other parameters of the simulation:
# Using 50 ks
pset mysim.par ExposureTime=50000
pset mysim.par OutputDir="hetg_plaw"
pset mysim.par DitherModel="INTERNAL"

# Use the HETG with ACIS-S:
pset mysim.par GratingType="HETG"
pset mysim.par DetectorType="ACIS-S"

# Some other parameters it can be useful to set:

# Date of observation (effects ACIS QE)
pset mysim.par TStart=2009.50
# Roll of the observation: 0 puts average dispersion along E -- W.
pset mysim.par Roll_Nom=0.0

# Source RA/DEC (degrees)
pset mysim.par SourceRA=250.000
pset mysim.par SourceDEC=-54.000
# Pointing RA/DEC (degrees)
pset mysim.par RA_Nom=250.000
pset mysim.par Dec_Nom=-54.000

# Finally,
# Run the simulation:
marx @@mysim.par

marx runs, ending with something similar to (since marx is a Monte-Carlo based simulation, the exact number of detected photons can vary):

Writing output to directory 'hetg_plaw' ...
Total photons: 3495031, Total Photons detected: 255172, (efficiency: 0.073010)
  (efficiency this iteration  0.073986)  Total time: 50000.000492

# Create the fits event file and the aspect solution file:
marx2fits hetg_plaw hetg_plaw_evt2.fits
marxasp MarxDir="hetg_plaw" OutputFile="hetg_plaw_asol1.fits"

Now, do another simulation keeping most things the same as above by starting with the mysim.par as it is left from above but changing a few things to use a DISK Source:

cp mysim.par mysim_disk.par
# New output dir for the Disk
pset mysim_disk.par OutputDir="hetg_pldisk"

# Change the SourceType:
pset mysim_disk.par SourceType="DISK"
# a thin disk with average radius ~ 2.0"
pset mysim_disk.par S-DiskTheta0=1.7
pset mysim_disk.par S-DiskTheta1=2.3

# Run the second simulation:
marx @@mysim_disk.par

Next, we create the fits event file and the aspect solution files for both simulations:

marx2fits hetg_plaw hetg_plaw_evt2.fits
marxasp MarxDir="hetg_plaw" OutputFile="hetg_plaw_asol1.fits"

marx2fits hetg_pldisk hetg_pldisk_evt2.fits
marxasp MarxDir="hetg_pldisk" OutputFile="hetg_pldisk_asol1.fits"

Note

Combining marx simulations

The tool marxcat allows simulations to be combined, e.g., we could do the following to make a combination of the point and disk events:

marxcat hetg_plaw hetg_pldisk hetg_plboth

and then create fits and asol files as above. This allows more complex spatial-spectral simulations to be done with marx.

We can look at the simulation output event files with ds9 to check that they are as expected before continuuing with ciao-processing. Figure Event files from both marx simulations shows both simulations.

Event files from both marx simulations

The simulation of the point source is shown on the top, the extended source on the bottom. The extended source has a much wider zeroth order, but the scaling of the image is chosen to bring out the faint features so this is hard to see. Above and below the zeroth order the read-out streak is visible. In both images the X-shape of the grating spectra can be seen. The spectra are much wider in the bottom image due to the source extension. Still, the grating extraction area (green outlines) is large enough to capture most of the signal.

We can also use ds9 to record the center of the disk (simulation 2) in X,Y coordinates (4096.5, 4096.5) for further processing.

Extract HETG spectra

We will extract the HETG spectra and then calculate the response matrix for the positive and negative first order in the MEG grating. There is very little signal in the higher orders, so they would not help to constrain the fit significantly. Extraction of the HEG grating works in a similar way, see the CIAO documentation for details.

#!/bin/sh

# Arguments: root name of files
# Call this with the root names of the files as first argument, e.g.
# > hetg_ciao.sh hetg_plaw
# to process the first example.

root=$1
evtfile="${root}_evt2.fits"
evt1afile="${root}_evt1a.fits"
reg1afile="${root}_reg1a.fits"
asol1file="${root}_asol1.fits"
pha2file="${root}_pha2.fits"

X=4096.5;    # determined by hand in ds9
Y=4096.5;    # determined by hand in ds9

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

tg_resolve_events infile="$evtfile" outfile="$evt1afile" \
  regionfile="$reg1afile" acaofffile="$asol1file" \
  osipfile=CALDB verbose=0 clobber=yes

tgextract infile="$evt1afile" outfile="$pha2file" mode=h clobber=yes tg_order_list="-1,1"

# make response matrices for the MEG +1 and -1 order
mkgrmf grating_arm="MEG" order=-1 outfile="${root}meg-1_rmf.fits" srcid=1 detsubsys=ACIS-S3 \
   threshold=1e-06 obsfile="$pha2file" regionfile="$pha2file" \
   wvgrid_arf="compute" wvgrid_chan="compute"

mkgrmf grating_arm="MEG" order=1 outfile="${root}meg1_rmf.fits" srcid=1 detsubsys=ACIS-S3 \
   threshold=1e-06 obsfile="$pha2file" regionfile="$pha2file" \
   wvgrid_arf="compute" wvgrid_chan="compute"

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

fullgarf phafile="$pha2file"  pharow=4 evtfile="$evtfile" \
   asol="$asol1file" engrid="grid(${root}meg1_rmf.fits[cols ENERG_LO,ENERG_HI])" \
   maskfile=NONE dafile=NONE dtffile=NONE badpix=NONE rootname="$root" ardlibqual=";UNIFORM;bpmask=0"

Perform spectral analysis

As an examples of the spectra analysis, we will fit an absorbed powerlaw to the point-like source spectrum. The spectral regions where the extra lines are located are ignored in the fit. Fitting the spectrum of the extended source works in the similar way. Naturally, a lot more analysis could be done on the simulated spectra, the most obvious point might be a fit of the spectral lines to determine the spectral resolution of Chandra in this example. (The lines from the extended source will appear much wider.) However, that is beyond the scope of this example. Please check the manual of your preferred X-ray spectral fitting software. In this example below, we use ISIS, but XSPEC and Sherpa work very similar.

This is ISIS code to fit an absorbed powerlaw to the simulated spectrum:

% Load the pha's and responses:

()=load_data("hetg_plaw_pha2.fits",[3,4]);

()=load_arf("hetg_plawMEG_-1_garf.fits");
()=load_arf("hetg_plawMEG_1_garf.fits");
assign_arf(1,1);
assign_arf(2,2);

()=load_rmf("hetg_plawmeg-1_rmf.fits");
()=load_rmf("hetg_plawmeg1_rmf.fits");
assign_rmf(1,1);
assign_rmf(2,2);

% Group and notice the data:
group_data(all_data, 8);
% 0.6 to 7.5 keV:
xnotice(all_data, Const_hc/7.5, Const_hc/0.6);

% List the data
list_data;

% open plot
plot_open("hetgmeg1.ps/CPS");

% Use keV for the X axis in plots:
plot_unit("keV");

% Use log x axis, linear y axis starting at 0.
xlog; xrange(0.5,10.0);
ylin; yrange(0.,);

% Can look at the data
title("MARX simulation of HETG observation");
plot_data_counts(1);
% show all of these in grey
oplot_data_counts(1,14);
oplot_data_counts(2,14);

% Ignore the line regions for the continuum fitting we'll do
%   Fe region
ignore(all_data, 0.95*Const_hc/6.4, 1.05*Const_hc/6.4);
%   Ne region
ignore(all_data, 0.97*Const_hc/0.8486, 1.03*Const_hc/0.8486);
%
% Replot the spectra again - use some colors:
oplot_data_counts(1,2);   %  MEG -1  red
oplot_data_counts(2,8);   %  MEG -1  orange

% Define the fitting function:
fit_fun("phabs(1)*powerlaw(1)");

% Set initial values and limits for the parameters;
% keep them all "thawed" with 0 as third argument.
set_par (1, 1.00,  0,  0.1,   10.0);
set_par (2, 0.01, 0,  1.e-4, 1.0);
set_par (3, 1.00,  0,  0.5,   2.5);

% Use Gehrels to help with low-counts bins
set_fit_statistic ("chisqr;sigma=gehrels");

set_fit_method ("lmdif;default");
()=fit_counts;
list_par;

% overplot one of the model fits
oplot_model_counts(1, 4);

plot_close();

The resulting fit parameters are similar, but not identical, to the parameters we put into the simulation above.

Simulated spectra for MEG order +1 and -1

The MEG orders are shown as the red and orange line. They differ from one another because the effective area is different, e.g. the gaps between the ACIS chips fall on different wavelengths. The blue line is the best fit model for the spectrum shown in red.