Data Simulation and Analysis Threads
Simulate an EPIC Image and Spectrum with SIXTE
SIXTE is a software package that can be used to simulate images and spectra from an array of X-ray telescopes, including XMM-Newton. In the walk-through below, we assume that it has been successfully installed, the environment variables have been set, and either "sixte-install.csh" or "sixte-install.sh" has been sourced according to the SIXTE documentation. It is also available on the joint NASA-JHU cloud computing service, SciServer.
We will demonstrate how to simulate a point source and an extended source in the PN with the thin filter. The procedures are very similar (indeed, largely identical!) for both, and begin with making a catalog file with SIMPUT. Once the catalog file is made, the event file can be simulated, and with the event file in hand, the simulated image and spectrum can be made.
Point Source
The first step is to make an Xspec command file (denoted by default with .xcm). This will be used by SIMPUT to make a catalog file. For this example, we will model an absorbed power law. In the text below, a pound sign (#) designates comments.
xspec abund wilm xsect vern model tbabs * powerlaw 0.1 1.5 0.001 flux 0.5 10 # we will need the flux for the region we are interested in save model abs_pl.xcm quit y
In addition to the .xcm file, SIMPUT needs the source's right ascension and declination. For this example, we'll set them to 0.0, though when modeling more seriously, it is a good idea to use a source's real RA and Dec to avoid confusion.
simputfile simput=abs_pl.fits src_name=example ra=0.0 dec=0.0 srcflux=7.4478e-12 \ elow=0.1 eup=15 nbins=1000 logegrid=yes emin=0.5 emax=10 xspecfile=abs_pl.xcm
where
-
simput - catalog file (output)
src_name - source name
ra - right ascension
dec - declination
srcflux - source flux (erg s-1 cm-2)
elow - lower bound of the generated spectrum
eup - upper bound of the generated spectrum
nbins - number of energy bins created between elow and eup
logegrid - use a logarithmic energy grid?
emin - minimum energy of flux range
emax - maximum energy of flux range
xspecfile - Xspec .xcm file
Now that we have a catalog file, we can get a simulated event file. We will use an exposure time of 10 ks in the full frame mode with the thin filter. Please note that this step can take several minutes to run. We can safely ignore warnings about keywords it can't find, and the spectrum not covering the full energy range of the ARF.
runsixt xmlfile=/full/path/to/fullframe_thinfilter.xml ra=0.0 dec=0.0 \ simput=abs_pl.fits evtfile=evt_abs_pl.fits exposure=10000
where the keywords are as above, and
-
xmlfile - the appropriate xml file, with its full path
evtfile - output event file
exposure - exposure time (s)
Now that we have an event file, we can make an image and display it with ds9. Note that while the keywords are not case-sensitive, the values they pass are. The keywords in are capitals below since how that is how they are in the header, so it is easier to pick them out by eye. The image is shown in Figure 1 (top).
imgev evtfile=evt_abs_pl.fits image=img_abs_pl.fits coordinatesystem=0 \ projection=TAN NAXIS1=512 NAXIS2=512 CUNIT1=deg CUNIT2=deg CRVAL1=0.0 \ CRVAL2=0.0 CRPIX1=256.5 CRPIX2=256.5 CDELT1=-6.207043e-04 CDELT2=6.207043e-04 \ history=true clobber=yes ds9 img_abs_pl.fits &
where the keywords are as above, and
-
image - output image file
coordinateSystem - coordinate system (0: equatorial, 1: galactic)
projection - projection type (tangential, i.e. "TAN", is recommended for most applications)
NAXIS1 - length of data axis 1
NAXIS2 - length of data axis 2
CUNIT1 - unit of coordinate increment and value
CUNIT2 - unit of coordinate increment and value
CRVAL1 - coordinate value at reference point
CRVAL2 - coordinate value at reference point
CRPIX1 - pixel coordinate of reference point
CRPIX2 - pixel coordinate of reference point
CDELT1 - coordinate increment at reference point
CDELT2 - coordinate increment at reference point
history - include the history with program parameters in the FITS header?
clobber - overwrite existing output file?
|
|
Now we will extract a spectrum and verify the simulated data. Again, it is safe to ignore warnings about keywords it can't find. We will use the eventfilter keyword to select the region with standard C Boolean operators (that is, || corresponds to logical OR, && corresponds to logical AND, and ! corresponds to logical NOT).
makespec evtfile=evt_abs_pl.fits spectrum=abs_pl.pha \ eventfilter="(ra>359.95 || ra<0.05) && (dec>-0.05 && dec<+0.05)" \ rsppath=/full/path/to/sixte/instrument/responses/ clobber=yes
where the keywords are as above, and
-
rsppath - full path to the response file (do not include the name of the response file)
spectrum - output spectrum
Note that region files from ds9 can also be used in the eventfilter parameter. To do that, we must first add X, Y sky coordinates and a WCS to the event file:
radec2xy evtfile=evt_abs_pl.fits projection=TAN refra=0 refdec=0
where the keywords are as above, and
-
refra and refdec - WCS reference points
Next, we make our region file, src.reg, making sure to save it in WCS coordinates, not the ds9 default physical coordinates. Our file looks like this:
$my_computer> more src.reg
# Region file format: DS9 version 4.1
global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 \
highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1
fk5
circle(0.0003167,-0.0003078,120")
And so the call to makespec would be:
makespec evtfile=evt_abs_pl.fits spectrum=abs_pl_reg.pha eventfilter="regfilter(\"src.reg\")" \ rsppath=/full/path/to/instrument/response/file/ clobber=yes
Now that we have a spectrum, we can analyze it in a convenient software package. SIXTE automatically sets the ANCRFILE and RESPFILE keywords in the header. Below, we display and fit the simulated spectrum using sherpa.
ciao
sherpa
set_xsxsect("vern") # use the atomic cross sections of Verner et al. 1996
set_xsabund("wilm") # use the elemental abundances of Wilms et al. 2000
set_stat("cstat") # use an implementation of Cash statistics
set_method("simplex") # use the Nelder-Mead Simplex method
load_pha("abs_pl.pha") # load the spectrum file and any response files and
# background spectrum (if there is one) that may be in the header
ignore() # ignore the full data range
notice("0.5:10") # use only the data between 0.5-10 keV
set_xlog() # set the x-axis to log scale
set_ylog() # ... and do the same for the y-axis
plot_data()
set_source(xstbabs.t1 * xspowerlaw.p1) # define the model
t1.nh=1 # ... and set parameter values
p1.phoindex=2
p1.norm=1
fit()
covar() # find the uncertainties
calc_energy_flux(0.5,10) # calculate the flux in erg s-1 cm-2
calc_data_sum(0.5,10) # find the number of counts between 0.5-10 keV
group_counts(30) # group the counts to a minimum of 30 per bin
plot_fit_resid() # plot the grouped spectrum, fit, and residuals
exit()
The binned spectrum and the fit are shown in Figure 1 (bottom).
Extended Source
The previous example was quite straightforward because it was a point source. But what if we want to simulate an extended source? In that case, in addition to an .xcm file from Xspec, we also need a source image that shows the distribution of the emission. This could be an actual image from an observation, or an image from modeled emission. Here, we will use a Chandra image of N49, a supernova remnant in the LMC. It can be found at the Chandra openFITS repository.
SIMPUT needs to know three things to make the catalog file: the source position, spectral shape, and source flux in the band in question. We can get the source position from the image itself. For the spectral shape and flux, we can use the fitting results of Hirayama et al. (2002):
FX(1-10 keV) = 3.2e-11 erg s-1 cm-2
Γ = 2
NH = 4e21 cm-2
To get the source location, we can display the image in ds9 and search the header for "CRVAL1" and "CRVAL2". Alternatively, we can use the ftool "fkeyprint" to print the parameter values. Note that we want CRVAL1 and CRVAL2, not CRVAL1P and CRVAL2P. Again, the text after the pound sign below is a comment.
fkeyprint n49_0.5-7.0_flux.fits CRVAL1 # we find CRVAL1 = 8.1491223145500E+01 fkeyprint n49_0.5-7.0_flux.fits CRVAL2 # we find CRVAL2 = -6.6080986499400E+01
Setting the normalization to 1e-2 produces a spectrum with 3.3e-11 erg s-1 cm-2 between 1-10 keV, which is reasonably close to Hirayama et al's result. The image covers a very slightly different energy range, 0.5-7 keV; the flux for that range is 3.0e-11 erg s-1 cm-2. Note that, if we were being careful, we would use an elemental abundance table appropriate for this source! For the sake of this example, our usual elemental abundances will suffice. So we now make our .xcm file:
xspec abund wilm xsect vern model phabs * powerlaw 0.4 2 1e-2 flux 0.5 7 save model n49.xcm quit y
We can make our catalog file, using the .xcm file, the CRVAL1 and CRVAL2 values, and the image file, as well as the other parameters as described above.
simputfile simput=n49.fits ra=8.1491223145500E+01 dec=-6.6080986499400E+01 \
srcname=n49 srcflux=3.0e-11 emin=0.5 emax=7 elow=0.3 eup=15. xspecfile=n49.xcm \
imagefile=n49_0.5-7.0_flux.fits clobber=yes
where parameters are as described above, and
-
imagefile - image that shows the emission distribution
The procedure for making an image and spectrum are the same as for a point source, so in the commands below, the keywords are as described above. For our extended source example, the command to simulate an event file for a 10 ks exposure would be:
runsixt xmlfile=/full/path/to/fullframe_thinfilter.xml ra=8.1491223145500E+01 \
dec=-6.6080986499400E+01 simput=n49.fits evtfile=evt_n49.fits exposure=10000
We can make an image with the simulated event file:
imgev evtfile=evt_n49.fits image=img_n49.fits coordinatesystem=0 projection=TAN \
NAXIS1=512 NAXIS2=512 CUNIT1=deg CUNIT2=deg CRVAL1=8.1491223145500E+01 \
CRVAL2=-6.6080986499400E+01 CRPIX1=256.5 CRPIX2=256.5 CDELT1=-6.207043e-04 \
CDELT2=6.207043e-04 history=true clobber=yes
Next, we will make a spectrum. We'll use a region file, src_n49.reg, that contains the following information, and remember to add X, Y sky coordinates and a WCS to the event file.
$my_computer> more src_n49.reg
# Region file format: DS9 version 4.1
global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 \
highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1
fk5
circle(81.5011768,-66.0831631,120.000")
radec2xy evtfile=evt_n49.fits projection=TAN refra=81.5011768 refdec=-66.0831631
makespec evtfile=evt_n49.fits spectrum=n49.pha eventfilter="regfilter(\"src_n49.reg\")" \
rsppath=/full/path/to/responsefile/ clobber=yes
Lastly, we can verify the fit, where the commands below are the same as those defined for the point source, and short explanantions of them can be found there.
sherpa
set_xsxsect("vern")
set_xsabund("wilm")
set_stat("cstat")
set_method("simplex")
load_pha("n49.pha")
ignore()
notice("0.5:7")
set_xlog()
set_ylog()
plot_data()
set_source(xsphabs.t1 * xspowerlaw.p1)
t1.nh=0.1
p1.phoindex=1
p1.norm=1
fit()
covar()
calc_energy_flux(0.5,7)
calc_data_sum(0.5,7)
group_counts(30)
plot_fit_resid()
exit()
The image and spectrum are shown in Figure 2.
|
|
If you have any questions concerning XMM-Newton send e-mail to xmmhelp@lists.nasa.gov