diff --git a/.gitignore b/.gitignore index 924b194..b13e084 100644 --- a/.gitignore +++ b/.gitignore @@ -54,3 +54,12 @@ venv.bak/ *.ref *.src *.dat +*.xfm +*.log +*.tmp +*.attcorr +*.aspsol +*.txt +*.png +.* +!.gitignore \ No newline at end of file diff --git a/scripts/01_init_events.py b/scripts/01_init_events.py index fe8c5ef..c8d9fc3 100755 --- a/scripts/01_init_events.py +++ b/scripts/01_init_events.py @@ -28,31 +28,31 @@ region="box({},{},4d,4d,0)".format(ra_cen,de_cen) """ # TM1 -do_evtool_esass(evlist=el+'tm1.txt', outfile=pr+'tm1_obs_1.fits', gti='621296896. 621304128.') +do_evtool_esass(evlist=el+'tm1.txt', outfile=pr+'tm1_obs_1.fits', gti='621296896. 621304128.', rmlock=True) # TM5 -do_evtool_esass(evlist=el+'tm5.txt', outfile=pr+'tm5_obs_1.fits', gti='620606016. 620614848.') -do_evtool_esass(evlist=el+'tm5.txt', outfile=pr+'tm5_obs_2.fits', gti='620676992. 620689792.') +do_evtool_esass(evlist=el+'tm5.txt', outfile=pr+'tm5_obs_1.fits', gti='620606016. 620614848.', rmlock=True) +do_evtool_esass(evlist=el+'tm5.txt', outfile=pr+'tm5_obs_2.fits', gti='620676992. 620689792.', rmlock=True) # TM6 -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_park_1.fits', gti='620174080. 620178002.620032') -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_scan_1.fits', gti='620178002.620032 620192246.62720') -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_park_2.fits', gti='620192448. 620194624.') -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_scan_2.fits', gti='620194666.606144 620208904.673408') -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_park_3.fits', gti='620209162.670976 620211316.650304') -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_scan_3.fits', gti='620211328. 620225600.') -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_park_4.fits', gti='620225853.609024 620227974.68832') -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_scan_4.fits', gti='620227904. 620242176.') -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_obs_1.fits', gti='620242432. 620258368.') +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_park_1.fits', gti='620174080. 620178002.620032', rmlock=True) +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_scan_1.fits', gti='620178002.620032 620192246.62720', rmlock=True) +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_park_2.fits', gti='620192448. 620194624.', rmlock=True) +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_scan_2.fits', gti='620194666.606144 620208904.673408', rmlock=True) +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_park_3.fits', gti='620209162.670976 620211316.650304', rmlock=True) +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_scan_3.fits', gti='620211328. 620225600.', rmlock=True) +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_park_4.fits', gti='620225853.609024 620227974.68832', rmlock=True) +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_scan_4.fits', gti='620227904. 620242176.', rmlock=True) +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_obs_1.fits', gti='620242432. 620258368.', rmlock=True) -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_obs_2.fits', gti='620607424. 620614656.') +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_obs_2.fits', gti='620607424. 620614656.', rmlock=True) do_badpix_tm6(filename=pr+'tm6_obs_2.fits') -do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_obs_3.fits', gti='620676992. 620690368.') +do_evtool_esass(evlist=el+'tm6.txt', outfile=pr+'tm6_obs_3.fits', gti='620676992. 620690368.', rmlock=True) 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, rmlock=True) -#do_evtool_esass(evlist=el+'tm7.txt', outfile=pr+'tm7_obs_2.fits', gti='621110272. 621117952.', region=region) +#do_evtool_esass(evlist=el+'tm7.txt', outfile=pr+'tm7_obs_2.fits', gti='621110272. 621117952.', region=region, rmlock=True) diff --git a/scripts/02_merge_events.py b/scripts/02_merge_events.py index 5c32cdf..aec6510 100755 --- a/scripts/02_merge_events.py +++ b/scripts/02_merge_events.py @@ -22,8 +22,12 @@ outkey="mosa_tm0" """ 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' @@ -33,18 +37,22 @@ create_folder(outfile_dir) index=5 # select energy band -do_init = True +do_init = False do_merge = True +do_rate = False do_adapt = False # requires CIAO -#do_rate = True 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]: + #if not ("scan" in datakey): + # continue print("--> {}".format(datakey)) """ Подготавливаем списки событий индивидуальных наблюдений """ outfile_evtool,outfile_expmap=init_events(key=datakey, eband_index=eband[index], @@ -61,26 +69,27 @@ for tmkey in keylist_tm.keys(): emax_kev=emax_kev[index]) events.append(outfile_evtool) expmaps.append(outfile_expmap) + bkgmaps.append("{}_BackMap3_en{}.fits".format(os.path.join(outfile_dir,datakey), eband[0])) """ Собираем общий список событий """ outfile_evtool="{}_EventList_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]) +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) if(do_merge==True): do_evtool_esass(events=events, outfile=outfile_evtool) - do_expmap_ftools(expmaps=expmaps, outfile=outfile_expmap) + do_fimgmerge_ftools(maps=expmaps, outfile=outfile_expmap) + do_fimgmerge_ftools(maps=bkgmaps, outfile=outfile_bkgmap) +outfile_rate="{}_RateMap_en{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], vignetting) +if(do_rate==True): + make_rate_map(cntmap=outfile_evtool, expmap=outfile_expmap, outfile=outfile_rate) + + function='gaussian' -outfile_adapt="{}_ImageAdapt_en{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], function) +outfile_adapt="{}_ImageAdapt_en{}.{}.{}.fits".format(os.path.join(outfile_dir,outkey), eband[index], function, vignetting) if(do_adapt==True): 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 2f630f8..9f6d9f6 100755 --- a/scripts/03_init_obs.py +++ b/scripts/03_init_obs.py @@ -1,6 +1,39 @@ #!/usr/bin/env python -"""Подготавливает списки событий в разных энергетических диапазонах. - Производит списки источников в наждом наблюдении и делает астрокоррекцию с помощью wcs_match +""" +НАЗВАНИЕ: + + 01_init_obs.py + + +НАЗНАЧЕНИЕ: + + Подготавливает списки событий в разных энергетических диапазонах. + Производит списки источников в наждом наблюдении и делает астрокоррекцию с помощью wcs_match/wcs_update + +ВЫЗОВ: + + conda activate ciao-4.15 + esass + ./01_init_obs.py + + +УПРАВЛЕНИЕ: + + Запуск отдельных команд управляется переменными, например: do_init = True + Выбранные энергетические диапазоны управляется массивом eband_selected + +ПАРАМЕТРЫ: + eband_selected : Выбранные энергетические диапазоны + + +ВЫВОД: + + Выходные файлы записываются в директорию outfile_dir + + +ИСТОРИЯ: + Роман Кривонос, ИКИ РАН, krivonos@cosmos.ru + Март 2023 """ @@ -26,25 +59,28 @@ infile_dir=root_path+'/data/processed' outfile_dir=root_path+'/products' create_folder(outfile_dir) -do_init = False -do_ermask = False +do_init = True +do_ermask = True -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_erbox1 = True # local mode +do_erbackmap1 = True # +do_erbox2 = True # map mode, with background map +do_erbackmap2 = True # +do_erbox3 = True # map mode, with background map +do_erbackmap3 = True # -do_ermldet = False -do_catprep = False -do_cross_match = True +do_ermldet = True +do_catprep = True +do_cross_match = False -do_wcs_match = False # Chandra task -do_wcs_update = False # Chandra task +do_wcs_match = True # Chandra task +do_wcs_update = True # Chandra task + +eband_selected=[0,1,2,3,4] + +vign=True +vignetting = 'vign' if (vign==True) else 'novign' -eband_selected=[0] - def runme(datakey): """ runs datakey over energy bands """ @@ -71,11 +107,17 @@ def runme(datakey): do_center=True, 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) + expmaps.append(outfile_expmap) + events.append(outfile_evtool) + + """ + After astrometry-corrected files (*.attcorr.fits) are obtained, one can take them as original, in order to check the full chain: + events.append(outfile_evtool.replace(".fits", ".attcorr.fits")) + """ """ Detmask """ @@ -305,6 +347,7 @@ def runme(datakey): #save_ermldet_ds9reg(mllist,scale=60*60) catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post) + catprep_en0="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,datakey), eband[0], outfile_post) if(do_catprep==True): cmd=["catprep", "infile={}".format(mllist), @@ -315,60 +358,27 @@ 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,dlmin=10,refimage=events[ii],crval=wcslist[datakey], catalog=root_path+"/data/Gaia_unWISE/Gaia_unWISE_UDS.fits.catalog") - """ - if(do_att_corr==True): - attcorr(mllist,events=outfile_evtool[ii]) - - - xfm=mllist.replace(".fits", ".xfm") - if(do_wcs_match==True): - src_list=mllist.replace(".fits", ".src") - if not (os.path.isfile(src_list)==True): - print("Not found {}".format(src_list)) - sys.exit() - - ref_list=mllist.replace(".fits", ".ref") - if not (os.path.isfile(ref_list)==True): - print("Not found {}".format(ref_list)) - sys.exit() - - #xfm=mllist.replace(".fits", ".xfm") - log=mllist.replace(".fits", ".xfm.log") - - cmd=["wcs_match", - "infile={}".format(src_list), - "refsrcfile={}".format(ref_list), - "outfile={}".format(xfm), - "wcsfile={}".format(src_list), - "logfile={}".format(log), - "radius=7", - "residlim=1", - "verbose=1", - "method=trans", - #"method=rst", - "clobber=yes", - ] - os.system((" ").join(cmd)) - print((" ").join(cmd)) - - + if(do_wcs_match==True and eband[index]==0): + """ run wcs_match for 0.3-2.3 keV only """ + wcs_match_ciao(catprep, method='rst',radius=12,residlim=5) + if(do_wcs_update==True): - wcs_update(outfile_evtool[ii],crval=wcslist[key],transformfile=xfm) - """ + """ use 0.3-2.3 keV transform matrix for all other bands """ + attcorr=wcs_update_ciao(events[ii],crval=wcslist[datakey],transformfile=catprep_en0.replace(".fits", ".xfm"),clean=False) + do_evtool_esass(evfile=attcorr,outfile=attcorr,rmlock=False, do_center=True, ra_cen=ra_cen, de_cen=de_cen) -runme("tm7_obs_1") """ +testing +runme("tm7_obs_1") +runme("tm5_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/04_mosaics.py b/scripts/04_mosaics.py new file mode 100755 index 0000000..21079df --- /dev/null +++ b/scripts/04_mosaics.py @@ -0,0 +1,308 @@ +#!/usr/bin/env python +""" +НАЗВАНИЕ: + + 04_mosaics.py + + +НАЗНАЧЕНИЕ: + + Собирает мозайки в разных энергетических диапазонах. + +ВЫЗОВ: + + esass + ./01_mosaics.py + + +УПРАВЛЕНИЕ: + + Запуск отдельных команд управляется переменными, например: do_init = True + Выбранный энергетический диапазон управляется переменной index + +ПАРАМЕТРЫ: + + 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 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) + +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 = True +do_catprep = True +do_filter_catalog = True + +#do_ermldet = False +#do_cross_match = False + + + +index=4 + +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=False, + do_expmap=False, + 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{}{}".format(os.path.join(outfile_dir,outkey), + eband[index], + 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], outfile_post) +cheese_mask="{}_CheeseMask3_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) +if(do_erbackmap3==True): + do_erbackmap_esass(outfile_evtool,outfile_expmap,outfile_boxlist3,detmask,emin_ev[index],emax_ev[index], + outfile_backmap3,cheese_mask) + +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) +cmd=["ermldet", + "mllist=%s" %(mllist), + "boxlist=%s" %(outfile_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=5.", + "extlikemin=6.", + "compress_flag=N", + "cutrad=10.", # was 15 + "multrad=20.", + "extmin=2.0", + "extmax=15.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=4", + "nmulsou=2", + "fitext_flag=yes", + "srcima_flag=yes", + "srcimages=\'{}\'".format(srcmap) + ] + +sensmap="{}_SensitivityMap_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) +area="{}_AreaTable_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) +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]) + + +if(do_ermldet==True): + remove_file(mllist) + remove_file(srcmap) + os.system((" ").join(cmd)) + print((" ").join(cmd)) + +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) + os.system((" ").join(cmd)) + print((" ").join(cmd)) + +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') + diff --git a/scripts/05_srctool.py b/scripts/05_srctool.py new file mode 100755 index 0000000..1a841c7 --- /dev/null +++ b/scripts/05_srctool.py @@ -0,0 +1,150 @@ +#!/usr/bin/env python +""" +НАЗВАНИЕ: + + 04_mosaics.py + + +НАЗНАЧЕНИЕ: + + Собирает мозайки в разных энергетических диапазонах. + +ВЫЗОВ: + + esass + ./01_mosaics.py + + +УПРАВЛЕНИЕ: + + Запуск отдельных команд управляется переменными, например: do_init = True + Выбранный энергетический диапазон управляется переменной index + +ПАРАМЕТРЫ: + + 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 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) + +srctool_dir="{}/{}".format(outfile_dir,"srctool-products") +create_folder(srctool_dir) + +outkey="tm0" + +outfile_srctool="{}_SrcTool_".format(outkey) + +do_init = False +do_merge = False +do_srctool = True +do_grppha = True + +index=4 +""" работаем именно в этом диапазоне, чтобы спектры покрывали все энергии """ + +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=False, + do_expmap=False, + 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) + + +suffix_srctool=".fits" +""" Output filename suffix - all output filenames appended with this string. +If suffix contains no filename extension (does not contain a "."), then ".fits" +is also appended to the filename. """ + +catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,outkey), eband[0], outfile_post) +""" take source catalog from 0.3-2.3 keV band """ +if not (os.path.isfile(catprep)==True): + print("{} not found, run 04_mosaics.py?".format(catprep)) + sys.exit() + + +if(do_srctool==True): + test_exe('srctool') + cmd=['srctool', + "insts=\'1 5 6 7\'", + "eventfiles={}".format(outfile_evtool), + "prefix=\'{}\'".format(os.path.join(srctool_dir,outfile_srctool)), + "suffix=\'{}\'".format(suffix_srctool), + "srccoord={}".format(catprep), + #"srcreg=\'fk5;circle * * 60s\'", + #"backreg=\'fk5;annulus * * 90s 120s\'", + "srcreg=AUTO", + "backreg=AUTO", + "clobber=yes",] + os.system((" ").join(cmd)) + print((" ").join(cmd)) + +if(do_grppha==True): + group_spectra("{}/*_SourceSpec_*.fits".format(srctool_dir)) + + + diff --git a/scripts/README.md b/scripts/README.md index 2d3f14d..7b9cbb2 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -25,9 +25,12 @@ source /eSASS4EDR/bin/esass-init.csh ### 03_init_obs.py - 1) Подготавливает списки событий в разных энергетических диапазонах. 2) Запускает ```erbox``` в три этапа, чтобы получить рабочий список источников для ```ermldet```. 3) Запускает ```ermldet``` 4) Делает кросс-корреляцию с каталогом Gaia-unWISE ```do_cross_match=True```, которая создает три файла: ```.cross``` -- все пересечения, и ```.ref``` / ```.src``` -- входные каталоги для последующей команды ```wcs_match``` +5) Делает матрицу преобразования координат и корректирует списки событий. Для запуска команд```wcs_match/wcs_update``` требуется запустить окружение ```ciao``` (см. выше) +### 04_mosaics.py + +Создает сборные изображения (мозайки) в разных энергетических диапазонах. \ No newline at end of file diff --git a/uds/uds/config.py b/uds/uds/config.py index 2b2b39e..1859eb3 100644 --- a/uds/uds/config.py +++ b/uds/uds/config.py @@ -43,14 +43,15 @@ wcslist={'tm1_obs_1':[34.7279760,-5.0680267], 'tm6_scan_3':[34.5405596,-4.8088748], 'tm6_scan_4':[34.5405596,-4.8088748]} -""" Диапазоны энергий. """ -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", "5"] +eband=[ "0", "1", "2", "3", "4", "5"] +""" Диапазоны энергий. """ +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] + outfile_post='.fits' diff --git a/uds/uds/utils.py b/uds/uds/utils.py index 36afd66..44a9bd5 100644 --- a/uds/uds/utils.py +++ b/uds/uds/utils.py @@ -47,17 +47,35 @@ 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, rmlock=False): +def do_evtool_esass(evfile=None,events=None,outfile=None,evlist=None, + gti=None,region=None,emin=None,emax=None, rmlock=False, + do_center=False, ra_cen=None, de_cen=None): eventfiles=None if(events): eventfiles="eventfiles=\'{}\'".format((" ").join(events)) if(evlist): eventfiles="eventfiles=@{}".format(evlist) + if(evfile): + eventfiles="eventfiles={}".format(evfile) if not (eventfiles): print("ERROR: Event files not provided") + if(do_center==True and evfile): + """ re-center event file """ + if not (ra_cen and de_cen): + print("Please provide center coordinates") + sys.exit() + test_exe('radec2xy') + cmd=["radec2xy", + "{}".format(evfile), + "{}".format(ra_cen), + "{}".format(de_cen)] + print((" ").join(cmd)) + os.system((" ").join(cmd)) + + emin="emin={}".format(emin) if(emin) else '' emax="emax={}".format(emax) if(emax) else '' gti="gti=\'{}\'".format(gti) if(gti) else '' @@ -83,6 +101,7 @@ def do_evtool_esass(events=None,outfile=None,evlist=None,gti=None,region=None,em if(os.path.isfile(filename)==True): os.remove(filename) pass + def set_bit(value, bit): return value | (1< expcut) + 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 @@ -588,3 +608,400 @@ def make_rate_map(cntmap=None, expmap=None, outfile=None): fits.writeto(outfile, rate, overwrite=True) +def wcs_match_ciao(catalog=None, method='trans',radius=15,residlim=7): + + xfm=catalog.replace(".fits", ".xfm") + + src_list=catalog.replace(".fits", ".shu2019.src.fits") + if not (os.path.isfile(src_list)==True): + print("Not found {}".format(src_list)) + print("Did you run cross-match?") + sys.exit() + + ref_list=catalog.replace(".fits", ".shu2019.ref.fits") + if not (os.path.isfile(ref_list)==True): + print("Not found {}".format(ref_list)) + print("Did you run cross-match?") + sys.exit() + + log=catalog.replace(".fits", ".xfm.log") + + cmd=["wcs_match", + "infile={}".format(src_list), + "refsrcfile={}".format(ref_list), + "outfile={}".format(xfm), + "wcsfile={}".format(src_list), + "logfile={}".format(log), + "radius={}".format(radius), + "residlim={}".format(residlim), + "verbose=1", + "method={}".format(method), + "clobber=yes", + ] + test_exe('wcs_match') + os.system((" ").join(cmd)) + print((" ").join(cmd)) + +def wcs_update_ciao(events,ext=2,crval=None,transformfile=None,clean=True): + """ Prepares eRosita attitude file for Chandra's wcs_update command """ + + if not (transformfile): + print("ERROR: Please provide transformfile") + return + + if not (os.path.isfile(transformfile)==True): + print("Not found {}".format(transformfile)) + print("Did you run wcs_update_ciao?") + sys.exit() + + fnaspsol=events.replace(".fits", ".aspsol") + fntmp=events.replace(".fits", ".tmp") + fnattcor=events.replace(".fits", ".attcorr.fits") + hdul = fits.open(events,mode='readonly') + corratt=hdul[ext] + naxis1=corratt.header['NAXIS1'] + naxis2=corratt.header['NAXIS2'] + corratt_extname=corratt.header['EXTNAME'] + + """ create new empty Q_ATT column """ + q_att=np.zeros((naxis2,4)) + cols = [] + cols.append(fits.Column(name='q_att', format='4D', array=q_att)) + + orig_cols = corratt.columns + new_cols = fits.ColDefs(cols) + + hdu0 = fits.PrimaryHDU() + + hdu0.header['EXTNAME'] = ('ASPSOL', 'Name of this binary table extension') + hdu0.header['HDUNAME'] = ('ASPSOL', 'ASCDM block name') + + if(crval): + hdu0.header['RA_NOM']=(crval[0], '[deg] Pointing RA') + hdu0.header['DEC_NOM']=(crval[1], '[deg] Pointing Dec') + hdu0.header['ROLL_NOM']=(0.0 , '[deg] Pointing Roll') + + hdu0.header['RA_PNT']=(crval[0], '[deg] Nominal RA') + hdu0.header['DEC_PNT']=(crval[1], '[deg] Nominal Dec') + hdu0.header['ROLL_PNT']=(0.0 , '[deg] Nominal Roll') + + hdu = fits.BinTableHDU.from_columns(orig_cols + new_cols, header=hdu0.header) + hdu.writeto(fnaspsol, overwrite=True) + + """ + First, update attitude information in ASPSOL (Chandra notation), + which corresponds to CORRATT (eRosita notation) + """ + cmd=["wcs_update", + "infile={}".format(fnaspsol), + "outfile={}".format(fntmp), + "transformfile={}".format(transformfile), + "clobber=yes", + ] + os.system((" ").join(cmd)) + print((" ").join(cmd)) + + """ + Replace attitude extension + """ + shutil.copyfile(events, fnattcor) + data, hdr = getdata(fntmp, 1, header=True) + hdr['EXTNAME']=corratt_extname + update(fnattcor, data, ext, header=hdr) + + """ + Second, calculate new RA/Dec for events using updated attitude, using evatt from eSASS + """ + cmd=["evatt", + "EVENTFILE={}".format(fnattcor), + "ATTFILE={}".format(fnattcor), + "GTIFILE={}".format(fnattcor), + ] + os.system((" ").join(cmd)) + print((" ").join(cmd)) + + if(clean==True): + remove_file(fnaspsol) + remove_file(fntmp) + return fnattcor + + +def create_detmask_merged(expmaps,outfile,minval=None): + """ + creates mask from each exposure and then stacks each masks. + Header is taken from first exposure file. + """ + tmlist={} + for expmap in expmaps: + + hdul = fits.open(expmap) + emap = hdul[0].data + ehdr = hdul[0].header + instrume = ehdr['INSTRUME'] + threshold = minval if(minval) else 0 + mask = np.where(emap > threshold, 1, 0) + + print("--- detmask {} -- {}".format(expmap,instrume)) + if not (instrume in tmlist.keys()): + tmlist[instrume]=1 + else: + pass + + if 'merged_mask' in locals(): + merged_mask = np.add(merged_mask, mask) + else: + merged_mask = mask + merged_hdr = ehdr + + if 'merged_hdr' in locals(): + index=1 + for tm in tmlist.keys(): + merged_hdr["INSTRUM{}".format(index)]=tm + index=index+1 + merged_hdr['INSTRUME']='merged' + merged_hdr['NINST']=len(tmlist) + merged_hdr['OBS_MODE']=' ' + if 'merged_mask' in locals() and 'merged_hdr' in locals(): + fits.writeto(outfile, np.where(merged_mask > 0, 1, 0), header=merged_hdr, overwrite=True) + +def create_expmap_merged(expmaps,outfile,scale=7.0): + """ + Adds exposure from TMs. Header is taken from first exposure file. + """ + tmlist={} + for expmap in expmaps: + + hdul = fits.open(expmap) + emap = np.divide(hdul[0].data, scale) + ehdr = hdul[0].header + instrume = ehdr['INSTRUME'] + + print("--- expmap {} -- {}".format(expmap,instrume)) + if not (instrume in tmlist.keys()): + tmlist[instrume]=1 + else: + pass + + if 'merged_map' in locals(): + merged_map = np.add(merged_map, emap) + else: + merged_map = emap + merged_hdr = ehdr + + if 'merged_hdr' in locals(): + index=1 + for tm in tmlist.keys(): + merged_hdr["INSTRUM{}".format(index)]=tm + index=index+1 + merged_hdr['INSTRUME']='merged' + merged_hdr['NINST']=len(tmlist) + merged_hdr['OBS_MODE']=' ' + if 'merged_map' in locals() and 'merged_hdr' in locals(): + fits.writeto(outfile, merged_map, header=merged_hdr, overwrite=True) + +def do_erbackmap_esass(image,expimage,boxlist,detmask,emin,emax,outfile_backmap,cheese_mask): + test_exe('erbackmap') + cmd=["erbackmap", + "image=%s" %(image), + "expimage=%s" %(expimage), + "boxlist={}".format(boxlist), + "detmask=%s" %(detmask), + "emin=%s" %(emin), + "emax=%s" %(emax), + "bkgimage=%s" %(outfile_backmap), + "cheesemask=%s" %(cheese_mask), + "idband=1", + "scut=0.001", + "mlmin=6", + "maxcut=0.5", + "fitmethod=smooth smoothval=15", + "snr=40.", + ] + + remove_file(cheese_mask) + remove_file(outfile_backmap) + os.system((" ").join(cmd)) + print((" ").join(cmd)) + +def filter_catprep(filename, expcut=100,dlmin=6.0,dlmax=10,scale=60*60,ext_like=0.0,outkey='main'): + if(os.path.isfile(filename)==False): + print("File not found {}".format(filename)) + print("Filter catprep {}".format(filename)) + fout_selected=filename.replace(".fits", ".{}.selected.reg".format(outkey)) + fout_skip=filename.replace(".fits", ".{}.skip.reg".format(outkey)) + fout_extended=filename.replace(".fits", ".extended.reg") + + hdul = fits.open(filename) + tbdata = hdul[1].data + + catsel=[] + catskip=[] + catext=[] + skip_count=0 + selected_count=0 + keepalive_count=0 + for rec in tbdata: + if (rec['ext_like'] > ext_like): + catext.append({'ra':rec['ra'],'dec':rec['dec'],'radec_err':rec['radec_err'], + 'det_like':rec['det_like_0'],'ext_like':rec['ext_like'], + 'ext':rec['ext'],'ext_err':rec['ext_err'], + 'src_id':rec['id_src']}) + if ((rec['det_like_0'] > dlmin and rec['det_like_0'] < dlmax)): + catsel.append({'ra':rec['ra'],'dec':rec['dec'],'radec_err':rec['radec_err'],'det_like':rec['det_like_0'], + 'src_id':rec['id_src']}) + selected_count=selected_count + 1 + else: + catskip.append({'ra':rec['ra'],'dec':rec['dec'],'radec_err':rec['radec_err'],'det_like':rec['det_like_0'], + 'src_id':rec['id_src']}) + skip_count = skip_count+1 + + print("total={} skip_count={} selected={}".format(len(tbdata),skip_count,selected_count)) + with open(fout_selected, 'w') as writer: + for rec in catsel: + writer.write("fk5;circle({}, {}, {}) # color=white text={{{} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale, + rec['src_id'],rec['det_like'])) + + with open(fout_skip, 'w') as writer: + for rec in catskip: + writer.write("fk5;circle({}, {}, {}) # color=red text={{{} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale, + rec['src_id'],rec['det_like'])) + + with open(fout_extended, 'w') as writer: + for rec in catext: + writer.write("fk5;circle({}, {}, {}) # color=magenta text={{{} {:.2f} {:.2f}}}\n".format(rec['ra'], + rec['dec'], + rec['ext']/scale, + rec['src_id'], + rec['det_like'], + rec['ext_like'])) +def filter_mllist(filename, expcut=100,dlcut=6.0,dlmin=6.0,dlmax=10,scale=60*60,ext_like=0.0): + if(os.path.isfile(filename)==False): + print("File not found {}".format(filename)) + print("Filter mllist {}".format(filename)) + fout_selected=filename.replace(".fits", ".selected.reg") + fout_skip=filename.replace(".fits", ".skip.reg") + fout_extended=filename.replace(".fits", ".extended.reg") + + hdul = fits.open(filename) + tbdata = hdul[1].data + + """ create list of unique source IDs """ + id_src_set = set(tbdata['id_src']) + src_id_list = list(id_src_set) + catsel=[] + catskip=[] + catext=[] + skip_count=0 + selected_count=0 + keepalive_count=0 + for src_id in src_id_list: + + tab_tm0=tbdata[(tbdata['id_src'] == src_id) & (tbdata['id_inst'] == 0) & (tbdata['id_band'] == 0)] + """ averaged over all instruments and bands """ + + ext_tm0=tbdata[(tbdata['id_src'] == src_id) & (tbdata['id_inst'] == 0) & (tbdata['id_band'] == 0) & (tbdata['ext_like'] > ext_like)] + """ averaged over all instruments and bands """ + if(len(ext_tm0)>0): + catext.append({'ra':ext_tm0['ra'][0],'dec':ext_tm0['dec'][0],'radec_err':ext_tm0['radec_err'][0], + 'det_like':ext_tm0['det_like'][0],'ext_like':ext_tm0['ext_like'][0], + 'ext':ext_tm0['ext'][0],'ext_err':ext_tm0['ext_err'][0], + 'src_id':ext_tm0['id_src'][0]}) + + mask_src_id = (tbdata['id_src'] == src_id) & (tbdata['id_inst'] == 0) + """ get all source records except merged one """ + + + tab_src=tbdata[(tbdata['id_src'] == src_id) & (tbdata['id_inst'] != 0) & (tbdata['id_band'] == 1) & (tbdata['ml_exp'] > expcut)] + """ only individual instruments, and only first energy band is selected here """ + + + keepalive=False + mask = tab_src['det_like'] > dlcut + if (any(mask) and tab_tm0['det_like'] < dlmin): + keepalive=True + keepalive_count=keepalive_count+1 + + + if(keepalive): + tab_src=tbdata[(tbdata['id_src'] == src_id) & (tbdata['id_inst'] != 0) & (tbdata['id_band'] == 1)] + print("KEEP ALIVE ID={} DL={:.2f} | radec: {:.4f} {:.4f} DL0={:.2f} TEXP={:.1f}".format(src_id, tab_tm0['det_like'][0], tab_tm0['ra'][0], tab_tm0['dec'][0], + tab_src['det_like'][0],tab_src['ml_exp'][0])) + + + + + + if ((tab_tm0['det_like'] > dlmin and tab_tm0['det_like'] < dlmax) or keepalive): + catsel.append({'ra':tab_tm0['ra'][0],'dec':tab_tm0['dec'][0],'radec_err':tab_tm0['radec_err'][0],'det_like':tab_tm0['det_like'][0], + 'src_id':tab_tm0['id_src'][0]}) + selected_count=selected_count + 1 + else: + catskip.append({'ra':tab_tm0['ra'][0],'dec':tab_tm0['dec'][0],'radec_err':tab_tm0['radec_err'][0],'det_like':tab_tm0['det_like'][0], + 'src_id':tab_tm0['id_src'][0]}) + skip_count = skip_count+1 + + + print("total={} skip_count={} keepalive_count={} selected={}".format(len(src_id_list),skip_count,keepalive_count,selected_count)) + with open(fout_selected, 'w') as writer: + for rec in catsel: + writer.write("fk5;circle({}, {}, {}) # color=white text={{{} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale, + rec['src_id'],rec['det_like'])) + + with open(fout_skip, 'w') as writer: + for rec in catskip: + writer.write("fk5;circle({}, {}, {}) # color=red text={{{} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale, + rec['src_id'],rec['det_like'])) + + with open(fout_extended, 'w') as writer: + for rec in catext: + writer.write("fk5;circle({}, {}, {}) # color=magenta text={{{} {:.2f} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['ext']/scale, + rec['src_id'],rec['det_like'],rec['ext_like'])) + + +def create_sensmap(sensmap=None,expmap=None,backmap=None,detmask=None,areatab=None, + emin=None,emax=None,ecf=None, detlike=10): + test_exe('ersensmap') + cmd=['ersensmap', + "expimages=\'{}\'".format(expmap), + "bkgimages=\'{}\'".format(backmap), + "detmasks=\'{}\'".format(detmask), + "sensimage={}".format(sensmap), + "emin=\'{}\'".format(emin), + "emax=\'{}\'".format(emax), + "ecf=\'{}\'".format(ecf), + "area_table={}".format(areatab), + "method=APER", + "aper_type=BOX", + "aper_size=4.5", + "likemin={}".format(detlike), + "extlikemin=6.0", + "ext=6.0", + "extentmodel=beta", + "detmask_flag=yes", + "shapelet_flag=no", + "photon_flag=no", + "ext_flag=no", + "extlike_flag=no", + "compress_flag=no", + "area_flag=no",] + remove_file(sensmap) + remove_file(areatab) + os.system((" ").join(cmd)) + print((" ").join(cmd)) + +def group_spectra(path,group_min=25): + test_exe('grppha') + flist=glob.glob(path) + for f in flist: + fout=f.replace(".fits", ".fits.grp") + print(f) + cmd=["grppha", + "infile={}".format(f), + "outfile={}".format(fout), + "comm=\'GROUP MIN {}\'".format(group_min), + "tempc=EXIT", + "clobber=yes", + "chatter=0", # with chatter <= 5 being very quite and chatter >= 20 very verbose. + ] + print((" ").join(cmd)) + os.system((" ").join(cmd))