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.

Caveats

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

Last 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.
 

 

    1. 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

    2. 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
       

    3. 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.

 

 

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 we 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 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}0466_0104860501_SCX00000SUM.SAS', f'workdir={work_dir}']
In [ ]:
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.
In [ ]:
eventfile = work_dir+'0466_0104860501_EPN_S003_ImagingEvts.ds'
In [ ]:
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.
In [ ]:
# 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
In [ ]:
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.
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)                               # 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")
In [ ]:
# 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
In [ ]:
# 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
In [ ]:
# Execute SAS task with parameters
              w(cmd, inargs).run()
In [ ]:
# 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.
In [ ]:
# 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
In [ ]:
# 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
In [ ]:
# Execute SAS task with parameters
              w(cmd, inargs).run()
In [ ]:
# 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
In [ ]:
# Execute SAS task with parameters
              w(cmd, inargs).run()
In [ ]:
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.

 

 

 

 

 

 

 

Caveats

None
</center>  

 

 

 

In [ ]: