diff --git a/.gitignore b/.gitignore index b13e084..eb5a742 100644 --- a/.gitignore +++ b/.gitignore @@ -61,5 +61,9 @@ venv.bak/ *.aspsol *.txt *.png +*.grp +*.csv .* -!.gitignore \ No newline at end of file +*.qdp +*.pickle +!.gitignore diff --git a/data/xspec/ML0001.xcm b/data/xspec/ML0001.xcm new file mode 100644 index 0000000..708ac16 --- /dev/null +++ b/data/xspec/ML0001.xcm @@ -0,0 +1,11 @@ +method leven 10 0.01 +abund wilm +xsect vern +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +model TBabs*powerlaw + 0.0123819 1 0 0 100000 1e+06 + 2 -1 -3 -2 9 10 + 0.000106473 0.01 0 0 1e+20 1e+24 +bayes off diff --git a/data/xspec/ML0003.xcm b/data/xspec/ML0003.xcm new file mode 100644 index 0000000..80ada6e --- /dev/null +++ b/data/xspec/ML0003.xcm @@ -0,0 +1,11 @@ +method leven 10 0.01 +abund wilm +xsect vern +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +model TBabs*powerlaw + 0.00327282 1 0 0 100000 1e+06 + 2 -1 -3 -2 9 10 + 0.000121635 0.01 0 0 1e+20 1e+24 +bayes off diff --git a/data/xspec/ML0004.xcm b/data/xspec/ML0004.xcm new file mode 100644 index 0000000..6789831 --- /dev/null +++ b/data/xspec/ML0004.xcm @@ -0,0 +1,11 @@ +method leven 10 0.01 +abund wilm +xsect vern +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +model TBabs*powerlaw + 0.00240283 1 0 0 100000 1e+06 + 2 -1 -3 -2 9 10 + 0.000297459 0.01 0 0 1e+20 1e+24 +bayes off diff --git a/data/xspec/ML0005.xcm b/data/xspec/ML0005.xcm new file mode 100644 index 0000000..6f0828e --- /dev/null +++ b/data/xspec/ML0005.xcm @@ -0,0 +1,11 @@ +method leven 10 0.01 +abund wilm +xsect vern +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +model TBabs*powerlaw + 0.0320991 1 0 0 100000 1e+06 + 2 -1 -3 -2 9 10 + 0.000123814 0.01 0 0 1e+20 1e+24 +bayes off diff --git a/data/xspec/ecf.xcm b/data/xspec/ecf.xcm new file mode 100644 index 0000000..5beb7b6 --- /dev/null +++ b/data/xspec/ecf.xcm @@ -0,0 +1,11 @@ +method leven 10 0.01 +abund wilm +xsect vern +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +model TBabs*powerlaw + 0.02 -1 0 0 100000 1e+06 + 2 -1 -3 -2 9 10 + 4.83094e-05 0.01 0 0 1e+20 1e+24 +bayes off diff --git a/data/xspec/load.xcm b/data/xspec/load.xcm new file mode 100644 index 0000000..0f0202d --- /dev/null +++ b/data/xspec/load.xcm @@ -0,0 +1,8 @@ +data tm0_SrcTool_020_SourceSpec_00005.fits.grp +ign bad +ign **-0.2 +ign 10.-** +cpd /xs +setpl en +pl lda + diff --git a/products/srctool-products/ecf.xcm b/products/srctool-products/ecf.xcm new file mode 100644 index 0000000..5beb7b6 --- /dev/null +++ b/products/srctool-products/ecf.xcm @@ -0,0 +1,11 @@ +method leven 10 0.01 +abund wilm +xsect vern +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +model TBabs*powerlaw + 0.02 -1 0 0 100000 1e+06 + 2 -1 -3 -2 9 10 + 4.83094e-05 0.01 0 0 1e+20 1e+24 +bayes off diff --git a/products/srctool-products/rate.pro b/products/srctool-products/rate.pro new file mode 100644 index 0000000..a9f8f8e --- /dev/null +++ b/products/srctool-products/rate.pro @@ -0,0 +1,13 @@ +pro rate + ar=[[1.8847140728347272e-13, 0.37631887197494507],$ + [1.521645609807299e-13, 0.40992528200149536],$ + [4.915687878270888e-13, 1.0217742919921875],$ + [1.7903951117830771e-13, 0.282895565032959],$ + [6.886249544234846e-14, 0.21058528125286102],$ + [9.49530897031588e-14, 0.2018834501504898]] + +for i=0L,5 do begin +print, i, (ar[1,i])/(ar[0,i]) +endfor + +end diff --git a/scripts/04_mosaics.py b/scripts/04_mosaics.py index 21079df..6f98997 100755 --- a/scripts/04_mosaics.py +++ b/scripts/04_mosaics.py @@ -19,6 +19,7 @@ Запуск отдельных команд управляется переменными, например: do_init = True Выбранный энергетический диапазон управляется переменной index + forced=True делает принудительную фотометрию ПАРАМЕТРЫ: @@ -66,30 +67,27 @@ create_folder(outfile_dir) outkey="tm0" -do_init = False -do_merge = False -do_detmask = False -do_expmap = False +do_init = True +do_merge = True +do_detmask = True +do_expmap = True +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_ersensmap = True +do_ermldet = False +do_catprep = False +do_filter_catalog = False +do_cross_match = False +index=1 -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 # +forced=True +""" If forced=True, take input catalog from energy range en0 """ -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' @@ -110,7 +108,7 @@ for tmkey in keylist_tm.keys(): do_obsmode=False, do_center=False, do_evtool=False, - do_expmap=False, + do_expmap=True, vign=vign, ra_cen=ra_cen, de_cen=de_cen, emin_kev=emin_kev[index], @@ -138,9 +136,11 @@ if(do_detmask==True): """ Собираем общую карту экспозиции, обратите внимание на коэффициент 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], +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) @@ -231,17 +231,38 @@ if(do_erbox3==True): 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) +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): - do_erbackmap_esass(outfile_evtool,outfile_expmap,outfile_boxlist3,detmask,emin_ev[index],emax_ev[index], - outfile_backmap3,cheese_mask) + 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], outfile_post) + srcmap="{}_SourceMap_en{}.forced{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) + """ give mllist from en0 as input """ + boxlist3="{}_MaxLikSourceList_en{}{}".format(os.path.join(outfile_dir,outkey), eband[0], outfile_post) + fitpos_flag="fitpos_flag=no" + fitext_flag="fitext_flag=no" + nmulsou = "nmulsou=1" + """ don't allow ermpldet to split sources """ + if(index == 3): + """ 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" + """ allow ermpldet to split sources (no more than two) """ -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), + "boxlist=%s" %(boxlist3), "images=\'{}\'".format(outfile_evtool), "expimages=\'{}\'".format(outfile_expmap), "detmasks=\'{}\'".format(detmask), @@ -250,7 +271,7 @@ cmd=["ermldet", "emax=\'{}\'".format(emax_ev[index]), "ecf=\'{}\'".format(ecf[index]), "hrdef=", - "likemin=5.", + "likemin=0.", "extlikemin=6.", "compress_flag=N", "cutrad=10.", # was 15 @@ -267,14 +288,13 @@ cmd=["ermldet", "thres_col=like", "thres_val=30.", "nmaxfit=4", - "nmulsou=2", - "fitext_flag=yes", + nmulsou, + fitpos_flag, + fitext_flag, "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, @@ -283,15 +303,32 @@ if(do_ersensmap==True): 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): - remove_file(mllist) - remove_file(srcmap) - os.system((" ").join(cmd)) - print((" ").join(cmd)) + if(vign==False): + print('Run ermldet with vignetted exposure!') + sys.exit() + remove_file(mllist) + remove_file(srcmap) + os.system((" ").join(cmd)) + print((" ").join(cmd)) + if(forced==True): + check_ermldet_forced(mllist) + +if(forced==True): + catprep="{}_SourceCatalog_en{}.forced{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) +else: + catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) -catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) if(do_catprep==True): cmd=["catprep", "infile={}".format(mllist), @@ -306,3 +343,6 @@ if(do_filter_catalog==True): 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") diff --git a/scripts/05_srctool.py b/scripts/05_srctool.py index 1a841c7..9560ba9 100755 --- a/scripts/05_srctool.py +++ b/scripts/05_srctool.py @@ -2,32 +2,31 @@ """ НАЗВАНИЕ: - 04_mosaics.py + 05_srctool.py НАЗНАЧЕНИЕ: - Собирает мозайки в разных энергетических диапазонах. - + Запускает scrtool для самого широкого канала 0.2-10 кэВ, чтобы спектры имели самое полное покрытие по энергиям. Список источников берется из 0.3-2.3 кэВ. + ВЫЗОВ: esass - ./01_mosaics.py + ./05_srctool.py УПРАВЛЕНИЕ: - Запуск отдельных команд управляется переменными, например: do_init = True - Выбранный энергетический диапазон управляется переменной index + Требуется запуск предыдущего скрипта 04_mosaics.py ПАРАМЕТРЫ: - index : Выбранный энергетический диапазон + index=4 : Выбранный энергетический диапазон ВЫВОД: - Выходные файлы записываются в директорию outfile_dir + Выходные файлы записываются в директорию outfile_dir/srctool_dir ИСТОРИЯ: @@ -49,6 +48,7 @@ import uds from uds.utils import * from uds.config import * +from uds.sherpa import * """ find UDS root dir """ @@ -73,10 +73,13 @@ outfile_srctool="{}_SrcTool_".format(outkey) do_init = False do_merge = False -do_srctool = True -do_grppha = True +do_srctool = False +do_grppha = False +do_ecf_calc = False +do_ecf_print = False +do_catalog = True -index=4 +index=1 """ работаем именно в этом диапазоне, чтобы спектры покрывали все энергии """ vign=True @@ -122,6 +125,7 @@ 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() @@ -146,5 +150,19 @@ if(do_srctool==True): if(do_grppha==True): group_spectra("{}/*_SourceSpec_*.fits".format(srctool_dir)) +ecfout="{}_SampleFlux.pickle".format(os.path.join(outfile_dir,outkey)) +if(do_ecf_calc==True): + calc_ecf("{}/tm0_SrcTool_020_ARF_?????.fits".format(srctool_dir), + catprep=catprep, emin=emin_kev, emax=emax_kev, eband=eband, outfile=ecfout, simnum=10000) +if(do_ecf_print==True): + print_ecf(infile=ecfout, emin=emin_kev, emax=emax_kev, eband=eband, skipfrac=10.0) + +index=0 +catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) +if(do_catalog==True): + make_catalog(infile=catprep, dlmin=10.0, dlmax=100000, ext_like=10, ecf=ecf[index], + emin=emin_kev[index], emax=emax_kev[index], eband=eband[index], + infile_en03cat='../products/tm0_SourceCatalog_en3.forced.fits', + infile_en03sens='../products/tm0_SensitivityMap_dl10_en3.fits') diff --git a/scripts/06_plot.py b/scripts/06_plot.py new file mode 100755 index 0000000..0c93158 --- /dev/null +++ b/scripts/06_plot.py @@ -0,0 +1,189 @@ +#!/usr/bin/env python +""" +НАЗВАНИЕ: + + 05_srctool.py + + +НАЗНАЧЕНИЕ: + + Запускает scrtool для самого широкого канала 0.2-10 кэВ, чтобы спектры имели самое полное покрытие по энергиям. Список источников берется из 0.3-2.3 кэВ. + +ВЫЗОВ: + + esass + ./05_srctool.py + + +УПРАВЛЕНИЕ: + + Требуется запуск предыдущего скрипта 04_mosaics.py + +ПАРАМЕТРЫ: + + index=4 : Выбранный энергетический диапазон + + +ВЫВОД: + + Выходные файлы записываются в директорию outfile_dir/srctool_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 + +import matplotlib.pyplot as plt + + + +from uds.utils import * +from uds.config import * +from uds.sherpa 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_flux_distr = False +do_sens_curve = True + +index=0 + +catalog = "{}_SourceCatalog_en{}.main.selected.csv".format(os.path.join(outfile_dir,outkey), eband[0]) + +if not (os.path.isfile(catalog)==True): + print("{} not found, run 05_srctool.py?".format(catalog)) + sys.exit() + +if(do_flux_distr==True): + + data, logbins = get_log_distr(infile=catalog, field='ml_rate', minval=1e-3, maxval=2) + fig, ax = plt.subplots() + for axis in ['top','bottom','left','right']: + ax.spines[axis].set_linewidth(1) + ax.tick_params(axis="both", width=1, labelsize=14) + plt.hist(data, bins=logbins, histtype='step', color='blue', linewidth=1, linestyle='solid') + plt.xlabel('Count rate (counts s$^{-1}$)',fontsize=14, fontweight='normal') + plt.ylabel('Number',fontsize=14, fontweight='normal') + plt.grid(visible=True) + plt.xscale('log') + plt.yscale('log') + plt.savefig(catalog.replace("main.selected.csv", "ml_rate.png"), bbox_inches='tight') + plt.close(fig) + + + data, logbins = get_log_distr(infile=catalog, field='ml_flux', minval=1e-15, maxval=2e-12) + fig, ax = plt.subplots() + for axis in ['top','bottom','left','right']: + ax.spines[axis].set_linewidth(1) + ax.tick_params(axis="both", width=1, labelsize=14) + plt.hist(data, bins=logbins, histtype='step', color='blue', linewidth=1, linestyle='solid') + plt.xlabel('Energy flux (erg s$^{-1}$ cm$^{-2}$)',fontsize=14, fontweight='normal') + plt.ylabel('Number',fontsize=14, fontweight='normal') + plt.grid(visible=True) + plt.xscale('log') + plt.yscale('log') + plt.savefig(catalog.replace("main.selected.csv", "ml_flux.png"), bbox_inches='tight') + plt.close(fig) + +if(do_sens_curve==True): + coeff = 2.4336e-13 / 3.4012e-13 # see below + + areatab="{}_AreaTable_dl10_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) + hdul = fits.open(areatab) + tbdata = hdul[1].data + """ convert limflux from 0.3-2.3 keV to 0.5-2 keV """ + limflux_dl10 = tbdata['LIMFLUX']*coeff + area_dl10 = tbdata['SKY_AREA'] + hdul.close() + + areatab="{}_AreaTable_dl6_en{}{}".format(os.path.join(outfile_dir,outkey), eband[index], outfile_post) + hdul = fits.open(areatab) + tbdata = hdul[1].data + """ convert limflux from 0.3-2.3 keV to 0.5-2 keV """ + limflux_dl6 = tbdata['LIMFLUX']*coeff + area_dl6 = tbdata['SKY_AREA'] + hdul.close() + + fig, ax = plt.subplots() + for axis in ['top','bottom','left','right']: + ax.spines[axis].set_linewidth(1) + ax.tick_params(axis="both", width=1, labelsize=14) + ax.set_xlim([3e-16,1e-13]) + ax.set_ylim([0.3,160]) + + df = pandas.read_csv("../data/surveys/eFEDS.dat",header=None,names=['limflux','area']) + plt.plot(df['limflux'],df['area'], color="brown", linestyle='solid',label="eFEDS (DL6)") + + plt.plot(limflux_dl10,area_dl10, color="black", linestyle='solid',label="eUDS (DL10)") + plt.plot(limflux_dl6,area_dl6, color="black", linestyle='dashed', label="eUDS (DL6)") + + df = pandas.read_csv("../data/surveys/cosmos-legacy.dat",header=None,names=['limflux','area']) + plt.plot(df['limflux'],df['area'], color="blue", linestyle='solid',label="COSMOS Legacy") + + df = pandas.read_csv("../data/surveys/CDWFS.dat",header=None,names=['limflux','area']) + plt.plot(df['limflux'],df['area'], color="green", linestyle='solid', label="CDWFS") + + df = pandas.read_csv("../data/surveys/XMM-RM.dat",header=None,names=['limflux','area']) + plt.plot(df['limflux'],df['area'], color="magenta", linestyle='solid',label="XMM-RM") + + df = pandas.read_csv("../data/surveys/XMM-XXL-N.dat",header=None,names=['limflux','area']) + plt.plot(df['limflux'],df['area'], color="red", linestyle='solid',label="XMM-XXL-N") + + + ax.legend() + + #plt.hist(data, bins=logbins, histtype='step', color='blue', linewidth=1, linestyle='solid') + plt.xlabel('Limiting flux (0.5-2 keV, erg s$^{-1}$ cm$^{-2}$)',fontsize=14, fontweight='normal') + plt.ylabel('Area (deg$^{2}$)',fontsize=14, fontweight='normal') + plt.grid(visible=True) + plt.xscale('log') + plt.yscale('log') + png="{}_AreaTable_en{}.png".format(os.path.join(outfile_dir,outkey), eband[index]) + plt.savefig(png, bbox_inches='tight') + plt.close(fig) + + """ + ======================================================================== + Model TBabs<1>*powerlaw<2> Source No.: 1 Active/On + Model Model Component Parameter Unit Value + par comp + 1 1 TBabs nH 10^22 2.00000E-02 frozen + 2 2 powerlaw PhoIndex 2.00000 frozen + 3 2 powerlaw norm 1.14851E-04 +/- 3.83988E-06 + ________________________________________________________________________ + + Model Flux 0.00028334 photons (3.4012e-13 ergs/cm^2/s) range (0.30000 - 2.3000 keV) + Model Flux 0.0001619 photons (2.4336e-13 ergs/cm^2/s) range (0.50000 - 2.0000 keV) + """ diff --git a/scripts/README.md b/scripts/README.md index 7b9cbb2..7d7e9d0 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -33,4 +33,8 @@ source /eSASS4EDR/bin/esass-init.csh ### 04_mosaics.py -Создает сборные изображения (мозайки) в разных энергетических диапазонах. \ No newline at end of file +Создает сборные изображения (мозайки) в разных энергетических диапазонах. + +### 05_scrtool.py + +Запускает scrtool для самого широкого канала 0.2-10 кэВ, чтобы спектры имели самое полное покрытие по энергиям. Список источников берется из 0.3-2.3 кэВ. diff --git a/uds/uds/config.py b/uds/uds/config.py index 1859eb3..14df81a 100644 --- a/uds/uds/config.py +++ b/uds/uds/config.py @@ -51,7 +51,25 @@ 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] +#ecf = [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] +ecf = [9.7817E+11, 3.2982E+12, 1.3903E+12, 2.3322E+12, 5.2022E+11, 5.8453E+11] + +""" +*** en0 ecf 9.7817E+11 +/- 2.4606E+10 2.52% N=17 +*** en1 ecf 3.2982E+12 +/- 8.2963E+10 2.52% N=17 +*** en2 ecf 1.3903E+12 +/- 3.5036E+10 2.52% N=17 +*** en3 ecf 2.3322E+12 +/- 5.8717E+10 2.52% N=17 +*** en4 ecf 5.2022E+11 +/- 1.3110E+10 2.52% N=17 +*** en5 ecf 5.8453E+11 +/- 1.4743E+10 2.52% N=17 +""" outfile_post='.fits' + +""" +Pavel Medvedev: +0.3-2.3: 9.135819435325375e-13 +0.3-0.6: 9.160477830652834e-13 +0.6-2.3: 9.201664167869427e-13 +2.3-5.0: 8.721504826794627e-12 +""" diff --git a/uds/uds/sherpa.py b/uds/uds/sherpa.py new file mode 100644 index 0000000..b805a75 --- /dev/null +++ b/uds/uds/sherpa.py @@ -0,0 +1,219 @@ +import glob +import os +import sys +import numpy as np +from astropy.io import fits +import pickle +from scipy.stats import sem + +#import matplotlib.pyplot as plt +#plt.style.use('ggplot') + +from sherpa.astro.ui import * + +def calc_ecf(path,catprep=None, emin=None, emax=None, eband=None, outfile=None, fitmin=0.2, fitmax=5.0, simnum=100): + """ + Do uncorrection of SPECRESP for CORRCOMP=CORRPSF*CORRVIGN, like this: + SPECRESP=SPECRESP/CORRCOMB, because because ermldet does not correct for PSF and vignetting. + """ + + data={} + + if (catprep == None): + print("catprep file is not provided, exit") + sys.exit() + else: + print("Reading {}".format(catprep)) + + catprep_hdul = fits.open(catprep,mode='readonly') + catprep_data = catprep_hdul[1].data + #print(catprep_data['detUID']) + mlrate = {} + mlrate_err = {} + dl = {} + for rec in catprep_data: + srcid=rec['detUID'] + + """ consider only non-extended sources """ + if(rec['EXT'] > 0.0): + continue + + if(srcid in mlrate.keys()): + print("Duplicate SrcID! {}".format(srcid)) + sys.exit() + else: + mlrate[srcid]=rec['ML_RATE_0'] + mlrate_err[srcid]=rec['ML_RATE_ERR_0'] + dl[srcid]=rec['DET_LIKE_0'] + + + flist=glob.glob(path) + for filename in flist: + + print("Reading {}".format(filename)) + + fout=filename.replace(".fits", ".uncorr.fits") + + pha=filename.replace("ARF", "SourceSpec") + pha=pha.replace(".fits", ".fits.grp") + pha_header = fits.getheader(pha,ext=1) + object = pha_header['OBJECT'] + + """ check a given source (testing...) """ + #if (int(object) > 10): + # continue + + srcid="ML{}".format(object) + if not (srcid in mlrate.keys()): + print("Skip {}".format(srcid)) + continue + + + hdul = fits.open(filename,mode='readonly') + specresp = np.divide(hdul[1].data['SPECRESP'], hdul[1].data['CORRCOMB']) + hdul[1].data['SPECRESP'] = specresp + hdul.writeto(fout, overwrite=True) + clean() # sherpa + set_xsxsect("vern") + set_xsabund('wilm') + load_pha(pha) # sherpa + + """ load uncorrected ARF """ + # load_arf(fout) # sherpa + + # Define the source model + #sh.set_source("powlaw1d.p1") + set_source("xstbabs.abs1 * powlaw1d.p1") # sherpa + abs1.nH = 0.02 + # Change reference energy of the model + p1.ref = 1 # 1 keV + p1.gamma = 2.0 + p1.ampl = 1e-4 # in cm**-2 s**-1 keV**-1 + freeze(abs1.nH, p1.gamma) # sherpa + + ### Define the statistic + set_stat("wstat") # sherpa + + ignore_bad() # sherpa + ### Define the fit range + print("Notice 0.3-5.0 keV:") + notice(fitmin, fitmax) # sherpa + + """ Check wether spectrum still has counts after filtering """ + dep = get_dep(filter=True) + if not (len(dep)): + print("*** Skip {} due to low counts".format(srcid)) + continue + + + + ### Do the fit + fit() + res_fit=get_fit_results() + rstat=res_fit.rstat + + set_analysis('energy') + + photon_flx=[] + energy_flx=[] + skipme=False + for ie in range(len(eband)): + ph_flx=calc_photon_flux(emin[ie], emax[ie]) + en_flx=calc_energy_flux(emin[ie], emax[ie]) + if not (ph_flx>0): + skipme=True + photon_flx.append(ph_flx) + energy_flx.append(en_flx) + print("en{} Photon flux: {:.4E} Energy flux: {:.4E}".format(eband[ie], ph_flx, en_flx)) + + if (skipme == True): + print("*** Skip zero flux for {}".format(srcid)) + continue + + print("----> Success: {} pow norm: {:.4E} reduced stat: {:.2f}".format(res_fit.succeeded,res_fit.parvals[0], rstat)) + + sample_flx=[] + sample_flx_lo=[] + sample_flx_hi=[] + if (res_fit.rstat < 3.0): + conf() + res_conf=get_conf_results() + #print(res_conf) + component=abs1+p1 + #print(component) + #covar() + #errs = get_covar_results().parmaxes + #ans = sample_flux(correlated=False, scales=errs, num=500) + for ie in range(len(eband)): + pars_flux = sample_flux(modelcomponent=component, lo=emin[ie], hi=emax[ie], num=simnum, correlated=True, numcores=10, confidence=68) + sample_flx.append(pars_flux[0][0]) + sample_flx_lo.append(pars_flux[0][2]) + sample_flx_hi.append(pars_flux[0][1]) + + else: + continue + + print("### SAMPLE FLUX: {}".format(sample_flx)) + data[srcid]={'sample_flx':sample_flx, + 'sample_flx_lo':sample_flx_lo, + 'sample_flx_hi':sample_flx_hi, + 'photon_flx':photon_flx, + 'energy_flx':energy_flx, + 'ml_rate':mlrate[srcid], + 'ml_rate_err':mlrate_err[srcid], + 'det_like':dl[srcid]} + + with open(outfile, 'wb') as f: + pickle.dump(data, f) + + +def print_ecf(infile=None, emin=None, emax=None, eband=None, skipfrac=5.0): + + + + with open(infile, 'rb') as f: + data = pickle.load(f) + + for ie in range(len(eband)): + print() + ecf=[] + + for key in data.keys(): + ar=data[key] + ml = ar['ml_rate'] + dml = ar['ml_rate_err'] + + flx = ar['sample_flx'][ie] + dflx1 = ar['sample_flx'][ie] - ar['sample_flx_lo'][ie] + dflx2 = ar['sample_flx_hi'][ie] - ar['sample_flx'][ie] + + # take maximum error as conservative approach + dflx = np.maximum(dflx1,dflx2) + + ecf0 = ml/flx + decf0 = abs(flx*dml - ml*dflx) / np.power(flx,2) + frac0 = 100 * abs(flx*dml - ml*dflx) / (flx*ml) + + skipme=False + if((100*dflx/flx) < skipfrac): + skipme=False + else: + skipme=True + + #if not (len(ar['sample_flx'])>0): + # continue + print("en{} {} DL={:.2f} FLUX: {:4E} (-{:.4E} +{:.4E} {:.2f}%) ML_RATE: {:.4f} +/- {:.4f} ({:.2f}%) --> {:.4E} +/- {:.4E} {:.2f}% skip={}".format(eband[ie],key,ar['det_like'], + flx,dflx1,dflx2,100*dflx/flx, + ml,dml,100*dml/ml, + ecf0,decf0,frac0,skipme)) + if(skipme==True): + continue + ecf.append(ecf0) + + + ecf_mean = np.mean(ecf) + ecf_err = np.std(ecf, ddof=1)/np.sqrt(np.size(ecf)) + print("\n*** en{} ecf {:.4E} +/- {:.4E} {:.2f}% N={}".format(eband[ie], + ecf_mean, + ecf_err,100*ecf_err/ecf_mean, + np.size(ecf))) diff --git a/uds/uds/utils.py b/uds/uds/utils.py index 44a9bd5..7386bc6 100644 --- a/uds/uds/utils.py +++ b/uds/uds/utils.py @@ -1,5 +1,6 @@ import os import sys +import csv import numpy as np from astropy.io import fits from astropy import wcs @@ -9,7 +10,7 @@ from astropy.io.fits import getdata import glob from fitsio import FITS from pathlib import Path - +import pandas from astropy.table import QTable, Table, Column from astropy import units as u @@ -983,7 +984,7 @@ def create_sensmap(sensmap=None,expmap=None,backmap=None,detmask=None,areatab=No "ext_flag=no", "extlike_flag=no", "compress_flag=no", - "area_flag=no",] + "area_flag=yes",] remove_file(sensmap) remove_file(areatab) os.system((" ").join(cmd)) @@ -1005,3 +1006,157 @@ def group_spectra(path,group_min=25): ] print((" ").join(cmd)) os.system((" ").join(cmd)) + +def check_ermldet_forced(infile): + + if(os.path.isfile(infile)==False): + print("File not found {}".format(infile)) + + print("reading mllist {}".format(infile)) + + hdul = fits.open(infile) + tbdata = hdul[1].data + hdul.close() + + for rec in tbdata: + idsrc=rec['ID_SRC'] + boxid = rec['BOX_ID_SRC'] + if(idsrc != boxid): + print("The ermldet catalog in forced mode should contain only unique sources.") + print("Use fitpos_flag=no, fitext_flag=no, nmulsou=1") + print("{} != {}".format(idsrc,boxid)) + sys.exit() + +def make_catalog(infile=None, dlmin=6.0,dlmax=10,scale=60*60,ext_like=0.0,outkey='main',emin=None, emax=None, eband=None, ecf=None, + infile_en03cat=None,infile_en03sens=None): + if(os.path.isfile(infile)==False): + print("File not found {}".format(infile)) + + print("reading catprep {}".format(infile)) + fout_selected=infile.replace(".fits", ".{}.selected.csv".format(outkey)) + fout_skip=infile.replace(".fits", ".{}.skip.reg".format(outkey)) + fout_extended=infile.replace(".fits", ".extended.reg") + + """ Read en03 energy band """ + en03cat={} + hdul = fits.open(infile_en03cat) + en03tab = hdul[1].data + hdul.close() + + for rec in en03tab: + key = rec['detUID'] + en03cat[key]={'ra':rec['RA'], 'dec':rec['DEC'], 'det_like':rec['DET_LIKE_0'], + 'ml_rate':rec['ML_RATE_0'], 'ml_rate_err':rec['ML_RATE_ERR_0'], + 'ml_flux':rec['ML_FLUX_0'], 'ml_flux_err':rec['ML_FLUX_ERR_0'],'sens':0.0} + + print("reading {}".format(infile_en03sens)) + hdul = fits.open(infile_en03sens) + en03sens = hdul[0].data + en03hdr = hdul[0].header + hdul.close() + en03wcs = WCS(en03hdr) + for key in en03cat.keys(): + ra=en03cat[key]['ra'] + dec=en03cat[key]['dec'] + crd = SkyCoord(ra, dec, frame=FK5(), unit="deg") + pix = en03wcs.wcs_world2pix([(ra, dec),], 1) + px=round(pix[0][0])-1 + py=round(pix[0][1])-1 + en03cat[key]['sens']=en03sens[py,px] + #print(key,ra,dec,px,py,value) + #(ra0,dec0), = en03wcs.wcs_pix2world(pix, 1) + #crd0 = SkyCoord(ra0, dec0, frame=FK5(), unit="deg") + #sep = crd.separation(crd0).arcsec + + hdul = fits.open(infile) + tbdata = hdul[1].data + hdul.close() + + catsel=[] + catskip=[] + catext=[] + skip_count=0 + selected_count=0 + keepalive_count=0 + for rec in tbdata: + src_crd = SkyCoord(rec['ra'], rec['dec'], frame=FK5(), unit="deg") + key = rec['detUID'] + if not (key in en03cat.keys()): + print("{} not found in {}".format(rec['detUID'],infile_en03cat)) + sys.exit() + else: + """ check coords """ + en03_crd = SkyCoord(en03cat[key]['ra'], en03cat[key]['dec'], frame=FK5(), unit="deg") + sep = src_crd.separation(en03_crd).arcsec + if(sep > 1.0): + print("{} significant offset: {.4f} arcsec".format(key,sep)) + sys.exit() + pass + + 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 and rec['ext_like'] < ext_like)): + + catsel.append({'ra':rec['ra'], + 'dec':rec['dec'], + 'radec_err':rec['radec_err'], + 'det_like':rec['det_like_0'], + 'src_id':rec['id_src'], + 'ml_rate':rec['ml_rate_0'], + 'ml_rate_err':rec['ml_rate_err_0'], + 'ml_flux':rec['ml_flux_0']/ecf, + 'ml_flux_err':rec['ml_flux_err_0']/ecf, + 'en03_dl':en03cat[key]['det_like'], + 'en03_flux':en03cat[key]['ml_flux'], + 'en03_flux_err':en03cat[key]['ml_flux_err'], + 'en03_sens':en03cat[key]['sens'], + }, + ) + + 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 csvfile: + fieldnames = ['src_id', 'ra', 'dec', 'radec_err', 'det_like','ml_rate','ml_rate_err','ml_flux','ml_flux_err','en03_dl','en03_flux','en03_flux_err','en03_sens'] + writer = csv.DictWriter(csvfile, fieldnames=fieldnames) + writer.writeheader() + for rec in catsel: + writer.writerow(rec) +# writer.writerow({'src_id':rec['src_id'], +# 'ra':rec['ra'], +# 'dec':rec['dec'], +# 'radec_err':rec['radec_err'], +# 'det_like':rec['det_like']}) + return + 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 get_log_distr(infile=None, index_col='src_id', field=None, minval=10, maxval=100, nbin=50): + logbins=np.logspace(np.log10(minval),np.log10(maxval), nbin) + df = pandas.read_csv(infile, index_col=index_col) + data = df[field] + #hist, edges = np.histogram(array, bins=logbins) + return data, logbins + + + + +