8. An RGS Data Processing and Analysis Primer

While a variety of analysis packages can be used for the following steps, the SAS was designed for the basic reduction and analysis of XMM-Newton data and will therefore be used here for demonstration purposes.

At this point, it is assumed that you have downloaded the data from the HEASARC archive onto a Hera server, standard or anonymous Hera is running (see § 4.2), you have prepared the data for processing with odfingest (see §6), and the working directory PROC has been made. Throughout this chapter, we will use the Mkn 421 dataset with ObsID 0153950701, available through links at the HEASARC archive.

8.1 Rerunning the pipeline

It is very likely that you will want to filter your data to some extent; in this case, you will need to reprocess it in order to determine the appropriate filters, regardless of the age of the observation. To do this, verify that the working directory PROC is highlighted in the GUI. Then,

1)

In the Available Tools panel, under ``XMM-SAS'', call rgsproc. (The task automatically processes all exposures, so if you want to process data only from RGS1, set ``enable explicit subset of exposures'' to ``yes'', and directly below, in the ``instrument exposure id list'' box, enter the six-digit string that corresponds to the RGS1 data, e.g., ``R1S001''. Verify that ``Dispersion binning method'' is set to ``lambda''.)

2)

Click ``Run'' in the task pop-up window.

This takes several minutes, and outputs 12 files per RGS, plus 3 general use FITS files. At this point, renaming files to something easy to type is a good idea. This is easily done by right-clicking on the event files. We will assume that the newly pipelined event files are named rgs1.fits and rgs2.fits.

8.2 Potentially useful tips for using the pipeline

The pipeline task, rgsproc, is very flexible and can address potential pitfalls for RGS users. In §8.1, we used the default parameter settings, and if this is sufficient for your data (and it should be for most), feel free to skip to §8.3. In the following sections, we will look at the cases of a nearby bright optical source, a nearby bright X-ray source, and a user-defined source.

8.2.1 A Nearby Bright Optical Source

With certain pointing angles, zeroth-order optical light may be reflected off the telescope optics and cast onto the RGS CCD detectors. If this falls on an extraction region, the current energy calibration will require a wavelength-dependent zero-offset. Stray light can be detected on RGS DIAGNOSTIC images taken before, during and after the observation. This test, and the offset correction, are not performed on the data before delivery.

To check for stray light and apply the appropriate offsets,

1)

Call rgsproc. Next to ``calculate offsets from diagnostic mode'' check ``yes''. Immediately below, next to ``Store histograms in a dataset'', check ``no''.

2)

Click ``Run''.

8.2.2 A Nearby Bright X-ray Source

In the example above, it is assumed that the field around the source contains sky only. Provided a bright background source is well-separated from the target in the cross-dispersion direction, a mask can be created that excludes it from the background region. Here the source has been identified in the EPIC images and its coordinates have been taken from the EPIC source list which is included among the pipeline products. The bright neighboring object is found to be the third source listed in the sources file. The first source is the target.

To exclude a nearby source from the background region,

1)

Call rgsproc. Next to ``add sources from an epic-derived source list'' click ``yes''. Immediately below, next to ``name of the epic- derived source list'', enter the name of the EPIC source list; it will have the format *EMSRLI*.FTZ. Next to ``select which sources to exclude from the background'', enter INDEX==1&&INDEX==3.

2)

Click ``Run''.

8.2.3 User-defined Source Coordinates

If the true coordinates of an object are not included in the EPIC source list or the science proposal, the user can define the coordinates of a new source.

1)

Call rgsproc. Next to ``add a source with user-defined position, rate, etc'' click ``yes''. Immediately below, next to ``label for the user-defined source'' enter the name of your source, Mkn421. Next to ``style for entering the position of the user-defined source'', enter how the source position is to be defined; the choices are radec (celestial coordinates) and wrtatt (RGS angular coordinates). For radec, in the two boxes just below, enter the RA and declination of the source in decimal degrees; for wrtatt, in the two boxes below RA and dec, enter the dispersion-direction offset from the pointing axis and the cross-dispersion-direction offset from the pointing axis, each in arminutes. For this example, lets use radec and set the RA to 166.113808 and the declination to +38.208833.

2)

Click ``Run''.

8.3 Examine and Analyze the Data

Since the event files are current, we can proceed with some simple analysis demonstrations, which will allow us to generate filters.

The heraXmmselect GUI, in the ``experimental'' folder of ``XMM-SAS'', provides a very simple method for producing and displaying images, spectra, and light curves, and is the recommended method for extracting data.

An event list from the RGS1 can be loaded by highlighting the event file ``rgs1.fits'' in our PROC directory, then clicking on heraXmmselect, and ``Run Tool...''. A Table Selection window will pop up to confirm your selection; clicking ``Go'' will start the task. Note that heraXmmselect always sees (and places) files in the directory that it was invoked in.

8.3.1 Create and Display an Image

Two commonly-made plots are those showing PI vs. BETA_CORR (also known as ``banana plots'') and XDSP_CORR vs. BETA_CORR.

To create such images,

1)

In the heraXmmselect window, check the square boxes to the left of the ``BETA_CORR'' and ``PI'' entries. Click the ``Image'' button near the bottom of the page. This brings up the evselect GUI.

2)

In the ``Image'' tab in the evselect GUI, confirm that the withimageset box is checked. In the imageset box, change the output image name from image.ds to something descriptive, in this case, pi_bc.fits.

3)

Click on the ``Run'' button on the lower left corner of the evselect GUI.

Different binnings and other selections can be invoked by accessing the ``Image'' tab at the top of the GUI. The default settings are reasonable, however, for a basic image. Similarly, plots can be made comparing BETA_CORR to XDSP_CORR. These two example plots can be seen in Figure 8.1.

Figure 8.1: Plots of XDSP_CORR vs. BETA_CORR (left) and PI vs. BETA_CORR (right). The gap is due to the missing CCD7. Similarly, CCD4 is missing in RGS2.

8.3.2 Create and Display a Light Curve

The background is assessed through examination of the light curve. We will extract a region, CCD9, that is most susceptible to proton events and generally records the least source events due to its location close to the optical axis. Also, to avoid confusing solar flares for source variability, a region filter that that removes the source from the final event list should be used. The region filters are kept in the source file product P*SRCLI_*.FIT.

To create a light curve,

1)

Enter the filtering criteria in the ``Selection expression'' box at the top of the xmmselect GUI:
(CCDNR==9)&&(REGION(P0153950701R1S001SRCLI_0000_PPS.FIT:RGS1_BACKGROUND,M_LAMBDA,XDSP_CORR))

2)

Check the circle to the left of the time entry.

3)

Click on the ``OGIP Rate Curve'' button near the bottom of the page. This brings up the evselect GUI (Figure 5.3). In the ``Lightcurve'' tab, confirm that the withrateset box is checked. Change the timebinsize to a reasonable amount, e.g. 10 or 100 s, and change the default output file name in the rateset box to something appropriate, in this case, r1_ltcrv.fits.

4)

Click on the ``Run'' button at the lower left corner of the evselect GUI.

The resultant light curve is displayed automatically using POWplot (see Figure 8.2).

Figure 8.2: Background event rate from the RGS1 CCD9 chip. The flares are solar events. The time units are elapsed mission time.

8.3.3 Generating the Good Time Interval (GTI) File

Examination of the lightcurve shows that there is a noisy section at the end of the observation, after 1.36975e8 seconds, where the count rate is well above the normal background count rate of $\sim$ 0.05 count/second. There are two procedures that make the GTI file (gtibuild and tabgtigen) that, when applied to the event file in another run of rgsproc, will excise these sections.

The first method, using gtibuild, requires a text file as input. This file can be made on your local machine and uploaded to your Hera account by right-clicking and dragging the file from your local directory to the remote directory. In the first two columns, refer to 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. In the example case, then, we would write in our ASCII file (named r1_gti.txt):

1.36958e8 1.36975e8 +

and proceed to the task gtibuild:

1)

Call gtibuild in the Hera GUI.

2)

Enter the name of the text file and the output GTI table (here, r1_gti.fits).

3)

Click ``Run''.

Alternatively, we can make the GTI file with tabgtigen and filter for RATE (though we could just as easily filter on TIME):

1)

Call the tabgtigen task.

2)

Next to ``Name of input data set'', enter the name of the lightcurve file, r1_ltcrv.fits. Directly below, enter the name of the output file, r1_gti.fits. Next to ``Booleam expression controlling the GTI creation'', enter the filtering expression. Since the nominal count rate is about 0.05 count/sec, we will set the upper limit to 0.2 count/sec: RATE$<$0.2

3)

Click ``Run''.

8.3.4 Applying the GTI

Now that we have GTI file, we can apply it to the event file by running rgsproc again. rgsproc is a complex task, running several steps, with five different entry and exit points. It is not necessary to rerun all the steps in the procedure, only the ones involving filtering.

To apply the GTI to the event file,

1)

Call rgsproc in the Available Tools panel.

2)

Next to ``reflection orders for spectrum generation'', make sure it is set for both orders, 1 2 . Next to ``(re)start processing at'' and ``conclude processing at'', use the pull-down menus to select 3:filter and 5:fluxing, respectively. Next to ``Additional GTI tables'', enter the name of the GTI file, r1_gti.fits. Verify that ``background correct the spectra'' is set to no. Verify that ``Produce column with m*lambda (source at requested position) is set to ``yes''.

3)

Click on ``Run''.

We will refer to the output event file as r1_filt.fits.

8.3.5 Creating the Response Matrices (RMFs)

Response matrices (RMFs) are not provided as part of the pipeline product package, so you must create your own before analyzing data. The task rgsproc generates a response matrix automatically, but as noted in §8.2.3, the source coordinates are under the observer's control. The source coordinates have a profound influence on the accuracy of the wavelength scale as recorded in the RMF that is produced automatically by rgsproc, and each RGS instrument and each order will have its own RMF.

Making the RMF is easily done with the package rgsrmfgen. Please note that, unlike with EPIC data, it is not necessary to make ancillary response files (ARFs).

To make the RMFs,

1)

Call the rgsrmfgen task.

2)

Next to ``input events and exposure map extensions'', enter the name of the filtered event list, rgs1_filt.fits. Next to ``output response file'', enter the name of the output file, r1_o1_rmf.fits. Next to ``extracted spectrum file'', enter the name of the spectrum file; it has the form *SRSPEC*. Next to ``maximum output energy'', enter 2.5; for the minimum output energy, enter 0.4. Next to ``number of energy bins'', enter 5000.

3)

Click ``Run''.

At this point, the spectra can be analyzed. If you you wish, skip the discussion on combining spectra (§8.3.6) and go straight to fitting the spectrum (§8.4.)

8.3.6 Combining Spectra

Spectra from the same order in RGS1 and RGS2 can be safely combined to create a spectrum with higher signal-to-noise if they were reprocessed using rgsproc with spectrumbinning=lambda, as we did in §8.1 (this also happens to be the default). The task rgscombine also merges response files and background spectra. When merging response files, be sure that they have the same number of bins. For this example, we assume that RMFs were made for order 1 in both RGS1 and RGS2.

To merge RGS1 and RGS2 spectra,

1)

Call the rgscombine task.

2)

Next to ``output combined spectrum file'', enter the name of the combined spectrum, r12_o1_srspec.fits. Next to ``output combined back spectrum file'', enter the name of the combined background spectrum, r12_o1_bgspec.fits. Next to ``output combined response file'', enter the name of the combined RMFs, r12_o1_rmf.fits. Next to ``spectra files'', enter the spectra to be merged, P0153950701R1S001SRSPEC1001.FIT P0153950701R2S002SRSPEC1001.FIT. Next to ``corresponding response files'' enter the response files to be merged, r1_o1_rmf.fits r2_o1_rmf.fits. Next to ``corresponding background files'', enter the background spectra to be combined, P0153950701R1S001BGSPEC1001.FIT P0153950701R2S002BGSPEC1001.FIT

3)

Click ``Run''.

The spectra are ready for analysis, so we can prepare the spectrum for fitting.

8.4 Approaches to Spectral Fitting and the Cash Statistic

For data sets of high signal-to-noise and low background, where counting statistics are within the Gaussian regime, the data products above are suitable for analysis using the default fitting scheme in XSPEC, $\chi^2$-minimization. However, for low count rates, in the Poisson regime, $\chi^2$-minimization is no longer suitable. With low count rates in individual channels, the error per channel can dominate over the count rate. Since channels are weighted by the inverse-square of the errors during $\chi^2$ model fitting, channels with the lowest count rates are given overly-large weights in the Poisson regime. Spectral continua are consequently often fit incorrectly, with the model lying underneath the true continuum level. This will be a common problem with most RGS sources. Even if count rates are large, much of the flux from these sources can be contained within emission lines, rather than the continuum. Consequently, even obtaining correct equivalent widths for such sources is non-trivial.

The traditional way to increase the signal-to-noise of a data set is to rebin or group the channels, since, if channels are grouped in sufficiently large numbers, the combined signal-to-noise of the groups will jump into the Gaussian regime. However, this results in the loss of information. For example, sharp features like an absorption edge or emission line can be completely washed out. Further, in the Poisson regime, the background spectrum cannot simply be subtracted, as is commonly done in the Gaussian regime, since this could result in negative counts. Therefore, rebinning should be reserved for fast, preliminary analysis of spectra without sharp features, or for making plots for publication. When working on the final analysis for a low-count data set, the (unbinned) background and source spectra should be fitted simultaneously using the Cash statistic. (If fitting with XSPEC, be sure you are running v11.1.0 or later. This is because RGS spectrum files have prompted a slight modification to the OGIP standard, since the RGS spatial extraction mask has a spatial-width which is a varying function of wavelength. Thus, it has become necessary to characterize the BACKSCL and AREASCL parameters as vectors (i.e., one number for each wavelength channel), rather than scalar keywords as they are for data from the EPIC cameras and past missions. These quantities map the size of the source extraction region to the size of the background extraction region and are essential for accurate fits. Only Xspec v11.1.0, or later versions, are capable of reading these vectors, so be certain that you have an up-to-date installation at your site.)

Finally, a caveat of using the Cash statistic in Xspec is that the scheme requires a ``total'' and ``background'' spectrum to be loaded into Xspec. This is in order to calculate parameter errors correctly. Consequently, be sure not to use the ``net'' spectra that were created as part of product packages by SAS v5.2 or earlier. To change schemes in Xspec before fitting the data, type:

XSPEC$>$ statistic cstat

For our sample spectrum, we will rebin and fit it with $\chi^2$ statistics.

8.4.1 Spectral Rebinning

SAS does not include fitting software, so HEASoft packages will be used, and all fitting tasks will be called from the command line.

There are two ways to rebin a spectrum: the FTOOL grppha, or the RGS pipeline. grppha can group channels using an algorithm which bins up consecutive channels until a count rate threshold is reached. This method conserves the resolution in emission lines above the threshold while improving statistics in the continuum.

To rebin the spectrum with grppha, type

grppha

and edit the parameters as needed:

   > Please enter PHA filename[] P0153950701R1S001SRSPEC1001.FIT
   > Please enter output filename[] P0153950701R1S001SRSPEC1001.bin30.FIT
   > GRPPHA[] group min 30 
   > GRPPHA[] exit

The disadvantage of using grppha is that, although channel errors are propagated through the binning process correctly, the errors column in the original spectrum product is not strictly accurate. The problem arises because there is no good way to treat the errors within channels containing no counts. To allow statistical fitting, these channels are arbitrarily given an error value of unity, which is subsequently propagated through the binning. Consequently, the errors are overestimated in the resulting spectra.

The other approach, which involves calling the RGS pipeline after it is complete, bins the data during spectral extraction. The following rebins the pipeline spectrum by a factor 3.

To rebin the spectrum with rgsproc, type

rgsproc orders='1 2' rebin=3 rmfbins=5000 entrystage=4:spectra
$$ finalstage=5:fluxing bkgcorrect=no

where

orders - dispersion orders to extract
rebin - wavelength rebinning factor
rmfbins - number of bins in the response file; this should be greater than 3000
entrystage - entry stage to the pipeline
finalstage - exit stage for the pipeline

One disadvantage of this approach is that you can only choose integer binning of the original channel size. To change the sampling of the events, the pipeline must be run from the second stage (``angles'') or earlier:

rgsproc orders='1 2' nbetabins=1133 rmfbins=5000 entrystage=2:angles
$$ finalstage=fluxing bkgcorrect=no

where the parameters are as defined previously, and

nbetabins - number of bins in the dispersion direction; the default is 3400

The disadvantage of using rgsproc, as opposed to grppha, is that the binning is linear across the dispersion direction. Velocity resolution is lost in the lines, so the accuracy of redshift determinations will be degraded, transition edges will be smoothed, and neighboring lines will become blended.

8.5 Fitting a Spectral Model

We can fit the spectrum using Xspec. This is easily done by entering

xspec

Enter the data, background, and response file at the prompts, and edit the fitting parameters as needed.

      XSPEC> data P0153950701R1S001SRSPEC1001.bin30.FIT ! input data 
      XSPEC> back P0153950701R1S001BGSPEC1001.bin30.FIT ! input background
      XSPEC> resp r1_o1_rmf.fits                        ! input response file
      XSPEC> ignore **-0.4                              ! set sensible limits
      XSPEC> model wabs*pow                             ! set spectral model to absorbed powerlaw
      1:wabs:nH> 0.01                                   ! enter reasonable initial values
      2:powerlaw:PhoIndex>2.0
      3:powerlaw:norm>1.0
      XSPEC> renorm
      XSPEC> fit
      XSPEC> cpd /xw
      XSPEC> setplot wave
      XSPEC> setplot command window all
      XSPEC> setplot command log x off
      XSPEC> setplot command wind 1
      XSPEC> setplot command r y 0 6e-17
      XSPEC> setplot command wind 2
      XSPEC> setplot command r y -10 10
      XSPEC> plot data chi
      XSPEC> exit

Figure 8.3 shows the fit to the spectrum.

Figure 8.3: 1st order RGS1 spectrum of Mkn 421. The fit is an absorbed power law model. The gap between 10-15Å is due to the absence of CCD7.

8.6 Analysis of Extended Sources

8.6.1 Region masks

The optics of the RGS allow spectroscopy of reasonably extended sources, up to a few arc minutes. The width of the spatial extraction mask is defined by the fraction of total events one wishes to extract. With the default pipeline parameter values, over 90% of events are extracted, assuming a point-like source.

Altering and optimizing the mask width for a spatially-extended source may take some trial and error, and, depending on the temperature distribution of the source, may depend on which lines one is currently interested in. While AB Dor is not an extended source, the following example increases the width of the extraction mask and ensures that the size of the background mask is reduced so that the two do not overlap.

To adjust the region mask with rgsproc in the GUI:

1)

Call rgsproc in the ``Available Tools'' panel of the Hera GUI.

3)

Next to ``(re)start processing at'', use the pull-down menu to select ``4:spectra''. Next to ``percent of cross-dispersion PSF included'' and ``percent of cross-dispersion PSF excluded from the background'', enter 99. Next to ``percent of pulse-height distribution included'', enter 95.

4)

Click ``Run''.

Observing extended sources effectively broadens the psf of the spectrum in the dispersion direction. Therefore, it is prudent to also increase the width of the PI masks using the pdistincl parameter in order to prevent event losses.

8.6.2 Making RMFs for extended sources

RGS response matrices as made in §8.3.5 are appropriate for use with point sources only. If we are interested in analyzing an extended source, the RMF must take into account the spatial degradation of the resolution. The most straight-forward way to do this is to modify the response matrix prior to spectral fitting. For sources extended up to about 1 arcminute, this can be done with the FTOOL rgsrmfsmooth. It requires three files: the point source RMF (as made in §8.3.5), an image of the source (from an EPIC camera, see §7.2.1, or different mission), and a text file. The better the resolution of the image, the more accurate the modified RMF will be, so if a Chandra image is available for a source, it should be used instead of an EPIC image. The text file must list the name of the image, the boresight, and the aperture size in the following format:

RGS_XSOURCE_IMAGE <name of source image>
RGS_XSOURCE_BORESIGHT <image boresight: RA (h:m:s), DEC(degrees:m:s), and PA (decimal
$$ degrees)>
RGS_XSOURCE_EXTRACTION <source radius in arcminutes>

For an example case, we will name our text file xsource.mod. We will assume that a RMF for the first order grating was made as in §8.3.5 and an MOS1 image was made as in §7.2.1; xsource.mod contains these lines:

RGS_XSOURCE_IMAGE my_MOS1_image.fits
RGS_XSOURCE_BORESIGHT 05:28:45 -65:26:55 219.769546508789
RGS_XSOURCE_EXTRACTION 1.2

This file can be made on the user's local machine and uploaded to the Hera server by right-clicking and dragging the file from the Local Directory panel to the Remote Directory.

To make an RMF for an extended source in the Command Window, type:

rgsrmfsmooth rmffil=my_ps_rmf.fits imgfil=xsource.mod order=1 outfil=my_es_rmf.fits

where

rmffil - the point source RGS RMF
imgfil - the text file with the name of the source image, boresight, and extraction region
order - grating order of RMF
outfil - output RMF name