Initial commit

This commit is contained in:
2024-04-13 14:35:10 +03:00
commit 9aac031a01
48 changed files with 65971 additions and 0 deletions

9
uds/setup.py Normal file
View File

@@ -0,0 +1,9 @@
from setuptools import setup, find_packages
setup(name="uds",
version="0.1",
description='eROSITA UDS field analysis',
author='Roman Krivonos',
author_email='krivonos@cosmos.ru',
url='https://www.srg.cosmos.ru',
packages=find_packages())

0
uds/uds/__init__.py Normal file
View File

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

109
uds/uds/config.py Normal file
View File

@@ -0,0 +1,109 @@
from pathlib import Path
work_dir = Path('/path/to/some/logical/parent/dir')
"""
Координаты сентрального кадра, к которому будут
приводиться изображения всех списков событий
"""
ra_cen=34.5342131
de_cen=-4.7956710
"""
Словарь камер со списком наблюдений каждой камеры.
Номера камер должны быть отсортированы
"""
keylist_tm={'1':['tm1_obs_1',],
'5':['tm5_obs_1','tm5_obs_2',],
'6':['tm6_obs_1','tm6_obs_2_badpix','tm6_obs_3_badpix',
'tm6_park_1','tm6_park_2','tm6_park_3','tm6_park_4',
'tm6_scan_1','tm6_scan_2','tm6_scan_3','tm6_scan_4'],
'7':['tm7_obs_1','tm7_obs_2',]}
"""
Примерные центры изображений каждого наблюдения.
Требуется для астрокоррекции. Будет расчитываться матрица сдвигов и поворотов,
так вот, повороты будут проводиться вокруг данных координат.
"""
wcslist={'tm1_obs_1':[34.7279760,-5.0680267],
'tm5_obs_1':[34.7351248,-4.4407314],
'tm5_obs_2':[34.8748997,-4.4871658],
'tm7_obs_1':[35.0015120,-4.7602124],
'tm7_obs_2':[34.9810029,-4.5915912],
'tm6_obs_1':[34.4227062,-4.7207170],
'tm6_obs_2_badpix':[34.7272339,-4.4294153],
'tm6_obs_3_badpix':[34.8750291,-4.4708468],
'tm6_park_1':[35.0544951,-4.0613619],
'tm6_park_2':[35.0558675,-4.0683084],
'tm6_park_3':[35.0565263,-4.0583538],
'tm6_park_4':[35.0602986,-4.0622220],
'tm6_scan_1':[34.5405596,-4.8088748],
'tm6_scan_2':[34.5405596,-4.8088748],
'tm6_scan_3':[34.5405596,-4.8088748],
'tm6_scan_4':[34.5405596,-4.8088748]}
""" like in the paper (Table 1) """
obslist={'tm1_obs_1':'N01',
'tm5_obs_1':'N02',
'tm5_obs_2':'N03',
'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'}
""" Это просто индекс диапазона для выходных файлов. """
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]
# 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
*** 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
"""
# 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'
"""
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
"""

359
uds/uds/sherpa.py Normal file
View File

@@ -0,0 +1,359 @@
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 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):
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)))
sys.exit()

3506
uds/uds/utils.py Normal file

File diff suppressed because it is too large Load Diff