uds/scripts/03_init_obs.py
2023-03-09 21:10:39 +03:00

375 lines
13 KiB
Python
Executable File
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#!/usr/bin/env python
"""Подготавливает списки событий в разных энергетических диапазонах.
Производит списки источников в наждом наблюдении и делает астрокоррекцию с помощью wcs_match
"""
from astropy.wcs import WCS
from astropy.io import fits
import sys, os, os.path, time, subprocess
from pathlib import Path
import numpy as np
import glob
from os.path import dirname
import inspect
import uds
from uds.utils import *
from uds.config import *
""" find UDS root dir """
root_path=dirname(dirname(dirname(inspect.getfile(uds))))
print("UDS root path: {}".format(root_path))
infile_dir=root_path+'/data/processed'
outfile_dir=root_path+'/products'
create_folder(outfile_dir)
do_init = False
do_ermask = False
do_erbox1 = False # local mode
do_erbackmap1 = False #
do_erbox2 = False # map mode, with background map
do_erbackmap2 = False #
do_erbox3 = False # map mode, with background map
do_erbackmap3 = False #
do_ermldet = False
do_catprep = False
do_cross_match = True
do_wcs_match = False # Chandra task
do_wcs_update = False # Chandra task
eband_selected=[0]
def runme(datakey):
""" runs datakey over energy bands """
events=[]
expmaps=[]
outfile_boxlist1=[]
outfile_boxlist2=[]
outfile_boxlist3=[]
outfile_backmap1=[]
outfile_backmap2=[]
outfile_backmap3=[]
cheesemask=[]
bkgimage=[]
srcmaps=[]
for ii in range(len(eband_selected)):
index=eband_selected[ii]
print("\t>>> Energy band en{} -- {}-{} keV".format(eband[index],emin_kev[index],emax_kev[index]))
outfile_evtool, outfile_expmap = init_events(key=datakey, eband_index=eband[index],
infile_dir=infile_dir,
outfile_dir=outfile_dir,
do_init=do_init,
do_obsmode=True,
do_center=True,
do_evtool=True,
do_expmap=True,
ra_cen=ra_cen, de_cen=de_cen,
emin_kev=emin_kev[index],
emax_kev=emax_kev[index])
events.append(outfile_evtool)
expmaps.append(outfile_expmap)
""" Detmask """
detmask="{}_DetectionMask{}".format(os.path.join(outfile_dir,datakey), outfile_post)
if(do_ermask==True):
cmd=["ermask",
"expimage=%s" %(expmaps[0]), # use the first exposure maps calculated for that skyfield, independent of the energy band
"detmask=%s" %(detmask),
"threshold1=0.01",
"threshold2=10.0",
"regionfile_flag=no"
]
remove_file(detmask)
print((" ").join(cmd))
os.system((" ").join(cmd))
for ii in range(len(eband_selected)):
index=eband_selected[ii]
print("\t>>> Energy band en{} -- {}-{} keV".format(eband[index],emin_kev[index],emax_kev[index]))
""" erbox in local mode """
outfile_boxlist1.append("{}_BoxList1_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
if(do_erbox1==True):
""" erbox in local mode """
cmd=["erbox",
"images=%s" %(events[ii]),
"boxlist=%s" %(outfile_boxlist1[ii]),
"expimages=%s" %(expmaps[ii]),
"detmasks=%s" %(detmask),
"emin=%s" %(emin_ev[index]),
"emax=%s" %(emax_ev[index]),
"ecf=1.0",
"nruns=2",
"likemin=6.0",
"boxsize=4",
"compress_flag=N",
"bkgima_flag=N",
"expima_flag=Y",
"detmask_flag=Y"
]
remove_file(outfile_boxlist1[ii])
print((" ").join(cmd))
os.system((" ").join(cmd))
save_ds9reg(outfile_boxlist1[ii])
outfile_backmap1.append("{}_BackMap1_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
cheese_mask="{}_CheeseMask1_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
if(do_erbackmap1==True):
""" back map 1 """
cmd=["erbackmap",
"image=%s" %(events[ii]),
"expimage=%s" %(expmaps[ii]),
"boxlist=%s" %(outfile_boxlist1[ii]),
"detmask=%s" %(detmask),
"emin=%s" %(emin_ev[index]),
"emax=%s" %(emax_ev[index]),
"bkgimage=%s" %(outfile_backmap1[ii]),
"cheesemask=%s" %(cheese_mask),
"idband=1",
"scut=0.001",
"mlmin=6",
"maxcut=0.5",
"fitmethod=smooth smoothval=15",
"snr=40.",
]
remove_file(cheese_mask)
remove_file(outfile_backmap1[ii])
os.system((" ").join(cmd))
print((" ").join(cmd))
outfile_boxlist2.append("{}_BoxList2_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
if(do_erbox2==True):
""" erbox in background mode """
cmd=["erbox",
"images=%s" %(events[ii]),
"boxlist=%s" %(outfile_boxlist2[ii]),
"expimages=%s" %(expmaps[ii]),
"detmasks=%s" %(detmask),
"emin=%s" %(emin_ev[index]),
"emax=%s" %(emax_ev[index]),
"ecf=1.0",
"nruns=2",
"likemin=4.0",
"boxsize=4",
"compress_flag=N",
"bkgima_flag=Y",
"bkgimages={}".format(outfile_backmap1[ii]),
"expima_flag=Y",
"detmask_flag=Y"
]
remove_file(outfile_boxlist2[ii])
print((" ").join(cmd))
os.system((" ").join(cmd))
save_ds9reg(outfile_boxlist2[ii])
outfile_backmap2.append("{}_BackMap2_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
cheese_mask="{}_CheeseMask2_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
if(do_erbackmap2==True):
""" back map 2 """
cmd=["erbackmap",
"image=%s" %(events[ii]),
"expimage=%s" %(expmaps[ii]),
"boxlist=%s" %(outfile_boxlist2[ii]),
"detmask=%s" %(detmask),
"emin=%s" %(emin_ev[index]),
"emax=%s" %(emax_ev[index]),
"bkgimage=%s" %(outfile_backmap2[ii]),
"cheesemask=%s" %(cheese_mask),
"idband=1",
"scut=0.001",
"mlmin=6",
"maxcut=0.5",
"fitmethod=smooth smoothval=15",
"snr=40.",
]
remove_file(cheese_mask)
remove_file(outfile_backmap2[ii])
os.system((" ").join(cmd))
print((" ").join(cmd))
outfile_boxlist3.append("{}_BoxList3_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
if(do_erbox3==True):
""" erbox in map mode FINAL """
cmd=["erbox",
"images=%s" %(events[ii]),
"boxlist=%s" %(outfile_boxlist3[ii]),
"expimages=%s" %(expmaps[ii]),
"detmasks=%s" %(detmask),
"emin=%s" %(emin_ev[index]),
"emax=%s" %(emax_ev[index]),
"ecf=1.0",
"nruns=2",
"likemin=4.0",
"boxsize=4",
"compress_flag=N",
"bkgima_flag=Y",
"bkgimages={}".format(outfile_backmap2[ii]),
"expima_flag=Y",
"detmask_flag=Y"
]
remove_file(outfile_boxlist3[ii])
print((" ").join(cmd))
os.system((" ").join(cmd))
save_ds9reg(outfile_boxlist3[ii])
outfile_backmap3.append("{}_BackMap3_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
cheese_mask="{}_CheeseMask3_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
if(do_erbackmap3==True):
""" back map 3 FINAL """
cmd=["erbackmap",
"image=%s" %(events[ii]),
"expimage=%s" %(expmaps[ii]),
"boxlist=%s" %(outfile_boxlist3[ii]),
"detmask=%s" %(detmask),
"emin=%s" %(emin_ev[index]),
"emax=%s" %(emax_ev[index]),
"bkgimage=%s" %(outfile_backmap3[ii]),
"cheesemask=%s" %(cheese_mask),
"idband=1",
"scut=0.001",
"mlmin=6",
"maxcut=0.5",
"fitmethod=smooth smoothval=15",
"snr=40.",
]
remove_file(cheese_mask)
remove_file(outfile_backmap3[ii])
os.system((" ").join(cmd))
print((" ").join(cmd))
mllist="{}_MaxLikSourceList_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
srcmap="{}_SourceMap_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
cmd=["ermldet",
"mllist=%s" %(mllist),
"boxlist=%s" %(outfile_boxlist3[ii]),
"images=%s" %(events[ii]),
"expimages=%s" %(expmaps[ii]),
"detmasks=%s" %(detmask),
"bkgimages=%s" %(outfile_backmap3[ii]),
"emin=%s" %(emin_ev[index]),
"emax=%s" %(emax_ev[index]),
"hrdef=",
"ecf={}".format(ecf[index]),
"likemin=5.",
"extlikemin=6.",
"compress_flag=N",
"cutrad=15.",
"multrad=20.",
"extmin=2.0",
"extmax=15.0",
#"bkgima_flag=Y", looks outdated
"expima_flag=Y",
"detmask_flag=Y",
"extentmodel=beta",
"thres_flag=N",
"thres_col=like",
"thres_val=30.",
"nmaxfit=4",
"nmulsou=2",
"fitext_flag=yes",
"srcima_flag=yes",
"srcimages=%s" %(srcmap)
]
if(do_ermldet==True):
remove_file(mllist)
remove_file(srcmap)
os.system((" ").join(cmd))
print((" ").join(cmd))
#save_ermldet_ds9reg(mllist,scale=60*60)
catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
if(do_catprep==True):
cmd=["catprep",
"infile={}".format(mllist),
"outfile={}".format(catprep),]
remove_file(catprep)
os.system((" ").join(cmd))
print((" ").join(cmd))
save_catprep_ds9reg(catprep,scale=60*60)
if(do_cross_match==True):
#hdulist = fits.open(detmask)
#w = WCS(hdulist[0].header)
#w.wcs.crval=wcslist[datakey]
crossmatch_shu2019(catprep,dlmin=10,refimage=events[ii],crval=wcslist[datakey],
catalog=root_path+"/data/Gaia_unWISE/Gaia_unWISE_UDS.fits.catalog")
"""
if(do_att_corr==True):
attcorr(mllist,events=outfile_evtool[ii])
xfm=mllist.replace(".fits", ".xfm")
if(do_wcs_match==True):
src_list=mllist.replace(".fits", ".src")
if not (os.path.isfile(src_list)==True):
print("Not found {}".format(src_list))
sys.exit()
ref_list=mllist.replace(".fits", ".ref")
if not (os.path.isfile(ref_list)==True):
print("Not found {}".format(ref_list))
sys.exit()
#xfm=mllist.replace(".fits", ".xfm")
log=mllist.replace(".fits", ".xfm.log")
cmd=["wcs_match",
"infile={}".format(src_list),
"refsrcfile={}".format(ref_list),
"outfile={}".format(xfm),
"wcsfile={}".format(src_list),
"logfile={}".format(log),
"radius=7",
"residlim=1",
"verbose=1",
"method=trans",
#"method=rst",
"clobber=yes",
]
os.system((" ").join(cmd))
print((" ").join(cmd))
if(do_wcs_update==True):
wcs_update(outfile_evtool[ii],crval=wcslist[key],transformfile=xfm)
"""
runme("tm7_obs_1")
"""
for tmkey in keylist_tm.keys():
print("TM{} in work... init events".format(tmkey))
for datakey in keylist_tm[tmkey]:
print("--> {}".format(datakey))
runme(datakey)
"""