#!/usr/bin/env python # -*- coding: utf-8 -*- """ НАЗВАНИЕ: 04_mosaics.py НАЗНАЧЕНИЕ: Собирает мозайки в разных энергетических диапазонах. ВЫЗОВ: esass ./01_mosaics.py УПРАВЛЕНИЕ: Запуск отдельных команд управляется переменными, например: do_init = True Выбранный энергетический диапазон управляется переменной index forced=True делает принудительную фотометрию ПАРАМЕТРЫ: index : Выбранный энергетический диапазон ВЫВОД: Выходные файлы записываются в директорию outfile_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 pickle import uds from uds.utils import * from uds.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("UDS root path: {}".format(root_path)) infile_dir=root_path+'/data/processed' outfile_dir=root_path+'/products' create_folder(outfile_dir) local_run = False outkey="tm0" do_init = False do_merge = False 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 = False do_fixcat = False # only for index=0 do_fixxmm = True # only for index=0 do_apetool = False do_catprep = False do_filter_catalog = False do_cross_match = False index=0 forced=False """ If forced=True, take input catalog from energy range en0 """ comm='-xmm' # for 4XMM-DR12 forced photometry use '-xmm' vign=True vignetting = 'vign' if (vign==True) else 'novign' events=[] expmaps=[] bkgmaps=[] 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=do_init, do_obsmode=False, do_center=False, do_evtool=True, do_expmap=True, vign=vign, 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) """ Собираем общий список событий """ outfile_evtool="{}_EventList_en{}.fits".format(os.path.join(outfile_dir,outkey), eband[index]) if(do_merge==True): do_evtool_esass(events=events, outfile=outfile_evtool) """ makes detmask from TM exposures """ detmask="{}/{}_DetectorMask_en{}{}".format(outfile_dir, outkey, eband[index], outfile_post) if(do_detmask==True): create_detmask_merged(expmaps,detmask,minval=100) """ Собираем общую карту экспозиции, обратите внимание на коэффициент 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{}.{}{}".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) if(do_erbox1==True): cmd=["erbox", "images=\'{}\'".format(outfile_evtool), "boxlist=%s" %(outfile_boxlist1), "expimages=\'{}\'".format(outfile_expmap), "detmasks=\'{}\'".format(detmask), "emin=\'{}\'".format(emin_ev[index]), "emax=\'{}\'".format(emax_ev[index]), "ecf=\'{}\'".format(ecf[index]), "nruns=2", "likemin=6.0", "boxsize=4", "compress_flag=N", "bkgima_flag=N", "expima_flag=Y", "detmask_flag=Y" ] remove_file(outfile_boxlist1) print((" ").join(cmd)) os.system((" ").join(cmd)) 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) 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) if(do_erbox2==True): cmd=["erbox", "images=\'{}\'".format(outfile_evtool), "boxlist=%s" %(outfile_boxlist2), "expimages=\'{}\'".format(outfile_expmap), "detmasks=\'{}\'".format(detmask), "emin=\'{}\'".format(emin_ev[index]), "emax=\'{}\'".format(emax_ev[index]), "ecf=\'{}\'".format(ecf[index]), "nruns=2", "likemin=6.0", "boxsize=4", "compress_flag=N", "bkgima_flag=Y", "bkgimages={}".format(outfile_backmap1), "expima_flag=Y", "detmask_flag=Y" ] remove_file(outfile_boxlist2) print((" ").join(cmd)) os.system((" ").join(cmd)) 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) 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) if(do_erbox3==True): cmd=["erbox", "images=\'{}\'".format(outfile_evtool), "boxlist=%s" %(outfile_boxlist3), "expimages=\'{}\'".format(outfile_expmap), "detmasks=\'{}\'".format(detmask), "emin=\'{}\'".format(emin_ev[index]), "emax=\'{}\'".format(emax_ev[index]), "ecf=\'{}\'".format(ecf[index]), "nruns=2", "likemin=6.0", "boxsize=4", "compress_flag=N", "bkgima_flag=Y", "bkgimages={}".format(outfile_backmap2), "expima_flag=Y", "detmask_flag=Y" ] remove_file(outfile_boxlist3) print((" ").join(cmd)) os.system((" ").join(cmd)) 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) if(do_erbackmap3==True): boxlist3 = outfile_boxlist3 if(forced == False) else "{}/{}_BoxList3_en{}{}".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) """ 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) if(index==0): boxlist3="{}_MaxLikSourceList_en{}.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() add_specific_columns(boxlist3) fitpos_flag="fitpos_flag=no" fitext_flag="fitext_flag=no" nmulsou = "nmulsou=1" nmaxfit="nmaxfit=1" """ don't allow ermpldet to split sources """ 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) 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) boxlist3 = outfile_boxlist3 fitpos_flag="fitpos_flag=yes" fitext_flag="fitext_flag=yes" nmulsou = "nmulsou=2" nmaxfit="nmaxfit=4" """ allow ermpldet to split sources (no more than two) """ cmd=["ermldet", "mllist=%s" %(mllist), "boxlist=%s" %(boxlist3), "images=\'{}\'".format(outfile_evtool), "expimages=\'{}\'".format(outfile_expmap), "detmasks=\'{}\'".format(detmask), "bkgimages=\'{}\'".format(outfile_backmap3), "emin=\'{}\'".format(emin_ev[index]), "emax=\'{}\'".format(emax_ev[index]), "ecf=\'{}\'".format(ecf[index]), "hrdef=", "likemin=0.", "extlikemin=5.", "compress_flag=N", "cutrad=15.", "multrad=20.", "extmin=2.0", "extmax=35.0", #"bkgima_flag=Y", looks outdated "expima_flag=Y", "detmask_flag=Y", "shapelet_flag=yes", "photon_flag=yes", "extentmodel=beta", "thres_flag=N", "thres_col=like", "thres_val=30.", nmaxfit, nmulsou, fitpos_flag, fitext_flag, "srcima_flag=yes", "srcimages=\'{}\'".format(srcmap) ] 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), 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), expmap=outfile_expmap, backmap=outfile_backmap3,detlike=detlike, detmask=detmask, emin=emin_ev[index], emax=emax_ev[index],ecf=ecf[index]) """ if(do_ermldet==True): if(vign==False): print('Run ermldet with vignetted exposure!') sys.exit() remove_file(mllist) remove_file(srcmap) runme(cmd, local_run=local_run) #correct_fluxerr_ermldet_forced(mllist) if(forced==True): result = check_ermldet_forced(mllist) """ for a some reason, for an arbitrary energy band, ermldet break order of sources. Do this forced correction. """ if(result == False): correct_srcid_ermldet_forced(mllist) if(do_fixcat==True): if not index == 0: print("ERROR: You can fix only reference catalog for en0.") sys.exit() if forced == True: print("ERROR: You can fix only non-forced catalog for en0.") sys.exit() srcs_remove=[174, 90, 299, 300, 504, 215, 401, 20,] # keep 671 as unclassified extended source srcs_add = {'4XMM J021925.4-042647':[34.8559099,-4.4465007, 1.366], # 147 '4XMM J021922.8-042655':[34.8451832,-4.4487901, 1.958], # 147 '4XMM J021929.4-043224':[34.8728586,-4.5400022, 0.660], # 90 '4XMM J021931.2-043222':[34.8801169,-4.5395495, 2.561], # 90 '4XMM J021911.2-050550':[34.7968110,-5.0972990, 0.732], # 504 '4XMM J021919.3-050511':[34.8307176,-5.0864242,4.988], # 504 '4XMM J021911.5-050501':[34.7981099,-5.0837146,6.834], # 504 '4XMM J021658.9-044900':[34.2455964,-4.8168126,4.449], # 300 '4XMM J021659.2-044945':[34.2468704,-4.8291892,1.548], # 300 '4XMM J021812.2-045814':[34.5510753,-4.9705972,0.550], # 215 '4XMM J021812.0-045813':[34.5502698,-4.9703004,0.497], # 215 '4XMM J021912.6-052756':[34.8028459,-5.4656239,0.579], # 401 '4XMM J021705.5-042254':[34.2730294,-4.3816810,0.288], # 20 '4XMM J021705.3-042314':[34.2720952,-4.3873162,0.587], # 20 '4XMM J021827.2-045456':[34.6134256,-4.9157208,0.252], '4XMM J021831.3-045504':[34.6306930,-4.9178676,0.242], '4XMM J021925.4-045201':[34.8558373,-4.8671200,0.529], } fix_catalog(mllist=mllist,refimage=outfile_evtool, srcs_remove=srcs_remove, srcs_add=srcs_add) """ Note that fix_catalog added ID_SRC to each XMM source. Next, we save forced XMM sources (with new ID_SRC!) for later catalog compilation """ with open(mllist.replace(".fits", ".xmm.pickle"), 'wb') as f: pickle.dump(srcs_add, f) if(do_fixxmm==True): if not index == 0: print("ERROR: You can fix only reference catalog for en0.") sys.exit() if forced == True: print("ERROR: You can fix only non-forced catalog for en0.") sys.exit() fix_xmm_sources(mllist=mllist,refimage=outfile_evtool, xmm_catalog='../data/4XMM-DR12/4XMM_DR12cat_slim_v1.0_UDS.fits.catalog') if(do_apetool==True): psfmap="{}_PsfMap{}".format(os.path.join(outfile_dir,outkey), outfile_post) #remove_file(psfmap) #cmd=["apetool", # "images=\'{}\'".format(outfile_evtool), # "psfmaps=\'{}\'".format(psfmap), # "psfmapflag=yes",] #runme(cmd, local_run=local_run) cmd=["apetool", "mllist={}".format(mllist), "apelistout={}".format(mllist), # give the same file "images=\'{}\'".format(outfile_evtool), "expimages=\'{}\'".format(outfile_expmap), "detmasks=\'{}\'".format(detmask), "bkgimages=\'{}\'".format(outfile_backmap3), "emin=\'{}\'".format(emin_ev[index]), "emax=\'{}\'".format(emax_ev[index]), "srcimages=\'{}\'".format(srcmap), "psfmaps={}".format(psfmap), "psfmapflag=no", "stackflag=no", "apexflag=yes", "apesenseflag=no", "eefextract=0.65", "cutrad=15", "eindex=1",] 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) else: catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) if(do_catprep==True): cmd=["catprep", "infile={}".format(mllist), "outfile={}".format(catprep),] remove_file(catprep) runme(cmd, local_run=local_run) if(do_filter_catalog==True): #filter_mllist(mllist,expcut=5000.0,dlcut=10.0,dlmin=10,dlmax=10000) """ works the same """ filter_catprep(catprep,expcut=5000.0,dlmin=10,dlmax=10000,outkey='bright') filter_catprep(catprep,expcut=5000.0,dlmin=6,dlmax=10,outkey='faint') if(do_cross_match==True): crossmatch_shu2019(catprep,dlmin=10,refimage=outfile_evtool,crval=[ra_cen, de_cen], catalog=root_path+"/data/Gaia_unWISE/Gaia_unWISE_UDS.fits.catalog")