From b94818554fe640c9d77d3342e9fd67e546547fcb Mon Sep 17 00:00:00 2001 From: Roman Krivonos Date: Wed, 8 Mar 2023 19:58:59 +0300 Subject: [PATCH] 8 March work --- .gitignore | 1 + data/evtlists/tm7.txt | 6 +- scripts/00_check.py | 58 +++++++++++++++++++ scripts/02_merge_events.py | 21 +++++-- scripts/03_init_obs.py | 14 ++--- scripts/README.md | 4 +- uds/uds/config.py | 12 ++-- uds/uds/utils.py | 111 +++++++++++++++++++++++++++++++++---- 8 files changed, 195 insertions(+), 32 deletions(-) create mode 100755 scripts/00_check.py diff --git a/.gitignore b/.gitignore index 041ca7e..924b194 100644 --- a/.gitignore +++ b/.gitignore @@ -53,3 +53,4 @@ venv.bak/ *.cross *.ref *.src +*.dat diff --git a/data/evtlists/tm7.txt b/data/evtlists/tm7.txt index 5b5d009..7d1b7fd 100644 --- a/data/evtlists/tm7.txt +++ b/data/evtlists/tm7.txt @@ -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 diff --git a/scripts/00_check.py b/scripts/00_check.py new file mode 100755 index 0000000..6302333 --- /dev/null +++ b/scripts/00_check.py @@ -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)) + + diff --git a/scripts/02_merge_events.py b/scripts/02_merge_events.py index 0edc318..9d8167a 100755 --- a/scripts/02_merge_events.py +++ b/scripts/02_merge_events.py @@ -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) + diff --git a/scripts/03_init_obs.py b/scripts/03_init_obs.py index 380f702..eab5f46 100755 --- a/scripts/03_init_obs.py +++ b/scripts/03_init_obs.py @@ -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) -""" + diff --git a/scripts/README.md b/scripts/README.md index 6a6c484..134b867 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -1,7 +1,7 @@ Подготовка рабочего окружения: ``` -source venv/bin/activate.csh -esass +source /venv/bin/activate.csh +source /eSASS4EDR/bin/esass-init.csh ``` diff --git a/uds/uds/config.py b/uds/uds/config.py index 9f34b07..2b2b39e 100644 --- a/uds/uds/config.py +++ b/uds/uds/config.py @@ -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' diff --git a/uds/uds/utils.py b/uds/uds/utils.py index 3eb6052..a90c667 100644 --- a/uds/uds/utils.py +++ b/uds/uds/utils.py @@ -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 +