From e7f6f9d257e45d2a4ecfc676fc674b8480f43fbd Mon Sep 17 00:00:00 2001 From: Roman Krivonos Date: Thu, 9 Mar 2023 21:10:39 +0300 Subject: [PATCH] evening job --- scripts/00_check.py | 7 +- scripts/01_init_events.py | 2 +- scripts/02_merge_events.py | 38 +++++++--- scripts/03_init_obs.py | 8 +- scripts/README.md | 8 ++ scripts/expmaps.txt | 15 ++++ uds/uds/utils.py | 146 +++++++++++++++++++++++++++++-------- 7 files changed, 178 insertions(+), 46 deletions(-) create mode 100644 scripts/expmaps.txt diff --git a/scripts/00_check.py b/scripts/00_check.py index 6302333..74474d7 100755 --- a/scripts/00_check.py +++ b/scripts/00_check.py @@ -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)) diff --git a/scripts/01_init_events.py b/scripts/01_init_events.py index d2c79be..fe8c5ef 100755 --- a/scripts/01_init_events.py +++ b/scripts/01_init_events.py @@ -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) diff --git a/scripts/02_merge_events.py b/scripts/02_merge_events.py index 9d8167a..5c32cdf 100755 --- a/scripts/02_merge_events.py +++ b/scripts/02_merge_events.py @@ -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) +""" diff --git a/scripts/03_init_obs.py b/scripts/03_init_obs.py index eab5f46..2f630f8 100755 --- a/scripts/03_init_obs.py +++ b/scripts/03_init_obs.py @@ -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) - +""" diff --git a/scripts/README.md b/scripts/README.md index 134b867..2d3f14d 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -16,6 +16,13 @@ source /eSASS4EDR/bin/esass-init.csh Попутно этот скрипт унифицирует оригинальные списки событий для последующей обработки. А именно, корректируются слова OBS_MODE=POINING/SURVEY в зависимости от типа наблюдения и производится центрирование на одни и те же координаты с помощью команды ```radec2xy```. +Для запуска адаптивного сглаживания ```do_adapt = True``` требуется запустить окружение ```ciao```, так как нужна команда ```dmimgadapt``` + +``` +conda activate ciao-4.15 +source /eSASS4EDR/bin/esass-init.csh +``` + ### 03_init_obs.py @@ -23,3 +30,4 @@ source /eSASS4EDR/bin/esass-init.csh 2) Запускает ```erbox``` в три этапа, чтобы получить рабочий список источников для ```ermldet```. 3) Запускает ```ermldet``` 4) Делает кросс-корреляцию с каталогом Gaia-unWISE ```do_cross_match=True```, которая создает три файла: ```.cross``` -- все пересечения, и ```.ref``` / ```.src``` -- входные каталоги для последующей команды ```wcs_match``` + diff --git a/scripts/expmaps.txt b/scripts/expmaps.txt new file mode 100644 index 0000000..c77f2cb --- /dev/null +++ b/scripts/expmaps.txt @@ -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 diff --git a/uds/uds/utils.py b/uds/uds/utils.py index a90c667..36afd66 100644 --- a/uds/uds/utils.py +++ b/uds/uds/utils.py @@ -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< 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) + +