forked from xmm/arches
325 lines
9.6 KiB
Python
325 lines
9.6 KiB
Python
import os
|
|
import re
|
|
import sys
|
|
import csv
|
|
import numpy as np
|
|
from astropy.io import fits
|
|
from astropy import wcs
|
|
from astropy.wcs import WCS
|
|
from astropy.io.fits import update
|
|
from astropy.io.fits import getdata
|
|
import glob
|
|
import json
|
|
from fitsio import FITS
|
|
from pathlib import Path
|
|
import pandas
|
|
import pickle
|
|
import pyds9
|
|
|
|
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
|
|
|
|
from arches.config import *
|
|
|
|
def test_exe(program):
|
|
""" Tests if executable exists in PATH """
|
|
import os
|
|
def is_exe(fpath):
|
|
return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
|
|
|
|
fpath, fname = os.path.split(program)
|
|
if fpath:
|
|
if is_exe(program):
|
|
return program
|
|
else:
|
|
for path in os.environ["PATH"].split(os.pathsep):
|
|
exe_file = os.path.join(path, program)
|
|
if is_exe(exe_file):
|
|
return exe_file
|
|
|
|
print("\n*** Command {} not found ***\n".format(program))
|
|
sys.exit()
|
|
|
|
def create_folder(folder):
|
|
if not (os.path.exists(folder)):
|
|
os.makedirs(folder)
|
|
|
|
def remove_file(filename):
|
|
if(os.path.isfile(filename)==True):
|
|
os.remove(filename)
|
|
|
|
# Function to plot Lightcurve
|
|
def plotLC(plt,threshold,fileName):
|
|
if fileName != "NOT FOUND":
|
|
fitsFile = fits.open(fileName)
|
|
prihdu = fitsFile[1].header
|
|
if ('CUTVAL' in prihdu):
|
|
threshold = prihdu['CUTVAL']
|
|
|
|
cols = fitsFile[1].columns
|
|
colName = None
|
|
for i,x in enumerate(cols.names):
|
|
if "RATE" in x:
|
|
colName = cols.names[i]
|
|
if "COUNTS" in x:
|
|
colName = cols.names[i]
|
|
|
|
data = fitsFile[1].data
|
|
|
|
xdata = data.field('TIME') - min(data.field('TIME')) # extract the x data column
|
|
ydata = data.field(colName)
|
|
|
|
xmax=np.amax(xdata)
|
|
xmin=np.amin(xdata)
|
|
|
|
plt.plot(xdata,ydata) # plot the data
|
|
|
|
if colName == 'RATE':
|
|
plt.title("Flaring particle background lightcurve")
|
|
plt.xlabel("Time (s)")
|
|
plt.ylabel("Cts/s")
|
|
else:
|
|
plt.title("Lightcurve")
|
|
plt.xlabel("Time (s)")
|
|
plt.ylabel("Counts")
|
|
|
|
if (threshold != 'None'):
|
|
if (colName == 'COUNTS'):
|
|
threshold=float(threshold)*100.
|
|
|
|
y2data = [threshold]*len(xdata)
|
|
plt.plot(xdata,y2data)
|
|
plt.text(xmin+0.1*(xmax-xmin), threshold+0.01*threshold, str(threshold)+" cts/sec", ha='center')
|
|
|
|
fitsFile.close()
|
|
|
|
else:
|
|
print("File not found "+fileName+"\n")
|
|
|
|
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
|
|
|
|
""" See also convregion task, for conversion of regions """
|
|
def get_det_coords(ra,dec,image,key):
|
|
logfile=f'{key}-ecoordconv.log'
|
|
remove_file(logfile)
|
|
""" Convert sky coordinates to correct DETX/Y """
|
|
inargs = [f'imageset={image}',
|
|
f'x={ra}',
|
|
f'y={dec}',
|
|
'coordtype=eqpos',]
|
|
w('ecoordconv', inargs, logFile=logfile).run()
|
|
|
|
detx=None
|
|
dety=None
|
|
with open(logfile, "r") as file:
|
|
lines = file.readlines()
|
|
for line in lines:
|
|
if re.search('DETX', line):
|
|
print(line, end='\n')
|
|
(tx,ty,dx,dy)=line.strip().split()
|
|
detx=float(dx)
|
|
dety=float(dy)
|
|
return detx,dety
|
|
|
|
def run_mosspectra(key,ccds,check=True):
|
|
|
|
eventfile = f'{key}-allevc.fits'
|
|
|
|
(detx,dety)=get_det_coords(arches_ra,arches_dec,f'{key}-fovimt.fits',key)
|
|
|
|
arches_Rmax_pix=int(np.round(arches_Rmax_as/det_pix_as))
|
|
region_fn=f'{key}.txt'
|
|
remove_file(region_fn)
|
|
freg = open(region_fn, "w")
|
|
expression=f"((DETX,DETY) IN circle({detx},{dety},{arches_Rmax_pix}))"
|
|
freg.write(f'&&{expression}\n')
|
|
freg.close()
|
|
|
|
# check selected region
|
|
if(check):
|
|
unfilt_image = eventfile.replace("allevc","region")
|
|
inargs = [f'table={eventfile}',
|
|
'withimageset=yes', f'imageset={unfilt_image}',
|
|
'filtertype=expression', f'expression={expression}',
|
|
'ignorelegallimits=yes', 'imagebinning=imageSize',
|
|
'xcolumn=DETX','ycolumn=DETY','ximagesize=780','yimagesize=780',
|
|
'ximagemax=19500','yimagemax=19500', 'ximagemin=-19499', 'yimagemin=-19499',]
|
|
w('evselect', inargs).run()
|
|
|
|
#d = pyds9.DS9()
|
|
#d.set("file "+unfilt_image)
|
|
#d.set('cmap bb')
|
|
#d.set('scale log')
|
|
#return
|
|
|
|
|
|
inargs = [f'eventfile={eventfile}',
|
|
'withsrcrem=yes',
|
|
f'cornerfile={key}-corevc.fits',
|
|
f'imagefile={key}-fovimt.fits',
|
|
f'expmap={key}-fovimtexp.fits',
|
|
f'maskdet={key}-bkgregtdet.fits',
|
|
f'masksky={key}-bkgregtsky.fits',
|
|
'withregion=yes', f'regionfile={region_fn}',
|
|
'keepinterfiles=yes',
|
|
'pattern=12',
|
|
'elow=350',
|
|
'ehigh=1100',
|
|
f'ccds="{ccds}"',]
|
|
|
|
w('mosspectra', inargs).run()
|
|
|
|
def run_pnspectra(key,quads,check=True):
|
|
eventfile = f'{key}-allevc.fits'
|
|
|
|
(detx,dety)=get_det_coords(arches_ra,arches_dec,f'{key}-fovimt.fits',key)
|
|
|
|
arches_Rmax_pix=int(np.round(arches_Rmax_as/det_pix_as))
|
|
region_fn=f'{key}.txt'
|
|
remove_file(region_fn)
|
|
freg = open(region_fn, "w")
|
|
expression=f"((DETX,DETY) IN circle({detx},{dety},{arches_Rmax_pix}))"
|
|
freg.write(f'&&{expression}\n')
|
|
freg.close()
|
|
|
|
# check selected region
|
|
if(check):
|
|
unfilt_image = eventfile.replace("allevc","region")
|
|
inargs = [f'table={eventfile}',
|
|
'withimageset=yes', f'imageset={unfilt_image}',
|
|
'filtertype=expression', f'expression={expression}',
|
|
'ignorelegallimits=yes', 'imagebinning=imageSize',
|
|
'xcolumn=DETX','ycolumn=DETY','ximagesize=780','yimagesize=780',
|
|
'ximagemax=19500','yimagemax=19500', 'ximagemin=-19499', 'yimagemin=-19499',]
|
|
w('evselect', inargs).run()
|
|
|
|
#d = pyds9.DS9()
|
|
#d.set("file "+unfilt_image)
|
|
#d.set('cmap bb')
|
|
#d.set('scale log')
|
|
#return
|
|
|
|
inargs = [f'eventfile={key}-allevc.fits',
|
|
'withsrcrem=yes',
|
|
f'ootevtfile={key}-allevcoot.fits',
|
|
f'cornerfile={key}-corevc.fits',
|
|
f'imagefile={key}-fovimt.fits',
|
|
f'expmap={key}-fovimtexp.fits',
|
|
f'maskdet={key}-bkgregtdet.fits',
|
|
f'masksky={key}-bkgregtsky.fits',
|
|
'withregion=yes',
|
|
f'regionfile={region_fn}',
|
|
'keepinterfiles=yes',
|
|
'pattern=0',
|
|
'elow=350',
|
|
'ehigh=1100',
|
|
f'quads="{quads}"',]
|
|
|
|
w('pnspectra', inargs, ).run()
|
|
|
|
def run_mosback(key,ccds):
|
|
|
|
inargs = [f'inspecfile={key}-fovt.pi',
|
|
'elow=350',
|
|
'ehigh=1100',
|
|
f'ccds="{ccds}"',]
|
|
|
|
w('mosback', inargs, ).run()
|
|
|
|
def run_pnback(key,quads):
|
|
|
|
inargs = [f'inspecfile={key}-fovt.pi',
|
|
'elow=350',
|
|
'ehigh=1100',
|
|
f'quads="{quads}"',]
|
|
|
|
w('pnback', inargs).run()
|
|
|
|
def group_spectrum(key=None, group_min=25, oot=False):
|
|
if not (key):
|
|
print("Provide key= parameter, e.g. mos1S001")
|
|
return
|
|
|
|
test_exe('grppha')
|
|
specfile=f'{key}-fovt.pi' if not (oot) else f'{key}-fovtootsub.pi'
|
|
backfile=f'{key}-bkg.pi'
|
|
respfile=f'{key}.rmf'
|
|
ancrfile=f'{key}.arf'
|
|
fout=f'{specfile}.grp'
|
|
|
|
cmd=["grppha",
|
|
"infile={}".format(specfile),
|
|
"outfile={}".format(fout),
|
|
f"comm=\'chkey BACKFILE {backfile} & chkey RESPFILE {respfile} & chkey ANCRFILE {ancrfile} & GROUP MIN {group_min}\'",
|
|
"tempc=EXIT",
|
|
"clobber=yes",
|
|
"chatter=0", # with chatter <= 5 being very quite and chatter >= 20 very verbose.
|
|
]
|
|
print((" ").join(cmd))
|
|
os.system((" ").join(cmd))
|
|
|
|
fxcm = open(fout.replace("grp","xcm"), "w")
|
|
fxcm.write(f'data {fout}\n')
|
|
fxcm.write('ign **-2.\n')
|
|
fxcm.write('ign 11.-**\n')
|
|
fxcm.write('ign bad\n')
|
|
fxcm.write('cpd /xs\n')
|
|
fxcm.write('setpl en\n')
|
|
fxcm.write('@../../data/xspec/wabs_apec_cflux.xcm\n')
|
|
fxcm.write('pl lda delchi\n')
|
|
fxcm.close()
|