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

evening job

Browse Source
This commit is contained in:
Roman Krivonos 2023-03-09 21:10:39 +03:00
parent b94818554f
commit e7f6f9d257
7 changed files with 178 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
scripts/00_check.py
View File

@ -34,7 +34,7 @@ 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))
#print("--> {}".format(datakey))
""" Запускаем полностью в холостом режиме, нам нужно получить только названия файлов """
outfile_evtool,outfile_expmap=init_events(key=datakey, eband_index=eband[index],
infile_dir=infile_dir,
@ -51,7 +51,10 @@ for tmkey in keylist_tm.keys():
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("{}, {} -- {}, {} -- {}, {:.2f} ks".format(outfile_evtool,tstart,tstop,
mission2date_utc(tstart),
mission2date_utc(tstop),
(tstop-tstart)/1000))
print("\n***\n*** Total exposure: {:.1f} ks\n***".format(totexp))

2
scripts/01_init_events.py
View File

@ -52,7 +52,7 @@ do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_obs_3.fits', gti='62067699
do_badpix_tm6(filename=pr+'tm6_obs_3.fits')
"""
# TM7
do_evtool_esass(evlist=el+'tm7.txt', outfile=pr+'tm7_obs_1.fits', gti='621043136. 621052416.', emin=0.2, emax=10.0, region=region)
do_evtool_esass(evlist=el+'tm7.txt', outfile=pr+'tm7_obs_1.fits', gti='621043136. 621052416.', emin=0.2, emax=10.0, region=region, rmlock=True)
#do_evtool_esass(evlist=el+'tm7.txt', outfile=pr+'tm7_obs_2.fits', gti='621110272. 621117952.', region=region)

38
scripts/02_merge_events.py
View File

@ -20,18 +20,25 @@ from uds.config import *
outkey="mosa_tm0"
""" find UDS root dir """
root_path=dirname(dirname(dirname(inspect.getfile(uds))))
#root_path=dirname(dirname(dirname(inspect.getfile(uds))))
root_path='..'
print("UDS root path: {}".format(root_path))
infile_dir=root_path+'/data/processed'
outfile_dir=root_path+'/products'
create_folder(outfile_dir)
index=5
index=5 # select energy band
do_init = True
do_merge = True
do_adapt = False # requires CIAO
#do_rate = True
vign=True
events=[]
expmaps=[]
@ -44,10 +51,11 @@ for tmkey in keylist_tm.keys():
infile_dir=infile_dir,
outfile_dir=outfile_dir,
do_init=do_init,
do_obsmode=True,
do_center=True,
do_evtool=True,
do_obsmode=False,
do_center=False,
do_evtool=False,
do_expmap=True,
vign=vign,
ra_cen=ra_cen, de_cen=de_cen,
emin_kev=emin_kev[index],
emax_kev=emax_kev[index])
@ -56,11 +64,23 @@ for tmkey in keylist_tm.keys():
""" Собираем общий список событий """
outfile_evtool="{}_EventList_en{}.fits".format(os.path.join(outfile_dir,outkey), eband[index])
outfile_expmap="{}_ExposureMap_en{}.fits".format(os.path.join(outfile_dir,outkey), eband[index])
if (vign==True):
outfile_expmap="{}_ExposureMap_en{}.vign.fits".format(os.path.join(outfile_dir,outkey), eband[index])
else:
outfile_expmap="{}_ExposureMap_en{}.novign.fits".format(os.path.join(outfile_dir,outkey), eband[index])
if(do_merge==True):
do_evtool_esass(events=events, outfile=outfile_evtool)
do_expmap_ftools(expmaps=expmaps, outfile=outfile_expmap)
outfile_adapt="{}_ImageAdapt_en{}.fits".format(os.path.join(outfile_dir,outkey), eband[index])
function='gaussian'
outfile_adapt="{}_ImageAdapt_en{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], function)
if(do_adapt==True):
do_adapt_ciao(infile=outfile_evtool, outfile=outfile_adapt)
do_adapt_ciao(infile=outfile_evtool, outfile=outfile_adapt, expmap=outfile_expmap, function=function, expcut=100)
"""
outfile_rate="{}_RateAdapt_en{}.fits".format(os.path.join(outfile_dir,outkey), eband[index])
if(do_rate==True):
make_rate_map(cntmap=outfile_adapt, expmap=outfile_expmap, outfile=outfile_rate)
"""

8
scripts/03_init_obs.py
View File

@ -318,7 +318,7 @@ def runme(datakey):
#hdulist = fits.open(detmask)
#w = WCS(hdulist[0].header)
#w.wcs.crval=wcslist[datakey]
crossmatch_shu2019(catprep,dlmin=10,refimage=detmask,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")
"""
@ -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)
"""

8
scripts/README.md
View File

@ -16,6 +16,13 @@ source <MY PATH>/eSASS4EDR/bin/esass-init.csh
Попутно этот скрипт унифицирует оригинальные списки событий для последующей обработки. А именно, корректируются слова OBS_MODE=POINING/SURVEY в зависимости от типа наблюдения и производится центрирование на одни и те же координаты с помощью команды ```radec2xy```.
Для запуска адаптивного сглаживания ```do_adapt = True``` требуется запустить окружение ```ciao```, так как нужна команда ```dmimgadapt```
```
conda activate ciao-4.15
source <MY PATH>/eSASS4EDR/bin/esass-init.csh
```
### 03_init_obs.py
@ -23,3 +30,4 @@ source <MY PATH>/eSASS4EDR/bin/esass-init.csh
2) Запускает ```erbox``` в три этапа, чтобы получить рабочий список источников для ```ermldet```.
3) Запускает ```ermldet```
4) Делает кросс-корреляцию с каталогом Gaia-unWISE ```do_cross_match=True```, которая создает три файла: ```.cross``` -- все пересечения, и ```.ref``` / ```.src``` -- входные каталоги для последующей команды ```wcs_match```

15
scripts/expmaps.txt Normal file
View File

@ -0,0 +1,15 @@
../products/tm5_obs_1_ExposureMap_en5.fits,0,0
../products/tm5_obs_2_ExposureMap_en5.fits,0,0
../products/tm6_obs_1_ExposureMap_en5.fits,0,0
../products/tm6_obs_2_badpix_ExposureMap_en5.fits,0,0
../products/tm6_obs_3_badpix_ExposureMap_en5.fits,0,0
../products/tm6_park_1_ExposureMap_en5.fits,0,0
../products/tm6_park_2_ExposureMap_en5.fits,0,0
../products/tm6_park_3_ExposureMap_en5.fits,0,0
../products/tm6_park_4_ExposureMap_en5.fits,0,0
../products/tm6_scan_1_ExposureMap_en5.fits,0,0
../products/tm6_scan_2_ExposureMap_en5.fits,0,0
../products/tm6_scan_3_ExposureMap_en5.fits,0,0
../products/tm6_scan_4_ExposureMap_en5.fits,0,0
../products/tm7_obs_1_ExposureMap_en5.fits,0,0
../products/tm7_obs_2_ExposureMap_en5.fits,0,0

146
uds/uds/utils.py
View File

@ -17,11 +17,27 @@ from astropy.coordinates import SkyCoord # High-level coordinates
from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames
from astropy.coordinates import Angle, Latitude, Longitude # Angles
from astropy.time import Time, TimeDelta, TimezoneInfo, TimeFromEpoch
import statistics
import shutil
MJDREF = 51543.875
TZ_UTC = TimezoneInfo(utc_offset=0*u.hour)
TZ_MSK = TimezoneInfo(utc_offset=3*u.hour)
TIMESTR_FORMAT = '%Y.%m.%d %H:%M:%S'
def mission2date_utc(timesec):
mjdref = Time(MJDREF, format='mjd')
delta = TimeDelta(timesec, format='sec')
dtime = mjdref + delta
return dtime.to_datetime(timezone=TZ_UTC)
def mission2date_msk(timesec):
mjdref = Time(MJDREF, format='mjd')
delta = TimeDelta(timesec, format='sec')
dtime = mjdref + delta + 3*u.hour
return dtime.to_datetime(timezone=TZ_MSK)
def create_folder(folder):
if not (os.path.exists(folder)):
@ -31,7 +47,7 @@ def remove_file(filename):
if(os.path.isfile(filename)==True):
os.remove(filename)
def do_evtool_esass(events=None,outfile=None,evlist=None,gti=None,region=None,emin=None,emax=None):
def do_evtool_esass(events=None,outfile=None,evlist=None,gti=None,region=None,emin=None,emax=None, rmlock=False):
eventfiles=None
if(events):
@ -60,6 +76,14 @@ def do_evtool_esass(events=None,outfile=None,evlist=None,gti=None,region=None,em
test_exe('evtool')
os.system((" ").join(cmd))
if(rmlock==True):
lockfiles=outfile.replace(".fits", "*.lock")
for filename in glob.glob(lockfiles):
print("Removing {}".format(filename))
if(os.path.isfile(filename)==True):
os.remove(filename)
pass
def set_bit(value, bit):
return value | (1<<bit)
@ -122,7 +146,7 @@ def do_badpix_tm6(filename):
f.close()
def init_events(key=None, eband_selected=[0], eband_index=None,
ra_cen=None, de_cen=None,do_init=True,
ra_cen=None, de_cen=None,do_init=True,vign=True,
emin_kev=None, emax_kev=None, infile_dir=None, outfile_dir=None,
do_obsmode=False,do_center=False,do_evtool=False,do_expmap=False):
expmaps=[]
@ -216,8 +240,13 @@ def init_events(key=None, eband_selected=[0], eband_index=None,
h.header['OBS_MODE'] = changeto
hdu.writeto(outfile_evtool,overwrite=True)
if(vign==True):
withvignetting='yes'
outfile_expmap="{}_ExposureMap_en{}.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{}{}".format(os.path.join(outfile_dir,key), eband_index, outfile_post)
if(do_expmap==True and do_init==True):
cmd=["expmap",
"inputdatasets=%s" %(outfile_evtool),
@ -225,7 +254,7 @@ def init_events(key=None, eband_selected=[0], eband_index=None,
"emax=%s" %(emax_kev),
"templateimage=%s" %(outfile_evtool),
"singlemaps=%s" %(outfile_expmap),
"withvignetting=yes",
"withvignetting={}".format(withvignetting),
"withsinglemaps=yes","withfilebadpix=yes",
"withmergedmaps=no",
]
@ -235,21 +264,37 @@ def init_events(key=None, eband_selected=[0], eband_index=None,
return outfile_evtool,outfile_expmap
def merge_expmaps_esass(expmaps,outfile):
def do_expmap_esass(expmaps=None,outfile=None,emin_kev=None,emax_kev=None,refimage=None):
cmd=["expmap",
"inputdatasets={}".format(outfile_evtool),
"inputdatasets=\'{}\'".format((" ").join(expmaps)),
"emin=%s" %(emin_kev),
"emax=%s" %(emax_kev),
"templateimage=%s" %(outfile_evtool),
"singlemaps=%s" %(outfile_expmap),
"templateimage=%s" %(refimage),
"singlemaps=%s" %(outfile),
"withvignetting=yes",
"withsinglemaps=yes","withfilebadpix=yes",
"withmergedmaps=no",
"withsinglemaps=no","withfilebadpix=yes",
"withmergedmaps=yes",
]
print((" ").join(cmd))
test_exe('expmap')
os.system((" ").join(cmd))
def do_expmap_ftools(expmaps=None,outfile=None):
""" takes first exposure as reference and merges the rest """
tmp="expmaps.txt"
f = open(tmp, "w")
for jj in range(1,len(expmaps)):
f.write("{},0,0\n".format(expmaps[jj]))
f.close()
cmd=["fimgmerge",
"infile={}".format(expmaps[0]),
"list=@{}".format(tmp),
"outfile={}".format(outfile),
"clobber=yes",
]
print((" ").join(cmd))
os.system((" ").join(cmd))
def test_exe(program):
""" Tests if executable exists in PATH """
import os
@ -326,7 +371,6 @@ def crossmatch_shu2019(filename,refimage=None,scale=1200,crval=None,devmax=30,dl
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 reference catalog: catalog=")
@ -428,14 +472,17 @@ def crossmatch_shu2019(filename,refimage=None,scale=1200,crval=None,devmax=30,dl
print("delta [arcsec] RA, Dec, sep {:.4f} {:.4f}".format((rec['src_ra']-rec['ref_ra'])*3600,
(rec['src_dec']-rec['ref_dec'])*3600,sep))
tstart, tstop = read_tstart_tstop(infile=refimage)
date_obs = mission2date_utc(tstart)
date_end = mission2date_utc(tstop)
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)))
writer.write("[DATE-OBS {}, DATE-END {}] nref, mean [arcsec] dRA, dDec, mean/median sep (mean/median): {} & {:.4f} & {:.4f} & {:.4f}/{:.4f}\n".format(
date_obs,date_end,len(delta_ra),
statistics.mean(delta_ra)*3600,
statistics.mean(delta_dec)*3600,
statistics.mean(delta_sep),
statistics.median(delta_sep)))
src_tab.write(src_fout, format='fits', overwrite=True)
with fits.open(src_fout) as f:
@ -456,31 +503,55 @@ def crossmatch_shu2019(filename,refimage=None,scale=1200,crval=None,devmax=30,dl
#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=")
def do_adapt_ciao(infile=None,outfile=None,expmap=None,function='tophat',expcut=None):
if not (infile and expmap and outfile):
print("ERROR: Please provide input and output files")
return
radfile=infile.replace(".fits", ".adapt.scale.fits")
sumfile=infile.replace(".fits", ".sum.scale.fits")
test_exe('dmimgadapt')
cmd=["dmimgadapt",
"infile={}".format(infile),
"outfile={}".format(outfile),
"function=tophat",
"function={}".format(function),
"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")),
"radfile={}".format(radfile),
"sumfile={}".format(sumfile),
"clobber=yes",]
remove_file(outfile)
print((" ").join(cmd))
os.system((" ").join(cmd))
cmd=["dmimgadapt",
"infile={}".format(infile),
"innormfile={}".format(expmap),
"outfile={}".format(outfile),
"function={}".format(function),
"minrad=1",
"maxrad=30",
"numrad=30",
"radscale=log",
"counts=30",
"inradfile={}".format(radfile),
"sumfile={}".format(sumfile),
"clobber=yes",]
remove_file(outfile)
print((" ").join(cmd))
os.system((" ").join(cmd))
if (expcut):
cntmap_data, cntmap_hdr = fits.getdata(outfile, ext=0, header=True)
expmap_data, expmap_hdr = fits.getdata(expmap, ext=0, header=True)
index = np.where(expmap_data > expcut)
cntmap_data[index]=0.0
fits.writeto(outfile, cntmap_data, cntmap_hdr, overwrite=True)
"""
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
@ -492,13 +563,28 @@ def read_tstart_tstop(infile=None):
if not infile:
print("ERROR: Please provide input file: infile=")
return
print("Reading {}".format(infile))
#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))
#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
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
cntmap_data, cntmap_hdr = fits.getdata(cntmap, ext=0, header=True)
expmap_data, expmap_hdr = fits.getdata(expmap, ext=0, header=True)
rate = np.zeros(expmap_data.shape)
index = np.where(expmap_data > 0.0)
rate[index]=cntmap_data[index]/expmap_data[index]
fits.writeto(outfile, rate, overwrite=True)