How to reduce RGS data and extract spectra of point-like sources

 

 

Introduction

This thread contains a step-by-step recipe to process RGS data of point-like sources.

The purpose of this jupyter notebook thread is to illustrate how to use SAS within python notebooks. For an in-depth explanation on how to reduce RGS data using SAS, please refer to the standard SAS threads,

Expected Outcome

Source and background RGS spectra and the associated response matrices as well as lightcurves.

SAS Tasks to be Used

Prerequisites

Useful Links

Caveats

Last Reviewed: 25 May 2023, for SAS v21.0

Last Updated: 25 May 2023
 

 

 

Procedure

This thread is intended to show the steps to follow to extract RGS spectra and lightcurves from a Jupiter notebook. It expects that the ODF or PPS have been downloaded into separate folders and that the SAS Startup Thread has been followed.

The task rgsproc can redo different stages of processing without starting from scratch. Allow enough time for it to finish. It may take between 5-10 minutes depending on the ODF size and your computer. The different entry and exit points are called processing stages, of which there are six:

  • 1:events preliminary tasks, source-independent calibrations
  • 2:angles aspect-drift corrections
  • 3:filter filters events and exposure
  • 4:spectra generates spectra
  • 5:fluxing generates response matrices and a combined flux-calibrated spectrum
  • 6:lightcurve generates lightcurves

If the task has been run only through spectra neither the response matrices nor the fluxed spectra or the lightcurves are created. The fluxing stage is the most time consuming, and therefore it is recommended to run it only once the results of the previous stages are satisfactory.

 

In order to run the following blocks, SAS needs to be inizialized in the terminal from which the Jupyter nootebook has been opened. This has to be done prior to launching the notebook. Follow the steps provided in the SAS Startup Thread in Python.

 

Import the following Python libraries,

In [ ]:
from pysas.wrapper import Wrapper as w
import os, sys,glob,shutil
import pathlib
from astropy.io import fits
from astropy import coordinates as coo
from astropy import units as u
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.table import Table
from astropy.convolution import Gaussian2DKernel, convolve
 

Before we start, lets see what we have already defined in terms of SAS variables,

In [ ]:
inargs = []
t = w('sasver', inargs)
t.run()
 

We still need to define the SAS_ODF and SAS_CCF variables (make sure the SAS_CCFPATH is also defined). Use the SAS task startsas as described in the SAS Startup Thread in Python to download an ODF. For this thread, we will assume that an ODF has already been downloaded. We will use the ODF with id 00096020101 as an example and assume that there is already a directory (sas_file in the block below) with an existing CIF and SAS Summary File. In the block below introduce the absoute path to the Working directory (finished with '/') and the location of the CIF and SAS Summary File,

Note: the path to the CIF and SAS Summary File must be an absolute path begining with '/'.

In [ ]:
work_dir = 'absolute_path_to_wrk_directory'
sas_file = 'absolute_path_to_cifSUMSAS_directory'
In [ ]:
inargs = [f'sas_ccf={sas_file}ccf.cif', f'sas_odf={sas_file}0082_0096020101_SCX00000SUM.SAS', f'workdir={work_dir}']
In [ ]:
w('startsas', inargs).run()
 

This process should have defined the SAS_CCF an SAS_ODF variables.

 

In the following blocks we will produce an RGS event file and associated files to extract a spectra and lightcurves of point-like source.

 
  • First, run the rgsproc task with no further parameters. The default entry and final stages are 1:events and 6:lightcurve, respectively.
In [ ]:
# SAS Command
cmd        = "rgsproc" # SAS task to be executed                  

# Arguments of SAS Command
inargs = [ ]           # do not provide any parameters 

# Execute the SAS task to obtain RGS products
w(cmd, inargs).run()
 

For most analysis, this is all is needed, however, there are a number of checks we can do.

 

Checking for periods of high background activity

 

 

In line with the other XMM-Newton instruments, RGS observations can be affected by high particle background periods caused by solar activity or by the Earth radiation belts. You should check your observations and, if necessary, filter out these periods before extracting any scientific products. In the blocks that follow will we use RGS1. The equivalent commands can be used for RGS2. We will use the ODF with id 0096020101 as an example.

 

  • Create a light curve of the pure background events. For this, only events in the background region of CCD number 9 are selected as it is the one closest to the optical axis of the telescope, and the most affected by background flares.
In [ ]:
# Define some RGS1 files needed in the next blocks

evt_list     = work_dir+'P0096020101R1S004EVENLI0000.FIT'  # event list
src_list     = work_dir+'P0096020101R1S004SRCLI_0000.FIT'  # source list
mer_list     = work_dir+'P0096020101R1S004merged0000.FIT'  # combined event list
In [ ]:
# Function to plot Lightcurve
def plotLC(plt,threshold,fileName):
        if fileName != "NOT FOUND":
            fitsFile = fits.open(fileName)
            prihdu   = fitsFile[1].header
            if ('CUTVAL' in prihdu):
                threshold = prihdu['CUTVAL']

            cols    = fitsFile[1].columns
            colName = None
            for i,x  in enumerate(cols.names):
                if "RATE" in x:
                    colName = cols.names[i]
                if "COUNTS" in x:
                    colName = cols.names[i]

            data = fitsFile[1].data

            xdata = data.field('TIME') - min(data.field('TIME')) # extract the x data column
            ydata = data.field(colName)

            xmax=np.amax(xdata)
            xmin=np.amin(xdata)

            plt.plot(xdata,ydata,label=fileName) # plot the data

            if colName == 'RATE':
                plt.title("Lightcurve")
                plt.xlabel("Time (s)")
                plt.ylabel("Cts/s")
            else:
                plt.title("Lightcurve")
                plt.xlabel("Time (s)")
                plt.ylabel("Counts")

            if (threshold != 'None'):
                if (colName == 'COUNTS'):
                    threshold=float(threshold)*100.

                y2data = [threshold]*len(xdata)
                plt.plot(xdata,y2data)
                plt.text(xmin+0.1*(xmax-xmin), threshold+0.01*threshold, str(threshold)+" cts/sec", ha='center')

            fitsFile.close()

        else:
            print("File not found "+fileName+"\n")
In [ ]:
# Function to plot RGS plots
def plotRGS(plt,fileName,rgsfile,srcfile,source):
    hdu=fits.open(fileName)
    hdr=hdu[1].header
    ev=Table(hdu[1].data)

    plt.subplot(2,1,1)

    img=plt.hist2d(ev['M_LAMBDA'],ev['PI'],[360,250],cmap='Oranges',vmin=0,vmax=5)
    hdu.close()

    hdu=fits.open(srcfile)
    inst=hdr['INSTRUME'].strip()
    reg=Table(hdu[inst+'_SRC'+str(source)+'_ORDER_1'].data)
    nreg=len(reg)

    for i in range(nreg):
        rr=reg[i]
        ii=np.nonzero(rr['LAMBDA'])
        if rr['SHAPE'][0] == '!':
            linest=':r'
        else:
            linest='-b'
        plt.plot(rr['LAMBDA'][ii],rr['PI'][ii],linest)

    plt.title(os.path.basename(fileName))
    plt.xlabel('Wavelength (A)')
    plt.ylabel('PI (eV)')

    hdu.close()

    plt.subplot(2,1,2)

    hdus=fits.open(rgsfile)
    r1=Table(hdus[1].data)
    xdsp=np.degrees(r1['XDSP_CORR'])*3600.
    plt.rcParams["figure.figsize"]=(12,6)
    img,xe,ye=np.histogram2d(xdsp,r1['M_LAMBDA'],[250,250])

    kernel=np.array([[1,1,1],[1.,3.,1.],[1.,9,1.],[1.,3.,1.],[1.,1,1.]])
    kernel=kernel/np.sum(kernel)
    imgs=convolve(img,kernel)

    plt.imshow(imgs,origin='lower',vmin=0,vmax=1.5,cmap='rainbow',
          extent=[ye[0],ye[-1],xe[0],xe[-1]],aspect='auto',)

    hdus.close()

    hdul=fits.open(srcfile)
    src=Table(hdul['SRCLIST'].data)
    dd=src['DELTA_XDSP']*(-60)
    reg=Table(hdul['RGS1_SRC'+str(source)+'_SPATIAL'].data)
    nreg=len(reg)

    for i in range(nreg):
        rr=reg[i]
        ii=np.nonzero(rr['LAMBDA'])
        plt.plot(rr['LAMBDA'][ii],np.degrees(rr['XDSP_CORR'][ii])*3600,'r')
    reg=Table(hdul['RGS1_BACKGROUND'].data)
    nreg=len(reg)

    for i in range(nreg):
        rr=reg[i]
        ii=np.nonzero(rr['LAMBDA'])
        if rr['SHAPE'][0] != '!':
            linest='-g'
        else:
            linest='--y'
        plt.plot(rr['LAMBDA'][ii],np.degrees(rr['XDSP_CORR'][ii])*3600,linest,clip_on=True)

    plt.xlabel('Wavelength')
    plt.ylabel('Cross-dispersion (arcsec)')
    plt.xlim([6.,32.])

    hdul.close()
In [ ]:
# Define a SAS filter expression to derive a background rate cut

lc_binsize   = 100      # lightcurve bin size in secs

out_LCFile   = work_dir+'RGS1_FlareBKGRate.fit'  # Name of the output BKG lightcurve
In [ ]:
# SAS Command
cmd        = "evselect" # SAS task to be executed                  

# Arguments of SAS Command
expression = f'(CCDNR==9)&&(REGION({src_list}:RGS1_BACKGROUND,M_LAMBDA,XDSP_CORR))'  # event filter expression
inargs     = [f'table={evt_list}','withrateset=Y',f'rateset={out_LCFile}','maketimecolumn=Y',f'timebinsize={lc_binsize}','makeratecolumn=Y',f'expression={expression}']

print("   Filter expression to use: "+expression+" \n")
print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs
In [ ]:
# Execute the SAS task with the parameters to produce the flaring background lightcurve.

w(cmd, inargs).run()
 
  • Check the created background light curve.
In [ ]:
# Inspect light curve

rgs_threshold = 0.1         # Rate cts/sec (only used here for display purposes. Draws a horizontal line at this value)
plt.figure(figsize=(20,8))  # Size of figure

plotLC(plt,rgs_threshold,out_LCFile) # Plot source region light curve

plt.legend()
 
  • If necessary, create an additional GTI for filtering periods of high background out of the event list, to be used in conjunction with the internal GTI tables. The expression in the tabgtigen task indicates that only periods with count rates less than rgs_rcut counts/sec should be selected. The value of rgs_rcut should be chosen after inspection of the background light curve; it goes typically from 0.1 to 2 counts/sec.
In [ ]:
# Define some parameters to produce the gti file

rgs_rcut    = 0.1     # Rate cts/sec 

gti_list    = work_dir+'RGS1_GTI.fit'  # Name of the output GTI file
In [ ]:
# SAS Command
cmd        = "tabgtigen" # SAS task to be executed                  

# Arguments of SAS Command
expression = f'(RATE<{rgs_rcut})'
inargs     = [f'table={out_LCFile}',f'gtiset={gti_list}',f'expression={expression}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs
In [ ]:
# Execute the SAS task with the parameters to produce GTI file

w(cmd, inargs).run()
 
  • Re-process the data with rgsfilter. The filter stage requires a combined event list named PxxxxxxyyyyRrzeeemerged0000.FIT as input. Since we have already run rgsproc already, such a file has been created and it is in the working directory.
In [ ]:
# SAS Command
cmd        = "rgsfilter" # SAS task to be executed                  

# Arguments of SAS Command
inargs     = [f'mergedset={mer_list}',f'evlist={evt_list}',f'auxgtitables={gti_list}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs
In [ ]:
# Execute the SAS task with the parameters to filter RGS file

w(cmd, inargs).run()
 
  • Run rgsproc from the spectra stage to obtain the spectral products after GTI filtering.
In [ ]:
# SAS Command
cmd        = "rgsproc" # SAS task to be executed                  

# Arguments of SAS Command
inargs     = [f'auxgtitables={gti_list}',f'entrystage=\'4:spectra\'']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs
In [ ]:
# Execute the SAS task with the parameters to produce RGS spectral products after gti filtering

w(cmd, inargs).run()
 
  • Display the dispersion vs cross dispersion and dispersion vs energy images and overlay the selected region masks.
In [ ]:
# Define some parameters for evselect

outfile   ='RGS1.fits' # Output event file name
In [ ]:
# SAS Command
cmd        = "evselect" # SAS task to be executed                  

# Arguments of SAS Command
expression = f'REGION({src_list}:RGS1_SRC1_ORDER_1,M_LAMBDA,PI)'
inargs=[f'table={evt_list}',f'expression={expression}',f'keepfilteroutput=yes',f'withfilteredset=yes',f'filteredset={outfile}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs
In [ ]:
# Execute the SAS task with the parameters to produce an image

w(cmd, inargs).run()
In [ ]:
# Inspect plots

plt.figure(figsize=(20,18))  # Size of figure

source = 1 # if user has defined the coordinate (1: default ; 3: user defined)
plotRGS(plt,evt_list,outfile,src_list,source) # Plot RGS plots

plt.legend()
 
  • In the figure above, red solid lines indicate the source spectrum extraction region; yellow dashed lines indicate excluded areas for background spectrum extraction; and the green solid lines correspond to the background spectrum extraction region.
  • It is very important to check that the selection regions are correctly centred on the prime source. Also, check that the region used for background extraction excludes other existing sources.
 

Changing the coordinates of the prime source

 
  • The coordinates of the prime source can be changed (in case for example that they are incorrect). The new coordinates can be entered and added to the source list. Add to the argument list (inargs) when calling rgsproc the following parameters: f'withsrc=yes','srclabel=USER',f'srcra={RA}',f'srcdec={DEC}' where the new coordinates are, RA = 184.6110 DEC =+29.81267 withscr=yes is used to define a new prime source by the parameters srclabel, srcra and srcdec. The label is the name that the new source will get in the source list, coordinates are J2000 in decimal degrees. See the thread rgsproc, Coordinates and Masks for further explanations.
In [ ]:
# Define some parameters for rgsproc

RA  = 184.6110  # source RA
DEC =+29.81267   # source DEC
In [ ]:
# SAS Command
cmd        = "rgsproc" # SAS task to be executed                  

# Arguments of SAS Command
inargs=[f'withsrc=yes','srclabel=USER',f'srcra={RA}',f'srcdec={DEC}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs
In [ ]:
# Execute the SAS task to run rgsproc changing the source coordinates

w(cmd, inargs).run()
 
  • Display again the dispersion vs cross dispersion and dispersion vs energy images and overlay the selected region masks.
In [ ]:
# Define some parameters for evselect

outfile   ='RGS1_User.fits' # Output event file name
 
  • In the expression defined in the block below, notice the change from, RGS1_SRC1_ORDER_1 to RGS1_SRC3_ORDER_1, which indicates that the coordinates of the source being used are those definde by the user.
In [ ]:
# SAS Command
cmd        = "evselect" # SAS task to be executed                  

# Arguments of SAS Command
expression = f'REGION({src_list}:RGS1_SRC3_ORDER_1,M_LAMBDA,PI)'
inargs=[f'table={evt_list}',f'expression={expression}',f'keepfilteroutput=yes',f'withfilteredset=yes',f'filteredset={outfile}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs
In [ ]:
# Execute the SAS task with the parameters to produce an image

w(cmd, inargs).run()
In [ ]:
# Inspect plots

plt.figure(figsize=(20,18))  # Size of figure

source = 3 # if user has defined the coordinate (1: default ; 3: user defined)
plotRGS(plt,evt_list,outfile,src_list,source) # Plot RGS plots

plt.legend()
 

 

 

 

 

 

 

Caveats

None