coma/scripts/04_mosaics_joint.py

#!/usr/bin/env python
# -*- coding:  utf-8  -*-
"""
НАЗВАНИЕ:

 04_mosaics.py
 
 
НАЗНАЧЕНИЕ:

   Собирает мозайки в разных энергетических диапазонах.
 
ВЫЗОВ:

   esass
   ./04_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 coma

from coma.utils import *
from coma.config import *


""" find Coma root dir """
root_path=dirname(dirname(dirname(inspect.getfile(coma))))

cwd = os.path.dirname(os.path.realpath(__file__))


"""
ftools does not like long file path names,
for this reason, we use relative path here
"""
#root_path='..'

print("Coma root path: {}".format(root_path))

infile_dir=root_path+'/data/processed'
outfile_dir=root_path+'/products'
create_folder(outfile_dir)

local_run = False

# check ../products/joint.evtool

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_ermldet_ds9reg  = True


do_resid           = False # residuals of data and source map
do_fixcat          = False # only for index=0
do_fixxmm          = False # prepare forced photometry, only for index=0
do_apetool         = False
do_catprep         = True
do_ext_bad         = True # mark extended sources as 'BAD'

do_filter_catalog  = True
do_cross_match     = True

index=0

forced=False
""" If forced=True, take input catalog from energy range en0 """

comm='' # for 4XMM-DR12 forced photometry use '-xmm'

vign=True
attcorr=True

outkey  = "mosa_joint_tm0{}".format('_attcorr' if (attcorr==True) else '')

vignetting = 'vign' if (vign==True) else 'novign'

outfile_evtool="{}/mosa_joint_tm0_attcorr_EventList_en{:02d}.fits".format(outfile_dir,eband[index]) # created by products/joint.evtool
outfile_expmap="{}/mosa_joint_tm0_attcorr_ExposureMap_en{:02d}.{}.fits".format(outfile_dir,eband[index],vignetting)
detmask="{}/mosa_parts_tm0_attcorr_DetectorMask_en{:02d}.fits".format(outfile_dir,0)






_emin_ev = "{}".format(emin_ev[index])
_emax_ev = "{}".format(emax_ev[index])
_ecf = "{}".format(ecf[index])
    
if(do_expmap==True):
    add_expmaps(["{}/mosa_parts_tm0_attcorr_ExposureMap_en{:02d}.{}.fits".format(outfile_dir,eband[index],vignetting),
                 "{}/mosa_survey_tm0_attcorr_ExposureMap_en{:02d}.{}.fits".format(outfile_dir,eband[index],vignetting)],
            "{}/mosa_joint_tm0_attcorr_ExposureMap_en{:02d}.{}.fits".format(outfile_dir,eband[index],vignetting))


outfile_boxlist1="{}/{}_BoxList1_en{:02d}{}".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),
        "emax=\'{}\'".format(_emax_ev),
        "ecf=\'{}\'".format(_ecf),
        "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{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
cheese_mask="{}_CheeseMask1_en{:02d}{}".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{:02d}{}".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),
        "emax=\'{}\'".format(_emax_ev),
        "ecf=\'{}\'".format(_ecf),
        "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{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
cheese_mask="{}_CheeseMask2_en{:02d}{}".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{:02d}{}".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),
        "emax=\'{}\'".format(_emax_ev),
        "ecf=\'{}\'".format(_ecf),
        "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{:02d}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], vignetting, outfile_post)
cheese_mask="{}_CheeseMask3_en{:02d}.{}{}".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{:02d}{}".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{:02d}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)
    srcmap="{}_SourceMap_en{:02d}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)
    residmap="{}_ResidMap_en{:02d}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)

    """ for en1,2,3,4,5,6,7 give mllist from en0 as input """
    boxlist3="{}_MaxLikSourceList_en{:02d}{}{}".format(os.path.join(outfile_dir,outkey), eband[0], comm, outfile_post)
    #if(index==0):
    boxlist3="{}_MaxLikSourceList_en{:02d}.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=10"
    multrad="multrad=15." 
    cutrad="cutrad=15."  
    if(index == 3 or index == 4):
        """ for hard band take unvignetted background """
        outfile_backmap3="{}_BackMap3_en{:02d}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], "novign", outfile_post)

else:
    mllist="{}_MaxLikSourceList_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
    srcmap="{}_SourceMap_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
    #mllist="{}_MaxLikSourceList_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
    #srcmap="{}_SourceMap_en{:02d}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post)
    residmap="{}_ResidMap_en{:02d}{}".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"
    multrad="multrad=20." 
    cutrad="cutrad=20." 
    """ allow ermpldet to split sources (no more than two) """ 

cmd=["ermldet",
     "mllist={}".format(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, 
     multrad, 
     "extmin=2.0",
     "extmax=35.0",
     #"bkgima_flag=Y", looks outdated
     "expima_flag=Y",
     "detmask_flag=Y",
     "shapelet_flag=no", # !!!
     "photon_flag=no",   # !!!
     "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):
   methods = ['nearest',]# 'linear', 'cubic']

   detlike=10
   sensmap="{}_SensitivityMap_dl{}_en{:02d}{}".format(os.path.join(outfile_dir,outkey), detlike, eband[index], outfile_post)
   create_sensmap(sensmap=sensmap,
                  areatab="{}_AreaTable_dl{}_en{:02d}{}".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], local_run=local_run, cwd=cwd)   
   
   for method in methods:
       print("Detlike {}, Method {}".format(detlike,method))
       corrmap="{}_SensitivityMap_{}_dl{}_en{:02d}{}".format(os.path.join(outfile_dir,outkey), method, detlike, eband[index], outfile_post)
       sensmap_corr(sensmap=sensmap, output=corrmap, method=method)
   
   """
   detlike=6
   sensmap="{}_SensitivityMap_dl{}_en{:02d}{}".format(os.path.join(outfile_dir,outkey), detlike, eband[index], outfile_post)
   create_sensmap(sensmap=sensmap,
                  areatab="{}_AreaTable_dl{}_en{:02d}{}".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], cwd=cwd)
   for method in methods:
       print("Detlike {}, Method {}".format(detlike,method))
       corrmap="{}_SensitivityMap_{}_dl{}_en{:02d}{}".format(os.path.join(outfile_dir,outkey), method, detlike, eband[index], outfile_post)
       sensmap_corr(sensmap=sensmap, output=corrmap, method=method)
   """

if(do_ermldet==True):
    
    test_exe('ermldet')
    if(vign==False):
        print('Run ermldet with vignetted exposure!')
        sys.exit()
    remove_file(mllist)
    remove_file(srcmap)
    print(cmd)
    runme(cmd, local_run=local_run, cwd=cwd)
    print(cmd)
    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_ermldet_ds9reg==True):   
    save_ermldet_ds9reg_dl(mllist,scale=60*60, dl=6, point=True)
    save_ermldet_ds9reg_id(mllist,scale=60*60, dl=6, point=True)
    
    correct_fluxerr_ermldet_forced(mllist)
    
    

if(do_resid==True):
    do_resid_map(data=outfile_evtool, model=srcmap, outfile=residmap)

    
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=[3200,1872,1134,2171,1219,448,357,973,4423,5195,3215,1134,119,
                 622,2916,2915,2315,6870,824,4027,1463,518,3726,393,4049,1876,2396,1569,
                 4837,4286,1870,3311,6752,3691,4273,1500,5303,6207,3327,133] 
    srcs_add = {'4XMM J130526.4+285519':[196.3603538, 28.9221712, 0.853],# 3200
                '4XMM J130526.8+285452':[196.3617587,28.9147139, 0.938],# 3200 
                '4XMM J130123.8+284744':[195.3494948, 28.7956932, 0.692], # 1872
                '4XMM J125555.1+283110':[193.9797853,28.5196369,1.648], # 448
                '4XMM J125708.3+264926':[194.2847924,26.8239298,0.220], 
                #'4XMM J022026.3-050251':[35.1098619,-5.0476199,0.551],
                #'4XMM J021925.4-042647':[34.8559099,-4.4465007,1.366],
                #'4XMM J021910.9-045108':[34.7954311,-4.8522901,0.898],
                #'4XMM J021945.2-045331':[34.9383593,-4.8919843,1.538],
                #'4XMM J021733.8-051311':[34.3910215,-5.2199877,2.247],
                }
    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)
    #sys.exit()
    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, cwd=cwd)

if(forced==True):
    catprep="{}_SourceCatalog_en{:02d}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[index], comm, outfile_post)
else:
    catprep="{}_SourceCatalog_en{:02d}{}".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_ext_bad==True):
    make_ext_bad(infile=catprep, bkg_map=outfile_backmap3, bkg_cut=1.05, exception=[388,])

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=500.0,dlmin=10,dlmax=10000000,outkey='dl10')
    #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_Coma.fits.catalog", errlim=5.0)

    # confused sources according to XMM data
    # 194.2812310 27.2174300