SAS Thread - epic SP LC for a point source in Python - XMM-Newton

How to extract a light curve and spectrum for an EPIC point-like source

Introduction

This thread contains a step-by-step recipe to extract a light curve and a spectrum of a point-like source which is valid for all the EPIC X-ray cameras. In the case of the light curve, it includes subtracting the background and correcting the light curve for exposure losses.

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 extract light curves and spectra for an EPIC point-like source using SAS, please refer to the standard SAS threads,

Expected Outcome

Spectrum and corrected light curve of a point-like source extracted from the XMM-Newton EPIC instruments.

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.

Useful Links

Caveats

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

Last Updated: 25 May 2023

Procedure

As an example case, we will consider the extraction of a spectrum and a light curve from a pn event list (PN_evt.fits). The same recipe applies for MOS. We highlight here the general steps which are included in the notebook cells below.

Set up your SAS environment (see Prerequisites for this thread at the top of the page).
Extract an image (in sky coordinates in this example; extraction in detector - DET[XY] - coordinates is possible as well).

evselect table=PN_evt.fits imagebinning=binSize imageset=PNimage.img withimageset=yes \ xcolumn=X ycolumn=Y ximagebinsize=80 yimagebinsize=80
Select the region in the image from which the spectrum and the light curve shall be accumulated. Lets assume its a circular region where the Centre coordinates are (27389.1,37183.40) and the radius 800. Units of sky coordinates (X,Y) are 0.05 arcsec, hence the radius in our example is 40 arcsec.
Repeat the previous step to determine the region from which the background for the spectrum and light curve is to be extracted. Lets assume its a circle, centered in (28855.1,40501.3) and radius also 800 pixels (have a look at the "EPIC status of calibration and data analysis" document (XMM-SOC-CAL-TN-0018) for the latest recommendations on how to select source and background regions).
Now you can extract a source+background and background spectrum and light curve, using the selection regions and including a quality selection appropriate for the extraction of a spectrum and light curve.
Light curve generation
1. Extract a source light curve, for PN, taking good events, singles and doubles with an energy range between 200 and 10000 eV (#XMMEA_EP && (PATTERN<=4) && (PI in [200:10000])); for MOS, taking good events, singles, doubles, triples and quadruples with an energy range between 200 and 10000 eV (#XMMEA_EM && (PATTERN<=12) && (PI in [200:10000])). In the example below, the bin size is 100 seconds,
  
  evselect table=PN_evt.fits energycolumn=PI \ expression='#XMMEA_EP&&(PATTERN<=4)&&((X,Y) IN circle(27389.1,37183.40,800))&&(PI in [200:10000])' \ withrateset=yes rateset="EPIC_PN_source_lightcurve.lc" timebinsize=100 \ maketimecolumn=yes makeratecolumn=yes
  
  The parameter makeratecolumn=yes produces a light curve in count rates (with errors). Otherwise the light curve is produced in counts (without errors).
2. Extract a background light curve, using all the selection expressions defined so far, and the same bin size (100 seconds) and energy range as for the source region,
  
  evselect table=PN_evt.fits energycolumn=PI \ expression='#XMMEA_EP&&(PATTERN<=4)&&((X,Y) IN circle(28855.1,40501.3,800))&&(PI in [200:10000])' \ withrateset=yes rateset="EPIC_PN_background_lightcurve.lc" timebinsize=100 \ maketimecolumn=yes makeratecolumn=yes
  
  The light curves are OGIP-complaint, and therefore can be processed by standard XRONOS-like LHEASOFT packages.
3. However, light curves obtained in such a way should be corrected for various effects affecting the detection efficiency, such as vignetting, bad pixels, PSF variations and quantum efficiency, as well as for variations affecting the stability of the detection within the exposure, like dead time and GTIs. Since all these effects can affect in a different manner source and background light curves, the background subtraction has to be done accordingly. A SAS task, epiclccorr, performs all of these corrections at once. It requires as input both light curves (which are used to establish the binning of the final corrected background subtracted light curve) and the event file. This is done with a simple command line call:
  
  epiclccorr srctslist=EPIC_PN_source_lightcurve.lc eventlist=PN_evt.fits \ outset=EPIC_PN_corrected_lightcurve.lc \ bkgtslist=EPIC_PN_background_lightcurve.lc withbkgset=yes applyabsolutecorrections=yes
Spectrum generation
1. Extract a source spectrum, for PN, taking good events, singles and doubles (#XMMEA_EP && (PATTERN<=4)); for MOS, taking good events, singles, doubles, triples and quadruples (#XMMEA_EM && (PATTERN<=12)). It is important to check withspecranges and use the following spectral bin range, specchannelmin=0 and specchannelmax=20479 (for MOS specchannelmin=0 and specchannelmax=11999), to accumulate the spectrum.
  
  evselect table=PN_evt.fits withspectrumset=yes \ spectrumset=EPIC_PN_source_spectrum.fits \ energycolumn=PI spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479 \ expression='(FLAG==0)&&(PATTERN<=4)&&((X,Y) IN circle(27389.1,37183.40,800))'
2. Extract a background spectrum, using all the selection expressions defined so far,
  
  evselect table=PN_evt.fits withspectrumset=yes \ spectrumset=EPIC_PN_background_spectrum.fits \ energycolumn=PI spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479 \ expression='(FLAG==0)&&(PATTERN<=4)&&((X,Y) IN circle(28855.1,40501.3,800))'
3. Calculate the area of source and background region used to make the spectral files. The area is written into the header of the SPECTRUM table of the file as keyword BACKSCAL.
  
  backscale spectrumset=EPIC_PN_source_spectrum.fits badpixlocation=PN_evt.fits
  backscale spectrumset=EPIC_PN_background_spectrum.fits badpixlocation=PN_evt.fits
4. Generate a redistribution matrix,
  
  rmfgen spectrumset=EPIC_PN_source_spectrum.fits rmfset=EPIC_PN.rmf
5. Generate an ancillary file,
  
  arfgen spectrumset=EPIC_PN_source_spectrum.fits arfset=EPIC_PN.arf \ withrmfset=yes rmfset=EPIC_PN.rmf \ badpixlocation=PN_evt.fits detmaptype=psf
6. Rebin the spectrum and link associated files. In the following example, the spectrum is rebinned in order to have at least 25 counts for each background-subtracted spectral channel and not to oversample the intrinsic energy resolution by a factor larger then 3.
  
  specgroup spectrumset=EPIC_PN_source_spectrum.fits \ mincounts=25 oversample=3 \ rmfset=EPIC_PN.rmf arfset=EPIC_PN.arf backgndset=EPIC_PN_background_spectrum.fits \ groupedset=EPIC_PN_spectrum_grp.fits
7. Fit the grouped spectrum in xspec.

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 matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
import pyds9
from xspec import *
import re

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

In [ ]:

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 '/'.

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.

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

Produce an image from this event file. The image is produced in sky coordinates (X,Y), but it can also be produced in detector coordinates (DETX,DETY).

In [ ]:

# Define some parameters to produce the image and the name of the output file

xbin=80     # xbin size
ybin=80     # ybin size
xcoord='X'  # coordinate system
ycoord='Y'  # coordinate system

out_IMFile   = work_dir+'EPIC_PN_Image.fit'  # Name of the output Image file

In [ ]:

# SAS Command
cmd        = "evselect" # SAS task to be executed                  

# Arguments of SAS Command
inargs     = [f'table={eventfile}','imagebinning=binSize',f'imageset={out_IMFile}','withimageset=yes',f'xcolumn={xcoord}',f'ycolumn={ycoord}',f'ximagebinsize={xbin}',f'yimagebinsize={ybin}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task with the parameters to produce an image

w(cmd, inargs).run()

The following block will open a ds9 window with the selected image, set the color map to bb and the scale to log.
This image can be used to select a source and background regions from which to extract a light curve for the source with its corresponding background.

In [ ]:

# Visualize the image with ds9

d = pyds9.DS9()
d.set("file "+out_IMFile)
d.set('cmap bb')
d.set('scale log')

Go to the ds9 window and select a source region using the tools provided within ds9. In this example, a circular region must be chosen (this can be controlled or changed in the following block). Do the same for the background region, using either circular or annulus shapes. The color of the source region must be changed to white in ds9. For the background region, use any other color except white, for example red, but the color must be also changed explicitly to include the #color identifier in the region information. The blocks below also accept an annulus for the background region. We will use this information to distinguish the source and background regions while parsing the information to Python. After this, execute the block below.
The following block should output the coordinates that have been extracted for the source and background region from the ds9 window and that will be used in the following blocks to extract the light curve and spectrum from the chosen source. Before proceeding with the notebook, double check that both regions have a #color identifier, and that the one from the source region is set to white.

In [ ]:

# Check the information that ds9 has about the regions defined in the ds9 window.

print(d.get("regions"))

In [ ]:

# Extract the relevant information from the ds9 regions.

region1=(re.split("circle|annulus",d.get("regions").partition("physical")[2]))[1].replace('(','').replace(')','').replace('#',',')
region2=(re.split("circle|annulus",d.get("regions").partition("physical")[2]))[2].replace('(','').replace(')','').replace('#',',')

print("Identified source region: ", region1)
print("Identified source region: ", region2)

In [ ]:

# Identify source and background regions using the 'white' color.

c1=region1.partition("color")[2].replace('=','').replace('\n','')
if(c1=='white'):
    regionsrc = region1
    regionbkg = region2
else:
    regionsrc = region2
    regionbkg = region1

# Save and print selected region coordinates.

x_source = regionsrc.split(",")[0].replace('\n','')
y_source = regionsrc.split(",")[1].replace('\n','')
r_source = regionsrc.split(",")[2].replace('\n','')
print("The coordinates of the selected source region are: \n")
print("  x_source = ", x_source, "(physical)")
print("  y_source = ", y_source, "(physical)")
print("  r_source = ", r_source, "(physical) \n")

x_bkg = regionbkg.split(",")[0].replace('\n','')
y_bkg = regionbkg.split(",")[1].replace('\n','')
r_bkg = regionbkg.split(",")[2].replace('\n','')
print("The coordinates of the selected background region are: \n")
print("  x_bkg = ", x_bkg, "(physical)")
print("  y_bkg = ", y_bkg, "(physical)")
print("  r_bkg = ", r_bkg, "(physical) \n")

# If the background is an annulus, save and print R2.

if "annulus" in str(d.get("regions")):
    r2_bkg = regionbkg.split(",")[3].replace('\n','')
    print("  r2_bkg = ", r2_bkg, "(physical)")

With both source and background regions defined, you can now start to extract the corresponding light curves and spectra.

Lightcurve Extraction

First, produce a light curve for the source region. In this example, we will produce a light curve for pn with a bin size of 100 seconds.
Use the source region and include a quality selection appropriate for a light curve extraction:
- For PN, take good events, singles and doubles with an energy range between 200 and 10000 eV (#XMMEA_EP && (PATTERN<=4) && (PI in [200:10000]));
- For MOS, take good events, singles, doubles, triples and quadruples with an energy range between 200 and 10000 eV (#XMMEA_EM && (PATTERN<=12) && (PI in [200:10000])).
NOTE: The parameter makeratecolumn=yes in the call to the SAS task evselect below produces a light curve in count rates (with errors). Otherwise, the light curve is produced in counts (without errors).

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,label=fileName) # plot the data

            if colName == 'RATE':
                plt.title("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 [ ]:

# Extract the source region lightcurve

# Define some parameters for filtering the event file and define the lightcurve binning

q_flag       = "#XMMEA_EP" # Quality flag for EPIC pn
n_pattern    = 4           # Pattern selection
pn_pi_min    = 200.        # Low energy range eV
pn_pi_max    = 10000.      # High energy range eV
lc_bin       = 100         # Lightcurve bin in secs

# Define the output ligthcurve file name

in_LCSRCFile = work_dir+'EPIC_PN_source_lightcurve.lc'   # Name of the output source lightcurve

In [ ]:

# SAS Command
cmd        = "evselect" # SAS task to be executed                  

# Arguments of SAS Command
expression = f'{q_flag}&&(PATTERN<={n_pattern})&&((X,Y) IN circle({x_source},{y_source},{r_source}))&&(PI in [{pn_pi_min}:{pn_pi_max}])'  # event filter expression
inargs     = [f'table={eventfile}','energycolumn=PI','withrateset=yes',f'rateset={in_LCSRCFile}',
              f'timebinsize={lc_bin}','maketimecolumn=yes','makeratecolumn=yes',f'expression={expression}']

print("   Filter expression to use: "+expression+" \n")
print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task with the parameters to produce the source region light curve

w(cmd, inargs).run()

Run the following cell to inspect the created source light curve.

In [ ]:

# Inspect light curve

pn_threshold = 1.0          # cts/sec (only used here for display purposes. Draws a horizontal line at this value)
plt.figure(figsize=(20,8))  # Size of figure

plotLC(plt,pn_threshold,in_LCSRCFile) # Plot source region light curve

plt.legend()

Second, produce a lightcurve for the background region, in the same energy range and using the same filtering criteria.

In [ ]:

# Extract the background region light curve

# Define some parameters for filtering the event file and define the lightcurve binning

q_flag       = "#XMMEA_EP" # Quality flag for EPIC pn
n_pattern    = 4           # Pattern selection
pn_pi_min    = 200.        # Low energy range eV
pn_pi_max    = 10000.      # High energy range eV
lc_bin       = 100         # Lightcurve bin in secs

# Define the output ligthcurve file name

in_LCBKGFile = work_dir+'EPIC_PN_background_lightcurve.lc'   # Name of the output background lightcurve

In [ ]:

# SAS Command
cmd        = "evselect" # SAS task to be executed                  

# Arguments of SAS Command
expression = f'{q_flag}&&(PATTERN<={n_pattern})&&((X,Y) IN circle({x_bkg},{y_bkg},{r_bkg}))&&(PI in [{pn_pi_min}:{pn_pi_max}])'  # event filter expression
inargs     = [f'table={eventfile}','energycolumn=PI','withrateset=yes',f'rateset={in_LCBKGFile}',
              f'timebinsize={lc_bin}','maketimecolumn=yes','makeratecolumn=yes',f'expression={expression}']

print("   Filter expression to use: "+expression+" \n")
print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task with the parameters to produce the background region light curve

w(cmd, inargs).run()

Run the following cell to inspect the created background light curve.

In [ ]:

# Inspect light curve

pn_threshold = 0.01         # cts/sec (only used here for display purposes. Draws a horizontal line at this value)
plt.figure(figsize=(20,8))  # Size of figure

plotLC(plt,pn_threshold,in_LCBKGFile) # Plot background region light curve

plt.legend()

The produced light curves can be corrected for a number of effects affecting the detection efficiency, such as vignetting, bad pixels, PSF variations and quantum efficiency, as well as for variations affecting the stability of the detection within the exposure, like dead time and GTIs. Since all these effects can affect in a different manner source and background light curves, the background subtraction has to be done accordingly.
A SAS task, epiclccorr, performs all of these corrections at once. It requires as input both light curves (which are used to establish the binning of the final corrected background subtracted light curve) and the event file. This is done with a simple command line call:

In [ ]:

# Define the output corrected ligthcurve file name

in_LCFile = work_dir+'EPIC_PN_corrected_lightcurve.lc'   # Name of the output corrected lightcurve

In [ ]:

# Correct the light curve with the SAS task epiclccorr

# SAS Command
cmd        = "epiclccorr" # SAS task to be executed                  

# Arguments of SAS Command
inargs     = [f'eventlist={eventfile}',f'srctslist={in_LCSRCFile}',f'outset={in_LCFile}',
              f'bkgtslist={in_LCBKGFile}','withbkgset=yes','applyabsolutecorrections=yes']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task with the parameters to produce the corrected light curve

w(cmd, inargs).run()

Run the following cell to inspect the created corrected light curve.

NOTE: make sure the previous block has finished running before proceeding with the next blocks. The task epiclccorr can take some time to run.

In [ ]:

# Inspect light curve

pn_threshold = 1.0          # cts/sec (only used here for display purposes. Draws a horizontal line at this value)
plt.figure(figsize=(20,8))  # Size of figure

plotLC(plt,pn_threshold,in_LCFile)          # Plot corrected light curve
plotLC(plt,pn_threshold,in_LCSRCFile)       # Plot source region light curve
plotLC(plt,pn_threshold,in_LCBKGFile)       # Plot background region light curve

plt.legend()

Spectra Extraction

First, produce a spectrum for the source region.
Use the source region and include a quality selection appropriate for spectrum extraction:
- For PN, take good events, singles and doubles (#XMMEA_EP&&(PATTERN<=4));
- For MOS, take good events, singles, doubles, triples and quadruples (#XMMEA_EM&&(PATTERN<=12)).

In [ ]:

# Extract the source region spectrum

# Define some parameters for filtering the event file

q_flag       = "#XMMEA_EP" # Quality flag for EPIC pn
n_pattern    = 4           # Pattern selection

# Define the output ligthcurve file name

in_SPSRCFile = work_dir+'EPIC_PN_source_spectrum.fits'   # Name of the output source spectrum

In [ ]:

# SAS Command
cmd        = "evselect" # SAS task to be executed                  

# Arguments of SAS Command
expression = f'{q_flag}&&(PATTERN<={n_pattern})&&((X,Y) IN circle({x_source},{y_source},{r_source}))'  # event filter expression
inargs     = [f'table={eventfile}','withspectrumset=yes',f'spectrumset={in_SPSRCFile}',
              'energycolumn=PI','spectralbinsize=5','withspecranges=yes','specchannelmin=0',
              'specchannelmax=20479',f'expression={expression}']
print("   Filter expression to use: "+expression+" \n")
print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task with the parameters to produce the source spectrum

w(cmd, inargs).run()

Second, produce a spectrum for the background region, using the same filtering criteria.

In [ ]:

# Extract the source region spectrum

# Define some parameters for filtering the event file

q_flag       = "#XMMEA_EP" # Quality flag for EPIC pn
n_pattern    = 4           # Pattern selection

# Define the output ligthcurve file name

in_SPBKGFile = work_dir+'EPIC_PN_background_spectrum.fits'   # Name of the output background spectrum

In [ ]:

# SAS Command
cmd        = "evselect" # SAS task to be executed                  

# Arguments of SAS Command
expression = f'{q_flag}&&(PATTERN<={n_pattern})&&((X,Y) IN circle({x_bkg},{y_bkg},{r_bkg}))'  # event filter expression
inargs     = [f'table={eventfile}','withspectrumset=yes',f'spectrumset={in_SPBKGFile}',
              'energycolumn=PI','spectralbinsize=5','withspecranges=yes','specchannelmin=0',
              'specchannelmax=20479',f'expression={expression}']
print("   Filter expression to use: "+expression+" \n")
print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task with the parameters to produce the background spectrum

w(cmd, inargs).run()

Third, run backscale on both spectra.

In [ ]:

# SAS Command
cmd        = "backscale" # SAS task to be executed                  

# Arguments of SAS Command
inargs     = [f'spectrumset={in_SPSRCFile}',f'badpixlocation={eventfile}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task for the source spectrum

w(cmd, inargs).run()

In [ ]:

# SAS Command
cmd        = "backscale" # SAS task to be executed                  

# Arguments of SAS Command
inargs     = [f'spectrumset={in_SPBKGFile}',f'badpixlocation={eventfile}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task for the background spectrum

w(cmd, inargs).run()

Fourth, generate a redistribution matrix. This task could take some time to run. Check the terminal or log to see when it finishes before executing the next cells.

In [ ]:

# Define some parameters for rmfgen

# Define the output redistribution matrix file name

in_RESPFile = work_dir+'EPIC_PN.rmf'   # Name of the output redistribution

In [ ]:

# SAS Command
cmd        = "rmfgen" # SAS task to be executed                  

# Arguments of SAS Command
inargs     = [f'spectrumset={in_SPSRCFile}',f'rmfset={in_RESPFile}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task for generation of redistribution matrix

w(cmd, inargs).run()

Fifth, generate an ancillary file. This task could take some time to run. Check the terminal or log to see when it finishes before executing the next cells.

In [ ]:

# Define some parameters for rmfgen

# Define the output ancillary file name

in_ARFFile = work_dir+'EPIC_PN.arf'   # Name of the output ancillary

In [ ]:

# SAS Command
cmd        = "arfgen" # SAS task to be executed                  

print("   Checking for Response File ..... \n")
# Check if RESP file is available.
if os.path.isfile(in_RESPFile):
    print ("File "+in_RESPFile+" exists. \n")
else:
    print ("File "+in_RESPFile+" does not exist, please check. \n")

# Arguments of SAS Command
inargs     = [f'spectrumset={in_SPSRCFile}',f'arfset={in_ARFFile}',
              'withrmfset=yes',f'rmfset={in_RESPFile}',f'badpixlocation={eventfile}','detmaptype=psf']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task for generation of the ancillary file 

w(cmd, inargs).run()

Sixth, rebin the spectra and link associated files.

In [ ]:

# Define some parameters for specgruop

# Define the output ancillary file name

in_GRPFile = work_dir+'EPIC_PN_spectrum_grp.fits'   # Name of the output specgruop

In [ ]:

# SAS Command
cmd        = "specgroup" # SAS task to be executed                  

# Arguments of SAS Command
inargs     = [f'spectrumset={in_SPSRCFile}','mincounts=25','oversample=3',
              f'rmfset={in_RESPFile}',f'arfset={in_ARFFile}',
              f'backgndset={in_SPBKGFile}',f'groupedset={in_GRPFile}']

print("   SAS command to be executed: "+cmd+", with arguments; \n")
inargs

In [ ]:

# Execute the SAS task for generation of the ancillary file 

w(cmd, inargs).run()

Last, load the group spectral file in xspec and fit it with a given model.
In the following cells we show some basic commands that can be used as a guideline to plot, fit and interact with xspec. xspec will output the log to the terminal from which the notebook has been launched. The following blocks are all illustrative.

In [ ]:

# Before you start, clear all data and models

AllModels.clear()  # Clear all models
AllData.clear()    # Clear all data

In [ ]:

# Load the group spectral file

# The group file already contains in the header the names of all the necessary files needed by xspec
# (source and background spectra, response and ancillery files)

data1 = Spectrum(in_GRPFile)  # load spectra groupped file
AllData.show()                # inspect loaded data

In [ ]:

# Plot results by opening a window directly on the notebook

Plot.device = "/svg"

In [ ]:

# First inspection of the spectrum

Plot.xAxis = "keV"        # set X axis to energy units
Plot.xLog  = True         # log scale
Plot.yLog  = True         # log scale
Plot("data background")   # plot source and background spectra

In the plot above, you can see the source region spectrum (top panel) and the background region spectra (bottom panel).

In [ ]:

# Define a model and its parameters

model1 = Model("phabs*pow")   # powerlaw with absorption

In [ ]:

Xset.chatter
# Fit model

data1.ignore("**-0.2 10.-**") # For the fit, ignore the range ourside (0.2-10.)
data1.ignore("bad")           # For the fit, ignore bad channels

Fit.nIterations = 100         # maximum number of iterations during fitting
Fit.statMethod  = "cstat"     # fit statistics to be used
Fit.perform()                 # perform fit
cl              = 2.706       # confidence level (90%)
Fit.error(str(cl)+" 1 2 3")   # derive errors at given confidence range for all parameters

fene_ini = 0.2                # Initial energy to derive flux (keV)
fene_end = 10.                # End energy to derive flux (keV)
AllModels.calcFlux(str(fene_ini)+" "+str(fene_end)+" err")  # derive the flux and error in the 0.2-10 keV energy range according to fitted model

In [ ]:

# Extract information from fitted model

nh       = model1(1).values[1]         # Model parameter 1
nherr    = AllModels(1)(1).error       # Error model parameter 1
alpha    = model1(2).values[0]         # Model parameter 2
alphaerr = AllModels(1)(2).error       # Error model parameter 2
norm     = model1(3).values[0]         # Model parameter 3
normerr  = AllModels(1)(3).error       # Error model parameter 3

dof      = Fit.dof                     # get fit dof
statv    = Fit.statMethod              # fit statistics method
stat     = Fit.statistic               # get value of fit statistics
tstatv   = Fit.statTest                # test statistics method
tstat    = Fit.testStatistic           # get value of test statistics
flux     = data1.flux                  # get flux

print(" Fit Statistics: \n")
print("   Value of fit stat: "+str(stat)+" ("+statv+")\n")
print("   Value of test stat: "+str(tstat)+" ("+tstatv+"\n")
print("   Fit dof: "+str(dof)+"\n")
print(" Model Parameters and their error: \n")
print("   Nh Value (10^22) : "+str(100.*nh)+" , "+str(nherr)+"\n")
print("   Alpha Value: "+str(alpha)+" , "+str(alphaerr)+"\n")
print("   Norm Value: "+str(norm)+" , "+str(normerr)+"\n")
print("   Flux Value ("+ str(fene_ini)+"-"+str(fene_end)+" keV) (ergs/cm^2/s): "+str(flux[0])+" , "+"("+str(flux[1])+", "+str(flux[2])+")\n")

In [ ]:

# Plot results by opening a window directly on the notebook

Plot.device = "/svg"

In [ ]:

# Plot data and folded model with residuals

Plot.xAxis = "keV"        # set X axis to energy units
Plot.xLog  = True         # log scale
Plot.yLog  = True         # log scale
Plot("ldata","residuals") # plot data and folded model, with residuals

In the plot above, you can see the source spectrum (top panel) and the residuals against the best fit model (bottom panel).

Caveats

None

xmm-newton

SAS Thread - epic SP LC for a point source in Python - XMM-Newton

How to extract a light curve and spectrum for an EPIC point-like source

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

Last Updated: 25 May 2023

Procedure

Lightcurve Extraction

Spectra Extraction

SCI Directorate

Science Resources

Science Faculty

Visiting us

Support

SAS Thread - epic SP LC for a point source in Python - XMM-Newton

How to extract a light curve and spectrum for an EPIC point-like source

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

Last Updated: 25 May 2023

Procedure

Lightcurve Extraction

Spectra Extraction

SCI Directorate

Science Resources

Science Faculty

Visiting us

Support

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