maps

arches/arches/config.py

@@ -1,109 +1,3 @@
 from pathlib import Path
 
-work_dir = Path('/path/to/some/logical/parent/dir')
+good_pn=['0862470801','0862470601']
-
-"""
-Координаты сентрального кадра, к которому будут
-приводиться изображения всех списков событий
-"""
-ra_cen=34.5342131
-de_cen=-4.7956710
-
-
-"""
-Словарь камер со списком наблюдений каждой камеры.
-Номера камер должны быть отсортированы
-"""
-keylist_tm={'1':['tm1_obs_1',],
-            '5':['tm5_obs_1','tm5_obs_2',],
-            '6':['tm6_obs_1','tm6_obs_2_badpix','tm6_obs_3_badpix',
-                 'tm6_park_1','tm6_park_2','tm6_park_3','tm6_park_4',
-                 'tm6_scan_1','tm6_scan_2','tm6_scan_3','tm6_scan_4'],
-            '7':['tm7_obs_1','tm7_obs_2',]}
-
-"""
-Примерные центры изображений каждого наблюдения.
-Требуется для астрокоррекции. Будет расчитываться матрица сдвигов и поворотов,
-так вот, повороты будут проводиться вокруг данных координат.
-"""
-wcslist={'tm1_obs_1':[34.7279760,-5.0680267],
-         'tm5_obs_1':[34.7351248,-4.4407314],
-         'tm5_obs_2':[34.8748997,-4.4871658],
-         'tm7_obs_1':[35.0015120,-4.7602124],
-         'tm7_obs_2':[34.9810029,-4.5915912],
-         'tm6_obs_1':[34.4227062,-4.7207170],
-         'tm6_obs_2_badpix':[34.7272339,-4.4294153],
-         'tm6_obs_3_badpix':[34.8750291,-4.4708468],
-         'tm6_park_1':[35.0544951,-4.0613619],
-         'tm6_park_2':[35.0558675,-4.0683084],
-         'tm6_park_3':[35.0565263,-4.0583538],
-         'tm6_park_4':[35.0602986,-4.0622220],
-         'tm6_scan_1':[34.5405596,-4.8088748],
-         'tm6_scan_2':[34.5405596,-4.8088748],
-         'tm6_scan_3':[34.5405596,-4.8088748],
-         'tm6_scan_4':[34.5405596,-4.8088748]}
-
-""" like in the paper (Table 1) """
-obslist={'tm1_obs_1':'N01',
-         'tm5_obs_1':'N02',
-         'tm5_obs_2':'N03',
-         'tm7_obs_1':'N15',
-         'tm7_obs_2':'N16',
-         'tm6_obs_1':'N12',
-         'tm6_obs_2_badpix':'N13',
-         'tm6_obs_3_badpix':'N14',
-         'tm6_park_1':'N04',
-         'tm6_park_2':'N06',
-         'tm6_park_3':'N08',
-         'tm6_park_4':'N10',
-         'tm6_scan_1':'N05',
-         'tm6_scan_2':'N07',
-         'tm6_scan_3':'N09',
-         'tm6_scan_4':'N11'}
-
-""" Это просто индекс диапазона для выходных файлов. """
-eband=[   "0", "1",  "2",  "3",   "4", "5",  "6", "7", "8"]
-""" Диапазоны энергий. """
-emin_ev=[ 300, 300,  600, 2300,   200, 300, 5000,  500, 1000]
-emax_ev=[2300, 600, 2300, 5000, 10000,8000, 8000, 1000, 2000]
-
-emin_kev=[0.3, 0.3,  0.6,  2.3,   0.2, 0.3, 5.0, 0.5, 1.0]
-emax_kev=[2.3, 0.6,  2.3,  5.0,  10.0, 8.0, 8.0, 1.0, 2.0]
-
-#ecf = [   1.0, 1.0,  1.0,  1.0,  1.0,  1.0]
-# something is wrong here
-#ecf = [9.7817E+11, 3.2982E+12,  1.3903E+12,  2.3322E+12,  5.2022E+11,  5.8453E+11, 3.8468E+12]
-"""
-*** en0 ecf 9.7817E+11 +/- 2.4606E+10 2.52% N=17
-*** en1 ecf 3.2982E+12 +/- 8.2963E+10 2.52% N=17
-*** en2 ecf 1.3903E+12 +/- 3.5036E+10 2.52% N=17
-*** en3 ecf 2.3322E+12 +/- 5.8717E+10 2.52% N=17
-*** en4 ecf 5.2022E+11 +/- 1.3110E+10 2.52% N=17
-*** en5 ecf 5.8453E+11 +/- 1.4743E+10 2.52% N=17
-"""
-
-# finally use Pavel's file ../data/ECF/ecf_tbabspow_g2nh0.02.pkl
-"""
-for e in [(0.3, 2.3), (0.3, 0.6), (0.6, 2.3), (2.3, 5.0), (5.0, 8.0)]:
-    print(f'{ecf[e]:g}')
-"""
-ecf = [1.0911e+12,  # (0.3, 2.3)
-       1.07252e+12, # (0.3, 0.6)
-       1.08963e+12, # (0.6, 2.3)
-       1.14674e+11, # (2.3, 5.0)
-       1.0,
-       1.0,
-       2.77581e+10, # (5.0, 8.0)
-       1354632916123.6191, # (0.5, 1.0) 4XMM-DR12 EP2 band
-       1014764099304.4968] # (1.0, 2.0) 4XMM-DR12 EP3 band
-
-
-outfile_post='.fits'
-
-"""
-Pavel Medvedev:
-0.3-2.3: 9.135819435325375e-13
-0.3-0.6: 9.160477830652834e-13
-0.6-2.3: 9.201664167869427e-13
-2.3-5.0: 8.721504826794627e-12
-"""

scripts/01_init_events.py

@@ -2,37 +2,74 @@
 
 from pysas.wrapper import Wrapper as w
 import os, sys
+from os.path import dirname
+import inspect
+import glob
 
+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))
 
-sys.exit()
+archive_dir=root_path+'/data/archive'
+events_dir=root_path+'/data/processed'
+products_dir=root_path+'/products'
 
-root='/data/xmm'
+create_folder(events_dir)
 
 inargs = ['--version']
 t = w('sasver', inargs)
 t.run()
 
-obsid='0862470501'
+files = glob.glob(archive_dir+'/0862*')
-local_ccf=f'{root}/work/{obsid}/ccf.cif'
-work_dir=f'{root}/work/{obsid}'
 
-os.environ['SAS_ODF'] = f'{root}/arc/{obsid}/odf/'
+for obsid in files:
+    obsid = os.path.basename(obsid)
     
-if not os.path.exists(work_dir):
+    local_ccf=f'{events_dir}/{obsid}/ccf.cif'
-    os.makedirs(work_dir)
+    work_dir=f'{events_dir}/{obsid}'
 
-os.chdir(work_dir)
+    #inargs = [f'odfid={obsid}',f'workdir={work_dir}']
+    #w('startsas', inargs).run()
 
-if not os.path.exists(local_ccf):
+    os.environ['SAS_ODF'] = f'{archive_dir}/{obsid}/odf'
-    w('cifbuild', []).run()
-else:
-    print("Skip cifbuild, SAS_CCF = {}".format(local_ccf))
 
-os.environ['SAS_CCF'] = local_ccf
+    create_folder(work_dir)
+
+    os.chdir(work_dir)
+
+    if not os.path.exists(local_ccf):
+        w('cifbuild', []).run()
+    else:
+        print("*** Skip cifbuild, SAS_CCF = {}".format(local_ccf))
+    
+    os.environ['SAS_CCF'] = local_ccf
+
+    print("Set SAS_CCF = {}".format(os.environ['SAS_CCF']))
+
+    w('sasver', []).run() # print info
+
+    sasfiles = glob.glob(work_dir+'/*SUM.SAS')
+    if(sasfiles):
+        print("*** Skip odfingest ***")
+    else:
+        w('odfingest', []).run()
+        sasfiles = glob.glob(work_dir+'/*SUM.SAS') # search files again
+
+    os.environ['SAS_ODF'] = sasfiles[0]
+
+    w('sasver', []).run() # print info
+
+    print(f'{work_dir}/????_{obsid}_EPN_S???_ImagingEvts.ds')
+    epfiles = glob.glob(f'{work_dir}/????_{obsid}_EPN_S???_ImagingEvts.ds')
+    if not (epfiles):
+        print('Running epproc')
+        w('epproc', []).run()    
+
+    emfiles = glob.glob(f'{work_dir}/????_{obsid}_EMOS?_S???_ImagingEvts.ds')
+    if not (emfiles):
+        print('Running emproc')
+        w('emproc', []).run()
     
-print(os.environ['SAS_CCF'])

scripts/02_filter_flares.py

@@ -0,0 +1,395 @@
+#!/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'
+
+create_folder(products_dir)
+
+inargs = ['--version']
+t = w('sasver', inargs)
+t.run()
+
+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}'
+
+    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)
+    if not (epfiles):
+        print("*** run 01_init_events.py ***")
+        #w('epproc', []).run()    
+    
+    eventfile = epfiles[0]
+    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")
+
+    #############################################################
+    # For display purposes only define the following event cuts #
+    #############################################################
+    pn_pattern   = 4        # pattern selection
+    pn_pi_min    = 300.     # Low energy range eV
+    pn_pi_max    = 12000.   # High energy range eV
+    pn_flag      = 0        # FLAG
+
+    plt.figure(figsize=(20,8))
+    pl=1
+
+    hdu_list = fits.open(eventfile, memmap=True)
+    evt_data = Table(hdu_list[1].data)
+
+    mask = ((evt_data['PATTERN'] <= pn_pattern) &
+            (evt_data['FLAG'] == pn_flag) &
+            (evt_data['PI'] >= pn_pi_min) &
+            (evt_data['PI'] <= pn_pi_max))
+
+    print("Events in event file" + " " + eventfile + ": " + str(len(evt_data)) + "\n")
+    print("Events in filtered event file" + " " + eventfile + ": " + str(np.sum(mask)) + "\n")
+    print("   Filter: PATTERN <= " + str(pn_pattern) +
+          " : " + str(pn_pi_min) + " <= E(eV) <= " + str(pn_pi_max) +
+          " : " + " FLAG == " + str(pn_flag)+ "\n")
+
+
+    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(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()
+    plt.show()
+
+    ##################################################################
+    # 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.75      # cts/sec (only used here for display purposes)
+
+    out_LCFile   = work_dir+'/EPIC_pn_FlareBKGRate.fit'  # Name of the output BKG lightcurve
+
+    # SAS Command
+    cmd        = "evselect" # SAS task to be executed                  
+    
+    # Arguments of SAS Command
+    expression = f'#XMMEA_EP&&(PI>={pn_pi_min}&&PI<={pn_pi_max})&&(PATTERN=={pn_pattern})'  # event filter expression
+    inargs     = [f'table={eventfile}','withrateset=Y',f'rateset={out_LCFile}','maketimecolumn=Y','timebinsize=100','makeratecolumn=Y',f'expression={expression}']
+
+    print("   Filter expression to use: "+expression+" \n")
+    print("   SAS command to be executed: "+cmd+", with arguments; \n")
+
+    # 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'] <= pn_pattern) &
+            (evt_data['FLAG'] == pn_flag) &
+            (evt_data['PI'] >= pn_pi_min) &
+            (evt_data['PI'] <= pn_pi_max))
+
+    # 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(pn_pattern) +
+          " : " + str(pn_pi_min) + " <= E(eV) <= " + str(pn_pi_max) +
+          " : " + " FLAG == " + str(pn_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,pn_threshold,out_LCFile)
+
+    pl=pl+1
+    plt.show()
+    hdu_list.close()
+
+    ############################################
+    # 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.75      # 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 input BKG lightcurve
+    out_gti_set       = work_dir+'/EPIC_pn_gti.fit'            # Name of the output file containing GTI intervals
+    out_clean_evtFile = work_dir+'/EPIC_pn_gtiFilteredEvts.ds' # Name of the output Event file filtered by GTI
+
+    # SAS Command
+    cmd        = "tabgtigen"
+    
+    # Arguments of SAS Command
+    expression = 'RATE<='+str(pn_threshold)  # event filter expression
+    inargs     = [f'table={in_LCFile}',f'gtiset={out_gti_set}',f'expression={expression}']
+    
+    print("   Filter expression to use: "+expression+" \n")
+    print("   SAS command to be executed: "+cmd+", with arguments; \n")
+    # Execute SAS task with parameters
+    w(cmd, inargs).run()
+
+    # SAS Command
+    cmd        = "evselect"
+    
+    # Arguments of SAS Command
+    expression = ('#XMMEA_EP&&FLAG==0&&(PI>='+str(pn_pi_min)+'&&PI<='+str(pn_pi_max)+
+                  ')&&(gti('+str(out_gti_set)+',TIME))')
+    inargs     = [f'table={eventfile}','withfilteredset=Y',f'filteredset={out_clean_evtFile}',
+                  'destruct=Y','keepfilteroutput=T',f'expression={expression}']
+    
+    print("   Filter expression to use: "+expression+" \n")
+    print("   SAS command to be executed: "+cmd+", with arguments; \n")
+    # Execute SAS task with parameters
+    w(cmd, inargs).run()
+
+
+    plt.figure(figsize=(20,8))
+
+    pl=1
+
+    hdu_list = fits.open(eventfile, memmap=True)
+    evt_data = Table(hdu_list[1].data)
+    prihdu   = hdu_list[1].header
+    print("Events in event file" + " " + eventfile + ": " + str(len(evt_data)) + "\n")
+
+    gti_hdu_list = fits.open(out_clean_evtFile, memmap=True)
+    gti_evt_data = Table(gti_hdu_list[1].data)
+    print("Events in GTI clean event file" + " " + out_clean_evtFile + ": " + str(len(gti_evt_data)) + "\n")
+
+    # Create Events image
+    
+    xmax=np.amax(evt_data['X'])
+    xmin=np.amin(evt_data['X'])
+    xmid=(xmax-xmin)/2.+xmin
+    ymax=np.amax(evt_data['Y'])
+    ymin=np.amin(evt_data['Y'])
+    xbin_size=80
+    ybin_size=80
+    NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))
+
+    plt.subplot(1, 2, pl)
+
+    img_zero_mpl = plt.hist2d(evt_data['X'], evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm())
+
+    cbar = plt.colorbar(ticks=[10.,100.,1000.])
+    cbar.ax.set_yticklabels(['10','100','1000'])
+
+    plt.title(eventfile)
+    plt.xlabel('x')
+    plt.ylabel('y')
+
+    pl=pl+1
+
+    # Create Clean Events image
+
+    xmax=np.amax(gti_evt_data['X'])
+    xmin=np.amin(gti_evt_data['X'])
+    xmid=(xmax-xmin)/2.+xmin
+    ymax=np.amax(gti_evt_data['Y'])
+    ymin=np.amin(gti_evt_data['Y'])
+    xbin_size=80
+    ybin_size=80
+    NBINS = (int((xmax-xmin)/xbin_size),int((ymax-ymin)/ybin_size))
+
+    plt.subplot(1, 2, pl)
+
+    img_zero_mpl = plt.hist2d(gti_evt_data['X'], gti_evt_data['Y'], NBINS, cmap='GnBu', norm=LogNorm())
+
+    cbar = plt.colorbar(ticks=[10.,100.,1000.])
+    cbar.ax.set_yticklabels(['10','100','1000'])
+
+    plt.title(out_clean_evtFile)
+    plt.xlabel('x')
+    plt.ylabel('y')
+
+    plt.text(xmid, ymin-0.1*(ymax-ymin), expression, ha='center')
+
+    pl=pl+1
+
+    gti_hdu_list.close()
+    hdu_list.close()
+    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+'/EPIC_PN_Image.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')
+    
+    continue
+    
+    emfiles = glob.glob(f'{work_dir}/????_{obsid}_EMOS?_S???_ImagingEvts.ds')
+    if not (emfiles):
+        print('Running emproc')
+        w('emproc', []).run()
+