{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

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

\n", "
\n", "
\n", "\n", "\t\n", "\t\t\n", "\t\t\t\n", "\t\t\n", "\t\n", "
\n", "\t\t\t

Introduction

\n", "\t\t\t

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

\n", "

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,

\n", " \n", "\t\t\t

Expected Outcome

\n", "\t\t\t

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

\n", "\t\t\t

SAS Tasks to be Used

\n", "\t\t\t\n", "\t\t\t

Prerequisites

\n", "\t\t\t\n", "\t\t\t

Useful Links

\n", "\t\t\t\n", "\t\t\t

Caveats

\n", "\t\t\t

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

\n", "\t\t\t

Last Updated: 25 May 2023

\n", "\t\t\t
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Procedure

\n", "

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.

\n", "\n", "

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:

\n", " \n", "\n", "

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.

\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

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.

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Import the following Python libraries,

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from pysas.wrapper import Wrapper as w\n", "import os, sys,glob,shutil\n", "import pathlib\n", "from astropy.io import fits\n", "from astropy import coordinates as coo\n", "from astropy import units as u\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from astropy.io import fits\n", "from astropy.table import Table\n", "from astropy.convolution import Gaussian2DKernel, convolve" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

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

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inargs = []\n", "t = w('sasver', inargs)\n", "t.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

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,

\n", "\n", "Note: the path to the CIF and SAS Summary File must be an absolute path begining with '/'." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "work_dir = 'absolute_path_to_wrk_directory'\n", "sas_file = 'absolute_path_to_cifSUMSAS_directory'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inargs = [f'sas_ccf={sas_file}ccf.cif', f'sas_odf={sas_file}0082_0096020101_SCX00000SUM.SAS', f'workdir={work_dir}']" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "w('startsas', inargs).run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

This process should have defined the SAS_CCF an SAS_ODF variables.

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

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

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"rgsproc\" # SAS task to be executed \n", "\n", "# Arguments of SAS Command\n", "inargs = [ ] # do not provide any parameters \n", "\n", "# Execute the SAS task to obtain RGS products\n", "w(cmd, inargs).run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

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

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Checking for periods of high background activity

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

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.

\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define some RGS1 files needed in the next blocks\n", "\n", "evt_list = work_dir+'P0096020101R1S004EVENLI0000.FIT' # event list\n", "src_list = work_dir+'P0096020101R1S004SRCLI_0000.FIT' # source list\n", "mer_list = work_dir+'P0096020101R1S004merged0000.FIT' # combined event list" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Function to plot Lightcurve\n", "def plotLC(plt,threshold,fileName): \n", " if fileName != \"NOT FOUND\":\n", " fitsFile = fits.open(fileName) \n", " prihdu = fitsFile[1].header\n", " if ('CUTVAL' in prihdu):\n", " threshold = prihdu['CUTVAL']\n", "\n", " cols = fitsFile[1].columns\n", " colName = None\n", " for i,x in enumerate(cols.names):\n", " if \"RATE\" in x:\n", " colName = cols.names[i]\n", " if \"COUNTS\" in x:\n", " colName = cols.names[i] \n", "\n", " data = fitsFile[1].data \n", " \n", " xdata = data.field('TIME') - min(data.field('TIME')) # extract the x data column\n", " ydata = data.field(colName)\n", "\n", " xmax=np.amax(xdata) \n", " xmin=np.amin(xdata) \n", "\n", " plt.plot(xdata,ydata,label=fileName) # plot the data\n", "\n", " if colName == 'RATE':\n", " plt.title(\"Lightcurve\")\n", " plt.xlabel(\"Time (s)\")\n", " plt.ylabel(\"Cts/s\")\n", " else:\n", " plt.title(\"Lightcurve\")\n", " plt.xlabel(\"Time (s)\")\n", " plt.ylabel(\"Counts\")\n", "\n", " if (threshold != 'None'):\n", " if (colName == 'COUNTS'):\n", " threshold=float(threshold)*100.\n", "\n", " y2data = [threshold]*len(xdata) \n", " plt.plot(xdata,y2data)\n", " plt.text(xmin+0.1*(xmax-xmin), threshold+0.01*threshold, str(threshold)+\" cts/sec\", ha='center') \n", " \n", " fitsFile.close() \n", " \n", " else:\n", " print(\"File not found \"+fileName+\"\\n\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Function to plot RGS plots\n", "def plotRGS(plt,fileName,rgsfile,srcfile,source): \n", " hdu=fits.open(fileName)\n", " hdr=hdu[1].header\n", " ev=Table(hdu[1].data)\n", "\n", " plt.subplot(2,1,1)\n", " \n", " img=plt.hist2d(ev['M_LAMBDA'],ev['PI'],[360,250],cmap='Oranges',vmin=0,vmax=5)\n", " hdu.close()\n", " \n", " hdu=fits.open(srcfile)\n", " inst=hdr['INSTRUME'].strip()\n", " reg=Table(hdu[inst+'_SRC'+str(source)+'_ORDER_1'].data)\n", " nreg=len(reg)\n", " \n", " for i in range(nreg):\n", " rr=reg[i]\n", " ii=np.nonzero(rr['LAMBDA'])\n", " if rr['SHAPE'][0] == '!':\n", " linest=':r'\n", " else:\n", " linest='-b'\n", " plt.plot(rr['LAMBDA'][ii],rr['PI'][ii],linest)\n", " \n", " plt.title(os.path.basename(fileName))\n", " plt.xlabel('Wavelength (A)')\n", " plt.ylabel('PI (eV)')\n", " \n", " hdu.close()\n", " \n", " plt.subplot(2,1,2)\n", "\n", " hdus=fits.open(rgsfile)\n", " r1=Table(hdus[1].data)\n", " xdsp=np.degrees(r1['XDSP_CORR'])*3600.\n", " plt.rcParams[\"figure.figsize\"]=(12,6)\n", " img,xe,ye=np.histogram2d(xdsp,r1['M_LAMBDA'],[250,250]) \n", "\n", " kernel=np.array([[1,1,1],[1.,3.,1.],[1.,9,1.],[1.,3.,1.],[1.,1,1.]])\n", " kernel=kernel/np.sum(kernel)\n", " imgs=convolve(img,kernel)\n", "\n", " plt.imshow(imgs,origin='lower',vmin=0,vmax=1.5,cmap='rainbow',\n", " extent=[ye[0],ye[-1],xe[0],xe[-1]],aspect='auto',)\n", "\n", " hdus.close()\n", " \n", " hdul=fits.open(srcfile)\n", " src=Table(hdul['SRCLIST'].data)\n", " dd=src['DELTA_XDSP']*(-60) \n", " reg=Table(hdul['RGS1_SRC'+str(source)+'_SPATIAL'].data)\n", " nreg=len(reg)\n", "\n", " for i in range(nreg):\n", " rr=reg[i]\n", " ii=np.nonzero(rr['LAMBDA'])\n", " plt.plot(rr['LAMBDA'][ii],np.degrees(rr['XDSP_CORR'][ii])*3600,'r')\n", " reg=Table(hdul['RGS1_BACKGROUND'].data)\n", " nreg=len(reg)\n", " \n", " for i in range(nreg):\n", " rr=reg[i]\n", " ii=np.nonzero(rr['LAMBDA'])\n", " if rr['SHAPE'][0] != '!':\n", " linest='-g'\n", " else:\n", " linest='--y' \n", " plt.plot(rr['LAMBDA'][ii],np.degrees(rr['XDSP_CORR'][ii])*3600,linest,clip_on=True)\n", " \n", " plt.xlabel('Wavelength')\n", " plt.ylabel('Cross-dispersion (arcsec)')\n", " plt.xlim([6.,32.])\n", "\n", " hdul.close()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define a SAS filter expression to derive a background rate cut\n", "\n", "lc_binsize = 100 # lightcurve bin size in secs\n", "\n", "out_LCFile = work_dir+'RGS1_FlareBKGRate.fit' # Name of the output BKG lightcurve" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"evselect\" # SAS task to be executed \n", "\n", "# Arguments of SAS Command\n", "expression = f'(CCDNR==9)&&(REGION({src_list}:RGS1_BACKGROUND,M_LAMBDA,XDSP_CORR))' # event filter expression\n", "inargs = [f'table={evt_list}','withrateset=Y',f'rateset={out_LCFile}','maketimecolumn=Y',f'timebinsize={lc_binsize}','makeratecolumn=Y',f'expression={expression}'] \n", "\n", "print(\" Filter expression to use: \"+expression+\" \\n\")\n", "print(\" SAS command to be executed: \"+cmd+\", with arguments; \\n\")\n", "inargs" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the SAS task with the parameters to produce the flaring background lightcurve.\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Inspect light curve\n", "\n", "rgs_threshold = 0.1 # Rate cts/sec (only used here for display purposes. Draws a horizontal line at this value)\n", "plt.figure(figsize=(20,8)) # Size of figure\n", "\n", "plotLC(plt,rgs_threshold,out_LCFile) # Plot source region light curve\n", "\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define some parameters to produce the gti file\n", "\n", "rgs_rcut = 0.1 # Rate cts/sec \n", "\n", "gti_list = work_dir+'RGS1_GTI.fit' # Name of the output GTI file " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"tabgtigen\" # SAS task to be executed \n", "\n", "# Arguments of SAS Command\n", "expression = f'(RATE<{rgs_rcut})'\n", "inargs = [f'table={out_LCFile}',f'gtiset={gti_list}',f'expression={expression}']\n", "\n", "print(\" SAS command to be executed: \"+cmd+\", with arguments; \\n\")\n", "inargs " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the SAS task with the parameters to produce GTI file\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"rgsfilter\" # SAS task to be executed \n", "\n", "# Arguments of SAS Command\n", "inargs = [f'mergedset={mer_list}',f'evlist={evt_list}',f'auxgtitables={gti_list}']\n", "\n", "print(\" SAS command to be executed: \"+cmd+\", with arguments; \\n\")\n", "inargs " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the SAS task with the parameters to filter RGS file\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"rgsproc\" # SAS task to be executed \n", "\n", "# Arguments of SAS Command\n", "inargs = [f'auxgtitables={gti_list}',f'entrystage=\\'4:spectra\\'']\n", "\n", "print(\" SAS command to be executed: \"+cmd+\", with arguments; \\n\")\n", "inargs " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the SAS task with the parameters to produce RGS spectral products after gti filtering\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define some parameters for evselect\n", "\n", "outfile ='RGS1.fits' # Output event file name" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"evselect\" # SAS task to be executed \n", "\n", "# Arguments of SAS Command\n", "expression = f'REGION({src_list}:RGS1_SRC1_ORDER_1,M_LAMBDA,PI)'\n", "inargs=[f'table={evt_list}',f'expression={expression}',f'keepfilteroutput=yes',f'withfilteredset=yes',f'filteredset={outfile}']\n", "\n", "print(\" SAS command to be executed: \"+cmd+\", with arguments; \\n\")\n", "inargs " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the SAS task with the parameters to produce an image\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Inspect plots\n", "\n", "plt.figure(figsize=(20,18)) # Size of figure\n", "\n", "source = 1 # if user has defined the coordinate (1: default ; 3: user defined)\n", "plotRGS(plt,evt_list,outfile,src_list,source) # Plot RGS plots\n", "\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Changing the coordinates of the prime source

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define some parameters for rgsproc\n", "\n", "RA = 184.6110 # source RA\n", "DEC =+29.81267 # source DEC" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"rgsproc\" # SAS task to be executed \n", "\n", "# Arguments of SAS Command\n", "inargs=[f'withsrc=yes','srclabel=USER',f'srcra={RA}',f'srcdec={DEC}']\n", "\n", "print(\" SAS command to be executed: \"+cmd+\", with arguments; \\n\")\n", "inargs " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the SAS task to run rgsproc changing the source coordinates\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define some parameters for evselect\n", "\n", "outfile ='RGS1_User.fits' # Output event file name\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"evselect\" # SAS task to be executed \n", "\n", "# Arguments of SAS Command\n", "expression = f'REGION({src_list}:RGS1_SRC3_ORDER_1,M_LAMBDA,PI)'\n", "inargs=[f'table={evt_list}',f'expression={expression}',f'keepfilteroutput=yes',f'withfilteredset=yes',f'filteredset={outfile}']\n", "\n", "print(\" SAS command to be executed: \"+cmd+\", with arguments; \\n\")\n", "inargs " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the SAS task with the parameters to produce an image\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Inspect plots\n", "\n", "plt.figure(figsize=(20,18)) # Size of figure\n", "\n", "source = 3 # if user has defined the coordinate (1: default ; 3: user defined)\n", "plotRGS(plt,evt_list,outfile,src_list,source) # Plot RGS plots\n", "\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
 \n", "
\n", "\n", "\t\n", "\t\t\n", "\t\t\t\n", "\t\t\n", "\t\n", "
\n", "\t\t\t

Caveats
\n", "\t\t\t
\n", "\t\t\tNone

\n", "\t\t\t
\n", "
\n", " \n", "\n", "
\n", "

 

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.2" } }, "nbformat": 4, "nbformat_minor": 4 }