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

How to filter EPIC event lists for flaring particle background

\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 describes how to filter an EPIC event list for periods of high background flaring activity.

\n", "\t\t\t

Expected Outcome

\n", "\t\t\t

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.

\n", "\t\t\t

SAS Tasks to be Used

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

Prerequisites

\n", "\t\t\t

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.

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

Useful Links

\n", "\t\t\t

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.

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

Caveats

\n", "\t\t\t

Last Reviewed: 31 January 2025, for SAS v22.0

\n", "\t\t\t

Last Updated: 31 January 2025

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

Procedure

\n", "

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.
\n", " 

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

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

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import the following Python libraries" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from pysas.wrapper import Wrapper as w\n", "import os\n", "import os.path\n", "from os import path\n", "import subprocess\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 matplotlib.colors import LogNorm" ] }, { "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 = []" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = w('sasver', inargs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "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 0104860501 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}0466_0104860501_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 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.

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eventfile = work_dir+'0466_0104860501_EPN_S003_ImagingEvts.ds'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Checking for EPIC-pn Event Files ..... \\n\")\n", "# Check if epproc has already run.\n", "if os.path.isfile(eventfile):\n", " print (\"File \"+eventfile+\" exists. \\n\")\n", "else:\n", " print (\"File \"+eventfile+\" does not exist, please check. \\n\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# For display purposes only define the following event cuts\n", "\n", "pn_pattern = 4 # pattern selection\n", "pn_pi_min = 300. # Low energy range eV\n", "pn_pi_max = 12000. # High energy range eV\n", "pn_flag = 0 # FLAG" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,8))\n", "\n", "pl=1\n", " \n", "hdu_list = fits.open(eventfile, memmap=True)\n", "evt_data = Table(hdu_list[1].data)\n", " \n", "mask = ((evt_data['PATTERN'] <= pn_pattern) & \n", " (evt_data['FLAG'] == pn_flag) & \n", " (evt_data['PI'] >= pn_pi_min) & \n", " (evt_data['PI'] <= pn_pi_max))\n", "\n", "print(\"Events in event file\" + \" \" + eventfile + \": \" + str(len(evt_data)) + \"\\n\")\n", "print(\"Events in filtered event file\" + \" \" + eventfile + \": \" + str(np.sum(mask)) + \"\\n\")\n", "print(\" Filter: PATTERN <= \" + str(pn_pattern) +\n", " \" : \" + str(pn_pi_min) + \" <= E(eV) <= \" + str(pn_pi_max) +\n", " \" : \" + \" FLAG == \" + str(pn_flag)+ \"\\n\")\n", "\n", "# Create Events image \n", "\n", "xmax=np.amax(evt_data['X']) \n", "xmin=np.amin(evt_data['X']) \n", "xmid=(xmax-xmin)/2.+xmin\n", "ymax=np.amax(evt_data['Y']) \n", "ymin=np.amin(evt_data['Y'])\n", "xbin_size=80\n", "ybin_size=80\n", "NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))\n", " \n", "plt.subplot(1, 2, pl)\n", "\n", "img_zero_mpl = plt.hist2d(evt_data['X'], evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm())\n", "\n", "cbar = plt.colorbar(ticks=[10.,100.,1000.])\n", "cbar.ax.set_yticklabels(['10','100','1000'])\n", " \n", "plt.title(eventfile)\n", "plt.xlabel('x')\n", "plt.ylabel('y')\n", " \n", "pl=pl+1\n", " \n", "# Create Filtered Events image\n", "\n", "xmax=np.amax(evt_data['X'][mask]) \n", "xmin=np.amin(evt_data['X'][mask]) \n", "xmid=(xmax-xmin)/2.+xmin\n", "ymax=np.amax(evt_data['Y'][mask]) \n", "ymin=np.amin(evt_data['Y'][mask])\n", "xbin_size=80\n", "ybin_size=80\n", "NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))\n", " \n", "plt.subplot(1, 2, pl)\n", "\n", "img_zero_mpl = plt.hist2d(evt_data['X'][mask], evt_data['Y'][mask], NBINS, cmap='GnBu', norm=LogNorm())\n", "\n", "cbar = plt.colorbar(ticks=[10.,100.,1000.])\n", "cbar.ax.set_yticklabels(['10','100','1000'])\n", " \n", "plt.title(eventfile)\n", "plt.xlabel('x')\n", "plt.ylabel('y')\n", "\n", "txt=(\"PATTERN <= \" + str(pn_pattern) +\n", " \" : \" + str(pn_pi_min) + \" <= E(eV) <= \" + str(pn_pi_max) +\n", " \" : \" + \" FLAG == \" + str(pn_flag))\n", "plt.text(xmid, ymin+0.1*(ymax-ymin), txt, ha='center')\n", " \n", "pl=pl+1\n", "\n", "hdu_list.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "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) # plot the data\n", "\n", " if colName == 'RATE':\n", " plt.title(\"Flaring particle background 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": [ "# Define a SAS filter expression to derive a background rate cut\n", "\n", "pn_pattern = 0 # pattern selection\n", "pn_pi_min = 10000. # Low energy range eV\n", "pn_pi_max = 12000. # High energy range eV\n", "pn_threshold = 0.4 # cts/sec (only used here for display purposes)\n", "\n", "out_LCFile = work_dir+'EPIC_pn_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'#XMMEA_EP&&(PI>={pn_pi_min}&&PI<={pn_pi_max})&&(PATTERN=={pn_pattern})' # event filter expression\n", "inargs = [f'table={eventfile}','withrateset=Y',f'rateset={out_LCFile}','maketimecolumn=Y','timebinsize=100','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 SAS task with parameters\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Open event file\n", "\n", "hdu_list = fits.open(eventfile, memmap=True)\n", "evt_data = Table(hdu_list[1].data) \n", "prihdu = hdu_list[1].header\n", "\n", "mask = ((evt_data['PATTERN'] <= pn_pattern) & \n", " (evt_data['FLAG'] == pn_flag) & \n", " (evt_data['PI'] >= pn_pi_min) & \n", " (evt_data['PI'] <= pn_pi_max))\n", " \n", "# Read some information from keywords to be used later on\n", "\n", "if ('INSTRUME' in prihdu):\n", " ins = prihdu['INSTRUME']\n", " print(\"Looking into instrument: \"+ins+\" \\n\")\n", "if ('EXPIDSTR' in prihdu):\n", " expid = prihdu['EXPIDSTR']\n", " print(\"Looking at exposure: \"+expid+\" \\n\") \n", "\n", "# Check number of event in initial event file\n", "\n", "print(\"Events in event file\" + \" \" + eventfile + \": \" + str(len(evt_data)) + \"\\n\")\n", "print(\"Events in filtered event file\" + \" \" + eventfile + \": \" + str(np.sum(mask)) + \"\\n\")\n", "print(\" Filter: PATTERN <= \" + str(pn_pattern) +\n", " \" : \" + str(pn_pi_min) + \" <= E(eV) <= \" + str(pn_pi_max) +\n", " \" : \" + \" FLAG == \" + str(pn_flag)+ \"\\n\")\n", "\n", "# Create events image and background lightcurve\n", "\n", "plt.figure(figsize=(20,8))\n", "\n", "pl=1\n", " \n", "xmax=np.amax(evt_data['X'][mask]) \n", "xmin=np.amin(evt_data['X'][mask]) \n", "xmid=(xmax-xmin)/2.+xmin\n", "ymax=np.amax(evt_data['Y'][mask]) \n", "ymin=np.amin(evt_data['Y'][mask])\n", "xbin_size=80\n", "ybin_size=80\n", "NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))\n", " \n", "# Plot image\n", " \n", "plt.subplot(1, 2, pl)\n", "\n", "img_zero_mpl = plt.hist2d(evt_data['X'][mask], evt_data['Y'][mask], NBINS, cmap='GnBu', norm=LogNorm())\n", "\n", "cbar = plt.colorbar(ticks=[10.,100.,1000.])\n", "cbar.ax.set_yticklabels(['10','100','1000'])\n", " \n", "plt.title(eventfile)\n", "plt.xlabel('x')\n", "plt.ylabel('y')\n", "\n", "plt.text(xmid, ymin+0.1*(ymax-ymin), expression, ha='center') \n", " \n", "pl=pl+1\n", "plt.subplot(1, 2, pl)\n", " \n", "# Plot BKG lightcurve\n", " \n", "plotLC(plt,pn_threshold,out_LCFile)\n", " \n", "pl=pl+1\n", " \n", "hdu_list.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define energy range to filter event file\n", "\n", "pn_pattern = 4 # pattern selection\n", "pn_pi_min = 200. # Low energy range eV\n", "pn_pi_max = 10000. # High energy range eV\n", "pn_threshold = 0.4 # Cut to be applied to filter event file (cts/sec)\n", "\n", "# Define the input and output file names\n", "\n", "in_LCFile = work_dir+'EPIC_pn_FlareBKGRate.fit' # Name of the input BKG lightcurve\n", "out_gti_set = work_dir+'EPIC_pn_gti.fit' # Name of the output file containing GTI intervals\n", "out_clean_evtFile = work_dir+'EPIC_pn_gtiFilteredEvts.ds' # Name of the output Event file filtered by GTI" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"tabgtigen\" \n", "\n", "# Arguments of SAS Command\n", "expression = 'RATE<='+str(pn_threshold) # event filter expression\n", "inargs = [f'table={in_LCFile}',f'gtiset={out_gti_set}',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 SAS task with parameters\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SAS Command\n", "cmd = \"evselect\" \n", "\n", "# Arguments of SAS Command\n", "expression = ('#XMMEA_EP&&FLAG==0&&(PI>='+str(pn_pi_min)+'&&PI<='+str(pn_pi_max)+\n", " ')&&(gti('+str(out_gti_set)+',TIME))')\n", "inargs = [f'table={eventfile}','withfilteredset=Y',f'filteredset={out_clean_evtFile}',\n", " 'destruct=Y','keepfilteroutput=T',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 SAS task with parameters\n", "\n", "w(cmd, inargs).run()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,8))\n", "\n", "pl=1\n", "\n", "hdu_list = fits.open(eventfile, memmap=True)\n", "evt_data = Table(hdu_list[1].data) \n", "prihdu = hdu_list[1].header\n", "print(\"Events in event file\" + \" \" + eventfile + \": \" + str(len(evt_data)) + \"\\n\")\n", " \n", "gti_hdu_list = fits.open(out_clean_evtFile, memmap=True)\n", "gti_evt_data = Table(gti_hdu_list[1].data) \n", "print(\"Events in GTI clean event file\" + \" \" + out_clean_evtFile + \": \" + str(len(gti_evt_data)) + \"\\n\")\n", " \n", "# Create Events image\n", " \n", "xmax=np.amax(evt_data['X']) \n", "xmin=np.amin(evt_data['X']) \n", "xmid=(xmax-xmin)/2.+xmin\n", "ymax=np.amax(evt_data['Y']) \n", "ymin=np.amin(evt_data['Y'])\n", "xbin_size=80\n", "ybin_size=80\n", "NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))\n", " \n", "plt.subplot(1, 2, pl)\n", "\n", "img_zero_mpl = plt.hist2d(evt_data['X'], evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm())\n", "\n", "cbar = plt.colorbar(ticks=[10.,100.,1000.])\n", "cbar.ax.set_yticklabels(['10','100','1000'])\n", " \n", "plt.title(eventfile)\n", "plt.xlabel('x')\n", "plt.ylabel('y')\n", " \n", "pl=pl+1\n", " \n", "# Create Clean Events image\n", " \n", "xmax=np.amax(gti_evt_data['X']) \n", "xmin=np.amin(gti_evt_data['X']) \n", "xmid=(xmax-xmin)/2.+xmin\n", "ymax=np.amax(gti_evt_data['Y']) \n", "ymin=np.amin(gti_evt_data['Y'])\n", "xbin_size=80\n", "ybin_size=80\n", "NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))\n", " \n", "plt.subplot(1, 2, pl)\n", "\n", "img_zero_mpl = plt.hist2d(gti_evt_data['X'], gti_evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm())\n", "\n", "cbar = plt.colorbar(ticks=[10.,100.,1000.])\n", "cbar.ax.set_yticklabels(['10','100','1000'])\n", " \n", "plt.title(out_clean_evtFile)\n", "plt.xlabel('x')\n", "plt.ylabel('y')\n", " \n", "plt.text(xmid, ymin-0.1*(ymax-ymin), expression, ha='center')\n", " \n", "pl=pl+1\n", " \n", "gti_hdu_list.close()\n", "hdu_list.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After running this thread three files are produced,\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

An alternative way to produce an EPIC filtered event list and corresponding GTI file

\n", "\n", "

 

\n", "\n", "

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

\n", "\n", "

In its simplest form, the task can be ran by providing the name of the EPIC event list in the following way:
\n", "
\n", "   espfilt eventfile=EPIC.fits

\n", "\n", "

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:
\n", "
\n", "   ****-allevc.FIT, contains the EPIC clean event list
\n", "   ****-gti.FIT, contains the definition of GTI
\n", "   ****-hist.qdp, contains QDP text file of histograms and lightcurves
\n", "
\n", "where **** is a combination of instrument and exposure IDs. More files are produced by espfilt and the user is referred to the documentation of the SAS task. 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
\n", "
\n", "Keep in mind that for the case of pn Small Window exposures, there are no corner events, hence, the task does not support this method. For further information please read the documentation of espfilt.

\n" ] }, { "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\tAs of SAS V22.0, the SAS task espfilt, does not apply any filtering criteria, like FLAG or PATTERN, to the output EPIC clean event list, other than filtering by the created GTI file. It is up to the users to filter further the event list according to their needs using, for example, the SAS task evselect For further information please check the documentation of espfilt.\n", "\t\t\t

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

 

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "PythonSAS", "language": "python", "name": "localpython" }, "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.10.6" } }, "nbformat": 4, "nbformat_minor": 4 }