spectra

arches/arches/config.py

@@ -1,3 +1,17 @@
 from pathlib import Path
 
 good_pn=['0862470801','0862470601']
+
+pn_threshold = 0.2
+mos_threshold=[0.2,0.3]
+
+"""
+Ядро скопления
+RA=17h 45m 50.3s
+Dec=-28d 49m 19s
+Radius=15"
+
+Фоновая скорость счёта для области ядра Арки оценивалась в кольце 70′′ < R < 130′′.
+"""
+
+skip=['0862471401','0862471501','0862470501']

arches/arches/utils.py

@@ -13,18 +13,24 @@ from fitsio import FITS
 from pathlib import Path
 import pandas
 import pickle
+
+from os.path import dirname
+import inspect
+from pysas.wrapper import Wrapper as w
+
+
 from astropy.table import QTable, Table, Column
 from astropy import units as u
-
 from astropy.coordinates import SkyCoord # High-level coordinates
 from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames
 from astropy.coordinates import Angle, Latitude, Longitude # Angles
-
 from astropy.time import Time, TimeDelta, TimezoneInfo, TimeFromEpoch
 
 import statistics
 import shutil
 
+
+
 def create_folder(folder):
     if not (os.path.exists(folder)):
         os.makedirs(folder)
@@ -81,3 +87,46 @@ def plotLC(plt,threshold,fileName):
     else:
         print("File not found "+fileName+"\n")
         
+def init_work_dir(obsid):
+
+    import arches
+
+    root_path=dirname(dirname(dirname(inspect.getfile(arches))))
+    print("Arches root path: {}".format(root_path))
+
+    archive_dir=root_path+'/data/archive'
+    events_dir=root_path+'/data/processed'
+    products_dir=root_path+'/products'
+    ds9reg_dir=root_path+'/data/ds9reg'
+
+    create_folder(products_dir)
+
+    #inargs = ['--version']
+    #t = w('sasver', inargs)
+    #t.run()
+
+    print(f'*** ObsID: {obsid} ***')
+    
+    local_ccf=f'{events_dir}/{obsid}/ccf.cif'
+
+    work_dir=f'{products_dir}/{obsid}'
+
+    create_folder(work_dir)
+
+    if not os.path.exists(local_ccf):
+        print("*** Not found SAS_CCF = {}".format(local_ccf))
+        sys.exit()
+
+    sasfiles = glob.glob(events_dir+f'/{obsid}/*SUM.SAS')
+    if not (sasfiles):
+        print("*** run 01_init_events.py ***")
+        sys.exit()
+
+    os.environ['SAS_ODF'] = sasfiles[0]
+    os.environ['SAS_CCF'] = local_ccf
+
+    w('sasver', []).run() # print info
+
+    inargs = [f'workdir={work_dir}']
+    w('startsas', inargs).run()
+    return work_dir

scripts/02_filter_flares_mos.py

@@ -0,0 +1,370 @@
+#!/usr/bin/env python
+
+from pysas.wrapper import Wrapper as w
+import os, sys
+from os.path import dirname
+import inspect
+import glob
+
+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
+
+import pyds9
+
+import arches
+from arches.utils import *
+from arches.config import *
+
+root_path=dirname(dirname(dirname(inspect.getfile(arches))))
+print("Arches root path: {}".format(root_path))
+
+archive_dir=root_path+'/data/archive'
+events_dir=root_path+'/data/processed'
+products_dir=root_path+'/products'
+ds9reg_dir=root_path+'/data/ds9reg'
+
+imos=2
+
+
+create_folder(products_dir)
+
+inargs = ['--version']
+t = w('sasver', inargs)
+t.run()
+
+files = glob.glob(archive_dir+'/*')
+
+for obsid in files:
+    obsid = os.path.basename(obsid)
+
+    work_dir = init_work_dir(obsid)
+
+    os.chdir(work_dir)
+
+    search_str = f'{events_dir}/{obsid}/????_{obsid}_EMOS{imos}_S???_ImagingEvts.ds'
+    print(search_str)
+    emfiles = glob.glob(search_str)
+    if not (emfiles):
+        print("*** run 01_init_events.py ***")
+        sys.exit()
+
+    eventfile = emfiles[0]
+    print(f"Checking for EPIC-MOS{imos} 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")
+    
+    #############################################################
+    # For display purposes only define the following event cuts #
+    #############################################################
+    mos_pattern   = 4        # pattern selection
+    mos_pi_min    = 300.     # Low energy range eV
+    mos_pi_max    = 12000.   # High energy range eV
+    mos_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'] <= mos_pattern) &
+            (evt_data['FLAG'] == mos_flag) &
+            (evt_data['PI'] >= mos_pi_min) &
+            (evt_data['PI'] <= mos_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(mos_pattern) +
+          " : " + str(mos_pi_min) + " <= E(eV) <= " + str(mos_pi_max) +
+          " : " + " FLAG == " + str(mos_flag)+ "\n")
+
+
+    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(obsid)
+    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(obsid)
+    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()
+    plt.show()
+    
+    ##################################################################
+    # Define a SAS filter expression to derive a background rate cut #
+    ##################################################################
+    mos_pattern   = 0        # pattern selection
+    mos_pi_min    = 10000.   # Low energy range eV
+    #mos_threshold = threshold[imos-1]      # cts/sec (only used here for display purposes)
+
+    out_LCFile   = work_dir+f'/EPIC_MOS{imos}_FlareBKGRate.fit'  # Name of the output BKG lightcurve
+
+    # SAS Command
+    cmd        = "evselect" # SAS task to be executed                  
+    
+    # Arguments of SAS Command
+    expression = f'#XMMEA_EM&&(PI>={mos_pi_min})&&(PATTERN=={mos_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")
+
+    # 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
+
+    mask = ((evt_data['PATTERN'] <= mos_pattern) &
+            (evt_data['FLAG'] == mos_flag) &
+            (evt_data['PI'] >= mos_pi_min))
+
+    # 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")
+    print("Events in filtered event file" + " " + eventfile + ": " + str(np.sum(mask)) + "\n")
+    print("   Filter: PATTERN <= " + str(mos_pattern) +
+          " : " + str(mos_pi_min) + " <= E(eV) : FLAG == " + str(mos_flag)+ "\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(obsid)
+    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,mos_threshold[imos-1],out_LCFile)
+
+    pl=pl+1
+    plt.show()
+    hdu_list.close()
+    
+    ############################################
+    # Define energy range to filter event file #
+    ############################################
+    mos_pattern   = 4        # pattern selection
+    mos_pi_min    = 200.     # Low energy range eV
+    mos_pi_max    = 10000.   # High energy range eV
+    #mos_threshold = 0.75      # Cut to be applied to filter event file (cts/sec)
+
+    # Define the input and output file names
+
+    in_LCFile         = work_dir+f'/EPIC_MOS{imos}_FlareBKGRate.fit'   # Name of the input BKG lightcurve
+    out_gti_set       = work_dir+f'/EPIC_MOS{imos}_gti.fit'            # Name of the output file containing GTI intervals
+    out_clean_evtFile = work_dir+f'/EPIC_MOS{imos}_gtiFilteredEvts.ds' # Name of the output Event file filtered by GTI
+
+    # SAS Command
+    cmd        = "tabgtigen"
+    
+    # Arguments of SAS Command
+    expression = 'RATE<='+str(mos_threshold[imos-1])  # 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")
+    # Execute SAS task with parameters
+    w(cmd, inargs).run()
+
+    # SAS Command
+    cmd        = "evselect"
+    
+    # Arguments of SAS Command
+    expression = ('#XMMEA_EM&&FLAG==0&&(PI>='+str(mos_pi_min)+'&&PI<='+str(mos_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")
+    # 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(obsid)
+    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()
+    plt.show()
+
+    ###############################################################################
+    # 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+f'/EPIC_MOS{imos}_Image_{obsid}.fit'  # Name of the output Image file 
+    # SAS Command
+    cmd        = "evselect" # SAS task to be executed                  
+
+    # Arguments of SAS Command
+    
+    inargs     = [f'table={out_clean_evtFile}','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")
+    # Execute the SAS task with the parameters to produce an image
+    w(cmd, inargs).run()
+    
+    # Visualize the image with ds9
+
+    d = pyds9.DS9()
+    d.set("file "+out_IMFile)
+    d.set('cmap bb')
+    d.set('scale log')
+    d.set(f"region {ds9reg_dir}/arches.reg")
+
+    
+    

scripts/02_filter_flares_pn.py

@@ -27,6 +27,7 @@ print("Arches root path: {}".format(root_path))
 archive_dir=root_path+'/data/archive'
 events_dir=root_path+'/data/processed'
 products_dir=root_path+'/products'
+ds9reg_dir=root_path+'/data/ds9reg'
 
 create_folder(products_dir)
 
@@ -38,32 +39,11 @@ files = glob.glob(archive_dir+'/0862*')
 
 for obsid in files:
     obsid = os.path.basename(obsid)
-    print(f'*** ObsID: {obsid} ***')
-    local_ccf=f'{events_dir}/{obsid}/ccf.cif'
 
-    work_dir=f'{products_dir}/{obsid}'
+    work_dir = init_work_dir(obsid)
-
-    create_folder(work_dir)
 
     os.chdir(work_dir)
 
-    if not os.path.exists(local_ccf):
-        print("*** Not found SAS_CCF = {}".format(local_ccf))
-        sys.exit()
-
-    sasfiles = glob.glob(events_dir+f'/{obsid}/*SUM.SAS')
-    if not (sasfiles):
-        print("*** run 01_init_events.py ***")
-        sys.exit()
-
-    os.environ['SAS_ODF'] = sasfiles[0]
-    os.environ['SAS_CCF'] = local_ccf
-
-    w('sasver', []).run() # print info
-
-    inargs = [f'workdir={work_dir}']
-    w('startsas', inargs).run()
-
     search_str = f'{events_dir}/{obsid}/????_{obsid}_EPN_S???_ImagingEvts.ds'
     print(search_str)
     epfiles = glob.glob(search_str)
@@ -323,7 +303,7 @@ for obsid in files:
     cbar = plt.colorbar(ticks=[10.,100.,1000.])
     cbar.ax.set_yticklabels(['10','100','1000'])
 
-    plt.title(eventfile)
+    plt.title(obsid)
     plt.xlabel('x')
     plt.ylabel('y')
 
@@ -367,7 +347,7 @@ for obsid in files:
     xcoord='X'  # coordinate system
     ycoord='Y'  # coordinate system
 
-    out_IMFile   = work_dir+'/EPIC_PN_Image.fit'  # Name of the output Image file 
+    out_IMFile   = work_dir+f'/EPIC_PN_Image_{obsid}.fit'  # Name of the output Image file 
     # SAS Command
     cmd        = "evselect" # SAS task to be executed                  
 
@@ -385,11 +365,7 @@ for obsid in files:
     d.set("file "+out_IMFile)
     d.set('cmap bb')
     d.set('scale log')
-    
+    d.set(f"region {ds9reg_dir}/arches.reg")
-    continue
+
-    
+
-    emfiles = glob.glob(f'{work_dir}/????_{obsid}_EMOS?_S???_ImagingEvts.ds')
-    if not (emfiles):
-        print('Running emproc')
-        w('emproc', []).run()
     

scripts/03_show_image_mos.py

@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+
+from pysas.wrapper import Wrapper as w
+import os, sys
+from os.path import dirname
+import inspect
+import glob
+
+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
+
+import pyds9
+
+import arches
+from arches.utils import *
+from arches.config import *
+
+root_path=dirname(dirname(dirname(inspect.getfile(arches))))
+print("Arches root path: {}".format(root_path))
+
+archive_dir=root_path+'/data/archive'
+events_dir=root_path+'/data/processed'
+products_dir=root_path+'/products'
+ds9reg_dir=root_path+'/data/ds9reg'
+
+imos=1
+
+create_folder(products_dir)
+
+inargs = ['--version']
+t = w('sasver', inargs)
+t.run()
+
+files = glob.glob(archive_dir+'/*')
+
+for obsid in files:
+    obsid = os.path.basename(obsid)
+
+    work_dir = init_work_dir(obsid)
+
+    os.chdir(work_dir)
+    # EPIC_pn_gtiFilteredEvts.ds EPIC_MOS2_Image_0862470501.fit
+    search_str = f'{products_dir}/{obsid}/EPIC_MOS{imos}_Image_{obsid}.fit'
+    print(search_str)
+    epfiles = glob.glob(search_str)
+    if not (epfiles):
+        print("*** run 02_filter_flares_pn.py ***")
+        sys.exit()
+        
+    out_IMFile   = epfiles[0]
+
+    # Visualize the image with ds9
+
+    d = pyds9.DS9()
+    d.set("file "+out_IMFile)
+    d.set('cmap bb')
+    d.set('scale log')
+    d.set(f"region {ds9reg_dir}/arches.reg")
+
+    print(f'{obsid} MOS{imos}')
+    user_input = input()
+    

scripts/03_show_image_pn.py

@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+
+from pysas.wrapper import Wrapper as w
+import os, sys
+from os.path import dirname
+import inspect
+import glob
+
+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
+
+import pyds9
+
+import arches
+from arches.utils import *
+from arches.config import *
+
+root_path=dirname(dirname(dirname(inspect.getfile(arches))))
+print("Arches root path: {}".format(root_path))
+
+archive_dir=root_path+'/data/archive'
+events_dir=root_path+'/data/processed'
+products_dir=root_path+'/products'
+ds9reg_dir=root_path+'/data/ds9reg'
+
+create_folder(products_dir)
+
+inargs = ['--version']
+t = w('sasver', inargs)
+t.run()
+
+files = glob.glob(archive_dir+'/*')
+
+for obsid in files:
+    obsid = os.path.basename(obsid)
+
+    work_dir = init_work_dir(obsid)
+
+    os.chdir(work_dir)
+    # EPIC_pn_gtiFilteredEvts.ds
+    search_str = f'{products_dir}/{obsid}/EPIC_PN_Image_{obsid}.fit'
+    print(search_str)
+    epfiles = glob.glob(search_str)
+    if not (epfiles):
+        print("*** run 02_filter_flares_pn.py ***")
+        sys.exit()
+        
+    out_IMFile   = epfiles[0]
+
+    # Visualize the image with ds9
+
+    d = pyds9.DS9()
+    d.set("file "+out_IMFile)
+    d.set('cmap bb')
+    d.set('scale log')
+    d.set(f"region {ds9reg_dir}/arches.reg")
+
+    print(obsid)
+    user_input = input()
+    

scripts/04_spectrum_mos.py

@@ -0,0 +1,242 @@
+#!/usr/bin/env python
+
+from pysas.wrapper import Wrapper as w
+import os, sys
+from os.path import dirname
+import inspect
+import glob
+
+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
+
+import re
+import pyds9
+
+import arches
+from arches.utils import *
+from arches.config import *
+
+root_path=dirname(dirname(dirname(inspect.getfile(arches))))
+print("Arches root path: {}".format(root_path))
+
+archive_dir=root_path+'/data/archive'
+events_dir=root_path+'/data/processed'
+products_dir=root_path+'/products'
+ds9reg_dir=root_path+'/data/ds9reg'
+
+imos=2
+
+create_folder(products_dir)
+
+inargs = ['--version']
+t = w('sasver', inargs)
+t.run()
+
+files = glob.glob(archive_dir+'/*')
+
+for obsid in files:
+    obsid = os.path.basename(obsid)
+
+    work_dir = init_work_dir(obsid)
+
+    os.chdir(work_dir)
+
+    eventfile=f'{products_dir}/{obsid}/EPIC_MOS{imos}_gtiFilteredEvts.ds'
+    if(os.path.isfile(eventfile)==False):
+        print("*** run 02_filter_flares_mos.py ***")
+        sys.exit()
+
+    
+    # EPIC_pn_gtiFilteredEvts.ds
+    search_str = f'{products_dir}/{obsid}/EPIC_MOS{imos}_Image_{obsid}.fit'
+    print(search_str)
+    emfiles = glob.glob(search_str)
+    if not (emfiles):
+        print("*** run 02_filter_flares_mos.py ***")
+        sys.exit()
+        
+    out_IMFile   = emfiles[0]
+
+    # Visualize the image with ds9
+
+    d = pyds9.DS9()
+    d.set("file "+out_IMFile)
+    d.set('cmap bb')
+    d.set('scale log')
+    #d.set(f"region {ds9reg_dir}/arches.reg")
+    src_fn=f'{ds9reg_dir}/{obsid}/src.reg'
+    if(os.path.isfile(src_fn)==True):
+        d.set(f"region {src_fn}")
+    bkg_fn=f'{ds9reg_dir}/{obsid}/bkg.reg'
+    if(os.path.isfile(bkg_fn)==True):
+        d.set(f"region {bkg_fn}")
+        
+    reply = input(f"{obsid}: Proceed to make spectrum/lightcurve? [y/[n]] ") 
+    
+    if reply!='y':
+        continue
+        
+    print(d.get("regions"))
+    
+    # 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 first region: ", region1)
+    print("Identified second region: ", region2)
+    
+    #print("region1 "+region1.partition("color")[2].replace('=','').replace('\n',''))
+    #print("region2 "+region2.partition("color")[2].replace('=','').replace('\n',''))
+    
+    # Identify source and background regions using the 'white' color.
+    c1=region1.partition("color")[2].replace('=','').replace('\n','')
+    print("Region1 color: "+c1)
+    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)")
+
+    #########################################
+    # Extract the source region light curve #
+    #########################################
+    
+    # Define some parameters for filtering the event file and define the light curve binning
+
+    q_flag       = "#XMMEA_EM" # Quality flag for EPIC pn
+    n_pattern    = 12           # Pattern selection
+    mos_pi_min    = 200.        # Low energy range eV
+    mos_pi_max    = 10000.      # High energy range eV
+    lc_bin       = 100         # Light curve bin in secs
+
+    # Define the output ligthcurve file name
+    in_LCSRCFile = work_dir+f'/EPIC_MOS{imos}_source_lightcurve.lc'   # Name of the output source light curve
+
+    # 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 [{mos_pi_min}:{mos_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")
+    # Execute the SAS task with the parameters to produce the source region light curve
+    w(cmd, inargs).run()
+
+    # Inspect light curve
+
+    plt.figure(figsize=(20,8))  # Size of figure
+
+    plotLC(plt,mos_threshold[imos-1],in_LCSRCFile) # Plot source region light curve
+
+    plt.legend()
+    plt.show()
+    
+    ######################################
+    # Extract the source region spectrum #
+    ######################################
+
+    # Define some parameters for filtering the event file
+    q_flag       = "#XMMEA_EM" # Quality flag for EPIC pn
+    n_pattern    = 12           # Pattern selection
+
+    # Define the output ligthcurve file name
+    in_SPSRCFile = work_dir+f'/EPIC_MOS{imos}_source_spectrum.fits'   # Name of the output source spectrum
+
+    # 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=11999',f'expression={expression}']
+
+    print("   Filter expression to use: "+expression+" \n")
+    print("   SAS command to be executed: "+cmd+", with arguments; \n")
+
+    w(cmd, inargs).run()
+
+    # Extract the background region spectrum with the same criteria
+
+    # Define the output ligthcurve file name    
+    in_SPBKGFile = work_dir+f'/EPIC_MOS{imos}_background_spectrum.fits'   # Name of the output background spectrum
+
+    # 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=11999',f'expression={expression}']
+
+    print("   Filter expression to use: "+expression+" \n")
+    print("   SAS command to be executed: "+cmd+", with arguments; \n")
+
+    w(cmd, inargs).run()
+
+    w("backscale", [f'spectrumset={in_SPSRCFile}',f'badpixlocation={eventfile}']).run()
+    w("backscale", [f'spectrumset={in_SPBKGFile}',f'badpixlocation={eventfile}']).run()
+
+    in_RESPFile = work_dir+f'/EPIC_MOS{imos}.rmf'   # Name of the output redistribution
+    w("rmfgen", [f'spectrumset={in_SPSRCFile}',f'rmfset={in_RESPFile}']).run()
+
+    in_ARFFile = work_dir+f'/EPIC_MOS{imos}.arf'   # Name of the output ancillary
+    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")
+        sys.exit()
+
+    # Arguments of SAS Command
+    inargs     = [f'spectrumset={in_SPSRCFile}',f'arfset={in_ARFFile}',
+                  'withrmfset=yes',f'rmfset={in_RESPFile}',f'badpixlocation={eventfile}','detmaptype=psf']
+    w(cmd, inargs).run()
+
+    # rebin the spectra and link associated files
+    in_GRPFile = work_dir+f'/EPIC_MOS{imos}_spectrum_grp.fits'   # Name of the output specgruop
+    inargs     = [f'spectrumset={in_SPSRCFile}','mincounts=30','oversample=3',
+                  f'rmfset={in_RESPFile}',f'arfset={in_ARFFile}',
+                  f'backgndset={in_SPBKGFile}',f'groupedset={in_GRPFile}']
+    w("specgroup", inargs).run()

scripts/04_spectrum_pn.py

@@ -0,0 +1,240 @@
+#!/usr/bin/env python
+
+from pysas.wrapper import Wrapper as w
+import os, sys
+from os.path import dirname
+import inspect
+import glob
+
+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
+
+import re
+import pyds9
+
+import arches
+from arches.utils import *
+from arches.config import *
+
+root_path=dirname(dirname(dirname(inspect.getfile(arches))))
+print("Arches root path: {}".format(root_path))
+
+archive_dir=root_path+'/data/archive'
+events_dir=root_path+'/data/processed'
+products_dir=root_path+'/products'
+ds9reg_dir=root_path+'/data/ds9reg'
+
+create_folder(products_dir)
+
+inargs = ['--version']
+t = w('sasver', inargs)
+t.run()
+
+files = glob.glob(archive_dir+'/*')
+
+for obsid in files:
+    obsid = os.path.basename(obsid)
+
+    work_dir = init_work_dir(obsid)
+
+    os.chdir(work_dir)
+
+    eventfile=f'{products_dir}/{obsid}/EPIC_pn_gtiFilteredEvts.ds'
+    if(os.path.isfile(eventfile)==False):
+        print("*** run 02_filter_flares_pn.py ***")
+        sys.exit()
+
+    
+    # EPIC_pn_gtiFilteredEvts.ds
+    search_str = f'{products_dir}/{obsid}/EPIC_PN_Image_{obsid}.fit'
+    print(search_str)
+    epfiles = glob.glob(search_str)
+    if not (epfiles):
+        print("*** run 02_filter_flares_pn.py ***")
+        sys.exit()
+        
+    out_IMFile   = epfiles[0]
+
+    # Visualize the image with ds9
+
+    d = pyds9.DS9()
+    d.set("file "+out_IMFile)
+    d.set('cmap bb')
+    d.set('scale log')
+    #d.set(f"region {ds9reg_dir}/arches.reg")
+    src_fn=f'{ds9reg_dir}/{obsid}/src.reg'
+    if(os.path.isfile(src_fn)==True):
+        d.set(f"region {src_fn}")
+    bkg_fn=f'{ds9reg_dir}/{obsid}/bkg.reg'
+    if(os.path.isfile(bkg_fn)==True):
+        d.set(f"region {bkg_fn}")
+        
+    reply = input(f"{obsid}: Proceed to make spectrum/lightcurve? [y/[n]] ") 
+    
+    if reply!='y':
+        continue
+        
+    print(d.get("regions"))
+    
+    # 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 first region: ", region1)
+    print("Identified second region: ", region2)
+    
+    #print("region1 "+region1.partition("color")[2].replace('=','').replace('\n',''))
+    #print("region2 "+region2.partition("color")[2].replace('=','').replace('\n',''))
+    
+    # Identify source and background regions using the 'white' color.
+    c1=region1.partition("color")[2].replace('=','').replace('\n','')
+    print("Region1 color: "+c1)
+    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)")
+
+    #########################################
+    # Extract the source region light curve #
+    #########################################
+    
+    # Define some parameters for filtering the event file and define the light curve 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         # Light curve bin in secs
+
+    # Define the output ligthcurve file name
+    in_LCSRCFile = work_dir+'/EPIC_PN_source_lightcurve.lc'   # Name of the output source light curve
+
+    # 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")
+    # Execute the SAS task with the parameters to produce the source region light curve
+    w(cmd, inargs).run()
+
+    # Inspect light curve
+
+    plt.figure(figsize=(20,8))  # Size of figure
+
+    plotLC(plt,pn_threshold,in_LCSRCFile) # Plot source region light curve
+
+    plt.legend()
+    plt.show()
+    
+    ######################################
+    # 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
+
+    # 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")
+
+    w(cmd, inargs).run()
+
+    # Extract the background region spectrum with the same criteria
+
+    # Define the output ligthcurve file name    
+    in_SPBKGFile = work_dir+'/EPIC_PN_background_spectrum.fits'   # Name of the output background spectrum
+
+    # 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")
+
+    w(cmd, inargs).run()
+
+    w("backscale", [f'spectrumset={in_SPSRCFile}',f'badpixlocation={eventfile}']).run()
+    w("backscale", [f'spectrumset={in_SPBKGFile}',f'badpixlocation={eventfile}']).run()
+
+    in_RESPFile = work_dir+'/EPIC_PN.rmf'   # Name of the output redistribution
+    w("rmfgen", [f'spectrumset={in_SPSRCFile}',f'rmfset={in_RESPFile}']).run()
+
+    in_ARFFile = work_dir+'/EPIC_PN.arf'   # Name of the output ancillary
+    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")
+        sys.exit()
+
+    # Arguments of SAS Command
+    inargs     = [f'spectrumset={in_SPSRCFile}',f'arfset={in_ARFFile}',
+                  'withrmfset=yes',f'rmfset={in_RESPFile}',f'badpixlocation={eventfile}','detmaptype=psf']
+    w(cmd, inargs).run()
+
+    # rebin the spectra and link associated files
+    in_GRPFile = work_dir+'/EPIC_PN_spectrum_grp.fits'   # Name of the output specgruop
+    inargs     = [f'spectrumset={in_SPSRCFile}','mincounts=30','oversample=3',
+                  f'rmfset={in_RESPFile}',f'arfset={in_ARFFile}',
+                  f'backgndset={in_SPBKGFile}',f'groupedset={in_GRPFile}']
+    w("specgroup", inargs).run()