SAS Thread - epic reprocessing in Python - XMM-Newton

How to reprocess ODFs to generate calibrated and concatenated EPIC event lists

Introduction

This thread illustrates how to reprocess Observation Data Files (ODFs) to obtain calibrated and concatenated event lists.

Expected Outcome

The user will obtain calibrated and concatenated event lists which can be directly used to generate scientific products (images, spectra, light curves) through the SAS tasks evselect or xmmselect.

SAS Tasks to be Used

Prerequisites

SAS Startup Thread

Useful Links

Alternatively to following this thread, the pipeline products contain calibrated and concatenated event lists. Both ODFs and pipeline products can be downloaded from the XMM-Newton Science Archive.

Caveats

Last Reviewed: 25 May 2023, for SAS v21

Last Updated: 15 March 2021

Procedure

Run the EPIC reduction meta-tasks.

For EPIC-MOS:
emproc
and for EPIC-pn:
epproc

That's it! The default values of these meta-tasks are appropriate for most practical cases. You may have a look at the next section in this thread to learn how to perform specific reduction sub-tasks using emproc or epproc.

The files produced by epproc are the following:

????_??????????_AttHk.ds, the reconstructed attitude file
????_??????????_EPN_????_01_Badpixels.ds, one table per reduced CCD containing the bad pixels
????_??????????_EPN_????_ImagingEvts.ds, the calibrated and concatenated event list, which shall be used as an input to extract scientific products via evselect or xmmselect.

The files produced by emproc are conceptually the same. The main difference in the naming convention is that the string EPN is replaced by EMOS1 and EMOS2 for each EPIC-MOS camera, respectively.

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
              import os.path
              from os import path
              import subprocess
              import numpy as np
              import matplotlib.pyplot as plt
              from astropy.io import fits
              from astropy.table import Table
              from matplotlib.colors import LogNorm

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

In [ ]:

inargs = []

In [ ]:

t = w('sasver', inargs)

In [ ]:

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 010486050 as an example and assume that there is already a directory (sas_file in the block below) with an existing CIF and 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}0466_0104860501_SCX00000SUM.SAS', f'workdir={work_dir}']

In [ ]:

w('startsas', inargs).run()

This process should have defined the SAS_CCF an SAS_ODF variables.

The following blocks will produce EPIC-pn and EPIC-MOS event files.

Run the SAS command epproc. The following blocks are written in such a way that if EPIC-pn event files already exist in the Working directory, epproc wont run again. To force epproc to run, simply delete the corresponding event files.

In [ ]:

# SAS Command
              cmd    = 'epproc'  # SAS task to be executed
              # Arguments of SAS Command
              inargs = []        # comma separated arguments for SAS task
              print("   SAS command to be executed: "+cmd+", with arguments; \n")
              inargs

In [ ]:

print("Running epproc ..... \n")
              # Check if epproc has already run. If it has, do not run again 
              exists = 0
              pnevt_list = []
              for root, dirs, files in os.walk("."):
              for filename in files:
              if (filename.find('EPN') != -1) and filename.endswith('ImagingEvts.ds'):
              pnevt_list.append(filename)
              exists = 1
              if exists:
              print(" > " + str(len(pnevt_list)) + " EPIC-pn event list found. Not running epproc again.\n")
              for x in pnevt_list:
              print("    " + x + "\n")
              print("..... OK")
              else:
              w(cmd,inargs).run()      # <<<<< Execute SAS task
              exists = 0
              pnevt_list = []
              for root, dirs, files in os.walk("."):
              for filename in files:
              if (filename.find('EPN') != -1) and filename.endswith('ImagingEvts.ds'):
              pnevt_list.append(filename)
              exists = 1
              if exists:
              print(" > " + str(len(pnevt_list)) + " EPIC-pn event list found after running epproc.\n")
              for x in pnevt_list:
              print("    " + x + "\n")
              print("..... OK")
              else:
              print("Something has gone wrong with epproc. I cant find any event list files after running. \n")

Run the SAS command emproc. The following blocks are written in such a way that if EPIC-MOS event files already exist in the Working directory, emproc wont run again. To force emproc to run, simply delete the event files.

In [ ]:

# SAS Command
              cmd    = 'emproc'  # SAS task to be executed
              # Arguments of SAS Command
              inargs = []        # comma separated arguments for SAS task
              print("   SAS command to be executed: "+cmd+", with arguments; \n")
              inargs

In [ ]:

print("Running emproc ..... \n")
              # Check if emproc has already run. If it has, do not run again 
              exists = 0
              m1evt_list = []
              m2evt_list = []
              for root, dirs, files in os.walk("."):
              for filename in files:
              if (filename.find('EMOS1') != -1) and filename.endswith('ImagingEvts.ds'):
              m1evt_list.append(filename)
              exists = 1
              if (filename.find('EMOS2') != -1) and filename.endswith('ImagingEvts.ds'):
              m2evt_list.append(filename)
              exists = 1
              if exists:
              print(" > " + str(len(m1evt_list)) + " EPIC-MOS1 event list found. Not running emproc again.\n")
              for x in m1evt_list:
              print("    " + x + "\n")
              print(" > " + str(len(m2evt_list)) + " EPIC-MOS2 event list found. Not running emproc again.\n")
              for x in m2evt_list:
              print("    " + x + "\n")
              print("..... OK")
              else:
              w(cmd,inargs).run()      # <<<<< Execute SAS task
              exists = 0
              m1evt_list = []
              m2evt_list = []
              for root, dirs, files in os.walk("."):
              for filename in files:
              if (filename.find('EMOS1') != -1) and filename.endswith('ImagingEvts.ds'):
              m1evt_list.append(filename)
              exists = 1
              if (filename.find('EMOS2') != -1) and filename.endswith('ImagingEvts.ds'):
              m2evt_list.append(filename)
              exists = 1
              if exists:
              print(" > " + str(len(m1evt_list)) + " EPIC-MOS1 event list found. Not running emproc again.\n")
              for x in m1evt_list:
              print("    " + x + "\n")
              print(" > " + str(len(m2evt_list)) + " EPIC-MOS2 event list found. Not running emproc again.\n")
              for x in m2evt_list:
              print("    " + x + "\n")
              print("..... OK")
              else:
              print("Something has gone wrong with emproc. I cant find any event list file. \n")

Visualize the contents of the event files just created. The next blocks can be used to visualize the contents of the newly created event files. A simple selection criteria is impossed for illustration purposes. The thread How to filter EPIC event lists for flaring particle background will show how to use more complex filtering expressions by using the SAS task evselect.

In [ ]:

# For display purposes only, define a minimum filtering criteria for EPIC-pn
              pn_pattern   = 4        # pattern selection
              pn_pi_min    = 300.     # Low energy range eV
              pn_pi_max    = 12000.   # High energy range eV
              pn_flag      = 0        # FLAG
              # For display purposes only, define a minimum filtering criteria for EPIC-MOS
              mos_pattern   = 12      # pattern selection
              mos_pi_min    = 300.    # Low energy range eV
              mos_pi_max    = 20000.  # High energy range eV
              mos_flag      = 0       # FLAG

In [ ]:

plt.figure(figsize=(15,20))
              pl=1
              evts=len(pnevt_list)+len(m1evt_list)+len(m2evt_list)
              if len(pnevt_list) >0:
              for x in pnevt_list:
              hdu_list = fits.open(x, memmap=True)
              evt_data = Table(hdu_list[1].data)
              mask = ((evt_data['PATTERN'] <= pn_pattern) &
              (evt_data['FLAG'] == pn_flag) &
              (evt_data['PI'] >= pn_pi_min) &
              (evt_data['PI'] <= pn_pi_max))
              print("Events in event file" + " " + x + ": " + str(len(evt_data)) + "\n")
              print("Events in filtered event file" + " " + x + ": " + str(np.sum(mask)) + "\n")
              # Create Events image        
              xmax=np.amax(evt_data['X'])
              xmin=np.amin(evt_data['X'])
              xmid=(xmax-xmin)/2.+xmin
              ymax=np.amax(evt_data['Y'])
              ymin=np.amin(evt_data['Y'])
              xbin_size=80
              ybin_size=80
              NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))
              plt.subplot(evts, 2, pl)
              img_zero_mpl = plt.hist2d(evt_data['X'], evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm())
              cbar = plt.colorbar(ticks=[10.,100.,1000.])
              cbar.ax.set_yticklabels(['10','100','1000'])
              plt.title(x)
              plt.xlabel('x')
              plt.ylabel('y')
              pl=pl+1
              # Create Filtered Events image
              xmax=np.amax(evt_data['X'][mask])
              xmin=np.amin(evt_data['X'][mask])
              xmid=(xmax-xmin)/2.+xmin
              ymax=np.amax(evt_data['Y'][mask])
              ymin=np.amin(evt_data['Y'][mask])
              xbin_size=80
              ybin_size=80
              NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))
              plt.subplot(evts, 2, pl)
              img_zero_mpl = plt.hist2d(evt_data['X'][mask], evt_data['Y'][mask], NBINS, cmap='GnBu', norm=LogNorm())
              cbar = plt.colorbar(ticks=[10.,100.,1000.])
              cbar.ax.set_yticklabels(['10','100','1000'])
              plt.title(x)
              plt.xlabel('x')
              plt.ylabel('y')
              txt=("PATTERN <= " + str(pn_pattern) +
              " : " + str(pn_pi_min) + " <= E(eV) <= " + str(pn_pi_max) +
              " : " + " FLAG = " + str(pn_flag))
              plt.text(xmid, ymin+0.1*(ymax-ymin), txt, ha='center')
              pl=pl+1
              hdu_list.close()
              if len(m1evt_list) >0:
              for x in m1evt_list:
              hdu_list = fits.open(x, memmap=True)
              evt_data = Table(hdu_list[1].data)
              mask = ((evt_data['PATTERN'] <= mos_pattern) &
              (evt_data['FLAG'] == mos_flag) &
              (evt_data['PI'] >= mos_pi_min) &
              (evt_data['PI'] <= mos_pi_max))
              print("Events in event file" + " " + x + ": " + str(len(evt_data)) + "\n")
              print("Events in filtered event file" + " " + x + ": " + str(np.sum(mask)) + "\n")
              # Create Events image            
              xmax=np.amax(evt_data['X'])
              xmin=np.amin(evt_data['X'])
              xmid=(xmax-xmin)/2.+xmin
              ymax=np.amax(evt_data['Y'])
              ymin=np.amin(evt_data['Y'])
              xbin_size=80
              ybin_size=80
              NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))
              plt.subplot(evts, 2, pl)
              img_zero_mpl = plt.hist2d(evt_data['X'], evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm())
              cbar = plt.colorbar(ticks=[10.,100.,1000.])
              cbar.ax.set_yticklabels(['10','100','1000'])
              plt.title(x)
              plt.xlabel('x')
              plt.ylabel('y')
              pl=pl+1
              # Create Filtered Events image
              xmax=np.amax(evt_data['X'][mask])
              xmin=np.amin(evt_data['X'][mask])
              xmid=(xmax-xmin)/2.+xmin
              ymax=np.amax(evt_data['Y'][mask])
              ymin=np.amin(evt_data['Y'][mask])
              xbin_size=80
              ybin_size=80
              NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))
              plt.subplot(evts, 2, pl)
              img_zero_mpl = plt.hist2d(evt_data['X'][mask], evt_data['Y'][mask], NBINS, cmap='GnBu', norm=LogNorm())
              cbar = plt.colorbar(ticks=[10.,100.,1000.])
              cbar.ax.set_yticklabels(['10','100','1000'])
              plt.title(x)
              plt.xlabel('x')
              plt.ylabel('y')
              txt=("PATTERN <= " + str(mos_pattern) +
              " : " + str(mos_pi_min) + " <= E(eV) <= " + str(mos_pi_max) +
              " : " + " FLAG = " + str(mos_flag))
              plt.text(xmid, ymin+0.1*(ymax-ymin), txt, ha='center')
              pl=pl+1
              hdu_list.close()
              if len(m2evt_list) >0:
              for x in m2evt_list:
              hdu_list = fits.open(x, memmap=True)
              evt_data = Table(hdu_list[1].data)
              mask = ((evt_data['PATTERN'] <= mos_pattern) &
              (evt_data['FLAG'] == mos_flag) &
              (evt_data['PI'] >= mos_pi_min) &
              (evt_data['PI'] <= mos_pi_max))
              print("Events in event file" + " " + x + ": " + str(len(evt_data)) + "\n")
              print("Events in filtered event file" + " " + x + ": " + str(np.sum(mask)) + "\n")
              # Create Events image
              xmax=np.amax(evt_data['X'])
              xmin=np.amin(evt_data['X'])
              xmid=(xmax-xmin)/2.+xmin
              ymax=np.amax(evt_data['Y'])
              ymin=np.amin(evt_data['Y'])
              xbin_size=80
              ybin_size=80
              NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))
              plt.subplot(evts, 2, pl)
              img_zero_mpl = plt.hist2d(evt_data['X'], evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm())
              cbar = plt.colorbar(ticks=[10.,100.,1000.])
              cbar.ax.set_yticklabels(['10','100','1000'])
              plt.title(x)
              plt.xlabel('x')
              plt.ylabel('y')
              pl=pl+1
              # Create Filtered Events image    
              xmax=np.amax(evt_data['X'][mask])
              xmin=np.amin(evt_data['X'][mask])
              xmid=(xmax-xmin)/2.+xmin
              ymax=np.amax(evt_data['Y'][mask])
              ymin=np.amin(evt_data['Y'][mask])
              xbin_size=80
              ybin_size=80
              NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))
              plt.subplot(evts, 2, pl)
              img_zero_mpl = plt.hist2d(evt_data['X'][mask], evt_data['Y'][mask], NBINS, cmap='GnBu', norm=LogNorm())
              cbar = plt.colorbar(ticks=[10.,100.,1000.])
              cbar.ax.set_yticklabels(['10','100','1000'])
              plt.title(x)
              plt.xlabel('x')
              plt.ylabel('y')
              txt=("PATTERN <= " + str(mos_pattern) +
              " : " + str(mos_pi_min) + " <= E(eV) <= " + str(mos_pi_max) +
              " : " + " FLAG = " + str(mos_flag))
              plt.text(xmid, ymin+0.1*(ymax-ymin), txt, ha='center')
              pl=pl+1
              hdu_list.close()

How to accomplish specific reduction tasks

emproc and epproc are highly flexible tasks, which allow the user to perform a wide range of customized reduction tasks. Some emproc examples are listed below. The same customized reduction tasks can be performed for the EPIC-pn as well, just by substituting emproc with epproc in the commands.

if you want to reduce only one of the cameras (EPIC-MOS1 in the example):

emproc selectinstruments=yes emos1=yes

In [ ]:

# SAS Command
              cmd    = "emproc" # SAS task to be executed                  
              # Arguments of SAS Command
              inargs = ['selectinstruments=yes','emos1=yes']
              print("   SAS command to be executed: "+cmd+", with arguments; \n")
              inargs

In [ ]:

w(cmd, inargs).run()

if you want to reduce only a subsample of exposures:

emproc withinstexpids=yes instexpids="M1S001 M2S002"

In [ ]:

# SAS Command
              cmd    = "emproc" # SAS task to be executed                  
              # Arguments of SAS Command
              inargs = ['withinstexpids=yes','instexpids=M1S001 M2S002']
              print("   SAS command to be executed: "+cmd+", with arguments; \n")
              inargs

In [ ]:

w(cmd, inargs).run()

if you want to reduce data from 1 CCD only (#4 and #5 in the example):

emproc selectccds=yes ccd4=yes ccd5=yes

In [ ]:

# SAS Command
              cmd    = "emproc" # SAS task to be executed                  
              # Arguments of SAS Command
              inargs = ['selectccds=yes','ccd4=yes','ccd5=yes']
              print("   SAS command to be executed: "+cmd+", with arguments; \n")
              inargs

In [ ]:

w(cmd, inargs).run()

if you want to change the reference pointing for the calculation of the sky coordinates to a value of your choice:

emproc referencepointing=user ra=34.65646 dec=-12.876546

In [ ]:

# SAS Command
              cmd    = "emproc" # SAS task to be executed                  
              # Arguments of SAS Command
              inargs = ['referencepointing=user','ra=34.65646','dec=-12.876546']
              print("   SAS command to be executed: "+cmd+", with arguments; \n")
              inargs

In [ ]:

w(cmd, inargs).run()

Please be aware that if you want to supply coordinates for the analysis of the EPIC-MOS Timing mode, the command is slightly different, e.g.:

emproc withsrccoords=yes srcra=34.65646 srcdec=-12.876546

if you want to filter the event list events, using an external Good Time Interval (GTI) file (see the corresponding thread on how to filter event files for flaring particle background by creating a GTI file):

emproc withgtiset=yes gtiset=mygti.gti filterevents=yes

In [ ]:

# SAS Command
              cmd    = "emproc" # SAS task to be executed                  
              # Arguments of SAS Command
              inargs = ['withgtiset=yes','gtiset=mygti.gti','filterevents=yes']
              print("   SAS command to be executed: "+cmd+", with arguments; \n")
              inargs

In [ ]:

w(cmd, inargs).run()

Parameters can be combined to accomplish simultaneously two or more of the above tasks during the same run.

The user is referred to the on-line documentation of emproc and epproc for a complete list of the available options.

Reduction of EPIC-pn Timing Mode exposures

Most exposures in EPIC-pn Timing Mode are affected by X-ray Loading (XRL; cf. Sect.3.1 in Guainazzi et al., 2013, XMM-SOC-CAL-TN-0083). Furthermore, a residual dependence of the energy scale on the total count rate is corrected through the "Rate-Dependent PHA" correction (Guainazzi, 2014, XMM-CCF-REL-312). In order to correct for these effects a set of default calibration settings have been identified. As of SAS v14.0, this is controlled by a single parameter within the tasks epproc and epchain. This parameter is called withdefaultcal and is set to yes by default. Setting withdefaultcal=yes implies runepreject=yes withxrlcorrection=yes runepfast=no withrdpha=yes. So one shall run the EPIC-pn reduction meta-tasks as follows:

epproc

or:

epchain datamode=TIMING

For more information please refer to the documentation of epproc and epchain.

Reduction of EPIC-pn Burst Mode exposures

Most exposures in EPIC-pn Burst Mode are affected by X-ray Loading (XRL; cf. Sect.3.1 in Guainazzi et al., 2013, XMM-SOC-CAL-TN-0083). Furthermore, a residual dependence of the energy scale on the total count rate is corrected through the "Rate-Dependent CTI" correction. In order to correct for these effects a set of default calibration settings have been identified. As of SAS v14.0, this is controlled by a single parameter within the tasks epproc and epchain. This parameter is called withdefaultcal and is set to yes by default. Setting withdefaultcal=yes implies runepreject=yes withxrlcorrection=yes runepfast=yes withrdpha=no. So one shall run the EPIC-pn reduction meta-tasks as follows:

epproc burst=yes

Notice the inclusion of the extra parameter burst=yes in the call to epproc. The meta-task epchain also needs an extra parameter:

epchain datamode=BURST

In [ ]:

# SAS Command
              cmd    = "epchain" # SAS task to be executed                  
              # Arguments of SAS Command
              inargs = ['datamode=BURST']
              print("   SAS command to be executed: "+cmd+", with arguments; \n")
              inargs

In [ ]:

w(cmd,inargs).run()

Should I indeed bother to regenerate my event lists?

Concatenated and calibrated EPIC event lists are already available in the PPS Pipeline Products. These are produced with the most updated software and calibrations available at the moment the observation was performed and - normally a couple of weeks later - the ODF are produced. There is therefore no need for you to regenerate the EPIC event lists yourself, unless substantial changes in the software and/or calibration occurred between the time when the Pipeline Products were generated and the moment when you are analyzing the data. Of course, there is no general answer or recipe one can apply to decide if/when this is the case. In order to collect all the available and necessary elements to formulate your judgment, follow these steps:

check in the XMM-Newton Science Archive the SAS version and the date when the Calibration Index File was generated to produce the Pipeline Products;
verify in the SAS release notes the software changes between the version employed to produce the Pipeline Products and the latest SAS version;
verify in the CCF release notes the changes in the calibration files between the moment when Pipeline Products were generated and the moment your analysis is being performed.

If you feel lost or unsure, and prefer to stay on the safe side, you had probably better regenerate your EPIC calibrated event lists.

Caveats

None