forked from xmm/arches
814 lines
28 KiB
Python
814 lines
28 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
|
|
import subprocess
|
|
|
|
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
|
|
from astropy.convolution import interpolate_replace_nans
|
|
|
|
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)
|
|
|
|
def move_files(pattern, destination):
|
|
for file in glob.glob(pattern):
|
|
shutil.move(file, os.path.join(destination,file))
|
|
def test_file(fn):
|
|
if not (os.path.isfile(fn)==True):
|
|
print(f'requested filename {fn} is not found')
|
|
cwd = os.getcwd()
|
|
print(f"Current directory: {cwd}")
|
|
sys.exit()
|
|
|
|
def get_first_file(pattern):
|
|
files = glob.glob(pattern)
|
|
if not (files):
|
|
print(f"Files not found from pattern {pattern}")
|
|
sys.exit()
|
|
filename = files[0]
|
|
return 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,products_dir=None):
|
|
|
|
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'
|
|
ds9reg_dir=root_path+'/data/ds9reg'
|
|
|
|
if (products_dir==None):
|
|
print('products_dir is not defined')
|
|
sys.exit()
|
|
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):
|
|
if(os.path.isfile(image)==False):
|
|
print(f"{image} is not found")
|
|
sys.exit()
|
|
|
|
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=False,elow="350",ehigh="1100"):
|
|
|
|
eventfile = f'{key}-allevc.fits'
|
|
if(os.path.isfile(eventfile)==False):
|
|
print(f"{eventfile} is not found, run 01_init_events*.py and above")
|
|
cwd = os.getcwd()
|
|
print(f"Current directory: {cwd}")
|
|
sys.exit()
|
|
|
|
unfilt_image = eventfile.replace("allevc","ffov")
|
|
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()
|
|
if(os.path.isfile(unfilt_image)==False):
|
|
print(f"{unfilt_image} is not found")
|
|
cwd = os.getcwd()
|
|
print(f"Current directory: {cwd}")
|
|
sys.exit()
|
|
|
|
if(check):
|
|
d = pyds9.DS9()
|
|
d.set("file "+unfilt_image)
|
|
d.set('cmap bb')
|
|
d.set('scale log')
|
|
#return
|
|
|
|
#(detx,dety)=get_det_coords(arches_ra,arches_dec,f'{key}-fovimt.fits',key)
|
|
(detx,dety)=get_det_coords(arches_ra,arches_dec,unfilt_image,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()
|
|
|
|
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()
|
|
|
|
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',
|
|
f'elow={elow}',
|
|
f'ehigh={ehigh}',
|
|
f'ccds="{ccds}"',]
|
|
|
|
w('mosspectra', inargs).run()
|
|
|
|
# check coordinate conversion
|
|
#(detx,dety)=get_det_coords(arches_ra,arches_dec,f'{key}-fovimt.fits',key) # image input in X,Y
|
|
#(detx,dety)=get_det_coords(arches_ra,arches_dec,unfilt_image,key) # image input in DETX,DETY
|
|
#sys.exit()
|
|
|
|
|
|
def run_bkgimsky(key,elow="350",ehigh="1100"):
|
|
# convert background image in det coordinates (e.g. mos1S001-bkgimdet-350-1100.fits) to sky coordinates
|
|
# runs after mosback!
|
|
bkgimsky = f'{key}-bkgimdet-{elow}-{ehigh}.fits'
|
|
inargs = [f'intemplate={key}-fovimsky-{elow}-{ehigh}.fits',
|
|
f'inimage={key}-bkgimdet-{elow}-{ehigh}.fits',
|
|
f'outimage={key}-bkgimsky-{elow}-{ehigh}.fits',
|
|
'withdetxy=false', 'withskyxy=false',]
|
|
w('rotdet2sky', inargs).run()
|
|
|
|
def run_pnspectra(key,quads,check=True,elow="350",ehigh="1100",pattern=0):
|
|
eventfile = f'{key}-allevc.fits'
|
|
if(os.path.isfile(eventfile)==False):
|
|
print(f"{eventfile} is not found, run 01_init_events*.py and above")
|
|
sys.exit()
|
|
|
|
# write image file with full fov
|
|
unfilt_image = eventfile.replace("allevc","ffov")
|
|
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()
|
|
|
|
if(check):
|
|
d = pyds9.DS9()
|
|
d.set("file "+unfilt_image)
|
|
d.set('cmap bb')
|
|
d.set('scale log')
|
|
#return
|
|
|
|
#(detx,dety)=get_det_coords(arches_ra,arches_dec,f'{key}-fovimt.fits',key)
|
|
(detx,dety)=get_det_coords(arches_ra,arches_dec,unfilt_image,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()
|
|
|
|
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()
|
|
|
|
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',
|
|
f'pattern={pattern}',
|
|
f'elow={elow}',
|
|
f'ehigh={ehigh}',
|
|
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 run_mossave(key,elow="350",ehigh="1100"):
|
|
dd=f'ffov_{key}'
|
|
create_folder(dd)
|
|
move_files(f'{key}*-{elow}-{ehigh}*',dd)
|
|
move_files(f'{key}*-{elow}-{ehigh}*',dd)
|
|
move_files(f'{key}*.pi',dd)
|
|
move_files(f'{key}*.qdp',dd)
|
|
move_files(f'{key}*.ps',dd)
|
|
move_files(f'{key}*imt*',dd)
|
|
move_files(f'{key}*.arf',dd)
|
|
move_files(f'{key}*.rmf',dd)
|
|
move_files(f'{key}*imspdet*',dd)
|
|
move_files(f'{key}*bkgimsky*',dd)
|
|
|
|
|
|
def run_pnsave(key,elow="350",ehigh="1100",pattern=0):
|
|
dd=f'ffov_{key}_{pattern}'
|
|
create_folder(dd)
|
|
move_files(f'pn*-{elow}-{ehigh}*',dd)
|
|
move_files(f'pn*-{elow}-{ehigh}*',dd)
|
|
move_files(f'pn*.pi',dd)
|
|
move_files(f'pn*.qdp',dd)
|
|
move_files(f'pn*.ps',dd)
|
|
move_files(f'pn*imt*',dd)
|
|
move_files(f'pn*.arf',dd)
|
|
move_files(f'pn*.rmf',dd)
|
|
move_files(f'pn*imspdet*',dd)
|
|
move_files(f'pn*bkgimsky*',dd)
|
|
|
|
|
|
def group_spectrum(key=None, group_min=25, oot=False, chdir=None, cpdir=None, pattern=None):
|
|
if not (key):
|
|
print("Provide key= parameter, e.g. mos1S001")
|
|
return
|
|
|
|
prefix=''
|
|
try:
|
|
if(chdir):
|
|
prefix='../'
|
|
os.chdir(chdir)
|
|
except:
|
|
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'
|
|
|
|
if(pattern!=None):
|
|
specfile0=specfile
|
|
backfile0=backfile
|
|
respfile0=respfile
|
|
ancrfile0=ancrfile
|
|
fout0=fout
|
|
|
|
specfile=f'{key}_{pattern}-fovt.pi' if not (oot) else f'{key}_{pattern}-fovtootsub.pi'
|
|
backfile=f'{key}_{pattern}-bkg.pi'
|
|
respfile=f'{key}_{pattern}.rmf'
|
|
ancrfile=f'{key}_{pattern}.arf'
|
|
fout=f'{specfile}.grp'
|
|
|
|
shutil.copyfile(specfile0, specfile)
|
|
shutil.copyfile(backfile0, backfile)
|
|
shutil.copyfile(respfile0, respfile)
|
|
shutil.copyfile(ancrfile0, ancrfile)
|
|
shutil.copyfile(fout0, fout)
|
|
|
|
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))
|
|
|
|
with fits.open(specfile) as hdul:
|
|
header = hdul['SPECTRUM'].header
|
|
backscal = header['BACKSCAL']
|
|
area_arcmin2 = float(backscal)*(det_pix_as/60)**2
|
|
print(f"Value of BACKSCAL: {backscal} area: {area_arcmin2:.4f} arcmin2")
|
|
|
|
fout_xcm=fout.replace("grp","xcm")
|
|
fxcm = open(fout_xcm, "w")
|
|
fxcm.write(f'# backscal {backscal} pix2\n')
|
|
fxcm.write(f'# area {area_arcmin2:.8f} arcmin2\n')
|
|
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(f'@{prefix}../../../data/xspec/wabs_apec_cflux.xcm\n')
|
|
fxcm.write('pl lda delchi\n')
|
|
fxcm.close()
|
|
|
|
if(cpdir):
|
|
shutil.copyfile(specfile, os.path.join('../',cpdir,specfile))
|
|
shutil.copyfile(backfile, os.path.join('../',cpdir,backfile))
|
|
shutil.copyfile(respfile, os.path.join('../',cpdir,respfile))
|
|
shutil.copyfile(ancrfile, os.path.join('../',cpdir,ancrfile))
|
|
shutil.copyfile(fout, os.path.join('../',cpdir,fout))
|
|
shutil.copyfile(fout_xcm, os.path.join('../',cpdir,fout_xcm))
|
|
|
|
if(chdir):
|
|
os.chdir('../')
|
|
|
|
|
|
|
|
def get_ds9_regions(image, src_fn, bkg_fn):
|
|
# Visualize the image with ds9
|
|
d = pyds9.DS9()
|
|
d.set("file "+image)
|
|
d.set('cmap bb')
|
|
d.set('scale log')
|
|
#d.set(f"region {ds9reg_dir}/arches.reg")
|
|
|
|
if(os.path.isfile(src_fn)==True):
|
|
d.set(f"region {src_fn}")
|
|
else:
|
|
print(f'src_fn: {src_fn} is not found')
|
|
sys.exit()
|
|
|
|
if(os.path.isfile(bkg_fn)==True):
|
|
d.set(f"region {bkg_fn}")
|
|
else:
|
|
print(f'bkg_fn: {bkg_fn} is not found')
|
|
sys.exit()
|
|
|
|
#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.
|
|
r2_bkg=None
|
|
if "annulus" in str(d.get("regions")):
|
|
r2_bkg = regionbkg.split(",")[3].replace('\n','')
|
|
print(" r2_bkg = ", r2_bkg, "(physical, annulus)")
|
|
return x_source,y_source,r_source,x_bkg,y_bkg,r_bkg,r2_bkg
|
|
|
|
|
|
def run_evqpb(key, work_dir=None, evtfile=None, attfile=None, pimin=[200,], pimax=[12000,], label=['',], exposurefactor=10.0):
|
|
|
|
test_file(evtfile)
|
|
test_file(attfile)
|
|
|
|
|
|
# Run the SAS task evqpb over each one of _original_ event files.
|
|
evqpb_outset=work_dir+f'/EPIC_{key}_QPB.fits' # Name of the output file
|
|
inargs = [f'table={evtfile}', f'attfile={attfile}',
|
|
f'outset={evqpb_outset}',
|
|
f'exposurefactor={exposurefactor}',]
|
|
w("evqpb", inargs).run()
|
|
|
|
gti=f'EPIC_{key}_gti.fit'
|
|
test_file(gti)
|
|
|
|
expression=f'gti({gti},TIME)'
|
|
# filter the EPIC event files to create cleaned and filtered for
|
|
# flaring particle background event files for your observation
|
|
sci_clean=work_dir+f'/EPIC_{key}_filtered.fits' # Name of the output file
|
|
inargs = [f'table={evtfile}', 'withfilteredset=yes',
|
|
f'filteredset={sci_clean}',
|
|
'destruct=yes','keepfilteroutput=true',
|
|
f'expression={expression}']
|
|
w("evselect", inargs).run()
|
|
|
|
# Use the same expression to filter the FWC event files:
|
|
qpb_clean=work_dir+f'/EPIC_{key}_QPB_clean.fits' # Name of the output file
|
|
inargs = [f'table={evqpb_outset}', 'withfilteredset=yes',
|
|
f'filteredset={qpb_clean}',
|
|
'destruct=yes','keepfilteroutput=true',
|
|
f'expression={expression}']
|
|
w("evselect", inargs).run()
|
|
|
|
for idx, lab in enumerate(label):
|
|
if(key=='PN'):
|
|
expression=f'(#XMMEA_EP)&&(PATTERN<=4)&&(PI>={pimin[idx]})&&(PI<={pimax[idx]})'
|
|
else:
|
|
expression=f'(#XMMEA_EM)&&(PATTERN<=12)&&(PI>={pimin[idx]})&&(PI<={pimax[idx]})'
|
|
print(lab,expression)
|
|
|
|
sci_image=work_dir+f'/EPIC_{key}_sci_image{lab}.fits' # Name of the output file
|
|
# Extract an image for the science exposure
|
|
inargs = [f'table={sci_clean}', f'expression={expression}',
|
|
'imagebinning=binSize', f'imageset={sci_image}',
|
|
'withimageset=yes','xcolumn=X', 'ycolumn=Y',
|
|
'ximagebinsize=80', 'yimagebinsize=80',
|
|
]
|
|
w("evselect", inargs).run()
|
|
|
|
qpb_image=work_dir+f'/EPIC_{key}_qpb_image{lab}.fits' # Name of the output file
|
|
# Extract an image for the FWC exposure
|
|
inargs = [f'table={qpb_clean}', f'expression={expression}',
|
|
'imagebinning=binSize', f'imageset={qpb_image}',
|
|
'withimageset=yes','xcolumn=X', 'ycolumn=Y',
|
|
'ximagebinsize=80', 'yimagebinsize=80','zcolumn=EWEIGHT',
|
|
]
|
|
w("evselect", inargs).run()
|
|
|
|
#hdul = fits.open(qpb_image)
|
|
#primary_hdu = hdul[0]
|
|
#qpb_data = hdul[0].data
|
|
#qpb_header = hdul[0].header
|
|
#hdul.close()
|
|
|
|
#result = interpolate_replace_nans(image, kernel)
|
|
|
|
# -- Subtract background --
|
|
# skip, due to low statistics
|
|
"""
|
|
cor_image=work_dir+f'/EPIC_{key}_cor_image{lab}.fits' # Name of the output file
|
|
cmd = ['farith', f'{sci_image}', f'{qpb_image}', f'{cor_image}', 'SUB',
|
|
'copyprime=yes',
|
|
'clobber=yes']
|
|
result = subprocess.run(cmd, capture_output=True, text=True, check=True)
|
|
"""
|
|
|
|
def run_mosaic(outfile_cts=None,outfile_qpb=None,outfile_exp=None,
|
|
outfile_sub=None,outfile_flx=None,outfile_pix=None,
|
|
cts=None,qpb=None,exp=None,nn=3,devmax=5.0,cutbox=100):
|
|
ref_crd = SkyCoord(arches_ra, arches_dec, frame=FK5(), unit="deg")
|
|
nn2=nn*nn
|
|
print(outfile_cts)
|
|
print(cts[0])
|
|
print(exp[0])
|
|
# take first image as reference
|
|
ref_hdul = fits.open(cts[0])
|
|
ref_data = ref_hdul[0].data
|
|
ref_header = ref_hdul[0].header
|
|
ref_hdul.close()
|
|
ref_wcs = WCS(ref_header)
|
|
|
|
(nx,ny) = ref_data.shape
|
|
map_flx = np.zeros(ref_data.shape,dtype=np.float64)
|
|
map_cts = np.zeros(ref_data.shape,dtype=np.float64)
|
|
map_qpb = np.zeros(ref_data.shape,dtype=np.float64)
|
|
map_exp = np.zeros(ref_data.shape,dtype=np.float64)
|
|
map_pix = np.zeros(ref_data.shape,dtype=np.float64) # number of subpixels inside
|
|
|
|
# take reference exposure
|
|
exp_hdul = fits.open(exp[0])
|
|
exp_data = exp_hdul[0].data
|
|
exp_header = exp_hdul[0].header
|
|
exp_hdul.close()
|
|
exp_wcs = WCS(exp_header)
|
|
|
|
"""
|
|
for row_index in range(len(ref_data)):
|
|
for col_index in range(len(ref_data[row_index])):
|
|
refx, refy = ref_wcs.world_to_pixel(ref_crd)
|
|
|
|
if(np.absolute(refx-row_index)>cutbox):
|
|
continue
|
|
if(np.absolute(refy-col_index)>cutbox):
|
|
continue
|
|
|
|
|
|
sky = ref_wcs.pixel_to_world(row_index, col_index)
|
|
|
|
exp_xx, exp_yy = exp_wcs.world_to_pixel(sky)
|
|
exp_x = int(np.round(exp_xx))
|
|
exp_y = int(np.round(exp_yy))
|
|
|
|
if not (exp_data[exp_y, exp_x]>0):
|
|
continue
|
|
map_cts[col_index][row_index]=1.0
|
|
|
|
sep_arcmin = sky.separation(ref_crd).arcmin
|
|
if(sep_arcmin > devmax):
|
|
continue
|
|
map_cts[col_index][row_index]=2.0
|
|
|
|
ref_hdul[0].data=map_cts
|
|
ref_hdul.writeto(outfile_cts,overwrite=True)
|
|
sys.exit()
|
|
"""
|
|
|
|
|
|
for idx,in_cts in enumerate(cts):
|
|
#if not (idx==2):
|
|
# continue
|
|
hdul0 = fits.open(in_cts)
|
|
c_data = hdul0[0].data
|
|
c_header = hdul0[0].header
|
|
c_wcs = WCS(c_header)
|
|
hdul0.close()
|
|
print(f"{in_cts} filter {c_header['FILTER']}")
|
|
|
|
hdul0 = fits.open(qpb[idx])
|
|
q_data = hdul0[0].data
|
|
q_header = hdul0[0].header
|
|
q_wcs = WCS(q_header)
|
|
hdul0.close()
|
|
|
|
hdul0 = fits.open(exp[idx])
|
|
e_data = hdul0[0].data
|
|
e_header = hdul0[0].header
|
|
e_wcs = WCS(e_header)
|
|
hdul0.close()
|
|
print(exp[idx])
|
|
|
|
|
|
xx, yy = c_wcs.world_to_pixel(ref_crd)
|
|
ref_x = int(np.rint(xx))
|
|
ref_y = int(np.rint(yy))
|
|
|
|
xmin = ref_x - cutbox
|
|
xmax = ref_x + cutbox
|
|
ymin = ref_y - cutbox
|
|
ymax = ref_y + cutbox
|
|
|
|
if(xmin<0):
|
|
xmin=0
|
|
if(xmax>=nx):
|
|
xmax=nx
|
|
if(ymin<0):
|
|
ymin=0
|
|
if(ymax>=ny):
|
|
ymax=ny
|
|
|
|
#for row_index in range(len(c_data)):
|
|
# for col_index in range(len(c_data[row_index])):
|
|
for row_index in range(xmin,xmax):
|
|
for col_index in range(ymin,ymax):
|
|
c_sky = c_wcs.pixel_to_world(row_index, col_index)
|
|
|
|
exp_xx, exp_yy = e_wcs.world_to_pixel(c_sky)
|
|
exp_x = int(np.rint(np.float64(exp_xx)))
|
|
exp_y = int(np.rint(np.float64(exp_yy)))
|
|
if not (e_data[exp_y, exp_x]>0):
|
|
continue
|
|
|
|
sep_arcmin = c_sky.separation(ref_crd).arcmin
|
|
if(sep_arcmin > devmax):
|
|
continue
|
|
|
|
#xx, yy = ref_wcs.world_to_pixel(c_sky)
|
|
#refx = int(np.round(xx))
|
|
#refy = int(np.round(yy))
|
|
|
|
# do subpixeling
|
|
for ii in range(nn):
|
|
for jj in range(nn):
|
|
xx0=np.float64(row_index)-0.5+(np.float64(ii+1)-0.5)/np.float64(nn)
|
|
yy0=np.float64(col_index)-0.5+(np.float64(jj+1)-0.5)/np.float64(nn)
|
|
subpix_sky = c_wcs.pixel_to_world(xx0, yy0)
|
|
|
|
sub_xx, sub_yy = ref_wcs.world_to_pixel(subpix_sky)
|
|
refx = np.int32(np.rint(np.float64(sub_xx)))
|
|
refy = np.int32(np.rint(np.float64(sub_yy)))
|
|
|
|
exp_xx, exp_yy = e_wcs.world_to_pixel(subpix_sky)
|
|
exp_x = np.int32(np.rint(np.float64(exp_xx)))
|
|
exp_y = np.int32(np.rint(np.float64(exp_yy)))
|
|
|
|
if (refx>=0 and refx<nx and refy>=0 and refy<ny):
|
|
map_cts[refy][refx]+=float(c_data[col_index][row_index])/nn2
|
|
map_qpb[refy][refx]+=float(q_data[col_index][row_index])/nn2
|
|
map_pix[refy][refx]+=1.0/nn2
|
|
map_exp[refy][refx]+=np.float64(e_data[exp_y][exp_x])/nn2
|
|
|
|
|
|
|
|
# Convert summed sky back to normal units
|
|
for x in range(nx):
|
|
for y in range(ny):
|
|
if(map_exp[y][x]>0):
|
|
map_flx[y][x]=(map_cts[y][x]-map_qpb[y][x])/map_exp[y][x]
|
|
map_cts[y][x]=map_cts[world_to_pixely][x]/map_pix[y][x]
|
|
map_qpb[y][x]=map_qpb[y][x]/map_pix[y][x]
|
|
map_exp[y][x]=map_exp[y][x]/map_pix[y][x]
|
|
else:
|
|
map_flx[y][x]=0.0
|
|
|
|
if(outfile_pix):
|
|
ref_hdul[0].data=map_pix
|
|
ref_hdul.writeto(outfile_pix,overwrite=True)
|
|
|
|
ref_hdul[0].data=map_cts
|
|
ref_hdul.writeto(outfile_cts,overwrite=True)
|
|
|
|
ref_hdul[0].data=map_exp
|
|
ref_hdul.writeto(outfile_exp,overwrite=True)
|
|
|
|
ref_hdul[0].data=map_qpb
|
|
ref_hdul.writeto(outfile_qpb,overwrite=True)
|
|
|
|
ref_hdul[0].data=map_flx
|
|
ref_hdul.writeto(outfile_flx,overwrite=True)
|
|
#sys.exit()
|
|
|
|
"""
|
|
cmd = ['farith', f'{outfile_cts}', f'{outfile_qpb}', f'{outfile_sub}', 'SUB',
|
|
'copyprime=yes', 'clobber=yes']
|
|
result = subprocess.run(cmd, capture_output=True, text=True, check=True)
|
|
|
|
cmd = ['farith', f'{outfile_sub}', f'{outfile_exp}', f'{outfile_rat}', 'DIV',
|
|
'copyprime=yes', 'clobber=yes']
|
|
result = subprocess.run(cmd, capture_output=True, text=True, check=True)
|
|
"""
|