8. An EPIC Data Processing and Analysis Primer (Timing Mode, Command Line)

So, you've received an XMM-Newton EPIC data set. What are you going to do with it? After checking what the observation consists of (see § 3.2), you should note when the observation was taken. If it is a recent observation, it was likely processed with the most recent calibrations and SAS, and you can immediately start to analyze the Pipeline Processed data. However, if it is more than a year old, it was probably processed with older versions of CCF and SAS prior to archiving, and the pipeline should be rerun to generate event files with the latest calibrations.

As noted in Chapter 4, a variety of analysis packages can be used for the following steps. However, as the SAS was designed for the basic reduction and analysis of XMM-Newton data (extraction of spatial, spectral, and temporal data), it will be used here for demonstration purposes. SAS will be required at any rate for the production of detector response files (RMFs and ARFs) and other observatory-specific requirements. (Although for the simple case of on-axis point sources the canned response files provided by the SOC can be used.)

It is strongly recommended that you keep all reprocessed data in its own directory! SAS places output files in whichever directory it is in when a task is called. Throughout this primer, it is assumed that the Pipleline Processed data are in the PPS directory, the ODF data (with upper case file names, and uncompressed) are in the directory ODF, the analysis is taking place in the PROC directory, and the CCF data are in the directory CCF.

If your data are recent, you need only to gunzip the files and prepare the data for processing (see §5. Feel free to skip the section on repipelining and proceed to the later discussions. In any case, for simplicity, it is recommended that you change the name of the unzipped event file to something easy to type. For example, a PN event list:

cp PPS/PiiiiiijjkkPNSlllTIEVLI0000.FTZ PROC/pn.fits

where

iiiiiijjkk - observation number
lll - exposure number within the observation

Various analysis procedures are demonstrated using the Cen X-3 dataset, ObsID 0400550201. The following procedures are applicable to all XMM-Newton datasets, so it is not required that you use this particular dataset; any Timing Mode observation should be sufficient.

For detailed descriptions of PP data nomenclature, file contents, and which tasks can be used to view them, see Tables 3.2 and 3.3. For detailed descriptions of ODF data nomenclature and file contents, see Table 3.1.

8.1 Rerun the Pipeline

We assume that the data was prepared and environment variables were set according to §5. In the window where SAS was initialized, in your ``processing directory'' PROC, run epproc to produce calibrated photon event files for the PN camera.

Note that epproc will automatically detect what mode the data were taken in.

On the command line, just type:

epproc

By default, the task does not keep any intermediate files that were generated. Epproc designates its output event files with ``Evts.ds'', so ``*TimingEvts.ds'' and ``*BurstEvts.ds'' denote the timing and burst mode event lists, respectively. In any case, you may want to name the new files something easy to type. For example, to make a soft link to the new PN event files output from epproc, type

ln -s 1206_0403530301_EPN_S003_TimingEvts.ds pn.fits

8.2 Create and Display an Image

To create an image of your timing data, type

evselect table=pn.fits withimageset=yes imageset=image.fits \
$$ xcolumn=RAWX ycolumn=RAWY imagebinning=binSize ximagebinsize=1 yimagebinsize=1

where

table - input event table
withimageset - make an image
imageset - name of output image
xcolumn - event column for X axis
ycolumn - event column for Y axis
imagebinning - form of binning, force entire image into a given size or bin by a specified number of pixels
ximagesize - output image pixels in X
yimagesize - output image pixels in Y

The output file image.fits can be viewed by using a standard FITS display, such as ds9 (see Figure 8.1) :

ds9 image.fits &

Figure 8.1: The PN Timing image, displayed in ds9.

8.3 Applying Standard Filters to the Data

The filtering expression for the PN in Timing mode is:

(PATTERN $<=$ 4)&&(PI in [200:15000])&&#XMMEA_EP

The first two expressions will select good events with PATTERN in the 0 to 4 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. Single pixel events have PATTERN == 0, while double pixel events have PATTERN in [1:4].

The second keyword in the expressions, PI, selects the preferred pulse height of the event; 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 or 400 eV) will eliminate much of the rest.

Finally, the #XMMEA_EP 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 PN data.

To filter the data, type

evselect table=pn.fits withfilteredset=yes \
$$ expression='(PATTERN $<=$ 4)&&(PI in [200:15000])&&#XMMEA_EP' \
$$ filteredset=pn_filt.fits filtertype=expression keepfilteroutput=yes \
$$ updateexposure=yes filterexposure=yes

where

table - input event table
filtertype - method of filtering
expression - filtering expression
withfilteredset - create a filtered set
filteredset - output file name
keepfilteroutput - save the filtered output
updateexposure - update exposure information in event list and in spectrum files
filterexposure - filter exposure extensions of event list with same time

8.4 Create and Display 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.

To determine if our observation is affected by background flaring, we can examine the light curve. For the time binning, we will set it to something reasonable (usually between 10 and 100 s):

evselect table=pn.fits withrateset=yes rateset=pn_ltcrv.fits \
$$ maketimecolumn=yes timecolumn=TIME timebinsize=50 makeratecolumn=yes

where

table - input event table
withrateset - make a light curve
rateset - name of output light curve file
maketimecolumn - control to create a time column
timecolumn - time column label
timebinsize - time binning (seconds)
makeratecolumn - control to create a count rate column, otherwise a count column will be created

The output file pn_ltcrv.fits can be viewed by using fv:

fv pn_ltcrv.fits &

In the fv pop-up window, the RATE extension will be available in the second row (index 1, as numbering begins with 0). Select ``PLOT'' from this row, and select the column name and axis on which to plot it. The light curve is shown in Fig. 8.2. No flares are evident, so we will continue to the next section. However, if a dataset does contain flares, they should be removed in the same way as shown for EPIC Imaging mode data in §6.5.

Figure 8.2: The light curve of Cen X-3 from 0.2-15 keV, displayed in fv.

8.5 Extract the Source and Background Spectra

The first step in extracting a spectrum from PN Timing data is to make an image of the event file over the energy range we are interested in; for this example, we'll say 0.5-15 keV. And since this is the PN, we need to remember to set (FLAG==0) to get a high-quality spectrum. Thus, our expression parameter would be set to (FLAG==0) && (PI in [500:15000]), and the entire command would be:

evselect table=pn_filt.fits withimageset=yes imageset=image.fits \
$$ xcolumn=RAWX ycolumn=RAWY imagebinning=binSize ximagebinsize=1 \
$$ yimagebinsize=1 expression='(FLAG==0) && (PI in [500:15000])'

where the parameters are as described in §8.2.

The image can be displayed with ds9 and is shown in Fig. 8.3 (top). The source is centered on RAWX=37; we will extract this and the 10 pixels on either side of it:

evselect table=pn_filt.fits spectrumset=source_pi_WithBore.fits \
$$ energycolumn=PI spectralbinsize=5 specchannelmin=0 specchannelmax=20479 \
$$ filteredset=pn_filt_source_WithBore.fits \
$$ expression='(FLAG==0) && (PI in [500:15000]) && (RAWX in [27:47])'

where the keywords are the same as in §8.3, and

energycolumn - energy column
spectrumset - name of output spectrum
spectralbinsize - size of bin, in eV
specchannelmin - minimum of spectral range
specchannelmax - maximum of spectral range

For the background, the extraction area should be as far from the source as possible. However, sources with $>$ 200 ct/s (like our example!) are so bright that they dominate the entire CCD area, and there is no source-free region from which to extract a background. (It goes without saying that this is highly energy-dependent.) In such a case, it may be best not to subtract a background. Users are referred to Ng et al. (2010, A&A, 522, 96) for an in-depth discussion. While this observation is too bright to have a good background extraction region, the process is shown below nonetheless for the sake of demonstration:

evselect table=pn_filt.fits withspectrumset=yes spectrumset=bkg_pi.fits \
$$ energycolumn=PI spectralbinsize=5 withspecranges=yes specchannelmin=0 \
$$ expression='(FLAG==0) && (PI in [500:15000]) && (RAWX in [3:5])' \
$$ specchannelmax=20479 withfilteredset=y filteredset=pn_filt_bkg.fits

Figure 8.3: TOP: The Cen X-3 Timing Mode image, from 0.5-15 keV. The green line is a projection cut. BOTTOM: The average counts in the cut across the CCD.

8.6 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. Briefly, we deal with it in PN Timing data essentially the same way as in Imaging data, that is, by using only single pixel events, and/or removing the regions with very high count rates, checking the amount of pile up, and repeating until it is no longer a problem. We recommend to always check for it.

Note that this procedure requires as input the event files created when the spectrum was made, not the usual time-filtered event file.

To check for pile up:

epatplot set=pn_filt_source_WithBore.fits plotfile=pn_epat.ps useplotfile=yes \
$$ withbackgroundset=yes backgroundset=pn_flt_bkg.fits

where

set - input events file
plotfile - output postscript file
useplotfile - flag to use file name from "plotfile"
withbackgroundset - use background event set for background subtraction?
backgroundset - name of background event file

The output of epatplot is a postscript file, pn_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; see Figure 8.4.

A few words about interpretting the plots are in order. The top 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-$\sigma$ 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 0.5-2.0 keV 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 should be the main factor in determining if an observation is affected by pile up or not.

Examining the plots, we see that there is a large difference between the modeled and observed single and double pattern events, and that the observed-to-model fraction for doubles is over 1.0, indicating that the observation is piled up.

Figure 8.4: The output of epatplot.

8.7 My Observation is Piled Up! Now What?

There are a couple ways to deal with pile up. First, you can use event file filtering procedures to include only single pixel events (PATTERN==0), as these events are less sensitive to pile up than other patterns.

You can also excise areas of high count rates, i.e., the boresight column and several columns to either side of it. (This is analogous to removing the inner-most regions of a source in Imaging data.) The spectrum can then be re-extracted and you can continue your analysis on the excised event file. As with Imaging data, 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 pattern distributions agree.

To extract only the columns to either side of the boresight:

evselect table=pn_filt.fits withspectrumset=yes spectrumset=source_pi_NoBore.fits \
$$ energycolumn=PI spectralbinsize=5 withspecranges=yes specchannelmin=0 \
$$ specchannelmax=20479 \
$$ expression='(FLAG ==0)&&(PI in [500:15000])&&(RAWX in [27:47])&&!(RAWX in [29:45]) \
$$ withfilteredset=yes filteredset=pn_filt_source_NoBore.fits

Be aware that if you do this and are using SAS v. 13.x or older, you will need to use a non-standard way to make the ancillary files (ARFs) for your spectrum! This is discussed further in §8.9.

8.8 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 process is identical to that used for IMAGING data.) 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.

To find the source and background extraction areas:

backscale spectrumset=source_pi_NoBore.fits badpixlocation=pn_filt.fits

backscale spectrumset=bkg_pi.fits badpixlocation=pn_filt.fits

where

spectrumset - spectrum file
badpixlocation - event file containing the bad pixels

8.9 Create the Photon Redistribution Matrix (RMF) and Ancillary File (ARF)

If you are using SAS v. 14 or higher, making the RMF and ARF for PN data in TIMING mode is exactly the same as in IMAGING mode, even if you had to excise piled up areas. This is a change from earlier SAS versions; if you are working with an older SAS, you will need to use the special recipe below to generate the ARF (the method to make a RMF file is the same as shown here.)

To make the RMF:

rmfgen rmfset=source_rmf_NoBore.fits spectrumset=source_pi_NoBore.fits

where

rmfset - output file
spectrumset - spectrum file

To make the ARF:

arfgen arfset=source_arf_NoBore.fits spectrumset=source_pi_NoBore.fits \
$$ detmaptype=psf withrmfset=yes rmfset=source_rmf_NoBore.fits badpixlocation=pn_filt.fits

where

arfset - output file
spectrumset - spectrum file
arfset - output file
detmaptype - origin of the detector map
withrmfset - use the RMF dataset to define the ARF energy grid?
rmfset - RMF file
badpixlocation - the file containing the bad pixel locations

If you excised regions to make a spectrum and are using SAS v. 13.x or older, you will need the spectra of the full extraction area and the excised area. We already have it for the full extraction area, so for the excised area,

evselect table=pn_filt.fits withspectrumset=yes spectrumset=source_pi_Excised.fits \
$$ energycolumn=PI spectralbinsize=5 withspecranges=yes specchannelmin=0 \
$$ specchannelmax=20479 \
$$ expression='(FLAG ==0)&&(PI in [500:15000])&&(RAWX in [27:47])&&!(RAWX in [29:45]) \
$$ withfilteredset=yes filteredset=pn_filt_source_Excised.fits

Now we can use the spectra to make the ARFs:

arfgen arfset=source_arf_WithBore.fits spectrumset=source_pi_WithBore.fits \
$$ detmaptype=psf badpixlocation=pn_filt.fits

arfgen arfset=source_arf_Excised.fits spectrumset=source_pi_Excised.fits \
$$ detmaptype=psf badpixlocation=pn_filt.fits

where the parameters are as previously described. Now we can subtract them:

addarf "source_arf_WithBore.fits source_arf_Excised.fits" "1.0 -1.0" \
$$ source_arf_NoBore.fits

At this point, the spectrum is ready to be analyzed, so skip ahead to prepare the spectrum for fitting (§13).

The timing data can also be examined in Xronos; this is discussed in §15.