1
0
forked from erosita/uds
This commit is contained in:
Roman Krivonos 2024-11-29 17:31:43 +03:00
parent 02a46712e3
commit 233967b1fd
31 changed files with 24620 additions and 61537 deletions

View File

@ -35,3 +35,8 @@ pip install --editable coma/
Непосредственная работа с данными происходит в директории scripts, где нужно последовательно запускать скрипты обработки. В данной директории находится подробное описание всех действий.
## Заметки
После фильтрации неба по флагу IKI/MPE перестают собираться общие списки событий с помощью evtool, появляется такое предупреждение:
evtool/compare_stdheader_eventfile: **WARNING0** Instruments in input files are different.

View File

@ -16,9 +16,9 @@ from astropy.table import QTable, Table, Column
from mpdaf.obj import WCS as mWCS
from sherpa.astro.ui import *
#from sherpa.astro.ui import *
from uds.config import *
from coma.config import *
from scipy.optimize import minimize
from random import uniform
@ -226,7 +226,7 @@ def do_astro_corr(src=None, ref=None, catalog=None):
fout=catalog.replace(".fits", ".shift-map.fits")
title=os.path.split(catalog)[1]
key=title.replace("_SourceCatalog_en0.fits","")
key=title.replace("_SourceCatalog_en00.fits","")
#png=catalog.replace(".fits", ".opti.png")
#png=png.replace(key,"{}_{}".format(obslist[key],key))
@ -331,7 +331,7 @@ def do_astro_corr(src=None, ref=None, catalog=None):
writer.write("{:.2f} {:.2f} {:.2f}\n".format(offset1[idx],sep1[idx],err1[idx]))
with open(log, 'w') as writer:
writer.write("{} {:+.2f} {:+.2f} {:+.2f} ({:.2f} - {:.2f}) = {:.2f} {:.2f} ({:.2f}) {:.2f} ({:.2f}) N={} -- {}\n".format(obslist[key],
writer.write("{} {:+.2f} {:+.2f} {:+.2f} ({:.2f} - {:.2f}) = {:.2f} {:.2f} ({:.2f}) {:.2f} ({:.2f}) N={} -- {}\n".format(key,
rota_opt,
shift_ra,
shift_dec,
@ -388,7 +388,7 @@ def do_astro_corr_rota(src=None, ref=None, catalog=None):
png=catalog.replace(".fits", ".rota.png")
log=catalog.replace(".fits", ".rota.log")
key=title.replace("_SourceCatalog_en0.fits","")
png=png.replace(key,"{}_{}".format(obslist[key],key))
#png=png.replace(key,"{}_{}".format(obslist[key],key))
print(src)
hdul = fits.open(src)
@ -441,14 +441,14 @@ def do_astro_corr_rota(src=None, ref=None, catalog=None):
rota_min=angle*60
with open(log, 'w') as writer:
writer.write("{} {} {:.2f} {:.2f} {:.2f} {:.2f}%\n".format(obslist[key],title,rota_min, rms_min, resid, (1.0-rms_min/resid)*100.0))
writer.write("{} {} {:.2f} {:.2f} {:.2f} {:.2f}%\n".format(key,title,rota_min, rms_min, resid, (1.0-rms_min/resid)*100.0))
fig, ax = plt.subplots( nrows=1, ncols=1 ) # create figure & 1 axis
ax.plot(x,y)
plt.ylabel('RMS/RMS0, RMS0={:.2f} arcsec'.format(resid))
plt.xlabel('Rotation, arcmin')
plt.title("{} {} {} pairs, {:.1f}%".format(obslist[key], title,len(dsrc),(1.0-rms_min/resid)*100.0))
plt.title("{} {} {} pairs, {:.1f}%".format(key, title,len(dsrc),(1.0-rms_min/resid)*100.0))
plt.grid(True)
plt.axvline(x = 0.0, color = 'b', label = 'zero')
fig.savefig(png)

View File

@ -1,5 +1,9 @@
#from pathlib import Path
#work_dir = Path('/path/to/some/logical/parent/dir')
"""
DEPRECATION: Legacy editable install of coma==0.1 from file:///data/erosita/work/coma/coma (setup.py develop) is deprecated. pip 25.0 will enforce this behaviour change. A possible replacement is to add a pyproject.toml or enable --use-pep517, and use setuptools >= 64. If the resulting installation is not behaving as expected, try using --config-settings editable_mode=compat. Please consult the setuptools documentation for more information. Discussion can be found at https://github.com/pypa/pip/issues/11457
"""
"""
Координаты сентрального кадра, к которому будут
@ -7,6 +11,7 @@
"""
ra_cen=194.9601418
de_cen=27.6857267
width=10000
"""
Словарь камер со списком наблюдений каждой камеры.
@ -20,6 +25,39 @@ keylist_tm={'1':['tm1_scan1','tm1_scan2','tm1_scan3','tm1_scan4','tm1_scan5','tm
'6':['tm6_scan1','tm6_scan2','tm6_scan3','tm6_scan4','tm6_scan5','tm6_scan6','tm6_scan7','tm6_survey'],
'7':['tm7_scan1','tm7_scan2','tm7_scan3','tm7_scan4','tm7_scan5','tm7_scan6','tm7_scan7','tm7_survey'],}
keylist_survey={'1':['tm1_survey'],
'2':['tm2_survey'],
'3':['tm3_survey'],
'4':['tm4_survey'],
'5':['tm5_survey'],
'6':['tm6_survey'],
'7':['tm7_survey'],}
keylist_parts={'1':['tm1_partI','tm1_partII',],
'2':['tm2_partI','tm2_partII',],
'3':['tm3_partI','tm3_partII',],
'4':['tm4_partI','tm4_partII',],
'5':['tm5_partI','tm5_partII',],
'6':['tm6_partI','tm6_partII',],
'7':['tm7_partI','tm7_partII',],}
keylist_partI={'1':['tm1_partI',],
'2':['tm2_partI',],
'3':['tm3_partI',],
'4':['tm4_partI',],
'5':['tm5_partI',],
'6':['tm6_partI',],
'7':['tm7_partI',],}
keylist_partII={'1':['tm1_partII',],
'2':['tm2_partII',],
'3':['tm3_partII',],
'4':['tm4_partII',],
'5':['tm5_partII',],
'6':['tm6_partII',],
'7':['tm7_partII',],}
"""
Примерные центры изображений каждого наблюдения.
Требуется для астрокоррекции. Будет расчитываться матрица сдвигов и поворотов,
@ -42,6 +80,8 @@ wcslist={'tm1_obs_1':[34.7279760,-5.0680267],
'tm6_scan_3':[34.5405596,-4.8088748],
'tm6_scan_4':[34.5405596,-4.8088748]}
wcslist={'tm1_scan1':[194.9601418,27.6857267],}
""" like in the paper (Table 1) """
obslist={'tm1_obs_1':'N01',
'tm5_obs_1':'N02',
@ -61,48 +101,26 @@ obslist={'tm1_obs_1':'N01',
'tm6_scan_4':'N11'}
""" Это просто индекс диапазона для выходных файлов. """
eband=[ "0", "1", "2", "3", "4", "5", "6", "7", "8"]
eband=[ 0, 1, 2, 3, 4, 5, 6, 7,]
""" Диапазоны энергий. """
emin_ev=[ 300, 300, 600, 2300, 200, 300, 5000, 500, 1000]
emax_ev=[2300, 600, 2300, 5000, 10000,8000, 8000, 1000, 2000]
emin_ev=[ 200, 200, 600, 2300, 5000, 500, 1000, 200,]
emax_ev=[2300, 600, 2300, 5000, 8000, 1000, 2000, 10000,]
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
emin_kev=[0.2, 0.2, 0.6, 2.3, 5.0, 0.5, 1.0, 0.2,]
emax_kev=[2.3, 0.6, 2.3, 5.0, 8.0, 1.0, 2.0, 10.0,]
# Use Pavel's file ../data/ECF/ecf_tbabspow_g2nh0.02.pkl
# see 00_print_ecf.py
ecf = [1.05717e+12, # E0 0.2 - 2.3 keV
9.80313e+11, # E1 0.2 - 0.6 keV
1.08963e+12, # E2 0.6 - 2.3 keV
1.14674e+11, # E3 2.3 - 5.0 keV
2.77581e+10, # E4 5.0 - 8.0 keV
1.35463e+12, # E5 0.5 - 1.0 keV 4XMM-DR12 EP2 band
1.01476e+12, # E6 1.0 - 2.0 keV 4XMM-DR12 EP3 band
6.40472e+11, # E7 0.2 - 10.0 keV
]
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
"""

View File

@ -16,16 +16,19 @@ import pickle
from astropy.table import QTable, Table, Column
from astropy import units as u
from astropy_healpix import HEALPix
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
import astropy.units as u
from astropy.time import Time, TimeDelta, TimezoneInfo, TimeFromEpoch
import statistics
import shutil
#from coma.astrometry import *
from coma.astrometry import *
MJDREF = 51543.875
TZ_UTC = TimezoneInfo(utc_offset=0*u.hour)
@ -93,7 +96,7 @@ def convert_dr12_to_erosita_flux(rec, field_prefix='SC_'):
return (flux, error)
def do_evtool_esass(evfile=None,events=None,outfile=None,evlist=None,
def do_evtool_esass(evfile=None,events=None,outfile=None,evlist=None,rusky=None,
gti=None,region=None,emin=None,emax=None, rmlock=False,
do_center=False, width=1024, ra_cen=None, de_cen=None):
@ -133,7 +136,7 @@ def do_evtool_esass(evfile=None,events=None,outfile=None,evlist=None,
"outfile={}".format(outfile),
"image=yes",
"size='{}'".format(width),
"flag=0x2000",
"flag=0x2002" if(rusky == True) else "flag=0x2000",
"pattern=15"
]
@ -149,6 +152,27 @@ def do_evtool_esass(evfile=None,events=None,outfile=None,evlist=None,
os.remove(filename)
pass
def do_check_events(evfile=None,events=None,outfile=None,evlist=None,rusky=None,
gti=None,region=None,emin=None,emax=None, rmlock=False,
do_center=False, width=1024, ra_cen=None, de_cen=None):
eventfiles=None
if(events):
eventfiles=events
if(evfile):
eventfiles=[evfile]
failed = False
for infile in eventfiles:
if os.path.isfile(infile):
with fits.open(infile) as hdu:
for h in hdu:
extname = h.header['EXTNAME'] if "EXTNAME" in h.header else "None"
#if not (h.header['OBS_MODE'] == 'SURVEY'):
print("{}: \t EXTNAME={} \t OBS_MODE={} \t INSTRUME={}".format(infile,extname,h.header['OBS_MODE'],h.header['INSTRUME'],))
#failed = True
#if(failed == True):
# print("All OBS_MODE's must be SURVEY")
# #sys.exit()
def set_bit(value, bit):
return value | (1<<bit)
@ -248,7 +272,7 @@ def init_events(key=None, eband_selected=[0], eband_index=None,
with fits.open(infile) as hdu:
for h in hdu:
extname = h.header['EXTNAME'] if "EXTNAME" in h.header else "None"
changeto = 'SURVEY' if ('scan' in infile) else 'POINTING'
changeto = 'SURVEY' #if ('scan' in infile) else 'POINTING'
if not (h.header['OBS_MODE'] == changeto):
print("*** {} {} --> {}".format(extname,h.header['OBS_MODE'],changeto))
h.header['OBS_MODE'] = changeto
@ -279,7 +303,7 @@ def init_events(key=None, eband_selected=[0], eband_index=None,
pass
outfile_evtool="{}/{}_EventList_en{}{}".format(outfile_dir, key, eband_index, outfile_post_events)
outfile_evtool="{}/{}_EventList_en{:02d}{}".format(outfile_dir, key, eband_index, outfile_post_events)
cmd=["evtool",
"eventfiles=%s" %(infile),
@ -302,7 +326,7 @@ def init_events(key=None, eband_selected=[0], eband_index=None,
with fits.open(outfile_evtool) as hdu:
for h in hdu:
extname = h.header['EXTNAME'] if "EXTNAME" in h.header else "None"
changeto = 'SURVEY' if ('scan' in infile) else 'POINTING'
changeto = 'SURVEY' #if ('scan' in infile) else 'POINTING'
if not (h.header['OBS_MODE'] == changeto):
print("*** {} {} --> {}".format(extname,h.header['OBS_MODE'],changeto))
h.header['OBS_MODE'] = changeto
@ -310,10 +334,10 @@ def init_events(key=None, eband_selected=[0], eband_index=None,
if(vign==True):
withvignetting='yes'
outfile_expmap="{}_ExposureMap_en{}.vign{}".format(os.path.join(outfile_dir,key), eband_index, outfile_post)
outfile_expmap="{}_ExposureMap_en{:02d}.vign{}".format(os.path.join(outfile_dir,key), eband_index, outfile_post)
else:
withvignetting='no'
outfile_expmap="{}_ExposureMap_en{}.novign{}".format(os.path.join(outfile_dir,key), eband_index, outfile_post)
outfile_expmap="{}_ExposureMap_en{:02d}.novign{}".format(os.path.join(outfile_dir,key), eband_index, outfile_post)
if(do_expmap==True and do_init==True):
cmd=["expmap",
@ -349,7 +373,9 @@ def do_expmap_esass(expmaps=None,outfile=None,emin_kev=None,emax_kev=None,refima
def do_fimgmerge_ftools(maps=None,outfile=None):
""" takes first map as reference and merges the rest """
tmp="maps.txt"
tmp="/tmp/{}.maps.txt".format(os.getpid())
f = open(tmp, "w")
for jj in range(1,len(maps)):
f.write("{},0,0\n".format(maps[jj]))
@ -363,6 +389,9 @@ def do_fimgmerge_ftools(maps=None,outfile=None):
print((" ").join(cmd))
test_exe('fimgmerge')
os.system((" ").join(cmd))
if(os.path.isfile(tmp)==True):
os.remove(tmp)
def runme(cmd, local_run=True, cwd='/srg/work/krivonos/erosita/UDS-paper/uds/scripts'):
@ -445,7 +474,7 @@ def save_catprep_ds9reg(filename,scale=60,id_instr=1,id_band=1,ext_like=0.0,labe
writer.write("fk5;circle({}, {}, {}) # text={{{:.1f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale,rec[label]))
def crossmatch_shu2019(filename,refimage=None,crval=None,devmax=30,dlmin=6.0,dlmax=10000,ext_like=0.0,outkey='shu2019', catalog=None, errlim=10.0):
def crossmatch_shu2019(filename,refimage=None,crval=None,devmax=30,dlmin=6.0,dlmax=10000,ext_like=0.0,outkey='shu2019', catalog=None, errlim=10.0, nside=1024, order='nested', datakey=None):
if(os.path.isfile(filename)==False):
print("File not found {}".format(filename))
print("Start cross-match with Gaia-unWISE")
@ -464,22 +493,28 @@ def crossmatch_shu2019(filename,refimage=None,crval=None,devmax=30,dlmin=6.0,dlm
w.wcs.crval=crval
naxis1=hdulist[0].header['NAXIS1']
naxis2=hdulist[0].header['NAXIS2']
print("Reading reference image {}x{}".format(naxis1,naxis2))
if not catalog:
print("ERROR: Please provide reference catalog: catalog=")
return
hp = HEALPix(nside=nside, order=order, frame=FK5())
hdul = fits.open(catalog)
tbdata_ref = hdul[1].data
cat=[]
for rec in tbdata_ref:
cat.append({"ra":rec['ra'],"dec":rec['dec'],
crd = SkyCoord(float(rec['ra']), float(rec['dec']), frame="fk5", unit="deg")
heal = hp.skycoord_to_healpix(crd)
cat.append({"ra":rec['ra'],"dec":rec['dec'],"healpix":heal,
"pmra_err":rec['pmra_err'] if(rec['pmra_err']>0) else 1.0,
"pmdec_err":rec['pmdec_err'] if(rec['pmdec_err']>0) else 1.0,
"gaia":rec["gaia_sourceid"],"unwise":rec["unwise_objid"]})
print("Reading {}".format(filename))
hdul = fits.open(filename)
tbdata_src = hdul[1].data
@ -517,9 +552,12 @@ def crossmatch_shu2019(filename,refimage=None,crval=None,devmax=30,dlmin=6.0,dlm
# do correlation
for src in srcs:
src_crd = SkyCoord(src['ra'], src['dec'], frame=FK5(), unit="deg")
src_heal = hp.skycoord_to_healpix(src_crd)
sep=0
for scat in cat:
if(scat['healpix'] != src_heal):
continue
scat_crd = SkyCoord(scat['ra'], scat['dec'], frame=FK5(), unit="deg")
sep = src_crd.separation(scat_crd).arcsec
if(sep < devmax):
@ -560,15 +598,20 @@ def crossmatch_shu2019(filename,refimage=None,crval=None,devmax=30,dlmin=6.0,dlm
print("Take ID={} DET_LIKE={:.2f}".format(rec['id'],rec['det_like']))
src_crd = SkyCoord(rec['src_ra'], rec['src_dec'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(rec['ref_ra'], rec['ref_dec'], frame=FK5(), unit="deg")
src_x, src_y = w.world_to_pixel(src_crd)
#print(rec['src_ra'],rec['src_dec'])
#print(rec['ref_ra'],rec['ref_dec'])
#print(src_crd)
src_x, src_y = w.world_to_pixel(src_crd)
src_map[int(np.round(src_y)), int(np.round(src_x))] = 1
#print(rec['src_ra'], rec['src_dec'],"-->",src_x, src_y)
# check
#sky = w.pixel_to_world(src_x, src_y)
#print(sky)
#print(rec['src_ra'], rec['src_dec'],"-->",src_x, src_y)
#sys.exit()
ref_x, ref_y = w.world_to_pixel(ref_crd)
@ -593,7 +636,8 @@ def crossmatch_shu2019(filename,refimage=None,crval=None,devmax=30,dlmin=6.0,dlm
date_end = mission2date_utc(tstop)
stat_fout=fout.replace(".cross.fits", ".dat")
with open(stat_fout, 'w') as writer:
writer.write("[DATE-OBS {}, DATE-END {}] nref, mean [arcsec] dRA, dDec, mean/median sep (mean/median): {} & {:.4f} & {:.4f} & {:.4f}/{:.4f}\n".format(
writer.write("{} [DATE-OBS {}, DATE-END {}] nref, mean [arcsec] dRA, dDec, mean/median sep (mean/median): {} & {:.4f} & {:.4f} & {:.4f}/{:.4f}\n".format(
datakey,
date_obs,date_end,len(delta_ra),
statistics.mean(delta_ra)*3600,
statistics.mean(delta_dec)*3600,
@ -601,6 +645,7 @@ def crossmatch_shu2019(filename,refimage=None,crval=None,devmax=30,dlmin=6.0,dlm
statistics.median(delta_sep)))
src_tab.write(src_fout, format='fits', overwrite=True)
with fits.open(src_fout) as f:
f[0].verify('fix')
f[1].verify('fix')
@ -617,7 +662,7 @@ def crossmatch_shu2019(filename,refimage=None,crval=None,devmax=30,dlmin=6.0,dlm
f[0].header=f[0].header+w.to_header()
f[0].data=ref_map
#f[1].header=f[1].header+w.to_header()
def crossmatch_dr12(filename,devmax=30,ext_like=0.0,outkey='dr12', catalog=None):
flux_scale=1e-15
if(os.path.isfile(filename)==False):
@ -810,9 +855,13 @@ def make_rate_map(cntmap=None, expmap=None, outfile=None):
if not (cntmap and expmap and outfile):
print("ERROR: Please provide counts, exposure ane output file")
return
hdu_list_cntmap = fits.open(cntmap)
#hdu_list_expmap = fits.open(expmap)
return
cntmap_data, cntmap_hdr = fits.getdata(cntmap, ext=0, header=True)
expmap_data, expmap_hdr = fits.getdata(expmap, ext=0, header=True)
return
rate = np.zeros(expmap_data.shape)
index = np.where(expmap_data > 0.0)
@ -945,11 +994,13 @@ def wcs_update_ciao(events,ext=2,crval=None,transformfile=None,clean=True):
"""
Second, calculate new RA/Dec for events using updated attitude, using evatt from eSASS
Below command runs with CAPITAL keywords in early versions of eSASS
"""
cmd=["evatt",
"EVENTFILE={}".format(fnattcor),
"ATTFILE={}".format(fnattcor),
"GTIFILE={}".format(fnattcor),
"eventfile={}".format(fnattcor),
"attfile={}".format(fnattcor),
"gtifile={}".format(fnattcor),
]
os.system((" ").join(cmd))
print((" ").join(cmd))
@ -959,7 +1010,7 @@ def wcs_update_ciao(events,ext=2,crval=None,transformfile=None,clean=True):
remove_file(fntmp)
return fnattcor
def wcs_update_shift(events,flog=None,ext=2,crval=None,clean=True):
def wcs_update_shift(events,flog=None,ext=2,crval=None,clean=True, rusky=None):
test_exe('evatt')
fnattcor=events.replace(".fits", ".attcorr.fits")
@ -1120,7 +1171,7 @@ def filter_catprep(filename, expcut=100,dlmin=6.0,dlmax=10,scale=60*60,ext_like=
'det_like':rec['det_like_0'],'ext_like':rec['ext_like'],
'ext':rec['ext'],'ext_err':rec['ext_err'],
'src_id':rec['id_src']})
print("{:.2f} {} {} & {:.4f} & {:.4f} & {:.2f} & {:.2f} & {:.2f} $\pm$ {:.2f} \\\\".format(rec['ra'],rec['id_src'],
print("{:.2f} {} {} & {:.4f} & {:.4f} & {:.2f} & {:.2f} & {:.2f} $\\pm$ {:.2f} \\\\".format(rec['ra'],rec['id_src'],
make_source_name(rec['ra'], rec['dec']), rec['ra'], rec['dec'],
rec['det_like_0'], rec['ext_like'],
rec['ext'], rec['ext_err']))
@ -1252,7 +1303,7 @@ def create_sensmap(sensmap=None,expmap=None,backmap=None,detmask=None,areatab=No
"area_table={}".format(areatab),
"method=APER",
"aper_type=BOX",
"aper_size=4.5",
"aper_size=10.5",
"likemin={}".format(detlike),
"extlikemin=6.0",
"ext=6.0",
@ -2180,7 +2231,7 @@ def make_extended(infile=None, elmin=5.0, outreg=None, scale=3600):
print("Forced: ", rec['forced_name'])
ext_like = rec['ext_like']
if(ext_like>elmin):
print("{:.4f} {} {} & {:.4f} & {:.4f} & {:.2f} & {:.2f} & {:.2f} $\pm$ {:.2f} \\\\".format(rec['ext_like'],rec['src_id'],
print("{:.4f} {} {} & {:.4f} & {:.4f} & {:.2f} & {:.2f} & {:.2f} $\\pm$ {:.2f} \\\\".format(rec['ext_like'],rec['src_id'],
rec['name'], rec['ra'], rec['dec'],
rec['det_like'], rec['ext_like'],
rec['ext'], rec['ext_err']))

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

Binary file not shown.

View File

@ -1,19 +0,0 @@
from astropy.io import fits
import sys
#filename='4XMM_DR12cat_slim_v1.0_UDS.fits.catalog'
filename='4XMM_slim_DR13cat_v1.0_DEMO.fits'
fout=filename.replace(".fits.catalog", ".names.reg")
hdul = fits.open(filename)
#hdul.info()
tbdata = hdul[1].data
hdul.close()
#with open("./{}".format(fout), 'w') as writer:
for rec in tbdata:
print("fk5;circle({}, {}, 0.008)".format(rec['sc_ra'],rec['sc_dec']))
#writer.write("fk5;circle({}, {}, {}) # text={{{}}}\n".format(rec['sc_ra'],rec['sc_dec'],rec['sc_poserr']/3600,rec['IAUNAME']))

View File

@ -5,7 +5,8 @@
Для вырезки источников для работы с полем UDS была использована эта команда:
```
ftselect '4XMM_DR12cat_slim_v1.0.fits[1]' 4XMM_DR12cat_slim_v1.0_UDS.fits 'SC_RA > 32.75 && SC_RA < 36.31 && SC_DEC > -6.55 && SC_DEC < -3.0' clobber=yes
ftselect '4XMM_slim_DR14cat_v1.0.fits[1]' 4XMM_slim_DR14cat_v1.0_Coma.fits 'SC_RA > 185 && SC_RA < 205 && SC_DEC > 20 && SC_DEC < 34.0' clobber=yes
ftselect '4XMM_DR12cat_v1.0.fits[1]' 4XMM_DR12cat_v1.0_UDS.fits 'SC_RA > 32.75 && SC_RA < 36.31 && SC_DEC > -6.55 && SC_DEC < -3.0' clobber=yes
Demo with DR13:

File diff suppressed because one or more lines are too long

19
data/4XMM/print_ds9reg.py Executable file
View File

@ -0,0 +1,19 @@
from astropy.io import fits
import sys
#filename='4XMM_DR12cat_slim_v1.0_UDS.fits.catalog'
filename='4XMM_slim_DR14cat_v1.0_Coma.fits'
fout=filename.replace(".fits", ".names.reg")
hdul = fits.open(filename)
#hdul.info()
tbdata = hdul[1].data
hdul.close()
with open("./{}".format(fout), 'w') as writer:
for rec in tbdata:
print("fk5;circle({}, {}, 0.008)".format(rec['sc_ra'],rec['sc_dec']))
writer.write("fk5;circle({}, {}, {}) # text={{{}}}\n".format(rec['sc_ra'],rec['sc_dec'],rec['sc_poserr']/3600,rec['IAUNAME']))

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View File

@ -2,9 +2,9 @@
Каталог был получен с этого ресурса: https://zenodo.org/record/6837642#.ZAcufYBByRR
Для вырезки источников для работы с полем UDS была использована эта команда:
Для вырезки источников для работы с полем Coma была использована эта команда:
```
ftselect 'Gaia_unWISE_AGNs.fits[1]' Gaia_unWISE_UDS.fits 'RA > 32.75 && RA < 36.31 && DEC > -6.55 && DEC < -3.0' clobber=yes
ftselect 'Gaia_unWISE_AGNs.fits[1]' Gaia_unWISE_Coma.fits 'RA > 185 && RA < 205 && DEC > 20 && DEC < 34.0' clobber=yes
```

File diff suppressed because one or more lines are too long

View File

@ -0,0 +1,18 @@
from astropy.io import fits
import sys
filename='Gaia_unWISE_Coma.fits.catalog'
fout=filename.replace(".fits.catalog", ".reg")
hdul = fits.open(filename)
#hdul.info()
tbdata = hdul[1].data
hdul.close()
with open("./{}".format(fout), 'w') as writer:
for rec in tbdata:
print("fk5;circle({}, {}, 0.008)".format(rec['RA'],rec['DEC']))
writer.write("fk5;circle({}, {}, {})\n".format(rec['RA'],rec['DEC'],0.0080000))

96
data/eRASS1/print_ds9reg.py Executable file
View File

@ -0,0 +1,96 @@
from astropy.io import fits
import sys, math
from astropy_healpix import HEALPix
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
import astropy.units as u
from coma.config import *
def load_erass1_fits(filename, names=False):
cen_crd = SkyCoord(ra_cen, de_cen, frame=FK5(), unit="deg")
rmax=7.0 # deg
#filename="erass1-m.fits.gz"
if(names==True):
fout=filename.replace(".fits.gz", ".names.reg")
else:
fout=filename.replace(".fits.gz", ".reg")
with open("./{}".format(fout), 'w') as writer:
with fits.open(filename, memmap=True) as hdul:
hdul.info()
data = hdul[1].data
header = hdul[1].header
cols = hdul[1].columns
print(cols)
#sys.exit()
for r in data:
#IdCluster = int(r['IdCluster'])
#UID = int(r['UId'])
#UIDhard = int(r['UIdHard'])
#RADEC_ERR = float(r['RADEErr'])
#EXT = float(r['EXT'])
#s_EXT = r['s_EXT']
#EXT_LIKE = r['EXTLike']
#DET_LIKE = r['DetLike0']
#ML_CTS = r['MLcts1']
#s_ML_CTS = r['s_MLcts1']
#ML_RATE = r['MLRate1']
#s_ML_RATE = r['s_MLRate1']
#ML_FLUX = r['MLFlux1']
#s_ML_FLUX = r['s_MLFlux1']
#ML_EXP = r['MLExp1']
name=r['IAUName']
poserr=float(r['posErr']) # 1-sigma positional uncertainty
ra=float(r['RAdeg'])
dec=float(r['DEdeg'])
P=0.98
koeff=math.sqrt(-2*math.log(1-P))
crd = SkyCoord(ra, dec, frame=FK5(), unit="deg")
sep=crd.separation(cen_crd).deg
if(sep < rmax):
print(ra,dec,sep)
if(names==True):
writer.write("fk5;circle({}, {}, {}) # text={{{}}}\n".format(ra,dec,poserr/60/60, name))
else:
writer.write("fk5;circle({}, {}, {})\n".format(ra,dec,poserr/60/60))
load_erass1_fits("erass1-m.fits.gz")
load_erass1_fits("erass1-s.fits.gz")
load_erass1_fits("erass1-h.fits.gz")
load_erass1_fits("erass1-m.fits.gz", names=True)
load_erass1_fits("erass1-s.fits.gz", names=True)
load_erass1_fits("erass1-h.fits.gz", names=True)
"""
#filename='4XMM_DR12cat_slim_v1.0_UDS.fits.catalog'
filename='4XMM_slim_DR14cat_v1.0_Coma.fits'
fout=filename.replace(".fits", ".names.reg")
hdul = fits.open(filename)
#hdul.info()
tbdata = hdul[1].data
hdul.close()
with open("./{}".format(fout), 'w') as writer:
for rec in tbdata:
print("fk5;circle({}, {}, 0.008)".format(rec['sc_ra'],rec['sc_dec']))
writer.write("fk5;circle({}, {}, {}) # text={{{}}}\n".format(rec['sc_ra'],rec['sc_dec'],rec['sc_poserr']/3600,rec['IAUNAME']))
"""

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

98
scripts/00_print_ecf.py Executable file
View File

@ -0,0 +1,98 @@
#!/usr/bin/env python
"""
НАЗВАНИЕ:
05_srctool.py
НАЗНАЧЕНИЕ:
Запускает scrtool для самого широкого канала 0.2-10 кэВ, чтобы спектры имели самое полное покрытие по энергиям. Список источников берется из 0.3-2.3 кэВ.
ВЫЗОВ:
esass
./05_srctool.py
УПРАВЛЕНИЕ:
Требуется запуск предыдущего скрипта 04_mosaics.py
ПАРАМЕТРЫ:
index=4 : Выбранный энергетический диапазон
ВЫВОД:
Выходные файлы записываются в директорию outfile_dir/srctool_dir
ИСТОРИЯ:
Роман Кривонос, ИКИ РАН, krivonos@cosmos.ru
Март 2023
"""
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 coma
from scipy.stats import norm
import matplotlib.pyplot as plt
import pandas as pd
from coma.utils import *
from coma.config import *
""" find UDS root dir """
#root_path=dirname(dirname(dirname(inspect.getfile(uds))))
"""
ftools does not like long file path names,
for this reason, we use relative path here
"""
root_path='..'
print("Coma root path: {}".format(root_path))
infile_dir=root_path+'/data/processed'
outfile_dir=root_path+'/products'
create_folder(outfile_dir)
outkey="tm0"
outfile_srctool="{}_SrcTool_".format(outkey)
do_print_ecf = True
if(do_print_ecf==True):
filename='../data/ECF/ecf_tbabspow_g2nh0.02.pkl'
with open(filename, 'rb') as f:
ecf_table = pickle.load(f)
"""
for key in table.keys():
print("{} --> {}".format(key,table[key]))
"""
print(ecf_table[(0.3,2.3)])
print(ecf_table[(0.3,0.6)])
print(ecf_table[(0.6,2.3)])
print(ecf_table[(2.3,5.0)])
print(ecf_table[(5.0,8.0)])
print()
print(ecf_table[(0.5,1.0)]) # 4XMM-DR12 EP2 band
print(ecf_table[(1.0,2.0)]) # 4XMM-DR12 EP3 band
print()
print(ecf_table[(0.2,2.3)]) #
print(ecf_table[(0.2,0.6)]) #
print()
print(ecf_table[(0.2,10.0)]) #

View File

@ -25,10 +25,19 @@ region="box({},{},12d,12d,0)".format(ra_cen,de_cen) # Selection region
tm_list=[1,2,3,4,5,6,7]
#tm_list=[1,]
for m in [1,2,3,4,5,6,7]:
print("tm{}_".format(m))
do_evtool_esass(evlist=el+'tm{}_partI.txt'.format(m), outfile=pr+'tm{}_partI.fits'.format(m), gti='628781744.631552 628981200.', emin=0.2, emax=10.0, region=region, rmlock=True)
do_partI = False
do_partI_scans = False
do_partII = False
do_partII_scans = False
do_survey = True
if(do_partI):
for m in tm_list:
print("tm{}_".format(m))
do_evtool_esass(evlist=el+'tm{}_partI.txt'.format(m), outfile=pr+'tm{}_partI.fits'.format(m), gti='628781744.631552 628981200.', emin=0.2, emax=10.0, region=region, rmlock=True)
### Part I ###
partI=[628781744.631552,
@ -39,32 +48,41 @@ partI=[628781744.631552,
628981202.60192,]
scan=1
for i in range(len(partI)-1):
print("partI scan={}".format(scan))
tstart=partI[i]
tstop=partI[i+1]
for m in [1,2,3,4,5,6,7]:
print("tm{}_partI_scan{}".format(m,i+1))
do_evtool_esass(evlist=el+'tm{}_partI.txt'.format(m), outfile=pr+'tm{}_scan{}.fits'.format(m,scan), gti='{} {}'.format(tstart,tstop), emin=0.2, emax=10.0, region=region, rmlock=True)
if(do_partI_scans):
for m in tm_list:
print("*** tm{}_partI_scan{}".format(m,scan))
do_evtool_esass(evlist=el+'tm{}_partI.txt'.format(m), outfile=pr+'tm{}_scan{}.fits'.format(m,scan), gti='{} {}'.format(tstart,tstop), emin=0.2, emax=10.0, region=region, rmlock=True)
scan=scan+1
### Part II ###
for m in [1,2,3,4,5,6,7]:
print("tm{}_".format(m))
do_evtool_esass(evlist=el+'tm{}_partII.txt'.format(m), outfile=pr+'tm{}_partII.fits'.format(m), gti='645681013.670848 645832002.60928', emin=0.2, emax=10.0, region=region, rmlock=True)
if(do_partII):
for m in tm_list:
print("tm{}_".format(m))
do_evtool_esass(evlist=el+'tm{}_partII.txt'.format(m), outfile=pr+'tm{}_partII.fits'.format(m), gti='645681013.670848 645832002.60928', emin=0.2, emax=10.0, region=region, rmlock=True)
partII=[645681013.670848,645746786.665408,645832002.60928]
partII=[645681013.670848,
645746786.665408,
645832002.60928]
for i in range(len(partII)-1):
print("partII scan={}".format(scan))
tstart=partII[i]
tstop=partII[i+1]
for m in [1,2,3,4,5,6,7]:
print("tm{}_partII_scan{}".format(m,i+1))
do_evtool_esass(evlist=el+'tm{}_partII.txt'.format(m), outfile=pr+'tm{}_scan{}.fits'.format(m,scan), gti='{} {}'.format(tstart,tstop), emin=0.2, emax=10.0, region=region, rmlock=True)
if(do_partII_scans):
for m in tm_list:
print("*** tm{}_partII_scan{}".format(m,scan))
do_evtool_esass(evlist=el+'tm{}_partII.txt'.format(m), outfile=pr+'tm{}_scan{}.fits'.format(m,scan), gti='{} {}'.format(tstart,tstop), emin=0.2, emax=10.0, region=region, rmlock=True)
scan=scan+1
"""
### Survey ###
for m in [1,2,3,4,5,6,7]:
do_evtool_esass(evlist=el+'tm{}_survey.txt'.format(m), outfile=pr+'tm{}_survey.fits'.format(m), emin=0.2, emax=10.0, rmlock=True)
if(do_survey):
for m in tm_list:
do_evtool_esass(evlist=el+'tm{}_survey.txt'.format(m), outfile=pr+'tm{}_survey.fits'.format(m), emin=0.2, emax=10.0, rmlock=True, rusky=True)
"""

View File

@ -17,7 +17,6 @@ import coma
from coma.utils import *
from coma.config import *
outkey="mosa_tm0"
""" find root dir """
@ -36,28 +35,33 @@ outfile_dir=root_path+'/products'
create_folder(outfile_dir)
do_init = False
do_merge = True
do_rate = False
do_init = True
do_merge = False
do_rate = False
do_adapt = False # requires CIAO
index=5 # select energy band
width=7000
index=0 # select energy band
width=10000
vign=True
vignetting = 'vign' if (vign==True) else 'novign'
#print(ra_cen, de_cen)
#sys.exit()
outkey="mosa_all_tm0"
events=[]
expmaps=[]
bkgmaps=[]
for tmkey in keylist_tm.keys():
keylist=keylist_survey
for tmkey in keylist.keys():
#if not (tmkey == '2'):
# continue
print("TM{} in work... init events".format(tmkey))
for datakey in keylist_tm[tmkey]:
if ("survey" in datakey):
continue
for datakey in keylist[tmkey]:
#if not ("survey" in datakey):
# continue
print("--> {}".format(datakey))
""" Подготавливаем списки событий индивидуальных наблюдений """
outfile_evtool,outfile_expmap=init_events(key=datakey, eband_index=eband[index],
@ -67,21 +71,21 @@ for tmkey in keylist_tm.keys():
do_obsmode=True, # also controlled by lock file!
do_center=True, # also controlled by lock file!
do_evtool=True,
do_expmap=True,
do_expmap=False,
vign=vign,
ra_cen=ra_cen, de_cen=de_cen, width=width,
emin_kev=emin_kev[index],
emax_kev=emax_kev[index])
events.append(outfile_evtool)
expmaps.append(outfile_expmap)
bkgmaps.append("{}_BackMap3_en{}.fits".format(os.path.join(outfile_dir,datakey), eband[0]))
bkgmaps.append("{}_BackMap3_en{:02d}.fits".format(os.path.join(outfile_dir,datakey), eband[0]))
""" Собираем общий список событий """
outfile_evtool="{}_EventList_en{}.fits".format(os.path.join(outfile_dir,outkey), eband[index])
outfile_evtool="{}_EventList_en{:02d}.fits".format(os.path.join(outfile_dir,outkey), eband[index])
outfile_expmap="{}_ExposureMap_en{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
outfile_bkgmap="{}_BackMap_en{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
outfile_expmap="{}_ExposureMap_en{:02d}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
outfile_bkgmap="{}_BackMap_en{:02d}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
if(do_merge==True):
evlist="{}.evlist.txt".format(os.getpid())
@ -91,22 +95,23 @@ if(do_merge==True):
f.close()
do_evtool_esass(evlist=evlist, outfile=outfile_evtool, width=width)
# old version, does not work with long list of files
#do_evtool_esass(events=events, outfile=outfile_evtool)
#do_fimgmerge_ftools(maps=expmaps, outfile=outfile_expmap)
do_fimgmerge_ftools(maps=expmaps, outfile=outfile_expmap)
#do_fimgmerge_ftools(maps=bkgmaps, outfile=outfile_bkgmap)
if(os.path.isfile(evlist)==True):
os.remove(evlist)
outfile_rate="{}_RateMap_en{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
outfile_rate="{}_RateMap_en{:02d}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
if(do_rate==True):
make_rate_map(cntmap=outfile_evtool, expmap=outfile_expmap, outfile=outfile_rate)
function='gaussian'
outfile_adapt="{}_ImageAdapt_en{}.{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], function, vignetting)
outfile_adapt="{}_ImageAdapt_en{:02d}.{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], function, vignetting)
if(do_adapt==True):
do_adapt_ciao(infile=outfile_evtool, outfile=outfile_adapt, expmap=outfile_expmap, function=function, expcut=100)

View File

@ -47,48 +47,50 @@ from multiprocessing import Pool
from os.path import dirname
import inspect
import uds
import coma
from uds.utils import *
from uds.config import *
from coma.utils import *
from coma.config import *
""" find UDS root dir """
root_path=dirname(dirname(dirname(inspect.getfile(uds))))
print("UDS root path: {}".format(root_path))
""" find Coma root dir """
root_path=dirname(dirname(dirname(inspect.getfile(coma))))
print("Coma root path: {}".format(root_path))
infile_dir=root_path+'/data/processed'
outfile_dir=root_path+'/products'
create_folder(outfile_dir)
run_Pool=False
keylist=keylist_survey
do_init = False
do_ermask = False
do_init = True
do_ermask = True
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_erbox1 = True # local mode
do_erbackmap1 = True #
do_erbox2 = True # map mode, with background map
do_erbackmap2 = True #
do_erbox3 = True # map mode, with background map
do_erbackmap3 = True #
do_ermldet = False
do_catprep = False
do_cross_match = False
do_ermldet = True
do_catprep = True
do_cross_match = True
do_astro_corr = False # search optimal shift
do_astro_update = True
do_astro_corr = True # search optimal shift
do_astro_update = True
do_wcs_match = False # Chandra task -- DEPRECATED
do_wcs_update = False # Chandra task -- DEPRECATED
eband_selected=[5]
eband_selected=[0]
vign=True
vign=False
vignetting = 'vign' if (vign==True) else 'novign'
width=10000
def runme(datakey):
""" runs datakey over energy bands """
events=[]
@ -120,7 +122,7 @@ def runme(datakey):
do_evtool=True,
do_expmap=True,
vign=vign,
ra_cen=ra_cen, de_cen=de_cen,
ra_cen=ra_cen, de_cen=de_cen, width=width,
emin_kev=emin_kev[index],
emax_kev=emax_kev[index])
expmaps.append(outfile_expmap)
@ -151,7 +153,7 @@ def runme(datakey):
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))
outfile_boxlist1.append("{}_BoxList1_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
if(do_erbox1==True):
""" erbox in local mode """
cmd=["erbox",
@ -178,8 +180,8 @@ def runme(datakey):
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)
outfile_backmap1.append("{}_BackMap1_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
cheese_mask="{}_CheeseMask1_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
if(do_erbackmap1==True):
""" back map 1 """
cmd=["erbackmap",
@ -205,7 +207,7 @@ def runme(datakey):
os.system((" ").join(cmd))
print((" ").join(cmd))
outfile_boxlist2.append("{}_BoxList2_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
outfile_boxlist2.append("{}_BoxList2_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
if(do_erbox2==True):
""" erbox in background mode """
cmd=["erbox",
@ -233,8 +235,8 @@ def runme(datakey):
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)
outfile_backmap2.append("{}_BackMap2_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
cheese_mask="{}_CheeseMask2_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
if(do_erbackmap2==True):
""" back map 2 """
cmd=["erbackmap",
@ -260,7 +262,7 @@ def runme(datakey):
os.system((" ").join(cmd))
print((" ").join(cmd))
outfile_boxlist3.append("{}_BoxList3_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
outfile_boxlist3.append("{}_BoxList3_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
if(do_erbox3==True):
""" erbox in map mode FINAL """
cmd=["erbox",
@ -287,8 +289,8 @@ def runme(datakey):
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)
outfile_backmap3.append("{}_BackMap3_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post))
cheese_mask="{}_CheeseMask3_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
if(do_erbackmap3==True):
""" back map 3 FINAL """
cmd=["erbackmap",
@ -314,8 +316,8 @@ def runme(datakey):
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)
mllist="{}_MaxLikSourceList_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
srcmap="{}_SourceMap_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
cmd=["ermldet",
"mllist=%s" %(mllist),
@ -359,8 +361,8 @@ def runme(datakey):
print((" ").join(cmd))
save_ermldet_ds9reg(mllist,scale=60*60)
catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
catprep_en0="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,datakey), eband[0], outfile_post)
catprep="{}_SourceCatalog_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post)
catprep_en0="{}_SourceCatalog_en{:02d}{}".format(os.path.join(outfile_dir,datakey), eband[0], outfile_post)
if(do_catprep==True):
cmd=["catprep",
"infile={}".format(mllist),
@ -371,21 +373,23 @@ def runme(datakey):
save_catprep_ds9reg(catprep,scale=60*60)
if(do_cross_match==True):
crossmatch_shu2019(catprep,dlmin=10,refimage=events[ii],crval=wcslist[datakey],
catalog=root_path+"/data/Gaia_unWISE/Gaia_unWISE_UDS.fits.catalog",errlim=5.0)
crossmatch_shu2019(catprep, dlmin=10,crval=[ra_cen, de_cen], refimage=events[ii],datakey=datakey,
catalog=root_path+"/data/Gaia_unWISE/Gaia_unWISE_Coma.fits.catalog",errlim=5.0)
if(do_astro_corr==True and eband[index]=='0'):
if(do_astro_corr==True and eband[index]==0):
""" run astro_corr for 0.3-2.3 keV only """
print("START")
wcs_astro_corr(catprep)
#wcs_match_ciao(catprep, method='rst',radius=12,residlim=0,residtype=0,residfac=1)
if(do_astro_update==True):
""" run astro_corr for 0.3-2.3 keV only """
attcorr=wcs_update_shift(events[ii],flog=catprep_en0.replace(".fits", ".shift.log"))
do_evtool_esass(evfile=attcorr,outfile=attcorr,rmlock=False, do_center=True, ra_cen=ra_cen, de_cen=de_cen)
do_evtool_esass(evfile=attcorr,outfile=attcorr,rmlock=False, do_center=True, ra_cen=ra_cen, de_cen=de_cen, width=width)
if(do_wcs_match==True and eband[index]=='0'):
if(do_wcs_match==True and eband[index] == 0):
""" run wcs_match for 0.3-2.3 keV only """
wcs_match_ciao(catprep, method='trans',radius=12,residlim=5)
#wcs_match_ciao(catprep, method='rst',radius=12,residlim=0,residtype=0,residfac=1)
@ -393,7 +397,7 @@ def runme(datakey):
if(do_wcs_update==True):
""" use 0.3-2.3 keV transform matrix for all other bands """
attcorr=wcs_update_ciao(events[ii],crval=wcslist[datakey],transformfile=catprep_en0.replace(".fits", ".xfm"),clean=False)
do_evtool_esass(evfile=attcorr,outfile=attcorr,rmlock=False, do_center=True, ra_cen=ra_cen, de_cen=de_cen)
do_evtool_esass(evfile=attcorr,outfile=attcorr,rmlock=False, do_center=True, ra_cen=ra_cen, de_cen=de_cen, width=width)
"""
@ -403,13 +407,16 @@ runme("tm5_obs_1")
runme("tm6_scan_1")
"""
#runme("tm1_scan1")
if(run_Pool==True):
# parallel run
items=[]
for tmkey in keylist_tm.keys():
for datakey in keylist_tm[tmkey]:
for tmkey in keylist.keys():
for datakey in keylist[tmkey]:
items.append(datakey)
with Pool() as pool:
@ -417,8 +424,9 @@ if(run_Pool==True):
else:
# conventional run
for tmkey in keylist_tm.keys():
for datakey in keylist_tm[tmkey]:
for tmkey in keylist.keys():
for datakey in keylist[tmkey]:
print("--> {}".format(datakey))
runme(datakey)
#sys.exit()

View File

@ -51,50 +51,51 @@ import pickle
import uds
import coma
from uds.utils import *
from uds.config import *
from coma.utils import *
from coma.config import *
""" find UDS root dir """
#root_path=dirname(dirname(dirname(inspect.getfile(uds))))
""" find Coma root dir """
root_path=dirname(dirname(dirname(inspect.getfile(coma))))
"""
ftools does not like long file path names,
for this reason, we use relative path here
"""
root_path='..'
print("UDS root path: {}".format(root_path))
#root_path='..'
print("Coma root path: {}".format(root_path))
infile_dir=root_path+'/data/processed'
outfile_dir=root_path+'/products'
create_folder(outfile_dir)
local_run = False
local_run = True
outkey="tm0"
do_init = True
do_merge = True
do_detmask = True
do_expmap = True
do_erbox1 = True # local mode
do_erbackmap1 = True #
do_erbox2 = True # map mode, with background map
do_erbackmap2 = True #
do_erbox3 = True # map mode, with background map
do_erbackmap3 = True #
do_detmask = False
do_expmap = 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_ersensmap = False
do_ermldet = True
do_ermldet = False
do_fixcat = False # only for index=0
do_fixxmm = False # prepare forced photometry, only for index=0
do_apetool = False
do_catprep = True
do_catprep = False
do_filter_catalog = False
do_cross_match = False
index=3
index=0
forced=False
""" If forced=True, take input catalog from energy range en0 """
@ -103,29 +104,39 @@ comm='' # for 4XMM-DR12 forced photometry use '-xmm'
vign=True
attcorr=False
keylist = keylist_tm
outkey = "mosa_all_tm0{}".format('_attcorr' if (attcorr==True) else '')
"""
mosa_scans_tm0 -- all scans individually (keylist_scans)
mosa_partI_tm0 -- partI only (as is) (keylist_partI)
mosa_partII_tm0 -- partII only (as is) (keylist_partII)
mosa_parts_tm0 -- all parts (partI and partII) (keylist_parts)
mosa_all_tm0 -- all scans (scan1-7) and survey data (keylist_tm)
"""
vignetting = 'vign' if (vign==True) else 'novign'
attcorr=True
events=[]
expmaps=[]
bkgmaps=[]
for tmkey in keylist_tm.keys():
for tmkey in keylist.keys():
print("TM{} in work... init events".format(tmkey))
for datakey in keylist_tm[tmkey]:
print("--> {}".format(datakey))
for datakey in keylist[tmkey]:
#print("--> {}".format(datakey))
""" Подготавливаем списки событий индивидуальных наблюдений """
outfile_evtool,outfile_expmap=init_events(key=datakey,attcorr=attcorr,
eband_index=eband[index],
infile_dir=infile_dir,
outfile_dir=outfile_dir,
do_init=do_init,
do_obsmode=False,
do_center=False,
do_obsmode=True,
do_center=True,
do_evtool=True,
do_expmap=False,
vign=vign,
ra_cen=ra_cen, de_cen=de_cen,
ra_cen=ra_cen, de_cen=de_cen, width=width,
emin_kev=emin_kev[index],
emax_kev=emax_kev[index])
events.append(outfile_evtool)
@ -133,15 +144,25 @@ for tmkey in keylist_tm.keys():
""" Собираем общий список событий """
outfile_evtool="{}_EventList_en{}.fits".format(os.path.join(outfile_dir,outkey),
outfile_evtool="{}_EventList_en{:02d}.fits".format(os.path.join(outfile_dir,outkey),
eband[index])
if(do_merge==True):
do_evtool_esass(events=events, outfile=outfile_evtool)
#do_evtool_esass(events=events, outfile=outfile_evtool)
evlist="{}.evlist.txt".format(os.getpid())
f = open(evlist, "w")
for s in events:
f.write("{}\n".format(s))
f.close()
print(outfile_evtool)
do_check_events(events=events)
do_evtool_esass(evlist=evlist, outfile=outfile_evtool, width=width)
#if(os.path.isfile(evlist)==True):
# os.remove(evlist)
""" makes detmask from TM exposures """
detmask="{}/{}_DetectorMask_en{}{}".format(outfile_dir,
detmask="{}/{}_DetectorMask_en{:02d}{}".format(outfile_dir,
outkey,
eband[index],
outfile_post)
@ -151,16 +172,16 @@ if(do_detmask==True):
"""
Собираем общую карту экспозиции, обратите внимание на коэффициент 7.
Экспозиция рассчитывается на 7 телескопов.
outfile_expmap="{}_ExposureMap_en{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
outfile_bkgmap="{}_BackMap_en{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
outfile_expmap="{}_ExposureMap_en{:02d}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
outfile_bkgmap="{}_BackMap_en{:02d}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting)
"""
outfile_expmap="{}_ExposureMap_en{}.{}{}".format(os.path.join(outfile_dir,outkey),
outfile_expmap="{}_ExposureMap_en{:02d}.{}{}".format(os.path.join(outfile_dir,outkey),
eband[index],vignetting,
outfile_post)
if(do_expmap==True):
create_expmap_merged(expmaps,outfile_expmap,scale=7.0)
outfile_boxlist1="{}/{}_BoxList1_en{}{}".format(outfile_dir,outkey, eband[index], outfile_post)
outfile_boxlist1="{}/{}_BoxList1_en{:02d}{}".format(outfile_dir,outkey, eband[index], outfile_post)
if(do_erbox1==True):
cmd=["erbox",
"images=\'{}\'".format(outfile_evtool),
@ -184,13 +205,13 @@ if(do_erbox1==True):
save_ds9reg(outfile_boxlist1)
""" Background map 1 """
outfile_backmap1="{}_BackMap1_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
cheese_mask="{}_CheeseMask1_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
outfile_backmap1="{}_BackMap1_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
cheese_mask="{}_CheeseMask1_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
if(do_erbackmap1==True):
do_erbackmap_esass(outfile_evtool,outfile_expmap,outfile_boxlist1,detmask,emin_ev[index],emax_ev[index],
outfile_backmap1,cheese_mask)
outfile_boxlist2="{}/{}_BoxList2_en{}{}".format(outfile_dir,outkey, eband[index], outfile_post)
outfile_boxlist2="{}/{}_BoxList2_en{:02d}{}".format(outfile_dir,outkey, eband[index], outfile_post)
if(do_erbox2==True):
cmd=["erbox",
"images=\'{}\'".format(outfile_evtool),
@ -215,13 +236,13 @@ if(do_erbox2==True):
save_ds9reg(outfile_boxlist2)
""" Background map 2 """
outfile_backmap2="{}_BackMap2_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
cheese_mask="{}_CheeseMask2_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
outfile_backmap2="{}_BackMap2_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
cheese_mask="{}_CheeseMask2_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
if(do_erbackmap2==True):
do_erbackmap_esass(outfile_evtool,outfile_expmap,outfile_boxlist2,detmask,emin_ev[index],emax_ev[index],
outfile_backmap2,cheese_mask)
outfile_boxlist3="{}/{}_BoxList3_en{}{}".format(outfile_dir,outkey, eband[index], outfile_post)
outfile_boxlist3="{}/{}_BoxList3_en{:02d}{}".format(outfile_dir,outkey, eband[index], outfile_post)
if(do_erbox3==True):
cmd=["erbox",
"images=\'{}\'".format(outfile_evtool),
@ -246,21 +267,21 @@ if(do_erbox3==True):
save_ds9reg(outfile_boxlist3)
""" Background map 3 """
outfile_backmap3="{}_BackMap3_en{}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], vignetting, outfile_post)
cheese_mask="{}_CheeseMask3_en{}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], vignetting, outfile_post)
outfile_backmap3="{}_BackMap3_en{:02d}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], vignetting, outfile_post)
cheese_mask="{}_CheeseMask3_en{:02d}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], vignetting, outfile_post)
if(do_erbackmap3==True):
boxlist3 = outfile_boxlist3 if(forced == False) else "{}/{}_BoxList3_en{}{}".format(outfile_dir,outkey, eband[0], outfile_post)
boxlist3 = outfile_boxlist3 if(forced == False) else "{}/{}_BoxList3_en{:02d}{}".format(outfile_dir,outkey, eband[0], outfile_post)
do_erbackmap_esass(outfile_evtool,outfile_expmap,boxlist3,detmask,emin_ev[index],emax_ev[index],
outfile_backmap3,cheese_mask)
if(forced==True):
mllist="{}_MaxLikSourceList_en{}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)
srcmap="{}_SourceMap_en{}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)
mllist="{}_MaxLikSourceList_en{:02d}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)
srcmap="{}_SourceMap_en{:02d}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)
""" for en1,2,3,6 give mllist from en0 as input """
#boxlist3="{}_MaxLikSourceList_en{}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[0], comm, outfile_post)
#boxlist3="{}_MaxLikSourceList_en{:02d}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[0], comm, outfile_post)
#if(index==0):
boxlist3="{}_MaxLikSourceList_en{}.fixed{}{}".format(os.path.join(outfile_dir,outkey), eband[0], comm, outfile_post)
boxlist3="{}_MaxLikSourceList_en{:02d}.fixed{}{}".format(os.path.join(outfile_dir,outkey), eband[0], comm, outfile_post)
if not (os.path.exists(boxlist3)):
print("{} not found. Run do_fixcat=True, index=0, forced=False".format(boxlist3))
sys.exit()
@ -274,11 +295,11 @@ if(forced==True):
cutrad="cutrad=15."
if(index == 3 or index == 6):
""" for hard band take unvignetted background """
outfile_backmap3="{}_BackMap3_en{}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], "novign", outfile_post)
outfile_backmap3="{}_BackMap3_en{:02d}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], "novign", outfile_post)
else:
mllist="{}_MaxLikSourceList_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
srcmap="{}_SourceMap_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
mllist="{}_MaxLikSourceList_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
srcmap="{}_SourceMap_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
boxlist3 = outfile_boxlist3
fitpos_flag="fitpos_flag=yes"
fitext_flag="fitext_flag=yes"
@ -325,18 +346,14 @@ cmd=["ermldet",
if(do_ersensmap==True):
detlike=10
create_sensmap(sensmap="{}_SensitivityMap_dl{}_en{}{}".format(os.path.join(outfile_dir,outkey), detlike,
eband[index], outfile_post),
areatab="{}_AreaTable_dl{}_en{}{}".format(os.path.join(outfile_dir,outkey),
detlike, eband[index], outfile_post),
create_sensmap(sensmap="{}_SensitivityMap_dl{}_en{:02d}{}".format(os.path.join(outfile_dir,outkey), detlike, eband[index], outfile_post),
areatab="{}_AreaTable_dl{}_en{:02d}{}".format(os.path.join(outfile_dir,outkey), detlike, eband[index], outfile_post),
expmap=outfile_expmap, backmap=outfile_backmap3,detlike=detlike,
detmask=detmask, emin=emin_ev[index], emax=emax_ev[index],ecf=ecf[index])
"""
detlike=6
create_sensmap(sensmap="{}_SensitivityMap_dl{}_en{}{}".format(os.path.join(outfile_dir,outkey), detlike,
eband[index], outfile_post),
areatab="{}_AreaTable_dl{}_en{}{}".format(os.path.join(outfile_dir,outkey),
detlike, eband[index], outfile_post),
create_sensmap(sensmap="{}_SensitivityMap_dl{}_en{:02d}{}".format(os.path.join(outfile_dir,outkey), detlike, eband[index], outfile_post),
areatab="{}_AreaTable_dl{}_en{:02d}{}".format(os.path.join(outfile_dir,outkey), detlike, eband[index], outfile_post),
expmap=outfile_expmap, backmap=outfile_backmap3,detlike=detlike,
detmask=detmask, emin=emin_ev[index], emax=emax_ev[index],ecf=ecf[index])
"""
@ -435,9 +452,9 @@ if(do_apetool==True):
runme(cmd, local_run=local_run)
if(forced==True):
catprep="{}_SourceCatalog_en{}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)
catprep="{}_SourceCatalog_en{:02d}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)
else:
catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
catprep="{}_SourceCatalog_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
if(do_catprep==True):
cmd=["catprep",