erosita/coma
1
0
Fork 0
forked from erosita/uds
Code Pull Requests Activity

Many changes

Browse Source
This commit is contained in:
Roman Krivonos 2024-03-18 14:46:38 +03:00
parent 4e5f88e166
commit afe3862cda
24 changed files with 55736 additions and 236 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
.gitignore vendored
View File

@ -16,6 +16,7 @@ lib64/
parts/ parts/
sdist/ sdist/
var/ var/
products/
wheels/ wheels/
share/python-wheels/ share/python-wheels/
*.egg-info/ *.egg-info/
@ -43,9 +44,15 @@ venv.bak/
*# *#
.#* .#*
# Astronomy # Astronomy
*.eps
*.cds
*.fits *.fits
*.fits.gz *.fits.gz
*.fits.fz
*.awk
*.tex
*.evt *.evt
*.evt.gz *.evt.gz
*.lock *.lock

52767
data/4XMM-DR12/4XMM_DR12cat_v1.0_UDS.fits.catalog Normal file
View File

File diff suppressed because one or more lines are too long

5
data/4XMM-DR12/README.md
View File

@ -6,5 +6,10 @@
``` ```
ftselect '4XMM_DR12cat_slim_v1.0.fits[1]' 4XMM_DR12cat_slim_v1.0_UDS.fits 'SC_RA > 32.75 && SC_RA < 36.31 && SC_DEC > -6.55 && SC_DEC < -3.0' clobber=yes ftselect '4XMM_DR12cat_slim_v1.0.fits[1]' 4XMM_DR12cat_slim_v1.0_UDS.fits 'SC_RA > 32.75 && SC_RA < 36.31 && SC_DEC > -6.55 && SC_DEC < -3.0' clobber=yes
ftselect '4XMM_DR12cat_v1.0.fits[1]' 4XMM_DR12cat_v1.0_UDS.fits 'SC_RA > 32.75 && SC_RA < 36.31 && SC_DEC > -6.55 && SC_DEC < -3.0' clobber=yes
Demo with DR13:
ftselect '4XMM_slim_DR13cat_v1.0.fits[1]' 4XMM_slim_DR13cat_v1.0_DEMO.fits 'SC_RA > 25.75 && SC_RA < 50.31 && SC_DEC > -16.55 && SC_DEC < 13.0' clobber=yes
``` ```

BIN
data/4XMM-DR12/Rplots.pdf Normal file
View File

Binary file not shown.

20
data/4XMM-DR12/barplot.r Normal file
View File

@ -0,0 +1,20 @@
library(ggplot2)
# see https://r-graph-gallery.com/4-barplot-with-error-bar.html
png(file = "barplot.png")
# create dummy data
data <- data.frame(
name=c('','yy','zz','ff','DR12'),
flux=c(1,2,3,4,5),
err=c(1,0.2,3,2,4)
)
# Most basic error bar
ggplot(data) +
geom_bar(aes(x=name, y=flux), stat="identity", fill="skyblue", alpha=0.7) +
geom_errorbar( aes(x=name, ymin=flux-err, ymax=flux+err), width=0.4, colour="black", alpha=0.9, linewidth=2)
dev.off()

40
data/4XMM-DR12/plot_corr.r Normal file
View File

@ -0,0 +1,40 @@
options(digits = 7) # knitr sometimes decreases number of digits shown
require(stats)
xgrid <- c(1e-17,1e-16,1e-15,1e-14,1e-13,1e-12);
postscript(paste("../../products/eUDS_4XMM-DR12.flux.eps",sep=""), horizontal = FALSE, onefile = FALSE, paper = "special",width = 9.0, height = 6.0)
par(oma=c(1,1,1,1))
#layout(matrix(c(1,2), 2, 1, byrow = TRUE), heights = c(1.5,1), TRUE)
par(mar=c(4.3, 4.3, 1, 1))
par(cex.lab=1.2)
par(cex.axis=1.2)
a <- read.table("../../products/eUDS_4XMM-DR12.flux.cvs", col.names=c("eflux","eerr","xflux","xerr"))
xlim <- c(1e-15,1e-12)
ylim <- c(1e-16,1e-12)
plot(NULL, pch=21, bg="gold",ylim=ylim, xlim=xlim, main="", ylab="eUDS, erg/s/cm2",type="p",xaxt="n",xlab="4XMM-DR12, erg/s/cm2", log="xy")
axis(side=1, at=xgrid, labels=xgrid)
segments(a$xflux,a$eflux-a$eerr,a$xflux,a$eflux+a$eerr)
segments(a$xflux-a$xerr,a$eflux,a$xflux+a$xerr,a$eflux)
#segments(a$xlo,d,a$xhi,d)
points(a$xflux,a$eflux,pch=21, bg="gold")
lines(c(1e-17,1e-12),c(1e-17,1e-12))
lines(c(1e-17,1e-12),10*c(1e-17,1e-12),lt=2)
lines(c(1e-17,1e-12),0.1*c(1e-17,1e-12),lt=2)
grid()
abline(h=0, col = "black",lty=2)
abline(v=xgrid, col="lightgray", lty="dotted")
dev.off()

14
data/4XMM-DR12/print_ds9reg.py
View File

@ -1,15 +1,19 @@
from astropy.io import fits from astropy.io import fits
import sys import sys
filename='4XMM_DR12cat_slim_v1.0_UDS.fits.catalog' #filename='4XMM_DR12cat_slim_v1.0_UDS.fits.catalog'
filename='4XMM_slim_DR13cat_v1.0_DEMO.fits'
fout=filename.replace(".fits.catalog", ".names.reg") fout=filename.replace(".fits.catalog", ".names.reg")
hdul = fits.open(filename) hdul = fits.open(filename)
#hdul.info() #hdul.info()
tbdata = hdul[1].data tbdata = hdul[1].data
hdul.close()
#with open("./{}".format(fout), 'w') as writer:
for rec in tbdata:
print("fk5;circle({}, {}, 0.008)".format(rec['sc_ra'],rec['sc_dec']))
#writer.write("fk5;circle({}, {}, {}) # text={{{}}}\n".format(rec['sc_ra'],rec['sc_dec'],rec['sc_poserr']/3600,rec['IAUNAME']))
with open("../../products/{}".format(fout), 'w') as writer:
for rec in tbdata:
#writer.write("fk5;point({}, {})\n".format(rec['sc_ra'],rec['sc_dec']))
writer.write("fk5;circle({}, {}, {}) # text={{{}}}\n".format(rec['sc_ra'],rec['sc_dec'],rec['sc_poserr']/3600,rec['IAUNAME']))

BIN
data/ECF/ecf_tbabspow_g2nh0.02.pkl Normal file
View File

Binary file not shown.

BIN
data/ECF/ecf_tbabspow_g2nh003.pkl Normal file
View File

Binary file not shown.

27
data/NuSTAR/nustar.fits.catalog Normal file
View File

File diff suppressed because one or more lines are too long

13
data/NuSTAR/print_ds9reg.py Executable file
View File

@ -0,0 +1,13 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from astropy.io import fits
hdul = fits.open('nustar.fits.catalog')
tbdata = hdul[1].data
hdul.close()
with open("../../products/nustar_UDS.reg", 'w') as writer:
for rec in tbdata:
writer.write("fk5;circle({}, {}, 0.0025) # color=magenta text={{{}}}\n".format(rec['RAJ2000'], rec['DEJ2000'], rec['NuSTAR']))

11
products/srctool-products/ecf.xcm
View File

@ -1,11 +0,0 @@
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

13
products/srctool-products/rate.pro
View File

@ -1,13 +0,0 @@
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

12
scripts/02_merge_events.py
View File

@ -35,14 +35,14 @@ outfile_dir=root_path+'/products'
create_folder(outfile_dir) create_folder(outfile_dir)
index=5 # select energy band index=3 # select energy band
do_init = False do_init = True
do_merge = True do_merge = True
do_rate = False do_rate = False
do_adapt = False # requires CIAO do_adapt = True # requires CIAO
vign=True vign=False
vignetting = 'vign' if (vign==True) else 'novign' vignetting = 'vign' if (vign==True) else 'novign'
events=[] events=[]
@ -61,8 +61,8 @@ for tmkey in keylist_tm.keys():
do_init=do_init, do_init=do_init,
do_obsmode=False, do_obsmode=False,
do_center=False, do_center=False,
do_evtool=False, do_evtool=True,
do_expmap=True, do_expmap=False,
vign=vign, vign=vign,
ra_cen=ra_cen, de_cen=de_cen, ra_cen=ra_cen, de_cen=de_cen,
emin_kev=emin_kev[index], emin_kev=emin_kev[index],

96
scripts/03_init_obs.py
View File

@ -43,6 +43,8 @@ import sys, os, os.path, time, subprocess
from pathlib import Path from pathlib import Path
import numpy as np import numpy as np
import glob import glob
from multiprocessing import Pool
from os.path import dirname from os.path import dirname
import inspect import inspect
import uds import uds
@ -59,24 +61,29 @@ infile_dir=root_path+'/data/processed'
outfile_dir=root_path+'/products' outfile_dir=root_path+'/products'
create_folder(outfile_dir) create_folder(outfile_dir)
do_init = True run_Pool=False
do_ermask = True
do_erbox1 = True # local mode do_init = False
do_erbackmap1 = True # do_ermask = False
do_erbox2 = True # map mode, with background map
do_erbackmap2 = True #
do_erbox3 = True # map mode, with background map
do_erbackmap3 = True #
do_ermldet = True do_erbox1 = False # local mode
do_catprep = True 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_ermldet = False
do_catprep = False
do_cross_match = False do_cross_match = False
do_wcs_match = True # Chandra task do_astro_corr = False # search optimal shift
do_wcs_update = True # Chandra task do_astro_update = True
eband_selected=[0,1,2,3,4] do_wcs_match = False # Chandra task -- DEPRECATED
do_wcs_update = False # Chandra task -- DEPRECATED
eband_selected=[5]
vign=True vign=True
vignetting = 'vign' if (vign==True) else 'novign' vignetting = 'vign' if (vign==True) else 'novign'
@ -95,6 +102,11 @@ def runme(datakey):
cheesemask=[] cheesemask=[]
bkgimage=[] bkgimage=[]
srcmaps=[] srcmaps=[]
print(datakey)
print('module name:', __name__)
print('parent process:', os.getppid())
print('process id:', os.getpid())
for ii in range(len(eband_selected)): for ii in range(len(eband_selected)):
index=eband_selected[ii] index=eband_selected[ii]
@ -114,10 +126,10 @@ def runme(datakey):
expmaps.append(outfile_expmap) expmaps.append(outfile_expmap)
events.append(outfile_evtool) 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: # 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")) #events.append(outfile_evtool.replace(".fits", ".attcorr.fits"))
"""
""" Detmask """ """ Detmask """
@ -339,12 +351,13 @@ def runme(datakey):
] ]
if(do_ermldet==True): if(do_ermldet==True):
test_exe('ermldet')
remove_file(mllist) remove_file(mllist)
remove_file(srcmap) remove_file(srcmap)
os.system((" ").join(cmd)) os.system((" ").join(cmd))
print((" ").join(cmd)) print((" ").join(cmd))
#save_ermldet_ds9reg(mllist,scale=60*60) save_ermldet_ds9reg(mllist,scale=60*60)
catprep="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,datakey), eband[index], outfile_post) 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) catprep_en0="{}_SourceCatalog_en{}{}".format(os.path.join(outfile_dir,datakey), eband[0], outfile_post)
@ -359,12 +372,24 @@ def runme(datakey):
if(do_cross_match==True): if(do_cross_match==True):
crossmatch_shu2019(catprep,dlmin=10,refimage=events[ii],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") catalog=root_path+"/data/Gaia_unWISE/Gaia_unWISE_UDS.fits.catalog",errlim=5.0)
if(do_wcs_match==True and eband[index]==0): if(do_astro_corr==True and eband[index]=='0'):
""" run astro_corr for 0.3-2.3 keV only """
wcs_astro_corr(catprep)
#wcs_match_ciao(catprep, method='rst',radius=12,residlim=0,residtype=0,residfac=1)
if(do_astro_update==True):
""" run astro_corr for 0.3-2.3 keV only """
attcorr=wcs_update_shift(events[ii],flog=catprep_en0.replace(".fits", ".shift.log"))
do_evtool_esass(evfile=attcorr,outfile=attcorr,rmlock=False, do_center=True, ra_cen=ra_cen, de_cen=de_cen)
if(do_wcs_match==True and eband[index]=='0'):
""" run wcs_match for 0.3-2.3 keV only """ """ run wcs_match for 0.3-2.3 keV only """
wcs_match_ciao(catprep, method='rst',radius=12,residlim=5) wcs_match_ciao(catprep, method='trans',radius=12,residlim=5)
#wcs_match_ciao(catprep, method='rst',radius=12,residlim=0,residtype=0,residfac=1)
if(do_wcs_update==True): if(do_wcs_update==True):
""" use 0.3-2.3 keV transform matrix for all other bands """ """ 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) attcorr=wcs_update_ciao(events[ii],crval=wcslist[datakey],transformfile=catprep_en0.replace(".fits", ".xfm"),clean=False)
@ -372,13 +397,28 @@ def runme(datakey):
""" """
testing # individual run, testing
runme("tm7_obs_1") runme("tm7_obs_1")
runme("tm5_obs_1") runme("tm5_obs_1")
runme("tm6_scan_1")
""" """
for tmkey in keylist_tm.keys():
print("TM{} in work... init events".format(tmkey))
for datakey in keylist_tm[tmkey]: if(run_Pool==True):
print("--> {}".format(datakey)) # parallel run
runme(datakey) items=[]
for tmkey in keylist_tm.keys():
for datakey in keylist_tm[tmkey]:
items.append(datakey)
with Pool() as pool:
pool.map(runme, items)
else:
# conventional run
for tmkey in keylist_tm.keys():
for datakey in keylist_tm[tmkey]:
print("--> {}".format(datakey))
runme(datakey)
#sys.exit()

120
scripts/04_mosaics.py
View File

@ -13,7 +13,7 @@
ВЫЗОВ: ВЫЗОВ:
esass esass
./01_mosaics.py ./04_mosaics.py
УПРАВЛЕНИЕ: УПРАВЛЕНИЕ:
@ -74,37 +74,39 @@ local_run = False
outkey="tm0" outkey="tm0"
do_init = False do_init = True
do_merge = False do_merge = True
do_detmask = False do_detmask = True
do_expmap = False do_expmap = True
do_erbox1 = False # local mode do_erbox1 = True # local mode
do_erbackmap1 = False # do_erbackmap1 = True #
do_erbox2 = False # map mode, with background map do_erbox2 = True # map mode, with background map
do_erbackmap2 = False # do_erbackmap2 = True #
do_erbox3 = False # map mode, with background map do_erbox3 = True # map mode, with background map
do_erbackmap3 = False # do_erbackmap3 = True #
do_ersensmap = False do_ersensmap = False
do_ermldet = False do_ermldet = True
do_fixcat = False # only for index=0 do_fixcat = False # only for index=0
do_fixxmm = True # only for index=0 do_fixxmm = False # prepare forced photometry, only for index=0
do_apetool = False do_apetool = False
do_catprep = False do_catprep = True
do_filter_catalog = False do_filter_catalog = False
do_cross_match = False do_cross_match = False
index=0 index=3
forced=False forced=False
""" If forced=True, take input catalog from energy range en0 """ """ If forced=True, take input catalog from energy range en0 """
comm='-xmm' # for 4XMM-DR12 forced photometry use '-xmm' comm='' # for 4XMM-DR12 forced photometry use '-xmm'
vign=True vign=True
vignetting = 'vign' if (vign==True) else 'novign' vignetting = 'vign' if (vign==True) else 'novign'
attcorr=True
events=[] events=[]
expmaps=[] expmaps=[]
bkgmaps=[] bkgmaps=[]
@ -113,7 +115,7 @@ for tmkey in keylist_tm.keys():
for datakey in keylist_tm[tmkey]: for datakey in keylist_tm[tmkey]:
print("--> {}".format(datakey)) print("--> {}".format(datakey))
""" Подготавливаем списки событий индивидуальных наблюдений """ """ Подготавливаем списки событий индивидуальных наблюдений """
outfile_evtool,outfile_expmap=init_events(key=datakey, outfile_evtool,outfile_expmap=init_events(key=datakey,attcorr=attcorr,
eband_index=eband[index], eband_index=eband[index],
infile_dir=infile_dir, infile_dir=infile_dir,
outfile_dir=outfile_dir, outfile_dir=outfile_dir,
@ -121,7 +123,7 @@ for tmkey in keylist_tm.keys():
do_obsmode=False, do_obsmode=False,
do_center=False, do_center=False,
do_evtool=True, do_evtool=True,
do_expmap=True, do_expmap=False,
vign=vign, vign=vign,
ra_cen=ra_cen, de_cen=de_cen, ra_cen=ra_cen, de_cen=de_cen,
emin_kev=emin_kev[index], emin_kev=emin_kev[index],
@ -138,7 +140,7 @@ if(do_merge==True):
do_evtool_esass(events=events, outfile=outfile_evtool) do_evtool_esass(events=events, outfile=outfile_evtool)
""" makes detmask from TM exposures """ """ makes detmask from TM exposures """
detmask="{}/{}_DetectorMask_en{}{}".format(outfile_dir, detmask="{}/{}_DetectorMask_en{}{}".format(outfile_dir,
outkey, outkey,
eband[index], eband[index],
@ -256,19 +258,20 @@ if(forced==True):
srcmap="{}_SourceMap_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 """ """ 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) #boxlist3="{}_MaxLikSourceList_en{}.forced{}{}".format(os.path.join(outfile_dir,outkey), eband[0], comm, outfile_post)
if(index==0): #if(index==0):
boxlist3="{}_MaxLikSourceList_en{}.fixed{}{}".format(os.path.join(outfile_dir,outkey), eband[0], comm, outfile_post) boxlist3="{}_MaxLikSourceList_en{}.fixed{}{}".format(os.path.join(outfile_dir,outkey), eband[0], comm, outfile_post)
if not (os.path.exists(boxlist3)): if not (os.path.exists(boxlist3)):
print("{} not found. Run do_fixcat=True, index=0, forced=False".format(boxlist3)) print("{} not found. Run do_fixcat=True, index=0, forced=False".format(boxlist3))
sys.exit() sys.exit()
add_specific_columns(boxlist3) add_specific_columns(boxlist3)
fitpos_flag="fitpos_flag=no" fitpos_flag="fitpos_flag=no"
fitext_flag="fitext_flag=no" fitext_flag="fitext_flag=no"
nmulsou = "nmulsou=1" nmulsou = "nmulsou=1"
nmaxfit="nmaxfit=1" nmaxfit="nmaxfit=10"
""" don't allow ermpldet to split sources """ multrad="multrad=15."
cutrad="cutrad=15."
if(index == 3 or index == 6): if(index == 3 or index == 6):
""" for hard band take unvignetted background """ """ for hard band take unvignetted background """
outfile_backmap3="{}_BackMap3_en{}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], "novign", outfile_post) outfile_backmap3="{}_BackMap3_en{}.{}{}".format(os.path.join(outfile_dir,outkey), eband[index], "novign", outfile_post)
@ -281,10 +284,12 @@ else:
fitext_flag="fitext_flag=yes" fitext_flag="fitext_flag=yes"
nmulsou = "nmulsou=2" nmulsou = "nmulsou=2"
nmaxfit="nmaxfit=4" nmaxfit="nmaxfit=4"
multrad="multrad=20."
cutrad="cutrad=20."
""" allow ermpldet to split sources (no more than two) """ """ allow ermpldet to split sources (no more than two) """
cmd=["ermldet", cmd=["ermldet",
"mllist=%s" %(mllist), "mllist={}".format(mllist),
"boxlist=%s" %(boxlist3), "boxlist=%s" %(boxlist3),
"images=\'{}\'".format(outfile_evtool), "images=\'{}\'".format(outfile_evtool),
"expimages=\'{}\'".format(outfile_expmap), "expimages=\'{}\'".format(outfile_expmap),
@ -297,8 +302,8 @@ cmd=["ermldet",
"likemin=0.", "likemin=0.",
"extlikemin=5.", "extlikemin=5.",
"compress_flag=N", "compress_flag=N",
"cutrad=15.", cutrad,
"multrad=20.", multrad,
"extmin=2.0", "extmin=2.0",
"extmax=35.0", "extmax=35.0",
#"bkgima_flag=Y", looks outdated #"bkgima_flag=Y", looks outdated
@ -337,17 +342,22 @@ if(do_ersensmap==True):
""" """
if(do_ermldet==True): if(do_ermldet==True):
test_exe('ermldet')
if(vign==False): if(vign==False):
print('Run ermldet with vignetted exposure!') print('Run ermldet with vignetted exposure!')
sys.exit() sys.exit()
remove_file(mllist) remove_file(mllist)
remove_file(srcmap) remove_file(srcmap)
print(cmd)
runme(cmd, local_run=local_run) runme(cmd, local_run=local_run)
print(cmd)
#correct_fluxerr_ermldet_forced(mllist) save_ermldet_ds9reg(mllist,scale=60*60,label='det_like')
save_ermldet_ds9reg(mllist,scale=60*60,label='id_src')
correct_fluxerr_ermldet_forced(mllist)
if(forced==True): if(forced==True):
result = check_ermldet_forced(mllist) result = check_ermldet_forced(mllist)
""" for a some reason, for an arbitrary energy band, ermldet break order of sources. Do this forced correction. """ # for a some reason, for an arbitrary energy band, ermldet break order of sources. Do this forced correction.
if(result == False): if(result == False):
correct_srcid_ermldet_forced(mllist) correct_srcid_ermldet_forced(mllist)
@ -359,32 +369,22 @@ if(do_fixcat==True):
print("ERROR: You can fix only non-forced catalog for en0.") print("ERROR: You can fix only non-forced catalog for en0.")
sys.exit() sys.exit()
srcs_remove=[174, srcs_remove=[341,446,346,96]
90, srcs_add = {'4XMM J021738.8-051257':[34.4117002, -5.2159135, 0.624],# 341
299, '4XMM J021733.8-051311':[34.3910215,-5.2199877,2.247],# 341
300, '4XMM J021929.4-051220':[34.8725460,-5.2056849,1.074],#446
504, '4XMM J021930.7-051225':[34.8782267,-5.2072112,0.624],# 446
215, '4XMM J021945.2-045331':[34.9383593,-4.8919843,1.538],#346
401, '4XMM J021929.4-043224':[34.8728586,-4.5400022,0.659555],#96
20,] # keep 671 as unclassified extended source #'4XMM J021929.4-043224':[34.8728586,-4.5400022, 0.660],
srcs_add = {'4XMM J021925.4-042647':[34.8559099,-4.4465007, 1.366], # 147 #'4XMM J021831.8-050059':[34.6328841,-5.0163909,1.529],
'4XMM J021922.8-042655':[34.8451832,-4.4487901, 1.958], # 147 #'4XMM J022131.1-050027':[35.3797879,-5.0075498,0.941],
'4XMM J021929.4-043224':[34.8728586,-4.5400022, 0.660], # 90 #'4XMM J022129.5-045914':[35.3732136,-4.9874025,0.332],
'4XMM J021931.2-043222':[34.8801169,-4.5395495, 2.561], # 90 #'4XMM J022026.3-050251':[35.1098619,-5.0476199,0.551],
'4XMM J021911.2-050550':[34.7968110,-5.0972990, 0.732], # 504 #'4XMM J021925.4-042647':[34.8559099,-4.4465007,1.366],
'4XMM J021919.3-050511':[34.8307176,-5.0864242,4.988], # 504 #'4XMM J021910.9-045108':[34.7954311,-4.8522901,0.898],
'4XMM J021911.5-050501':[34.7981099,-5.0837146,6.834], # 504 #'4XMM J021945.2-045331':[34.9383593,-4.8919843,1.538],
'4XMM J021658.9-044900':[34.2455964,-4.8168126,4.449], # 300 #'4XMM J021733.8-051311':[34.3910215,-5.2199877,2.247],
'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) fix_catalog(mllist=mllist,refimage=outfile_evtool, srcs_remove=srcs_remove, srcs_add=srcs_add)
""" """
@ -451,7 +451,7 @@ if(do_filter_catalog==True):
#filter_mllist(mllist,expcut=5000.0,dlcut=10.0,dlmin=10,dlmax=10000) #filter_mllist(mllist,expcut=5000.0,dlcut=10.0,dlmin=10,dlmax=10000)
""" works the same """ """ works the same """
filter_catprep(catprep,expcut=5000.0,dlmin=10,dlmax=10000,outkey='bright') filter_catprep(catprep,expcut=5000.0,dlmin=10,dlmax=10000,outkey='bright')
filter_catprep(catprep,expcut=5000.0,dlmin=6,dlmax=10,outkey='faint') #filter_catprep(catprep,expcut=5000.0,dlmin=6,dlmax=10,outkey='faint')
if(do_cross_match==True): if(do_cross_match==True):
crossmatch_shu2019(catprep,dlmin=10,refimage=outfile_evtool,crval=[ra_cen, de_cen], crossmatch_shu2019(catprep,dlmin=10,refimage=outfile_evtool,crval=[ra_cen, de_cen],

111
scripts/05_srctool.py
View File

@ -73,20 +73,34 @@ outkey="tm0"
outfile_srctool="{}_SrcTool_".format(outkey) outfile_srctool="{}_SrcTool_".format(outkey)
do_init = False do_init = False
do_merge = False do_merge = False
do_srctool = False do_srctool = False
do_grppha = False do_grppha = False
do_ecf_calc = False # for all bands do_ecf_calc = False # for all bands
do_ecf_print = False # for all bands do_ecf_print = False # for all bands
do_catalog = False do_flux_calc = False # for all bands
do_extended = False
do_ds9reg = False do_catalog = False
do_xmm_catalog = False do_extended = False
do_xmm_final = True do_ds9reg = False
do_euds_final = False
do_euds_dr12 = False # crossmatch eUDS with DR12
do_euds_stat = False
do_euds_cds = False
do_xmm_catalog = False
do_xmm_final = False
do_xmm_xmatch = False
do_xmm_ds9reg = False
do_xmm_cds = False
do_euds_cosmatch = True
do_cross_check = False # Check whether all E3,E6 sources are detected in E0
index=0 index=0
""" работаем именно в этом диапазоне, чтобы спектры покрывали все энергии """ """ чтобы спектры покрывали все энергии работаем в диапазоне 5 """
vign=True vign=True
vignetting = 'vign' if (vign==True) else 'novign' vignetting = 'vign' if (vign==True) else 'novign'
@ -99,7 +113,7 @@ for tmkey in keylist_tm.keys():
for datakey in keylist_tm[tmkey]: for datakey in keylist_tm[tmkey]:
print("--> {}".format(datakey)) print("--> {}".format(datakey))
""" Подготавливаем списки событий индивидуальных наблюдений """ """ Подготавливаем списки событий индивидуальных наблюдений """
outfile_evtool,outfile_expmap=init_events(key=datakey, outfile_evtool,outfile_expmap=init_events(key=datakey,attcorr=True,
eband_index=eband[index], eband_index=eband[index],
infile_dir=infile_dir, infile_dir=infile_dir,
outfile_dir=outfile_dir, outfile_dir=outfile_dir,
@ -144,19 +158,21 @@ if(do_srctool==True):
"eventfiles={}".format(outfile_evtool), "eventfiles={}".format(outfile_evtool),
"prefix=\'{}\'".format(os.path.join(srctool_dir,outfile_srctool)), "prefix=\'{}\'".format(os.path.join(srctool_dir,outfile_srctool)),
"suffix=\'{}\'".format(suffix_srctool), "suffix=\'{}\'".format(suffix_srctool),
"srccoord={}".format(catprep), "srccoord={}".format("../products/eUDS_for_srctool.fits"),
# the same as original file eUDS.fits, but changed ML_CTS --> ML_CTS_0, ML_BKG --> ML_BKG_0, ML_EXP --> ML_EXP_0
#"srcreg=\'fk5;circle * * 60s\'", #"srcreg=\'fk5;circle * * 60s\'",
#"backreg=\'fk5;annulus * * 90s 120s\'", #"backreg=\'fk5;annulus * * 90s 120s\'",
"srcreg=AUTO", "srcreg=AUTO",
"backreg=AUTO", "backreg=AUTO",
"clobber=yes",] "clobber=yes",]
os.system((" ").join(cmd))
print((" ").join(cmd)) print((" ").join(cmd))
#os.system((" ").join(cmd))
#print((" ").join(cmd))
if(do_grppha==True): if(do_grppha==True):
group_spectra("{}/*_SourceSpec_*.fits".format(srctool_dir)) group_spectra("{}/*020_SourceSpec_*.fits".format(srctool_dir))
ecfout="{}_SampleFlux.pickle".format(os.path.join(outfile_dir,outkey)) ecfout="{}_SampleFlux_v1.pickle".format(os.path.join(outfile_dir,outkey))
if(do_ecf_calc==True): if(do_ecf_calc==True):
calc_ecf("{}/tm0_SrcTool_020_ARF_?????.fits".format(srctool_dir), calc_ecf("{}/tm0_SrcTool_020_ARF_?????.fits".format(srctool_dir),
@ -165,16 +181,27 @@ if(do_ecf_calc==True):
if(do_ecf_print==True): if(do_ecf_print==True):
print_ecf(infile=ecfout, emin=emin_kev, emax=emax_kev, eband=eband, skipfrac=10.0) print_ecf(infile=ecfout, emin=emin_kev, emax=emax_kev, eband=eband, skipfrac=10.0)
index=0 fluxout="{}_SherpaFlux.pickle".format(os.path.join(outfile_dir,outkey))
if(do_flux_calc==True):
calc_flux("{}/tm0_SrcTool_020_ARF_?????.fits".format(srctool_dir),
catprep=catprep, emin=emin_kev, emax=emax_kev, eband=eband, outfile=ecfout, simnum=100)
#index=0
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)
rawcat="{}_SourceCatalog_en{}.pickle".format(os.path.join(outfile_dir,outkey), eband[index]) rawcat="{}_SourceCatalog_en{}.pickle".format(os.path.join(outfile_dir,outkey), eband[index])
if(do_catalog==True): if(do_catalog==True):
forced_xmm_sources="{}_MaxLikSourceList_en{}.xmm.pickle".format(os.path.join(outfile_dir,outkey), eband[index]) forced_xmm_sources="{}_MaxLikSourceList_en{}.xmm.pickle".format(os.path.join(outfile_dir,outkey), eband[index])
with open(forced_xmm_sources, 'rb') as f: with open(forced_xmm_sources, 'rb') as f:
print("Reading forced XMM sources from {}".format(forced_xmm_sources))
srcs_forced = pickle.load(f) srcs_forced = pickle.load(f)
print()
print(srcs_forced)
print()
make_catalog(infile='../products/tm0_SourceCatalog_en0.forced.fits', rawcat=rawcat, dlmin=10.0, dlmax=100000, ext_like=1000, make_euds_catalog(infile='../products/tm0_SourceCatalog_en0.forced.fits', rawcat=rawcat, dlmin=10.0, dlmax=100000, ext_like=1000,
emin=emin_kev[index], emax=emax_kev[index], eband=eband[index], emin=emin_kev[index], emax=emax_kev[index], eband=eband[index],
infile_en00cat=catprep, infile_en00cat=catprep,
infile_en01cat='../products/tm0_SourceCatalog_en1.forced.fits', infile_en01cat='../products/tm0_SourceCatalog_en1.forced.fits',
@ -194,7 +221,27 @@ if(do_extended==True):
make_extended(infile=rawcat,outreg="{}_ExtendedCat_en{}.reg".format(os.path.join(outfile_dir,outkey), eband[index])) make_extended(infile=rawcat,outreg="{}_ExtendedCat_en{}.reg".format(os.path.join(outfile_dir,outkey), eband[index]))
if(do_ds9reg==True): if(do_ds9reg==True):
make_final_ds9reg(infile=rawcat,outreg="{}_FinalCat_dl10.reg".format(os.path.join(outfile_dir,outkey))) #make_final_ds9reg(infile=rawcat,outreg="{}_FinalCat_dl10.reg".format(os.path.join(outfile_dir,outkey)))
make_final_ds9reg(infile=rawcat,scale=(60*60)/10,outreg="{}_FinalCat_dl10_talk.reg".format(os.path.join(outfile_dir,outkey)))
if(do_euds_final==True):
""" make final eUDS catalog """
final_euds_catalog(infile=rawcat, outfile_fits='../products/eUDS.fits')
if(do_euds_dr12==True):
crossmatch_dr12('../products/eUDS.fits', catalog=root_path+"/data/4XMM-DR12/4XMM_DR12cat_slim_v1.0_UDS.fits.catalog", devmax=15)
if(do_euds_stat==True):
make_euds_stat(infile="../products/eUDS.fits",fluxlim=5e-14)
if(do_euds_cds==True):
make_euds_cds(infile="../products/eUDS.fits",outfile='../products/eUDS.cds')
if(do_cross_check==True):
""" cross check final eUDS catalog """
cross_check_euds(infile=catprep, euds='../products/eUDS.fits', outkey="../products/en{}_FinalCat_dl10".format(index))
if(do_xmm_catalog==True): if(do_xmm_catalog==True):
""" complile raw forced XMM catalog """ """ complile raw forced XMM catalog """
@ -211,3 +258,27 @@ if(do_xmm_final==True):
""" make final XMM-forced catalog """ """ make final XMM-forced catalog """
final_xmm_catalog(infile='../products/tm0_4XMM-DR12.pickle', outfile_fits='../products/eUDS_4XMM-DR12.fits') final_xmm_catalog(infile='../products/tm0_4XMM-DR12.pickle', outfile_fits='../products/eUDS_4XMM-DR12.fits')
if(do_xmm_xmatch==True):
""" cross-match XMM-forced catalog
outfile_cvs contains 0.3-2.3 keV eUDS flux (col1=flux,col2=err) vs. 4XMM-DR12 flux (col3=flux, col4=err)
XMM flux was converted from 0.2-2.0 keV to 0.3-2.3 keV using wabs*powerlow with wabs=0.02 gamma=2.0
"""
final_xmm_xmatch(infile='../products/eUDS_4XMM-DR12.fits',
xmmslim='../data/4XMM-DR12/4XMM_DR12cat_slim_v1.0_UDS.fits.catalog',
xmmfull='../data/4XMM-DR12/4XMM_DR12cat_v1.0_UDS.fits.catalog',
xmmlim=2e-14,
outfile_flux="../products/eUDS_4XMM-DR12.flux.csv")
if(do_xmm_ds9reg==True):
""" show XMM-forced catalog """
make_xmm_ds9reg_confused(infile='../products/eUDS_4XMM-DR12.fits', outfile='../products/eUDS_4XMM-DR12.confused.reg')
""" obsolete """
if(do_xmm_cds==True):
make_xmm_cds(infile="../products/eUDS_4XMM-DR12.fits", outfile='../products/eUDS_4XMM-DR12.cds')
if(do_euds_cosmatch==True):
""" prepare eUDS catalog for CosMatch (Mescheryakov) """
make_euds_cosmatch(infile='../products/eUDS.fits', outfile='../products/eUDS-CosMatch.fits')

272
scripts/06_plot.py
View File

@ -45,10 +45,10 @@ import glob
from os.path import dirname from os.path import dirname
import inspect import inspect
import uds import uds
from scipy.stats import norm
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import pandas as pd
from uds.utils import * from uds.utils import *
from uds.config import * from uds.config import *
@ -76,7 +76,13 @@ outkey="tm0"
outfile_srctool="{}_SrcTool_".format(outkey) outfile_srctool="{}_SrcTool_".format(outkey)
do_flux_distr = False do_flux_distr = False
do_sens_curve = True do_sens_curve = False
do_4xmm_ratio = False
do_euds_radec_err = False
do_euds_dr12_diff = False
do_euds_dr12_stat = True
do_print_ecf = False
index=0 index=0
@ -87,8 +93,9 @@ if not (os.path.isfile(catalog)==True):
sys.exit() sys.exit()
if(do_flux_distr==True): if(do_flux_distr==True):
data, logbins = get_log_distr(infile=catalog, field='ml_rate', minval=1e-3, maxval=2) data, logbins, mean, median = get_log_distr(infile=catalog, field='ml_rate', minval=1e-3, maxval=2)
fig, ax = plt.subplots() fig, ax = plt.subplots()
for axis in ['top','bottom','left','right']: for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(1) ax.spines[axis].set_linewidth(1)
@ -101,9 +108,9 @@ if(do_flux_distr==True):
plt.yscale('log') plt.yscale('log')
plt.savefig(catalog.replace("main.selected.csv", "ml_rate.png"), bbox_inches='tight') plt.savefig(catalog.replace("main.selected.csv", "ml_rate.png"), bbox_inches='tight')
plt.close(fig) plt.close(fig)
data, logbins = get_log_distr(infile=catalog, field='ml_flux', minval=1e-15, maxval=2e-12)
data, logbins, mean, median = get_log_distr(infile=catalog, field='ml_flux', minval=1e-15, maxval=2e-12)
fig, ax = plt.subplots() fig, ax = plt.subplots()
for axis in ['top','bottom','left','right']: for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(1) ax.spines[axis].set_linewidth(1)
@ -143,11 +150,12 @@ if(do_sens_curve==True):
ax.set_xlim([3e-16,1e-13]) ax.set_xlim([3e-16,1e-13])
ax.set_ylim([0.3,160]) ax.set_ylim([0.3,160])
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/eFEDS.dat",header=None,names=['limflux','area']) 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(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']) 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") plt.plot(df['limflux'],df['area'], color="blue", linestyle='solid',label="COSMOS Legacy")
@ -187,3 +195,251 @@ if(do_sens_curve==True):
Model Flux 0.00028334 photons (3.4012e-13 ergs/cm^2/s) range (0.30000 - 2.3000 keV) 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) Model Flux 0.0001619 photons (2.4336e-13 ergs/cm^2/s) range (0.50000 - 2.0000 keV)
""" """
if(do_4xmm_ratio==True):
filename="../products/eUDS_4XMM-DR12.flux.csv"
data, logbins, mean, median = get_log_distr(infile=filename, field='ratio', minval=8e-2, maxval=60, nbin=60)
print("Median {}".format(median))
print(" Mean {}".format(mean))
print("Ntotal {}".format(len(data)))
fig, ax = plt.subplots()
#plt.figure(figsize=(5,5))
#plt.figure().set_figheight(3.6)
#plt.rcParams['figure.figsize'] = [4, 4]
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.axvline(x = median, color = 'black', label = 'Median value', linestyle='dashed')
plt.xlabel('Energy flux ratio',fontsize=14, fontweight='normal')
plt.ylabel('Number',fontsize=14, fontweight='normal')
plt.grid(visible=True)
plt.xscale('log')
plt.yscale('log')
plt.savefig(filename.replace(".csv", ".distr.png"), bbox_inches='tight')
plt.close(fig)
df = pandas.read_csv(filename)
fig, ax = plt.subplots()
#plt.rcParams['figure.figsize'] = [4, 4]
#plt.figure(figsize=(5,5))
#plt.figure().set_figheight(4)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(1)
ax.tick_params(axis="both", width=1, labelsize=14)
plt.plot(np.array([1e-17,1e-11]), np.array([1e-17,1e-11]), linestyle = 'solid', color='black')
plt.plot(np.array([1e-17,1e-11]), 10*np.array([1e-17,1e-11]), linestyle = 'dashed', color='black')
plt.plot(np.array([1e-17,1e-11]), 0.1*np.array([1e-17,1e-11]), linestyle = 'dotted', color='black')
#plt.plot(np.array([1e-17,1e-11]), median*np.array([1e-17,1e-11]), linestyle = 'dashed', color='blue')
plt.errorbar(df['dr12_flux'],df['euds_flux'], yerr=df['euds_flux_err'], xerr=df['dr12_flux_err'], linestyle='None', color='black', label='Errors')
plt.plot(df['dr12_flux'],df['euds_flux'], marker="o", linewidth=1, linestyle='None', markerfacecolor='Gold',markeredgecolor="black",)
#plt.axvline(x = median, color = 'black', label = 'Median value', linestyle='dashed')
plt.xlabel('4XMM-DR12 Energy flux (erg s$^{-1}$ cm$^{-2}$)',fontsize=14, fontweight='normal')
plt.ylabel('eUDS Energy flux (erg s$^{-1}$ cm$^{-2}$)',fontsize=14, fontweight='normal')
plt.grid(visible=True)
plt.xlim(2e-16,1e-12)
plt.ylim(1e-15,2e-12)
plt.xscale('log')
plt.yscale('log')
plt.savefig(filename.replace(".csv", ".png"), bbox_inches='tight')
plt.close(fig)
if(do_euds_radec_err==True):
filename="../products/eUDS.fits"
ecat={}
hdul = fits.open(filename)
etab = hdul[1].data
hdul.close()
x=[]
y=[]
for s in etab:
if ("XMM" in s['DR12_IAU_NAME']):
print("Skip ",s['DR12_IAU_NAME'])
continue
if (s['EXT_LIKE']>0.0):
print("Skip extended ",s['ID_SRC'])
continue
x.append(s['DET_LIKE'])
y.append(s['RADEC_ERR'])
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([10,10000])
ax.set_ylim([0.1,40])
plt.plot(x,y,".", color="gold", markersize=12, markeredgewidth=1.5, markeredgecolor="black")
plt.xlabel('DET_LIKE',fontsize=14, fontweight='normal')
plt.ylabel('RADEC_ERR (arcsec)',fontsize=14, fontweight='normal')
plt.grid(visible=True)
plt.xscale('log')
plt.yscale('log')
png="../products/en0_radec_err.png"
plt.savefig(png, bbox_inches='tight')
plt.close(fig)
stat={}
hdul = fits.open("../products/eUDS.dr12.cross.fits")
tab = hdul[1].data
hdul.close()
print(len(np.unique(tab['ID'])))
for rec in tab:
sid=rec['ID']
if sid in stat.keys():
stat[sid]=stat[sid]+1
else:
stat[sid]=1
unique=0
double=0
triple=0
for key in stat.keys():
if(stat[key]==1):
unique=unique+1
if(stat[key]==2):
double=double+1
if(stat[key]==3):
triple=triple+1
print("unique: {}, double: {}. triple: {}".format(unique, double, triple))
if(do_euds_dr12_diff==True):
# read forced first
# fieldnames = ['euds_det_like', 'euds_flux', 'euds_flux_err', 'dr12_flux', 'dr12_flux_err', 'ratio']
filename="../products/eUDS_4XMM-DR12.flux.csv"
df = pandas.read_csv(filename)
# calculate difference similar to Bruner
threshold=5e-14
# forced catalog contains CONF flag, take it to remove from histogram
hdul = fits.open("../products/eUDS_4XMM-DR12.fits")
forced = hdul[1].data
fcat={}
for rec in forced:
key=rec['DR12_SRCID']
fcat[key]={'conf':rec['CONF'],}
# read cross-match results (not forced)
hdul = fits.open('../products/eUDS.dr12.cross.fits')
tb = hdul[1].data
cross_diff=[]
cross_ratio=[]
mean=0.0
count=0
for ind,s in enumerate(tb):
key=rec['DR12_SRCID']
#if not (tb[ind]['dr12_flux']>threshold and tb[ind]['src_flux']>threshold):
if not (tb[ind]['dr12_flux']>threshold):
continue
ape_flux=(tb[ind]['ape_cts']-tb[ind]['ape_bkg'])/tb[ind]['ape_exp']/ecf[index]
cross_dr12_flux=(tb[ind]['dr12_flux'])
cross_dr12_flux_error=(tb[ind]['dr12_flux_error'])
cross_euds_flux=(tb[ind]['src_flux'])
#cross_euds_flux=ape_flux
cross_euds_flux_error=(tb[ind]['src_flux_error'])
print(cross_euds_flux,cross_euds_flux_error,tb[ind]['ape_cts'],tb[ind]['det_like'])
d=(cross_dr12_flux - cross_euds_flux) / np.sqrt(cross_dr12_flux_error**2 + cross_euds_flux_error**2)
#if(d < -15):
# continue
r=cross_euds_flux/cross_dr12_flux
cross_diff.append(d)
cross_ratio.append(r)
#if(d < -15.0):
#print("GGG",d,r,cross_euds_flux,cross_dr12_flux)
#print(tb[ind]['DR12_NAME'],tb[ind]['DR12_RA'],tb[ind]['DR12_DEC'])
mean=mean + (cross_euds_flux/cross_dr12_flux)
count=count+1
print("min={:.2f}, max={:.2f}, mean={:.2f}, median={:.2f}, std={:.2f} N={}".format(min(cross_diff),
max(cross_diff),
np.mean(cross_diff),
np.median(cross_diff),
np.std(cross_diff),
len(cross_diff)))
print("ratio mean",mean/count)
print("ratio median",np.median(cross_ratio))
nbin=50
bins=np.linspace(-30,30, nbin)
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(cross_diff, bins=bins,histtype='step', color='green', linewidth=1, linestyle='solid', density=True)
#x = np.arange(-10, 10, 0.001)
#plot normal distribution with mean 0 and standard deviation 1
#plt.plot(x, norm.pdf(x, 0, 1), color='red', linewidth=2)
plt.ylabel('Relative fraction',fontsize=14, fontweight='normal')
plt.xlabel('(F$_{XMM}$-F$_{eUDS}$)/$\sqrt{\Delta F_{XMM}^{2}+\Delta F_{eUDS}^{2}}$',fontsize=14, fontweight='normal')
plt.grid(visible=True)
plt.xscale('linear')
plt.yscale('linear')
ax.set_xlim([-31, 31])
ax.set_ylim([0, 0.2])
plt.savefig("../products/cross-match_diff.png", bbox_inches='tight')
plt.close(fig)
nbin=200
bins=np.linspace(-1,1.0, nbin)
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(cross_ratio, bins=bins, histtype='step', color='green', linewidth=1, linestyle='solid', density=False)
#x = np.arange(-10, 10, 0.001)
#plot normal distribution with mean 0 and standard deviation 1
#plt.plot(x, norm.pdf(x, 0, 1), color='red', linewidth=2)
plt.xlabel('eUDS and 4XMM-DR12 flux ratio',fontsize=14, fontweight='normal')
plt.ylabel('Relative fraction',fontsize=14, fontweight='normal')
plt.grid(visible=True)
plt.xscale('linear')
plt.yscale('linear')
plt.savefig("../products/cross-match_ratio.png", bbox_inches='tight')
plt.close(fig)
if(do_print_ecf==True):
filename='../data/ECF/ecf_tbabspow_g2nh0.02.pkl'
with open(filename, 'rb') as f:
ecf_table = pickle.load(f)
"""
for key in table.keys():
print("{} --> {}".format(key,table[key]))
"""
print(ecf_table[(0.3,2.3)])
print(ecf_table[(0.3,0.6)])
print(ecf_table[(0.6,2.3)])
print(ecf_table[(2.3,5.0)])
print(ecf_table[(5.0,8.0)])
print()
print(ecf_table[(0.5,1.0)]) # 4XMM-DR12 EP2 band
print(ecf_table[(1.0,2.0)]) # 4XMM-DR12 EP3 band
if(do_euds_dr12_stat==True):
dr12_stat(infile='../products/eUDS_4XMM-DR12.fits',
xmmslim='../data/4XMM-DR12/4XMM_DR12cat_slim_v1.0_UDS.fits.catalog',
xmmlim=None,
outfile_reg='../products/dr12_fluxlim.reg')

55
scripts/print_radec_err.py Normal file
View File

@ -0,0 +1,55 @@
import glob, sys
from astropy.io import fits
import statistics
import pandas as pd
import numpy as np
import csv
colnames=['Match', 'Ref', 'Dup', 'RA', 'Dec', 'Incl']
flist=glob.glob('../products/tm[1,5,6,7]*en0*src.fits')
total1=0
total2=0
for f in flist:
hdul = fits.open(f)
try:
table=hdul[1].data['RADEC_ERR']
except:
continue
hdul.close()
dref=f.replace("src.fits","ref.fits")
hdul = fits.open(dref)
try:
rtable=hdul[1].data['RA']
except:
continue
hdul.close()
err=[]
err10=[]
dfile=f.replace("shu2019.src.fits","xfm.log.dat.awk")
with open(dfile) as csvfile:
spamreader = csv.reader(csvfile, delimiter=' ')
total=0
for row in spamreader:
if(row[5] == 'N'):
continue
for ii,values in np.ndenumerate(rtable):
if(int(row[2]) == int(ii[0])):
total=total+1
#print(">>",row[5],rtable[ii],table[ii],ii[0])
if(table[ii]>0.0):
err.append(table[ii])
if(table[ii]>10.0):
err10.append(table[ii])
total1=total1+total
total2=total2+len(err)
#print("{} MIN {:.2f} MAX {:.2f} ({}/{})".format(f[12:15],min(err),max(err),len(err10),len(err)))
print("{} MIN {:.2f} MAX {:.2f} ({}/{})".format(f,min(err),max(err),len(err10),len(err)))
#sys.exit()
#print(total1,total2)

67
scripts/print_transform.py Normal file
View File

@ -0,0 +1,67 @@
import glob, sys
from astropy.io import fits
import statistics
import pandas as pd
import numpy as np
import csv
from uds.config import *
colnames=['Match', 'Ref', 'Dup', 'RA', 'Dec', 'Incl']
flist=glob.glob('../products/tm[1,5,6,7]*en0*.xfm')
total1=0
total2=0
for f in flist:
pos = f.find("_SourceCatalog")
key=f[12:pos]
hdul = fits.open(f)
try:
tab=hdul[1].data
except:
continue
hdul.close()
dsrc=f.replace(".xfm",".shu2019.src.fits")
hdul = fits.open(dsrc)
try:
hdr=hdul[0].header
cdelt=hdr['CDELT2']*3600
except:
continue
hdul.close()
a11=tab[0]['a11']
a12=tab[0]['a12']
a22=tab[0]['a22']
a21=tab[0]['a21']
x_scale=tab[0]['x_scale']
y_scale=tab[0]['y_scale']
t1=tab[0]['t1']*cdelt
t2=tab[0]['t2']*cdelt
sx1=np.sign(a11)*np.sqrt(a11**2 + a12**2)
sy1=np.sign(a22)*np.sqrt(a21**2 + a22**2)
sx2=np.sign(a11)*np.sqrt(a11**2 + a21**2)
sy2=np.sign(a22)*np.sqrt(a12**2 + a22**2)
a1 = np.arctan2(a12/sx1,a11/sy2) * 180 / np.pi * 60
a2 = np.arctan2(a12/sx1,a22/sy2) * 180 / np.pi * 60
b1 = np.arctan2(a21,a11) * 180 / np.pi * 60
b2 = np.arctan2(a21,a22) * 180 / np.pi * 60
det=a11*a22-a12*a21
alpha=355*np.pi/180
a3 = np.arctan2(np.sin(alpha),np.cos(alpha)) * 180 / np.pi
print("{:+.8} {:+.8}".format(a11,a12))
print("{:+.8} {:+.8}".format(a21,a22))
#print("{} d {:.10f} s {:.8f} {:.8f} {:.8f} {:.8f} | {:+.8f} {:+.8f} {:+.8f} {:+.8f} t {:+.2f} {:+.2f}".format(obslist[key],np.sqrt(det)*cdelt,sx1*cdelt,sx2*cdelt,sy1*cdelt,sy2*cdelt,a1,a2,b1,b2,t1,t2))
#print("{} d {:.10f} s {:.8f} {:.8f} {:.8f} {:.8f} | {:+.8f} {:+.8f} {:+.8f} {:+.8f} t {:+.2f} {:+.2f}".format(obslist[key],np.sqrt(det)*cdelt,sx1*cdelt,sx2*cdelt,sy1*cdelt,sy2*cdelt,a1,a2,b1,b2,t1,t2))

468
uds/uds/astrometry.py Normal file
View File

@ -0,0 +1,468 @@
import glob
import os
import sys
import numpy as np
from astropy.io import fits
import pickle
from scipy.stats import sem
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord # High-level coordinates
from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames
from astropy.coordinates import Angle, Latitude, Longitude # Angles
import matplotlib.pyplot as plt
from astropy.table import QTable, Table, Column
from mpdaf.obj import WCS as mWCS
from sherpa.astro.ui import *
from uds.config import *
from scipy.optimize import minimize
from random import uniform
def sky2pix(wcs,dsrc):
# converts ra,dec --> x,y
# checks
#sk = wcs.pix2sky([0.0,0.0])
#print(sk)
#px = wcs.sky2pix([-6.86525105,36.54071803],nearest=False)
#print(px)
#sk = wcs.pix2sky([px[0][1],px[0][0]])
#print(sk)
for idx, s in enumerate(dsrc):
ra0=s['ra']
dec0=s['dec']
px = wcs.sky2pix([dec0,ra0],nearest=False)
s['x']=px[0][1]
s['y']=px[0][0]
#sk = wcs.pix2sky([y0,x0])
#ra=sk[0][1]
#dec=sk[0][0]
#px = wcs.sky2pix([dec,ra],nearest=False)
#x=px[0][1]
#y=px[0][0]
#print("{:.4f} {:.4f} --> {} {} --> {:.4f} {:.4f} --> {} {}".format(ra0,dec0,x0,y0,ra,dec,x,y))
def pix2sky(wcs,dsrc):
# converts x,y --> ra1,dec1
for idx, s in enumerate(dsrc):
sk = wcs.pix2sky([s['y'],s['x']])
s['ra1']=sk[0][1]
s['dec1']=sk[0][0]
def plot_resid(dsrc,dref,crval1,crval2):
# calculates RMS of original ra,dec
center_crd = SkyCoord(crval1, crval2, frame=FK5(), unit="deg")
offset=[]
resid=[]
error=[]
for idx, s in enumerate(dsrc):
src_crd = SkyCoord(dsrc[idx]['ra'], dsrc[idx]['dec'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(dref[idx]['ra'], dref[idx]['dec'], frame=FK5(), unit="deg")
offset.append(src_crd.separation(center_crd).arcmin)
resid.append(src_crd.separation(ref_crd).arcsec)
error.append(s['radec_err'])
indices = sorted(
range(len(offset)),
key=lambda index: offset[index]
)
return ([offset[index] for index in indices],
[resid[index] for index in indices],
[error[index] for index in indices])
def plot_resid1(dsrc,dref,crval1,crval2):
# calculates RMS of original ra,dec
center_crd = SkyCoord(crval1, crval2, frame=FK5(), unit="deg")
offset=[]
resid=[]
error=[]
for idx, s in enumerate(dsrc):
src_crd = SkyCoord(dsrc[idx]['ra1'], dsrc[idx]['dec1'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(dref[idx]['ra'], dref[idx]['dec'], frame=FK5(), unit="deg")
offset.append(src_crd.separation(center_crd).arcmin)
resid.append(src_crd.separation(ref_crd).arcsec)
error.append(s['radec_err'])
indices = sorted(
range(len(offset)),
key=lambda index: offset[index]
)
return ([offset[index] for index in indices],
[resid[index] for index in indices],
[error[index] for index in indices])
def get_rms(dsrc,dref):
# calculates RMS of original ra,dec
resid=0.0
n=0
for idx, s in enumerate(dsrc):
radec_err=s['radec_err']
src_crd = SkyCoord(dsrc[idx]['ra'], dsrc[idx]['dec'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(dref[idx]['ra'], dref[idx]['dec'], frame=FK5(), unit="deg")
sep=src_crd.separation(ref_crd).arcsec
resid = resid+sep**2
n=n+1
return np.sqrt(resid/n)
def get_rms_w(dsrc,dref):
# calculates RMS of original ra,dec
resid=0.0
n=0
A=0.0
B=0.0
for idx, s in enumerate(dsrc):
radec_err=s['radec_err']
src_crd = SkyCoord(dsrc[idx]['ra'], dsrc[idx]['dec'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(dref[idx]['ra'], dref[idx]['dec'], frame=FK5(), unit="deg")
sep=src_crd.separation(ref_crd).arcsec
A=A+sep/radec_err**2
B=B+1.0/radec_err**2
n=n+1
return A/B
# A=Total(table_lc_sav.rate/(table_lc_sav.err)^2)
# B=Total(1.0/(table_lc_sav.err)^2)
def get_rms1_w(dsrc,dref):
# calculates RMS of original ra,dec
resid=0.0
n=0
A=0.0
B=0.0
for idx, s in enumerate(dsrc):
radec_err=s['radec_err']
src_crd = SkyCoord(dsrc[idx]['ra1'], dsrc[idx]['dec1'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(dref[idx]['ra'], dref[idx]['dec'], frame=FK5(), unit="deg")
sep=src_crd.separation(ref_crd).arcsec
A=A+sep/radec_err**2
B=B+1.0/radec_err**2
n=n+1
return A/B
def get_chi2_radec(dsrc,dref):
# calculates RMS of original ra,dec
resid=0.0
n=0
for idx, s in enumerate(dsrc):
radec_err=s['radec_err']
src_crd = SkyCoord(dsrc[idx]['ra'], dsrc[idx]['dec'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(dref[idx]['ra'], dref[idx]['dec'], frame=FK5(), unit="deg")
sep=src_crd.separation(ref_crd).arcsec
resid = resid+sep**2/radec_err**2
n=n+1
return resid
def get_rms1(dsrc,dref):
# calculates RMS of modified ra1,dec1
resid=0.0
n=0
for idx, s in enumerate(dsrc):
radec_err=s['radec_err']
src_crd = SkyCoord(dsrc[idx]['ra1'], dsrc[idx]['dec1'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(dref[idx]['ra'], dref[idx]['dec'], frame=FK5(), unit="deg")
sep=src_crd.separation(ref_crd).arcsec
resid = resid+sep*sep
n=n+1
return np.sqrt(resid/n)
def get_chi2_radec1(dsrc,dref):
# calculates RMS of modified ra1,dec1
resid=0.0
n=0
for idx, s in enumerate(dsrc):
radec_err=s['radec_err']
src_crd = SkyCoord(dsrc[idx]['ra1'], dsrc[idx]['dec1'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(dref[idx]['ra'], dref[idx]['dec'], frame=FK5(), unit="deg")
sep=src_crd.separation(ref_crd).arcsec
resid = resid+sep**2/radec_err**2
n=n+1
return resid
def do_transform_resid(params, wcs, dsrc, dref, crval1, crval2):
angle_in, crval1_in, crval2_in = params
# fills x,y
sky2pix(wcs,dsrc)
wcs1=wcs.copy()
# do transform here
#wcs1.rotate(angle_in/60)
wcs1.set_crval1(crval1_in)
wcs1.set_crval2(crval2_in)
# fills ra1,dec1
pix2sky(wcs1,dsrc)
# caclulate new statistics
resid1 = get_chi2_radec1(dsrc,dref)
#print("--> {:+.4f} {:+.4f} {:.4f}".format((crval1-wcs.get_crval1())*3600,(crval2-wcs.get_crval2())*3600,resid1))
return resid1
def do_astro_corr(src=None, ref=None, catalog=None):
""" calculates astronomy correction """
log=catalog.replace(".fits", ".shift.log")
log0=catalog.replace(".fits", ".orig.qdp")
log1=catalog.replace(".fits", ".shift.qdp")
fout=catalog.replace(".fits", ".shift-map.fits")
title=os.path.split(catalog)[1]
key=title.replace("_SourceCatalog_en0.fits","")
#png=catalog.replace(".fits", ".opti.png")
#png=png.replace(key,"{}_{}".format(obslist[key],key))
t = Table(names=('rota', 'shift_ra','shift_dec','chi2',), dtype=('f8','f8','f8','f8'))
print(src)
hdul = fits.open(src)
src_map_data = hdul[0].data
src_map_hdr = hdul[0].header
src_tab_data = hdul[1].data
src_tab_hdr = hdul[1].header
hdul.close()
dsrc=[]
for rec in src_tab_data:
d = dict(ra = rec['RA'], dec = rec['DEC'], ra_err = rec['RA_ERR'], dec_err = rec['DEC_ERR'], radec_err = rec['RADEC_ERR'])
dsrc.append(d)
print(ref)
hdul = fits.open(ref)
ref_map_data = hdul[0].data
ref_map_hdr = hdul[0].header
ref_tab_data = hdul[1].data
ref_tab_hdr = hdul[1].header
hdul.close()
dref=[]
for rec in ref_tab_data:
d = dict(ra = rec['RA'], dec = rec['DEC'], ra_err = rec['RA_ERR'], dec_err = rec['DEC_ERR'])
dref.append(d)
rms = get_rms(dsrc,dref)
chi2_radec = get_chi2_radec(dsrc,dref)
wcs_src=mWCS(src_map_hdr)
wcs_src.info()
crval1=wcs_src.get_crval1()
crval2=wcs_src.get_crval2()
#print(crval1,crval2)
offset0, sep0, err0 = plot_resid(dsrc,dref,crval1,crval2)
with open(log0, 'w') as writer:
writer.write("label Y Separation, arcsec\nlabel X Offset, arcmin\nr y -1 15\n")
writer.write("ma 9 on\n")
writer.write("read serr 2\n")
for idx, s in enumerate(offset0):
writer.write("{:.2f} {:.2f} {:.2f}\n".format(offset0[idx],sep0[idx],err0[idx]))
wcs_ref=mWCS(ref_map_hdr)
wcs_ref.info()
chi2_radec1_opt=1000
rota_opt=0.0
crval1_opt=crval1
crval2_opt=crval2
rota_min=-30
rota_max=30
Rsim=10.0
Nsim=20000
crd = SkyCoord(crval1, crval2, frame=FK5(), unit="deg")
i=0
while i < Nsim:
#rota_rand = uniform(rota_min, rota_max)
rota_rand=0.0
crval1_rand = uniform(crval1-Rsim/3600, crval1+Rsim/3600)
crval2_rand = uniform(crval2-Rsim/3600, crval2+Rsim/3600)
#crval1_rand = crval1
#crval2_rand = crval2
sim_crd = SkyCoord(crval1_rand, crval2_rand, frame=FK5(), unit="deg")
if(sim_crd.separation(crd).arcsec > Rsim):
continue
chi2_radec1 = do_transform_resid([rota_rand, crval1_rand, crval2_rand], wcs_src, dsrc, dref, crval1, crval2)
t.add_row((rota_rand,(crval1-crval1_rand)*3600,(crval2-crval2_rand)*3600,chi2_radec1))
if(chi2_radec1 < chi2_radec1_opt):
chi2_radec1_opt = chi2_radec1
rota_opt = rota_rand
crval1_opt = crval1_rand
crval2_opt = crval2_rand
i += 1
# set optimal parameters
r = do_transform_resid([rota_opt, crval1_opt, crval2_opt], wcs_src, dsrc, dref, crval1, crval2)
print(chi2_radec1_opt,'==',r)
shift_ra=(crval1-crval1_opt)*3600
shift_dec=(crval2-crval2_opt)*3600
t.write(fout, format='fits', overwrite=True)
rms1 = get_rms1(dsrc,dref)
offset1, sep1, err1 = plot_resid1(dsrc,dref,crval1,crval2)
with open(log1, 'w') as writer:
writer.write("label Y Separation, arcsec\nlabel X Offset, arcmin\nr y -1 15\n")
writer.write("ma 9 on\n")
writer.write("read serr 2\n")
for idx, s in enumerate(offset1):
writer.write("{:.2f} {:.2f} {:.2f}\n".format(offset1[idx],sep1[idx],err1[idx]))
with open(log, 'w') as writer:
writer.write("{} {:+.2f} {:+.2f} {:+.2f} ({:.2f} - {:.2f}) = {:.2f} {:.2f} ({:.2f}) {:.2f} ({:.2f}) N={} -- {}\n".format(obslist[key],
rota_opt,
shift_ra,
shift_dec,
get_chi2_radec(dsrc,dref),
get_chi2_radec1(dsrc,dref),
chi2_radec-chi2_radec1_opt,
rms,get_rms_w(dsrc,dref),
rms1,get_rms1_w(dsrc,dref),
len(dsrc),key))
"""
Nsim=20000 Shift only
N01 +0.00 +4.76 +0.07 (171.22 - 28.06) = 143.15 6.35 (5.11) 3.92 (1.40) N=22 -- tm1_obs_1
N02 +0.00 +2.77 -2.74 (354.15 - 225.27) = 128.88 10.10 (5.34) 9.38 (2.61) N=31 -- tm5_obs_1
N03 +0.00 +3.88 -3.65 (542.25 - 303.25) = 239.00 8.74 (6.87) 6.73 (3.07) N=36 -- tm5_obs_2
N04 +0.00 +0.28 +0.03 (43.16 - 43.04) = 0.12 6.04 (3.39) 6.07 (3.36) N=17 -- tm6_park_1
N05 +0.00 +0.65 -0.78 (175.10 - 169.55) = 5.56 6.77 (4.47) 6.54 (4.51) N=38 -- tm6_scan_1
N06 +0.00 -0.38 +0.27 (26.60 - 26.37) = 0.23 6.54 (4.27) 6.52 (4.23) N=11 -- tm6_park_2
N07 +0.00 +0.04 -0.35 (128.82 - 127.98) = 0.84 5.76 (4.21) 5.75 (4.18) N=35 -- tm6_scan_2
N08 +0.00 +0.58 -0.22 (10.07 - 9.79) = 0.29 3.94 (2.83) 3.89 (2.82) N= 9 -- tm6_park_3
N09 +0.00 +1.35 -0.56 (95.06 - 87.04) = 8.02 5.73 (4.45) 5.60 (4.28) N=38 -- tm6_scan_3
N10 +0.00 +0.56 +0.36 (15.71 - 15.41) = 0.30 4.22 (3.49) 4.22 (3.43) N=10 -- tm6_park_4
N11 +0.00 +1.08 -0.62 (116.34 - 108.60) = 7.73 5.38 (4.41) 5.11 (4.26) N=37 -- tm6_scan_4
N12 +0.00 -0.56 -1.56 (84.93 - 59.86) = 25.07 4.20 (2.45) 4.09 (1.88) N=30 -- tm6_obs_1
N13 +0.00 +0.19 -1.03 (17.58 - 13.28) = 4.30 3.38 (1.90) 3.23 (1.31) N=13 -- tm6_obs_2_badpix
N14 +0.00 +1.03 -1.99 (63.23 - 35.27) = 27.97 4.36 (2.73) 3.27 (2.05) N=23 -- tm6_obs_3_badpix
N15 +0.00 +0.88 -1.66 (187.51 - 164.26) = 23.25 6.93 (3.29) 6.62 (2.68) N=33 -- tm7_obs_1
N16 +0.00 +2.74 -0.52 (113.70 - 75.09) = 38.61 5.34 (4.23) 4.73 (3.42) N=26 -- tm7_obs_2
Nsim=20000 Shift+Rotation
N01 -1.50 +4.93 +0.00 (171.22 - 27.67) = 143.55 6.35 (5.11) 3.83 (1.38) N=22 -- tm1_obs_1
N02 -2.15 +2.47 -3.05 (354.15 - 223.27) = 130.88 10.10 (5.34) 9.29 (2.71) N=31 -- tm5_obs_1
N03 -5.08 +3.87 -4.04 (542.25 - 296.00) = 246.25 8.74 (6.87) 6.93 (2.82) N=36 -- tm5_obs_2
N04 +5.90 -0.19 +0.50 (43.16 - 39.56) = 3.60 6.04 (3.39) 5.57 (3.57) N=17 -- tm6_park_1
N05 -5.84 +1.37 -0.85 (175.10 - 111.02) = 64.08 6.77 (4.47) 5.53 (3.58) N=38 -- tm6_scan_1
N06 -0.30 -0.27 +0.31 (26.60 - 26.36) = 0.24 6.54 (4.27) 6.50 (4.27) N=11 -- tm6_park_2
N07 -5.49 +0.73 -0.98 (128.82 - 86.16) = 42.66 5.76 (4.21) 4.94 (3.62) N=35 -- tm6_scan_2
N08 -0.28 +0.51 -0.29 (10.07 - 9.78) = 0.29 3.94 (2.83) 3.90 (2.82) N= 9 -- tm6_park_3
N09 -4.56 +0.91 -0.40 (95.06 - 64.56) = 30.49 5.73 (4.45) 4.93 (3.56) N=38 -- tm6_scan_3
N10 +2.01 +0.54 +0.63 (15.71 - 15.17) = 0.54 4.22 (3.49) 4.03 (3.39) N=10 -- tm6_park_4
N11 -3.83 +1.01 -1.29 (116.34 - 82.25) = 34.08 5.38 (4.41) 4.75 (3.72) N=37 -- tm6_scan_4
N12 -1.35 -0.23 -1.27 (84.93 - 60.30) = 24.63 4.20 (2.45) 4.20 (1.82) N=30 -- tm6_obs_1
N13 -1.52 +0.17 -1.38 (17.58 - 13.30) = 4.29 3.38 (1.90) 3.21 (1.28) N=13 -- tm6_obs_2_badpix
N14 -1.24 +0.63 -2.18 (63.23 - 34.42) = 28.81 4.36 (2.73) 3.28 (1.92) N=23 -- tm6_obs_3_badpix
N15 -5.23 +0.82 -1.47 (187.51 - 161.20) = 26.31 6.93 (3.29) 6.31 (2.66) N=33 -- tm7_obs_1
N16 -5.41 +2.47 -0.97 (113.70 - 60.40) = 53.30 5.34 (4.23) 4.59 (2.89) N=26 -- tm7_obs_2
"""
def do_astro_corr_rota(src=None, ref=None, catalog=None):
""" calculates astronomy correction """
title=os.path.split(catalog)[1]
png=catalog.replace(".fits", ".rota.png")
log=catalog.replace(".fits", ".rota.log")
key=title.replace("_SourceCatalog_en0.fits","")
png=png.replace(key,"{}_{}".format(obslist[key],key))
print(src)
hdul = fits.open(src)
src_map_data = hdul[0].data
src_map_hdr = hdul[0].header
src_tab_data = hdul[1].data
src_tab_hdr = hdul[1].header
hdul.close()
dsrc=[]
for rec in src_tab_data:
d = dict(ra = rec['RA'], dec = rec['DEC'], ra_err = rec['RA_ERR'], dec_err = rec['DEC_ERR'], radec_err = rec['RADEC_ERR'])
dsrc.append(d)
print(ref)
hdul = fits.open(ref)
ref_map_data = hdul[0].data
ref_map_hdr = hdul[0].header
ref_tab_data = hdul[1].data
ref_tab_hdr = hdul[1].header
hdul.close()
dref=[]
for rec in ref_tab_data:
d = dict(ra = rec['RA'], dec = rec['DEC'], ra_err = rec['RA_ERR'], dec_err = rec['DEC_ERR'])
dref.append(d)
resid = get_rms(dsrc,dref)
print(resid)
wcs_src=mWCS(src_map_hdr)
wcs_src.info()
wcs_ref=mWCS(ref_map_hdr)
wcs_ref.info()
x=[]
y=[]
rms_min=1000.0
rota_min=1000.0
for angle in np.arange(-0.5, 0.5, 0.01):
resid1=do_transform(angle,wcs_src,dsrc,dref)
x.append(angle*60)
y.append(resid1/resid)
if(resid1<rms_min):
rms_min=resid1
rota_min=angle*60
with open(log, 'w') as writer:
writer.write("{} {} {:.2f} {:.2f} {:.2f} {:.2f}%\n".format(obslist[key],title,rota_min, rms_min, resid, (1.0-rms_min/resid)*100.0))
fig, ax = plt.subplots( nrows=1, ncols=1 ) # create figure & 1 axis
ax.plot(x,y)
plt.ylabel('RMS/RMS0, RMS0={:.2f} arcsec'.format(resid))
plt.xlabel('Rotation, arcmin')
plt.title("{} {} {} pairs, {:.1f}%".format(obslist[key], title,len(dsrc),(1.0-rms_min/resid)*100.0))
plt.grid(True)
plt.axvline(x = 0.0, color = 'b', label = 'zero')
fig.savefig(png)
plt.close(fig) # close the figure window
def do_transform(angle,wcs,dsrc,dref):
# fills x,y
sky2pix(wcs,dsrc)
# do transform here
wcs.rotate(angle)
# fills ra1,dec1
pix2sky(wcs,dsrc)
resid1 = get_rms1(dsrc,dref)
return resid1

52
uds/uds/config.py
View File

@ -43,18 +43,36 @@ wcslist={'tm1_obs_1':[34.7279760,-5.0680267],
'tm6_scan_3':[34.5405596,-4.8088748], 'tm6_scan_3':[34.5405596,-4.8088748],
'tm6_scan_4':[34.5405596,-4.8088748]} 'tm6_scan_4':[34.5405596,-4.8088748]}
""" Это просто индекс диапазона для выходных файлов. """ """ like in the paper (Table 1) """
eband=[ "0", "1", "2", "3", "4", "5", "6"] obslist={'tm1_obs_1':'N01',
""" Диапазоны энергий. """ 'tm5_obs_1':'N02',
emin_ev=[ 300, 300, 600, 2300, 200, 300, 5000] 'tm5_obs_2':'N03',
emax_ev=[2300, 600, 2300, 5000, 10000,8000, 8000] 'tm7_obs_1':'N15',
'tm7_obs_2':'N16',
'tm6_obs_1':'N12',
'tm6_obs_2_badpix':'N13',
'tm6_obs_3_badpix':'N14',
'tm6_park_1':'N04',
'tm6_park_2':'N06',
'tm6_park_3':'N08',
'tm6_park_4':'N10',
'tm6_scan_1':'N05',
'tm6_scan_2':'N07',
'tm6_scan_3':'N09',
'tm6_scan_4':'N11'}
emin_kev=[0.3, 0.3, 0.6, 2.3, 0.2, 0.3, 5.0] """ Это просто индекс диапазона для выходных файлов. """
emax_kev=[2.3, 0.6, 2.3, 5.0, 10.0, 8.0, 8.0] eband=[ "0", "1", "2", "3", "4", "5", "6", "7", "8"]
""" Диапазоны энергий. """
emin_ev=[ 300, 300, 600, 2300, 200, 300, 5000, 500, 1000]
emax_ev=[2300, 600, 2300, 5000, 10000,8000, 8000, 1000, 2000]
emin_kev=[0.3, 0.3, 0.6, 2.3, 0.2, 0.3, 5.0, 0.5, 1.0]
emax_kev=[2.3, 0.6, 2.3, 5.0, 10.0, 8.0, 8.0, 1.0, 2.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, 3.8468E+12] # something is wrong here
#ecf = [9.7817E+11, 3.2982E+12, 1.3903E+12, 2.3322E+12, 5.2022E+11, 5.8453E+11, 3.8468E+12]
""" """
*** en0 ecf 9.7817E+11 +/- 2.4606E+10 2.52% N=17 *** en0 ecf 9.7817E+11 +/- 2.4606E+10 2.52% N=17
*** en1 ecf 3.2982E+12 +/- 8.2963E+10 2.52% N=17 *** en1 ecf 3.2982E+12 +/- 8.2963E+10 2.52% N=17
@ -64,6 +82,22 @@ ecf = [9.7817E+11, 3.2982E+12, 1.3903E+12, 2.3322E+12, 5.2022E+11, 5.8453E+1
*** en5 ecf 5.8453E+11 +/- 1.4743E+10 2.52% N=17 *** en5 ecf 5.8453E+11 +/- 1.4743E+10 2.52% N=17
""" """
# finally use Pavel's file ../data/ECF/ecf_tbabspow_g2nh0.02.pkl
"""
for e in [(0.3, 2.3), (0.3, 0.6), (0.6, 2.3), (2.3, 5.0), (5.0, 8.0)]:
print(f'{ecf[e]:g}')
"""
ecf = [1.0911e+12, # (0.3, 2.3)
1.07252e+12, # (0.3, 0.6)
1.08963e+12, # (0.6, 2.3)
1.14674e+11, # (2.3, 5.0)
1.0,
1.0,
2.77581e+10, # (5.0, 8.0)
1354632916123.6191, # (0.5, 1.0) 4XMM-DR12 EP2 band
1014764099304.4968] # (1.0, 2.0) 4XMM-DR12 EP3 band
outfile_post='.fits' outfile_post='.fits'
""" """

154
uds/uds/sherpa.py
View File

@ -35,8 +35,8 @@ def calc_ecf(path,catprep=None, emin=None, emax=None, eband=None, outfile=None,
srcid=rec['detUID'] srcid=rec['detUID']
""" consider only non-extended sources """ """ consider only non-extended sources """
if(rec['EXT'] > 0.0): #if(rec['EXT'] > 0.0):
continue # continue
if(srcid in mlrate.keys()): if(srcid in mlrate.keys()):
print("Duplicate SrcID! {}".format(srcid)) print("Duplicate SrcID! {}".format(srcid))
@ -52,7 +52,7 @@ def calc_ecf(path,catprep=None, emin=None, emax=None, eband=None, outfile=None,
print("Reading {}".format(filename)) print("Reading {}".format(filename))
fout=filename.replace(".fits", ".uncorr.fits") #fout=filename.replace(".fits", ".uncorr.fits")
pha=filename.replace("ARF", "SourceSpec") pha=filename.replace("ARF", "SourceSpec")
pha=pha.replace(".fits", ".fits.grp") pha=pha.replace(".fits", ".fits.grp")
@ -69,10 +69,10 @@ def calc_ecf(path,catprep=None, emin=None, emax=None, eband=None, outfile=None,
continue continue
hdul = fits.open(filename,mode='readonly') #hdul = fits.open(filename,mode='readonly')
specresp = np.divide(hdul[1].data['SPECRESP'], hdul[1].data['CORRCOMB']) #specresp = np.divide(hdul[1].data['SPECRESP'], hdul[1].data['CORRCOMB'])
hdul[1].data['SPECRESP'] = specresp #hdul[1].data['SPECRESP'] = specresp
hdul.writeto(fout, overwrite=True) #hdul.writeto(fout, overwrite=True)
clean() # sherpa clean() # sherpa
set_xsxsect("vern") set_xsxsect("vern")
set_xsabund('wilm') set_xsabund('wilm')
@ -165,6 +165,144 @@ def calc_ecf(path,catprep=None, emin=None, emax=None, eband=None, outfile=None,
with open(outfile, 'wb') as f: with open(outfile, 'wb') as f:
pickle.dump(data, f) pickle.dump(data, f)
def calc_flux(path,catprep=None, emin=None, emax=None, eband=None, outfile=None, fitmin=0.2, fitmax=5.0, simnum=100):
""" """
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))
pha=filename.replace("ARF", "SourceSpec")
pha=pha.replace(".fits", ".fits.grp")
pha_header = fits.getheader(pha,ext=1)
object = pha_header['OBJECT']
srcid="ML{}".format(object)
if not (srcid in mlrate.keys()):
print("Skip {}".format(srcid))
continue
clean() # sherpa
set_xsxsect("vern")
set_xsabund('wilm')
load_pha(pha) # 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
photon_flx=[]
energy_flx=[]
skipme=False
for ie in range(len(eband)):
ignore_bad() # sherpa
### Define the fit range
print("Notice {}-{} keV:".format(emin[ie], emax[ie]))
notice(emin[ie], emax[ie]) # 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')
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))
print("----> Success: {} pow norm: {:.4E} reduced stat: {:.2f}".format(res_fit.succeeded,res_fit.parvals[0], rstat))
if (skipme == True):
print("*** Skip zero flux for {}".format(srcid))
continue
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): def print_ecf(infile=None, emin=None, emax=None, eband=None, skipfrac=5.0):
@ -210,6 +348,7 @@ def print_ecf(infile=None, emin=None, emax=None, eband=None, skipfrac=5.0):
continue continue
ecf.append(ecf0) ecf.append(ecf0)
ecf_mean = np.mean(ecf) ecf_mean = np.mean(ecf)
ecf_err = np.std(ecf, ddof=1)/np.sqrt(np.size(ecf)) ecf_err = np.std(ecf, ddof=1)/np.sqrt(np.size(ecf))
@ -217,3 +356,4 @@ def print_ecf(infile=None, emin=None, emax=None, eband=None, skipfrac=5.0):
ecf_mean, ecf_mean,
ecf_err,100*ecf_err/ecf_mean, ecf_err,100*ecf_err/ecf_mean,
np.size(ecf))) np.size(ecf)))
sys.exit()

1648
uds/uds/utils.py
View File

File diff suppressed because it is too large Load Diff