SAS Thread - epic filterbackground in Python - XMM-Newton
How to filter EPIC event lists for flaring particle background
Introduction This thread describes how to filter an EPIC event list for periods of high background flaring activity. Expected Outcome The outcome of this thread is an EPIC filtered event list clean of time intervals of high background activity, and a Good Time Interval (GTI) file containing the definition of GTI for a given observation. SAS Tasks to be Used Prerequisites It is assumed that the processed event lists, like the ones produced by the SAS tasks epproc and emproc, and corresponding attitude and summary files as well as the ccf.cif file are available. Calibrated event lists may also be obtained as a starting point for this thread from the archive.
You might have to install the python package, beautifultable. For that, run in one of the blocks at the very beginning of the notebook the following command, !pip install beautifultable Useful Links In the following links, a full description is given of the different background components of the EPIC cameras, including the external flaring background component which is treated in this thread.
Last Reviewed: 25 May 2023, for SAS v21.0Last Updated: 15 March 2021 |
Procedure
As an example case, we will consider cleaning an event list, EPIC.fits, higlighting the differences between MOS and pn. Below are the main general steps which are included in the notebook cells that follow.
-
- Extract a single event (i.e. pattern zero only), high energy light curve, from the event file to identify intervals of flaring particle background:
evselect table=EPIC.fits withrateset=Y rateset=rateEPIC.fits \
maketimecolumn=Y timebinsize=100 makeratecolumn=Y \
expression='Selection_Expression'
where Selection_Expression is:
#XMMEA_EM && (PI>10000) && (PATTERN==0) for EPIC-MOS
#XMMEA_EP && (PI>10000&&PI<12000) && (PATTERN==0) for EPIC-pn
Notice how the EPIC-pn energy range selected for producing the background light curve includes events only up to 12 keV. The reason for this is to avoid hot pixels being miss-identified as very high energy events. One can look at the produced light curve by running:
dsplot table=rateEPIC.fits x=TIME y=RATE
- Determine where the light curve is low and steady. Choose a threshold (count/second) just above the low steady background to define the "low background" intervals, to create the corresponding GTI file:
tabgtigen table=rateEPIC.fits expression='RATE<=Rate_Expression' gtiset=EPICgti.fits
where Rate_Expression is:
RATE<=0.35 for EPIC-MOS
RATE<=0.4 for EPIC-pn
This value will be different for each data set (the value given here for each instrument represents a good reference value for a standard observation). Use this expression to filter the event list:
evselect table=EPIC.fits withfilteredset=Y filteredset=EPICclean.fits \
destruct=Y keepfilteroutput=T \
expression='Selection_Expression'
where Selection_Expression is:
#XMMEA_EM && gti(EPICgti.fits,TIME) && (PI>150) for EPIC-MOS
#XMMEA_EP && gti(EPICgti.fits,TIME) && (PI>150) for EPIC-pn
- The file EPICclean.fits contains the filtered EPIC event list and can be now used to produce scientific products. The file EPICgti.fits contains the definition of Good Time Intervals.
- Extract a single event (i.e. pattern zero only), high energy light curve, from the event file to identify intervals of flaring particle background:
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
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 we have already defined in terms of SAS variables,
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 010486050 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 '/'.
work_dir = 'absolute_path_to_wrk_directory'
sas_file = 'absolute_path_to_cifSUMSAS_directory'
inargs = [f'sas_ccf={sas_file}ccf.cif', f'sas_odf={sas_file}0466_0104860501_SCX00000SUM.SAS', f'workdir={work_dir}']
w('startsas', inargs).run()
This process should have defined the SAS_CCF and SAS_ODF variables.
In the following blocks we will filter an EPIC-pn event file for periods of high background flaring activity as an example. The same procedure is valid as well for an EPIC-MOS event file.
- Introduce the EPIC-pn event file to process including the absolute path. If no event files are available, please follow the How to reprocess ODFs to generate calibrated and concatenated EPIC event lists thread first.
eventfile = work_dir+'0466_0104860501_EPN_S003_ImagingEvts.ds'
print("Checking for EPIC-pn Event Files ..... \n")
# Check if epproc has already run.
if os.path.isfile(eventfile):
print ("File "+eventfile+" exists. \n")
else:
print ("File "+eventfile+" does not exist, please check. \n")
- Events can be filtered according to different criteria. Cuts can be impossed in, for example, event pattern, energy (PI) and quality flag. The following two blocks allow to impose some event cuts and demonstrate the effects of applying these cuts by visualizing an image before and after filtering.
# For display purposes only define the following event cuts
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
plt.figure(figsize=(20,8))
pl=1
hdu_list = fits.open(eventfile, 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" + " " + eventfile + ": " + str(len(evt_data)) + "\n")
print("Events in filtered event file" + " " + eventfile + ": " + str(np.sum(mask)) + "\n")
print(" Filter: PATTERN <= " + str(pn_pattern) +
" : " + str(pn_pi_min) + " <= E(eV) <= " + str(pn_pi_max) +
" : " + " FLAG == " + str(pn_flag)+ "\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(1, 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(eventfile)
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(1, 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(eventfile)
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()
- Extract a single event (i.e. pattern zero only), high energy light curve, to identify intervals of flaring particle background.
- Notice how the EPIC-pn energy range selected for producing the background light curve includes events only up to 12 keV. The reason for this is to avoid hot pixels being miss-identified as very high energy events (for EPIC-MOS the cut can include events with energies up to 20 keV).
- The value for the threshold in rate will be different for each data set (for EPIC-pn a good reference value is 0.4 cts/sec and for EPIC-MOS a value of 0.35 cts/sec). The value given here represents a good reference value for a standard observation.
- One can look at the produced background light curve by running the next blocks. The next blocks also show how to use more complex filtering expressions by using the SAS task evselect.
# 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) # plot the data
if colName == 'RATE':
plt.title("Flaring particle background 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")
# Define a SAS filter expression to derive a background rate cut
pn_pattern = 0 # pattern selection
pn_pi_min = 10000. # Low energy range eV
pn_pi_max = 12000. # High energy range eV
pn_threshold = 0.4 # cts/sec (only used here for display purposes)
out_LCFile = work_dir+'EPIC_pn_FlareBKGRate.fit' # Name of the output BKG lightcurve
# SAS Command
cmd = "evselect" # SAS task to be executed
# Arguments of SAS Command
expression = f'#XMMEA_EP&&(PI>={pn_pi_min}&&PI<={pn_pi_max})&&(PATTERN=={pn_pattern})' # event filter expression
inargs = [f'table={eventfile}','withrateset=Y',f'rateset={out_LCFile}','maketimecolumn=Y','timebinsize=100','makeratecolumn=Y',f'expression={expression}']
print(" Filter expression to use: "+expression+" \n")
print(" SAS command to be executed: "+cmd+", with arguments; \n")
inargs
# Execute SAS task with parameters
w(cmd, inargs).run()
# Open event file
hdu_list = fits.open(eventfile, memmap=True)
evt_data = Table(hdu_list[1].data)
prihdu = hdu_list[1].header
# Read some information from keywords to be used later on
if ('INSTRUME' in prihdu):
ins = prihdu['INSTRUME']
print("Looking into instrument: "+ins+" \n")
if ('EXPIDSTR' in prihdu):
expid = prihdu['EXPIDSTR']
print("Looking at exposure: "+expid+" \n")
# Check number of event in initial event file
print("Events in event file" + " " + eventfile + ": " + str(len(evt_data)) + "\n")
# Create events image and background lightcurve
plt.figure(figsize=(20,8))
pl=1
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))
# Plot image
plt.subplot(1, 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(eventfile)
plt.xlabel('x')
plt.ylabel('y')
plt.text(xmid, ymin+0.1*(ymax-ymin), expression, ha='center')
pl=pl+1
plt.subplot(1, 2, pl)
# Plot BKG lightcurve
plotLC(plt,pn_threshold,out_LCFile)
pl=pl+1
hdu_list.close()
- Determine where the light curve is low and steady. Choose a threshold (count/second) just above the low steady background to define the "low background" intervals, to create the corresponding GTI file.
- As already mentioned, this value will be different for each data set.
# Define energy range to filter event file
pn_pattern = 4 # pattern selection
pn_pi_min = 200. # Low energy range eV
pn_pi_max = 10000. # High energy range eV
pn_threshold = 0.4 # Cut to be applied to filter event file (cts/sec)
# Define the input and output file names
in_LCFile = work_dir+'EPIC_pn_FlareBKGRate.fit' # Name of the output BKG lightcurve
out_gti_set = work_dir+'EPIC_pn_gti.fit' # Name of the output file containing GTI intervals
out_clean_evtFile = work_dir+'EPIC_pn_gtiFilteredEvts.ds' # Name of the output Event file filtered by GTI
# SAS Command
cmd = "tabgtigen"
# Arguments of SAS Command
expression = 'RATE<='+str(pn_threshold) # event filter expression
inargs = [f'table={in_LCFile}',f'gtiset={out_gti_set}',f'expression={expression}']
print(" Filter expression to use: "+expression+" \n")
print(" SAS command to be executed: "+cmd+", with arguments; \n")
inargs
# Execute SAS task with parameters
w(cmd, inargs).run()
# SAS Command
cmd = "evselect"
# Arguments of SAS Command
expression = ('#XMMEA_EP&&FLAG==0&&(PI>='+str(pn_pi_min)+'&&PI<='+str(pn_pi_max)+
')&&(gti('+str(out_gti_set)+',TIME))')
inargs = [f'table={eventfile}','withfilteredset=Y',f'filteredset={out_clean_evtFile}',
'destruct=Y','keepfilteroutput=T',f'expression={expression}']
print(" Filter expression to use: "+expression+" \n")
print(" SAS command to be executed: "+cmd+", with arguments; \n")
inargs
# Execute SAS task with parameters
w(cmd, inargs).run()
plt.figure(figsize=(20,8))
pl=1
hdu_list = fits.open(eventfile, memmap=True)
evt_data = Table(hdu_list[1].data)
prihdu = hdu_list[1].header
print("Events in event file" + " " + eventfile + ": " + str(len(evt_data)) + "\n")
gti_hdu_list = fits.open(out_clean_evtFile, memmap=True)
gti_evt_data = Table(gti_hdu_list[1].data)
print("Events in GTI clean event file" + " " + out_clean_evtFile + ": " + str(len(gti_evt_data)) + "\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(1, 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(eventfile)
plt.xlabel('x')
plt.ylabel('y')
pl=pl+1
# Create Clean Events image
xmax=np.amax(gti_evt_data['X'])
xmin=np.amin(gti_evt_data['X'])
xmid=(xmax-xmin)/2.+xmin
ymax=np.amax(gti_evt_data['Y'])
ymin=np.amin(gti_evt_data['Y'])
xbin_size=80
ybin_size=80
NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))
plt.subplot(1, 2, pl)
img_zero_mpl = plt.hist2d(gti_evt_data['X'], gti_evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm())
cbar = plt.colorbar(ticks=[10.,100.,1000.])
cbar.ax.set_yticklabels(['10','100','1000'])
plt.title(out_clean_evtFile)
plt.xlabel('x')
plt.ylabel('y')
plt.text(xmid, ymin-0.1*(ymax-ymin), expression, ha='center')
pl=pl+1
gti_hdu_list.close()
hdu_list.close()
After running this thread three files are produced,
- The file *FlareBKGRate.fit contain the background lightcurve used to define Good Time Intervals.
- The file *gtiFilteredEvts.ds contain the filtered EPIC event list and can be now used to produce scientific products.
- The file *gti.fit contain the definition of Good Time Intervals.
An alternative way to produce an EPIC filtered event list and corresponding GTI file
The SAS task espfilt, one of the tools within the XMM-ESAS package (available as of SAS v9.0), is able to perform GTI filtering of EPIC event files. The task generates soft proton contamination-filtered products, including clean event lists and GTI files (refer to the espfilt documentation for a description of all the products generated and information on the task input parameters).
In its simplest form, the task can be ran by providing the name of the EPIC event list in the following way:
espfilt eventset=EPIC.fits
Relevant to this thread, three files are produced which contain the EPIC clean event list, GTI file and QDP text file of histograms and lightcurves. These files use the following naming convention:
****-objevlifilt.FIT, contains the EPIC clean event list
****-gti.FIT, contains the definition of GTI
****-hist.qdp, contains QDP text file of histograms and lightcurves
where **** is a combination of the observation and exposure IDs. The QDP text file should be examined as it provides an indication of the quality of the data. It can be plotted using the command, for example, qdp ****-hist.qdp
Keep in mind that for the case of pn Small Window exposures, there are no corner events, hence, the products related to corner events will not be included as part of the output.