THE XMM-NEWTON ABC GUIDE, STREAMLINED
EPIC (IMAGING mode), SAS Command Line
Contents
Prepare the Data
Reprocess the Data
Apply Standard Filters
Make a Light Curve
Apply Time Filters
Making an Image and Source Detection
Make an Exposure Map
Extract the Source and Background Spectra
Check for Pile Up
My Observation is Piled Up! Now What?
Determine the Spectrum Extraction Areas
Create the Photon Redistribution Matrix (RMF) and Ancillary File (ARF)
Prepare the Data
Please note that the two tasks in this section (cifbuild and odfingest) must be run in the ODF directory. These are the only tasks with that requirement, and after this section, we will work exclusively in our reprocessing directory.Many SAS tasks require calibration information from the Calibration Access Layer (CAL). Relevant files are accessed from the set of Current Calibration File (CCF) data using a CCF Index File (CIF). To set the environment parameter and make the ccf.cif file, type
cd ODF setenv SAS_ODF /full/path/to/ODF/directory/ setenv SAS_ODFPATH /full/path/to/ODF/directory/ cifbuild
To use the updated CIF file in further processing, you will need to reset the environment variable SAS_CCF:
setenv SAS_CCF /full/path/to/ODF/ccf.cif
The task odfingest extends the Observation Data File (ODF) summary file with data extracted from the instrument housekeeping data files and the calibration database. It is only necessary to run it once on any dataset, and will cause problems if it is run a second time. If for some reason odfingest must be rerun, you must first delete the earlier file it produced. This file largely follows the standard XMM naming convention, but has SUM.SAS appended to it. After running odfingest, you will need to reset the environment variable SAS_ODF to its output file. To run odfingest and reset the environment variable, type
odfingest setenv SAS_ODF /full/path/to/ODF/full_name_of_*SUM.SAS
You will likely find it useful to alias these environment variable resets in your login shell (.cshrc, .bashrc, etc.).
Reprocess the Data
To reprocess the data, make a new working directory and call emproc for MOS data, and epproc for PN data.
- emproc
epproc
By default, these tasks do not keep any intermediate files they
generate.
Emproc and epproc designate
their output event files with "ImagingEvts.ds". It is
convenient to rename them something easy to type. Here,
we will make a soft link:
- ln -s 0070_0123700101_EMOS1_S001_ImagingEvts.ds mos1.fits
Apply Standard Filters
The filtering expressions for the MOS and PN are:
(PATTERN <= 12)&&(PI in [200:12000])&&#XMMEA_EM
and
(PATTERN <= 4)&&(PI in [200:15000])&&#XMMEA_EP
If you are going to extract a spectrum from the PN data, you will need the most stringent filters in order to get a high-quality spectrum:
(PATTERN <= 4)&&(PI in [200:15000])&&#XMMEA_EP&&(FLAG == 0)
For the MOS, the standard filters should be appropriate for most cases.
The first two expressions will select good events with PATTERN in the 0 to 12 range. The PATTERN value is similar the GRADE selection for ASCA data, and is related to the number and pattern of the CCD pixels triggered for a given event.The PATTERN assignments are: single pixel events: PATTERN == 0, double pixel events: PATTERN in [1:4], triple and quadruple events: PATTERN in [5:12].
The second keyword in the expressions, PI, selects the preferred pulse height of the event; for the MOS, this should be between 200 and 12000 eV. For the PN, this should be between 200 and 15000 eV. This should clean up the image significantly with most of the rest of the obvious contamination due to low pulse height events. Setting the lower PI channel limit somewhat higher (e.g., to 300 eV) will eliminate much of the rest.
Finally, the #XMMEA_EM (#XMMEA_EP for the PN) filter provides a canned screening set of FLAG values for the event. (The FLAG value provides a bit encoding of various event conditions, e.g., near hot pixels or outside of the field of view.) Setting FLAG == 0 in the selection expression provides the most conservative screening criteria and should always be used when serious spectral analysis is to be done on the PN. It typically is not necessary for the MOS.
To filter the data, type
evselect table=mos1.fits filtertype=expression filteredset=mos1_filt.fits \
expression='(PATTERN <= 12) && (PI in [200:12000]) && #XMMEA_EM'
where
-
table - input event table
filtertype - method of filtering
expression - filtering expression
filteredset - output file name
Make a Light Curve
Sometimes, it is necessary to use filters on time in addition to those mentioned above. This is because of soft proton background flaring, which can have count rates of 100 counts/sec or higher across the entire bandpass.
It should be noted that the amount of flaring that needs to be removed depends in
part on the object observed; a faint, extended object will be more affected than
a very bright X-ray source.
To determine if our observation is affected by background flaring, we can make a light curve and display it with, for instance, fv:
evselect table=mos1_filt.fits withrateset=Y rateset=mos1_ltcrv.fits \
maketimecolumn=Y timebinsize=100 makeratecolumn=yes
where
-
table - input event table
rateset - name of output light curve file
maketimecolumn - make a time column
timebinsize - time binning (seconds)
makeratecolumn - make a count rate column, otherwise a count column will be created
The light curve is shown below.
Apply Time Filters
Taking a look at the light curve, we can see that there is a very large flare toward the end of the observation and two much smaller ones in the middle of the exposure.
There are several ways to filter an event file: on TIME, with an explicit reference to the TIME parameter in the filtering expression or by creating a secondary Good Time Interval (GTI) file, or on RATE, which requires making a new GTI file. New GTI files are easily made with the task tabgtigen or gtibuild. These are discussed in detail below. Any of these methods will produce a cleaned file, so which one to use is a matter of the user's preference. An image of the MOS1 events file with both the standard and time filters applied is shown in Figure 2.
Filter on RATE With tabgtigen
Examining the light curve shows us that during non-flare times, the count rate is quite low, about 1.3 ct/s, with a small increase at 7.3223e7 seconds to about 6 ct/s. We can use that to generate the GTI file and apply it:
tabgtigen table=mos1_ltcrv.fits gtiset=gtiset.fits \
expression='RATE<=6'
evselect table=mos1_filt.fits filtertype=expression \
filteredset=mos1_filt_time.fits expression='GTI(gtiset.fits,TIME)'
Filter on TIME With tabgtigen
Alternatively, we could have chosen to make a new GTI file by noting the times of the flaring in the light curve and using that as a filtering parameter. The big flare starts around 7.32276e7 s, and the smaller ones are at 7.32119e7 s and 7.32205e7 s. The expression to remove these would be (TIME <= 73227600)&&!(TIME IN [7.32118e7:7.3212e7])&&!(TIME IN [7.32204e7:7.32206e7]). The syntax &&(TIME < 73227600) includes only events with times less than 73227600, and the "!" symbol stands for the logical "not", so use &&!(TIME in [7.32118e7:7.3212e7]) to exclude events in that time interval.
tabgtigen table=mos1_ltcrv.fits gtiset=gtiset.fits timecolumn=TIME \
expression='(TIME <= 73227600)&&!(TIME IN [7.32118e7:7.3212e7])&&!(TIME IN [7.32204e7:7.32206e7])'
evselect table=mos1_filt.fits filtertype=expression \
filteredset=mos1_filt_time.fits expression='GTI(gtiset.fits,TIME)'
Filter on TIME With gtibuild
This method requires a text file as input. In the first two columns, enter the start and end times (in seconds) that you are interested in, and in the third column, indicate with either a + or - sign whether that region should be kept or removed. Each good (or bad) time interval should get its own line, with any optional comments preceeded by a "#". In the example case, we would write in our ASCII file (named gti.txt):
0 73227600 + # Good time from the start of the observation... 7.32118e7 7.3212e7 - # but without a small flare here, 7.32204e7 7.32206e7 - # and here.
and proceed to gtibuild and evselect:
gtibuild file=gti.txt table=gtiset.fits
evselect table=mos1_filt.fits filtertype=expression \
filteredset=mos1_filt_time.fits expression='GTI(gtiset.fits,TIME)'
where
-
file - input text file
table - output GTI file
Filter on TIME By Explicit Reference
Finally, we could have chosen to forgo using tabgtigen or gtibuild altogether, and simply filtered on TIME with the standard filtering expression. In that case, the full filtering expression would be:
(PATTERN <= 12) && (PI in [200:12000]) && #XMMEA_EM &&
(TIME <= 73227600) &&!(TIME IN [7.32118e7:7.3212e7])&&!(TIME IN [7.32204e7:7.32206e7])
This expression can then be used to
filter the original event file, or only the time requirements can be used to filter
the file that has already had the standard filters applied:
evselect table=mos1_filt.fits filtertype=expression filteredset=mos1_filt_time.fits \
expression='(TIME <= 73227600)&&!(TIME IN [7.32118e7:7.3212e7])&&!(TIME IN [7.32204e7:7.32206e7])'
Making an Image and Source DetectionSource Detection
The edetect_chain task does nearly all the work involved with EPIC source detection. It can process up to three intruments (both MOS cameras and the PN) with up to five images in different energy bands simultaneously. All images must have identical binning and WCS keywords. For this example, we will perform source detection on MOS1 images in two bands (``soft'' X-rays with energies between 300 and 2000 eV, and ``hard'' X-rays, with energies between 2000 and 10000 eV) using the filtered event files produced here. Edetect_chain has many parameters that users might want to edit, such as detection likelihood threshold; more information about them can be found here.
We will start by generating some files that edetect_chain needs: an attitude file and images of the sources in the desired energy bands.
First, the attitude file:
atthkgen atthkset=attitude.fits
where
- atthkset - output file name
Next, the images. We can choose either to bin to a certain image size (by setting the imagebinning parameter to imageSize and specifying values for ximagesize and yimagesize; the default values makes a 600x600 pixel image), or we can set the bin size directly (by setting imagebinning to binSize and specifying values for ximagebinsize and yimagebinsize).
In order to select a reasonable binsize, we need to know a little about the detectors themselves. The pixels of the MOS and PN cameras correspond to 1.1'' and 4.1'' on the sky, respectively, but photons incident on the detector have their positions randomized within these pixels into "virtual pixels" of size 0.05''. So, if we want image pixels that correspond to 1.1'' (or 4.1'') we need a bin size of 22 (or 82).
In addition to making images in soft and hard bands, we'll also make an
image that includes both bands, for display purposes.
evselect table=mos1_filt_time.fits withimageset=yes imageset=mos1-s.fits \
imagebinning=binSize xcolumn=X ximagebinsize=22 ycolumn=Y yimagebinsize=22 \
filtertype=expression expression='(FLAG == 0)&&(PI in [300:2000])'
evselect table=mos1_filt_time.fits withimageset=yes imageset=mos1-h.fits \
imagebinning=binSize xcolumn=X ximagebinsize=22 ycolumn=Y yimagebinsize=22 \
filtertype=expression expression='(FLAG == 0)&&(PI in [2000:10000])'
evselect table=mos1_filt_time.fits withimageset=yes imageset=mos1-all.fits \
imagebinning=binSize xcolumn=X ximagebinsize=22 ycolumn=Y yimagebinsize=22 \
filtertype=expression expression='(FLAG == 0)&&(PI in [300:10000])'
The parameters for all these are as listed previously,
and
-
withimageset - flag to create an image
imageset - fits image name to be created
imagebinning - how to bin the image
xcolumn - table column to use for the X axis
ximagebinsize - binning in X axis
ycolumn - table column to use for the Y axis
yimagebinsize - binning in Y axis
Now we can run edetect_chain.
edetect_chain imagesets='mos1-s.fits mos1-h.fits' eventsets='mos1_filt_time.fits' \
attitudeset=attitude.fits pimin='300 2000' pimax='2000 10000' \
ecf='0.878 0.220'
where
-
imagesets - list of count images
eventsets - list of event files
attitudeset - attitude file name
pimin - list of minimum PI channels for the bands
pimax - list of maximum PI channels for the bands
ecf - energy conversion factors for the bands
The energy conversion factors (ECFs) convert the source count rates into fluxes. The ECFs for each detector and energy band depend on the pattern selection and filter used during the observation. For more information, please consult the 3XMM Catalogue User Guide. Those used here are derived from PIMMS using the flux in the 0.1-10.0 keV band, a source power-law index of 1.9, an absorption of 0.5×1020 cm-2.
We can display the results of eboxdetect using srcdisplay and produce a region file for the sources.
srcdisplay boxlistset=emllist.fits imageset=mos1-all.fits \
regionfile=regionfile.txt sourceradius=0.01 withregionfile=yes
where
-
boxlistset - eboxdetect source list
imageset - image file name over which the source circles are to be plotted
includesources - flag to include the source positions on the display
regionfile - file name of output file containing source regions
sourceradius - radius of circle plotted to locate sources
withregionfile - flag to create a region file
Figure 3 shows the MOS1 event overlayed with the detected sources.
Make an Exposure Map
Prior to making an exposure map, you will need to make an image in the energy band you are interested in. You will also need the attitude file from atthkgen. Let's say that we're interested in the 1-2 keV band; we can use the attitude file that we made before, and to make the image,
evselect table=mos1_filt_time.fits withimageset=yes imageset=mos1_1-2keV.fits \
imagebinning=binSize xcolumn=X ximagebinsize=22 ycolumn=Y yimagebinsize=22 \
filtertype=expression expression='(PI in [1000:2000])'
Then,
eexpmap imageset=mos1_1-2keV.fits attitudeset=attitude.fits eventset=mos1_filt_time.fits \
expimageset=mos1_expmap_1-2keV.fits pimin=1000 pimax=2000
Extract the Source and Background Spectra
Throughout the following, please keep in mind that some parameters are instrument-dependent. The parameter specchannelmax should be set to 11999 for the MOS, or 20479 for the PN. Also, for the PN, the most stringent filters, (FLAG==0)&&(PATTERN<=4), must be included in the expression to get a high-quality spectrum.For the MOS, as noted before, the standard filters should be appropriate for many cases. However, there are some instances where tightening the selection requirements might be needed. For example, if obtaining the best possible spectral resolution is critical to your work, and the corresponding loss of counts is not important, only the single pixel events should be selected (PATTERN==0). If your observation is of a bright source, you again might want to select only the single pixel events to mitigate pile up.
In any case, you'll need to know spatial information about the area over which you want to extract the spectrum, so display the filtered event file with ds9:
ds9 mos1_filt_time.fits &
Select the object whose spectrum you wish to extract and resize as desired. For this example,
we will choose the source at (26188.5,22816.5) and set the extraction radius to 300 (in physical units.)
To extract the source spectrum,
evselect table='mos1_filt_time.fits' energycolumn='PI' filteredset='mos1_filtered.fits' \
filtertype='expression' expression='((X,Y) in CIRCLE(26188.5,22816.5,300))' \
spectrumset='mos1_pi.fits' spectralbinsize=5 \
withspecranges=yes specchannelmin=0 specchannelmax=11999
where the parameters are as indicated previously,
and
-
spectrumset - name of output spectrum
spectralbinsize - size of bin, in eV
withspecranges - covering a certain spectral range
specchannelmin - minimum of spectral range
specchannelmax - maximum of spectral range
evselect table=mos1_filt_time.fits energycolumn='PI' filteredset='bkg_filtered.fits' \
filtertype='expression' spectrumset='bkg_pi.fits' spectralbinsize=5 \
expression='((X,Y) in CIRCLE(26188.5,22816.5,1500))&&!((X,Y) in CIRCLE(26188.5,22816.5,500))' \
withspecranges=yes specchannelmin=0 specchannelmax=11999
Check for Pile Up
Depending on how bright the source is and what modes the EPIC detectors are in, event pile up may be a problem. Pile up occurs when a source is so bright that incoming X-rays strike two neighboring pixels or the same pixel in the CCD more than once in a read-out cycle. In such cases the energies of the two events are in effect added together to form one event. If this happens sufficiently often, 1) the spectrum will appear to be harder than it actually is, and 2) the count rate will be underestimated, since multiple events will be undercounted. To check whether an observation is piled up, use the SAS task epatplot. Heavily piled sources will be immediately obvious, as they will have a "hole" in the center of their image, but pile up is not always so conspicuous. Therefore, we recommend to always check for it.
Note that this procedure requires as input the event file created when the spectrum was made, not the usual time-filtered event file.
To check for pile up in our Lockman Hole example,
epatplot set=mos1_filtered.fits plotfile=mos1_epat.ps useplotfile=yes \
withbackgroundset=yes backgroundset=bkg_filtered.fits
where
-
set - name of the event file
useplotfile - use file name indicated in plotfile for output?
plotfile - output file name
withbackgroundset - use the background event set for background subtraction?
backgroundset - name of background event file
The output of epatplot is a postscript file, mos1_epat.ps, which may be viewed with viewers such as gv, containing two graphs describing the distribution of counts as a function of PI channel, as seen in Figure 3 (left).
A few words about interpretting the plots are in order. The top plot is the distribution of counts versus PI channel for each pattern class (single, double, triple, quadruple), and the bottom is the expected pattern distribution (smooth lines) plotted over the observed distribution (histogram). The lower plot shows the model distributions for single and double events and the observed distributions. It also gives the ratio of observed-to-modeled events with 1-σ uncertainties for single and double pattern events over a given energy range. (The default is 0.5-2.0 keV; this can be changed with the pileupnumberenergyrange parameter.) If the data is not piled up, there will be good agreement between the modeled and observed single and double event pattern distributions. Also, the observed-to-modeled fractions for both singles and doubles in the specified energy range will be unity, within errors. In contrast, if the data is piled up, there will be clear divergence between the modeled and observed pattern distributions, and the observed-to-modeled fraction for singles will be less than 1.0, and for doubles, it will be greater than 1.0.
Finally, when examining the plots, it should noted that the observed-to-modeled fractions can be inaccurate. Therefore, the agreement between the modeled and observed single and double event pattern distributions in the lower plot should be the main factor in determining if an observation is affected by pile up or not.
The source used in our Lockman Hole example is too faint to provide reasonable statistics for epatplot and is far from being affected by pile up. For comparison, an example of a bright source (from a different observation) which is strongly affected by pileup is shown in Figure 3 (right). Note that the observed-to-model fraction for doubles is over 1.0, and there is severe divergence between the model and the observed pattern distribution.
|
My Observation is Piled Up! Now What?
If you're working with a different (much brighter) dataset that does show signs of pile up, there are a few ways to deal with it. First, using the region selection and event file filtering procedures demonstrated in earlier sections, you can excise the inner-most regions of a source (as they are the most heavily piled up), re-extract the spectrum, and continue your analysis on the excised event file. For this procedure, it is recommended that you take an iterative approach: remove an inner region, extract a spectrum, check with epatplot, and repeat, each time removing a slightly larger region, until the model and observed distribution functions agree. If you do this, be aware that removing too small a region with respect to the instrumental pixel size (1.1'' for the MOS, 4.1'' for the PN) can introduce systematic inaccuracies when calculating the source flux; these are less than 4%, and decrease to less than 1% when the excised region is more than 5 times the instrumental pixel half-size. In any case, be certain that the excised region is larger than the instrumental pixel size!
You can also use the event file filtering procedures to include only single pixel events (PATTERN==0), as these events are less sensitive to pile up than other patterns.
Determine the Spectrum Extraction Areas
Now that we are confident that our spectrum is not piled up, we can continue by finding the source and background region areas. This is done with the task backscale, which takes into account any bad pixels or chip gaps, and writes the result into the BACKSCAL keyword of the spectrum table. Alternatively, we can skip running backscale, and use a keyword in arfgen below. We will show both options for the curious.
To find the source and background extraction areas explicitly:
backscale spectrumset=mos1_pi.fits badpixlocation=mos1_filt_time.fits
backscale spectrumset=bkg_pi.fits badpixlocation=mos1_filt_time.fits
where
-
spectrumset - name of the spectrum file
badpixlocation - the file containing the bad pixel locations (the event file)
Create the Photon Redistribution Matrix (RMF) and Ancillary File (ARF)
Now that a source spectrum has been extracted, we need to reformat the detector response by making a redistribution matrix file (RMF) and ancillary response file (ARF). To make the RMF:
rmfgen rmfset=mos1_rmf.fits spectrumset=mos1_pi.fits
where
-
rmfset - output file
spectrumset - spectrum file
Please note that rmfgen can take several minutes (~10) to run.
Now use the RMF, spectrum, and event file to make the ancillary file.
arfgen arfset=mos1_arf.fits spectrumset=mos1_pi.fits withrmfset=yes \
rmfset=mos1_rmf.fits withbadpixcorr=yes badpixlocation=mos1_filt_time.fits
If we had not run backscale, we could set a keyword in arfgen to find the region area:
arfgen arfset=mos1_arf.fits spectrumset=mos1_pi.fits withrmfset=yes \
rmfset=mos1_rmf.fits withbadpixcorr=yes badpixlocation=mos1_filt_time.fits \
setbackscale=yes
where
-
arfset - output ARF file name
spectrumset - input spectrum file name
withrmfset - flag to use the RMF
rmfset - RMF file created by rmfgen
withbadpixcorr - flag to include the bad pixel correction
badpixlocation - file containing the bad pixel information; should be set to the event
file from which the spectrum was extracted.
setbackscale - flag to calculate the area of the source region and write it to the
BACKSCAL keyword in the spectrum header
The spectrum can be fit using HEASoft or CIAO packages, as SAS does not include fitting software.
If you have any questions concerning XMM-Newton send e-mail to xmmhelp@lists.nasa.gov