diff --git a/arches/arches/config.py b/arches/arches/config.py index c637416..325655b 100644 --- a/arches/arches/config.py +++ b/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'] diff --git a/arches/arches/utils.py b/arches/arches/utils.py index 48def2e..ef95b6a 100644 --- a/arches/arches/utils.py +++ b/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 diff --git a/scripts/02_filter_flares_mos.py b/scripts/02_filter_flares_mos.py new file mode 100755 index 0000000..388482c --- /dev/null +++ b/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") + + + diff --git a/scripts/02_filter_flares.py b/scripts/02_filter_flares_pn.py similarity index 92% rename from scripts/02_filter_flares.py rename to scripts/02_filter_flares_pn.py index 09b10b3..48efc7c 100755 --- a/scripts/02_filter_flares.py +++ b/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}' - - create_folder(work_dir) + work_dir = init_work_dir(obsid) 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') - - continue - - emfiles = glob.glob(f'{work_dir}/????_{obsid}_EMOS?_S???_ImagingEvts.ds') - if not (emfiles): - print('Running emproc') - w('emproc', []).run() + d.set(f"region {ds9reg_dir}/arches.reg") + + diff --git a/scripts/03_show_image_mos.py b/scripts/03_show_image_mos.py new file mode 100755 index 0000000..ef930c3 --- /dev/null +++ b/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() + diff --git a/scripts/03_show_image_pn.py b/scripts/03_show_image_pn.py new file mode 100755 index 0000000..15375fa --- /dev/null +++ b/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() + diff --git a/scripts/04_spectrum_mos.py b/scripts/04_spectrum_mos.py new file mode 100755 index 0000000..30a8a6d --- /dev/null +++ b/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() diff --git a/scripts/04_spectrum_pn.py b/scripts/04_spectrum_pn.py new file mode 100755 index 0000000..d113c09 --- /dev/null +++ b/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()