erosita/coma
1
0
Fork 0
forked from erosita/uds
Code Pull Requests Activity

8 March work

Browse Source
This commit is contained in:
Roman Krivonos 2023-03-08 19:58:59 +03:00
parent fd370a8209
commit b94818554f
8 changed files with 195 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -53,3 +53,4 @@ venv.bak/
*.cross
*.ref
*.src
*.dat

6
data/evtlists/tm7.txt
View File

@ -1,6 +1,6 @@
/srg/work/krivonos/erosita/work/events/cef_43122_7_P003_00.fits.gz
/srg/work/krivonos/erosita/work/events/cef_43123_7_P003_00.fits.gz
/srg/work/krivonos/erosita/work/events/cef_43127_7_P003_00.fits.gz
#/srg/work/krivonos/erosita/work/events/cef_43122_7_P003_00.fits.gz
#/srg/work/krivonos/erosita/work/events/cef_43123_7_P003_00.fits.gz
#/srg/work/krivonos/erosita/work/events/cef_43127_7_P003_00.fits.gz
/srg/work/krivonos/erosita/work/events/cef_43129_7_P003_00.fits.gz
/srg/work/krivonos/erosita/work/events/cef_43133_7_P003_00.fits.gz
/srg/work/krivonos/erosita/work/events/cef_43134_7_P003_00.fits.gz

58
scripts/00_check.py Executable file
View File

@ -0,0 +1,58 @@
#!/usr/bin/env python
"""Печатает информацию о наблюдениях
"""
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 *
outkey="mosa_tm0"
""" 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)
index=5
events=[]
expmaps=[]
totexp=0.0
for tmkey in keylist_tm.keys():
print("TM{} in work... init events".format(tmkey))
for datakey in keylist_tm[tmkey]:
print("--> {}".format(datakey))
""" Запускаем полностью в холостом режиме, нам нужно получить только названия файлов """
outfile_evtool,outfile_expmap=init_events(key=datakey, eband_index=eband[index],
infile_dir=infile_dir,
outfile_dir=outfile_dir,
do_init=False,
do_obsmode=False,
do_center=False,
do_evtool=False,
do_expmap=False,
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)
tstart, tstop = read_tstart_tstop(infile=outfile_evtool)
totexp=totexp+(tstop-tstart)
print("{} {} {} {}".format(outfile_evtool,tstart,tstop,(tstop-tstart)/1000))
print("\n***\n*** Total exposure: {:.1f} ks\n***".format(totexp))

21
scripts/02_merge_events.py
View File

@ -27,9 +27,14 @@ infile_dir=root_path+'/data/processed'
outfile_dir=root_path+'/products'
create_folder(outfile_dir)
index=4
index=5
do_init = True
do_merge = True
do_adapt = False # requires CIAO
events=[]
expmaps=[]
for tmkey in keylist_tm.keys():
print("TM{} in work... init events".format(tmkey))
for datakey in keylist_tm[tmkey]:
@ -38,6 +43,7 @@ for tmkey in keylist_tm.keys():
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,
@ -46,8 +52,15 @@ for tmkey in keylist_tm.keys():
emin_kev=emin_kev[index],
emax_kev=emax_kev[index])
events.append(outfile_evtool)
expmaps.append(outfile_expmap)
""" Собираем общий список событий """
outfile_evtool="{}_EventList_en{}.fits".format(os.path.join(outfile_dir,outkey), eband[index])
do_evtool_esass(events=events, outfile=outfile_evtool)
outfile_expmap="{}_ExposureMap_en{}.fits".format(os.path.join(outfile_dir,outkey), eband[index])
if(do_merge==True):
do_evtool_esass(events=events, outfile=outfile_evtool)
outfile_adapt="{}_ImageAdapt_en{}.fits".format(os.path.join(outfile_dir,outkey), eband[index])
if(do_adapt==True):
do_adapt_ciao(infile=outfile_evtool, outfile=outfile_adapt)

14
scripts/03_init_obs.py
View File

@ -315,10 +315,10 @@ def runme(datakey):
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,w=w,naxis1=hdulist[0].header['NAXIS1'],naxis2=hdulist[0].header['NAXIS2'],
#hdulist = fits.open(detmask)
#w = WCS(hdulist[0].header)
#w.wcs.crval=wcslist[datakey]
crossmatch_shu2019(catprep,dlmin=10,refimage=detmask,crval=wcslist[datakey],
catalog=root_path+"/data/Gaia_unWISE/Gaia_unWISE_UDS.fits.catalog")
"""
@ -362,13 +362,13 @@ def runme(datakey):
wcs_update(outfile_evtool[ii],crval=wcslist[key],transformfile=xfm)
"""
runme("tm7_obs_1")
#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)
"""

4
scripts/README.md
View File

@ -1,7 +1,7 @@
Подготовка рабочего окружения:
```
source venv/bin/activate.csh
esass
source <MY PATH>/venv/bin/activate.csh
source <MY PATH>/eSASS4EDR/bin/esass-init.csh
```

12
uds/uds/config.py
View File

@ -44,13 +44,13 @@ wcslist={'tm1_obs_1':[34.7279760,-5.0680267],
'tm6_scan_4':[34.5405596,-4.8088748]}
""" Диапазоны энергий. """
emin_ev=[300, 300, 600, 2300, 200]
emax_ev=[2300, 600, 2300, 5000,10000]
emin_kev=[0.3, 0.3, 0.6, 2.3, 0.2]
emax_kev=[2.3, 0.6, 2.3, 5.0, 10.0]
ecf=[1.0, 1.0, 1.0, 1.0, 1.0]
emin_ev=[300, 300, 600, 2300, 200,300]
emax_ev=[2300, 600, 2300, 5000,10000,8000]
emin_kev=[0.3, 0.3, 0.6, 2.3, 0.2, 0.3]
emax_kev=[2.3, 0.6, 2.3, 5.0, 10.0, 8.0]
ecf=[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
""" Это просто индекс диапазона для выходных файлов. """
eband=["0", "1", "2", "3", "4"]
eband=["0", "1", "2", "3", "4", "5"]
outfile_post='.fits'

111
uds/uds/utils.py
View File

@ -235,6 +235,20 @@ def init_events(key=None, eband_selected=[0], eband_index=None,
return outfile_evtool,outfile_expmap
def merge_expmaps_esass(expmaps,outfile):
cmd=["expmap",
"inputdatasets={}".format(outfile_evtool),
"emin=%s" %(emin_kev),
"emax=%s" %(emax_kev),
"templateimage=%s" %(outfile_evtool),
"singlemaps=%s" %(outfile_expmap),
"withvignetting=yes",
"withsinglemaps=yes","withfilebadpix=yes",
"withmergedmaps=no",
]
print((" ").join(cmd))
test_exe('expmap')
os.system((" ").join(cmd))
def test_exe(program):
""" Tests if executable exists in PATH """
@ -293,13 +307,29 @@ 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,scale=1200,devmax=30,w=None,naxis1=1,naxis2=1,dlmin=6.0,dlmax=10000,ext_like=0.0,outkey='shu2019', catalog=None):
def crossmatch_shu2019(filename,refimage=None,scale=1200,crval=None,devmax=30,dlmin=6.0,dlmax=10000,ext_like=0.0,outkey='shu2019', catalog=None):
if(os.path.isfile(filename)==False):
print("File not found {}".format(filename))
print("Start cross-match with Gaia-unWISE")
if not crval:
print("ERROR: Please provide center coordinates array: crval=")
print("NOTES: See uds.config.wcslist")
return
if not refimage:
print("ERROR: Please provide reference image: refimage=")
return
hdulist = fits.open(refimage)
w = WCS(hdulist[0].header)
w.wcs.crval=crval
naxis1=hdulist[0].header['NAXIS1']
naxis2=hdulist[0].header['NAXIS2']
date_obs=hdulist[0].header['DATE-OBS']
if not catalog:
print("ERROR: Please provide catalog.")
print("ERROR: Please provide reference catalog: catalog=")
return
@ -324,7 +354,8 @@ def crossmatch_shu2019(filename,scale=1200,devmax=30,w=None,naxis1=1,naxis2=1,dl
"dec":rec['dec'],
"radec_err":rec['radec_err']/scale})
fout=filename.replace(".fits", ".{}.cross".format(outkey))
fout=filename.replace(".fits", ".{}.cross.fits".format(outkey))
print("Output file: {}".format(fout))
if(os.path.isfile(fout)==True):
os.remove(fout)
@ -361,9 +392,9 @@ def crossmatch_shu2019(filename,scale=1200,devmax=30,w=None,naxis1=1,naxis2=1,dl
t.write(fout, format='fits')
src_tab = Table(names=('ID','RA', 'DEC'), dtype=('i4','f8','f8',))
src_fout=fout.replace(".cross", ".{}.src".format(outkey))
src_fout=fout.replace(".cross.fits", ".src.fits")
ref_fout=fout.replace(".cross", ".{}.ref".format(outkey))
ref_fout=fout.replace(".cross.fits", ".ref.fits")
ref_tab = Table(names=('GAIA','RA', 'DEC'), dtype=('i8','f8','f8',))
hdul = fits.open(fout)
@ -374,12 +405,14 @@ def crossmatch_shu2019(filename,scale=1200,devmax=30,w=None,naxis1=1,naxis2=1,dl
ref_map=np.zeros((naxis1, naxis2))
delta_ra = []
delta_dec = []
delta_sep = []
for rec in tbdata:
#if (rec['det_like'] < 10.0):
# continue
if not (rec['det_like'] > dlmin and rec['det_like'] < dlmax):
continue
if (rec['id'] in skip.keys()):
print("Skip ID={} DET_LIKE={:.2f}".format(rec['id'],rec['det_like']))
else:
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 = np.round(w.world_to_pixel(src_crd))
@ -390,10 +423,20 @@ def crossmatch_shu2019(filename,scale=1200,devmax=30,w=None,naxis1=1,naxis2=1,dl
ref_tab.add_row((rec['gaia'],rec['ref_ra'],rec['ref_dec']))
delta_ra.append(rec['src_ra']-rec['ref_ra'])
delta_dec.append(rec['src_dec']-rec['ref_dec'])
print("delta RA {:.6f} delta DEC {:.6f}".format(rec['src_ra']-rec['ref_ra'],
rec['src_dec']-rec['ref_dec']))
sep=src_crd.separation(ref_crd).arcsec
delta_sep.append(sep)
print("delta [arcsec] RA, Dec, sep {:.4f} {:.4f}".format((rec['src_ra']-rec['ref_ra'])*3600,
(rec['src_dec']-rec['ref_dec'])*3600,sep))
stat_fout=fout.replace(".cross.fits", ".dat")
with open(stat_fout, 'w') as writer:
writer.write("[DATE-OBS {}] nref, mean [arcsec] dRA, dDec, mean/median sep (mean/median): {} & {:.4f} & {:.4f} & {:.4f}/{:.4f} \\\\ \n".format(date_obs,
len(delta_ra),
statistics.mean(delta_ra)*3600,
statistics.mean(delta_dec)*3600,
statistics.mean(delta_sep),
statistics.median(delta_sep)))
print("Mean delta RA: {}, delta Dec: {} (arcsec)".format(statistics.mean(delta_ra)*3600,statistics.mean(delta_dec)*3600))
src_tab.write(src_fout, format='fits', overwrite=True)
with fits.open(src_fout) as f:
f[0].verify('fix')
@ -411,3 +454,51 @@ def crossmatch_shu2019(filename,scale=1200,devmax=30,w=None,naxis1=1,naxis2=1,dl
f[0].header=f[0].header+w.to_header()
f[0].data=ref_map
#f[1].header=f[1].header+w.to_header()
def do_adapt_ciao(infile=None,outfile=None):
if not infile:
print("ERROR: Please provide input file: infile=")
return
if not outfile:
print("ERROR: file for output? outfile=")
return
test_exe('dmimgadapt')
cmd=["dmimgadapt",
"infile={}".format(infile),
"outfile={}".format(outfile),
"function=tophat",
"minrad=1",
"maxrad=30",
"numrad=30",
"radscale=log",
"counts=30",
"radfile=map.scale.fits".format(infile.replace(".fits", ".adapt.scale.fits")),
"sumfile={}".format(infile.replace(".fits", ".sum.scale.fits")),
"clobber=yes",]
remove_file(outfile)
print((" ").join(cmd))
os.system((" ").join(cmd))
"""
dmimgadapt artmap_ait.fits adapt.fits tophat 1 30 30 log 30 radfile=map.scale.fits sumfile=sum.fits clobber=yes
#dmimgadapt artmap_ait.fits adapt.fits func=gaussian min=1.0 max=30 num=30 rad=linear counts=30 sumfile=sum.fits normfile=norm.fits radfile=scale.fits
"""
def read_tstart_tstop(infile=None):
"""
looks like evtool does not set START-TSTOP values correctly. Take these times from event list.
"""
if not infile:
print("ERROR: Please provide input file: infile=")
return
print("Reading {}".format(infile))
hdul = fits.open(infile)
tbdata_ref = hdul[1].data
tstart=tbdata_ref['TIME'][0]
tstop=tbdata_ref['TIME'][-1]
print("*** TSTART diff {:.1f} ks".format((tstart-hdul[1].header['TSTART'])/1000))
print("*** TSTOP diff {:.1f} ks".format((hdul[1].header['TSTOP']-tstop)/1000))
return tstart, tstop