kea/sgrb2
1
0
Fork 0
forked from xmm/arches
Code Pull Requests Activity

good night

Browse Source
This commit is contained in:
Roman Krivonos
2025-10-28 00:11:18 +03:00
parent 22c164cd41
commit 777c92ef2e
6 changed files with 526 additions and 70 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
arches/arches/config.py
View File

@@ -23,3 +23,19 @@ Radius=15"
"""
skip=['0862471401','0862471501','0862470501']
# energy bands taken from Tatischeff 2012
emin=[
'2000', # continuum, 2-4 keV
'4170', # continuum, 4.17-5.86 keV (T12)
'6300', # 6.4 keV line, 6.3-6.48 keV (T12)
'6564', # 6.7 keV line, 6.564-6.753 keV (T12)
'8000', # high-energy continuum, 8-12 keV
]
emax=[
'4000',
'5860',
'6480',
'6753',
'12000',
]

270
arches/arches/utils.py
View File

@@ -15,6 +15,7 @@ from pathlib import Path
import pandas
import pickle
import pyds9
import subprocess
from os.path import dirname
import inspect
@@ -27,6 +28,7 @@ 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
from astropy.time import Time, TimeDelta, TimezoneInfo, TimeFromEpoch
from astropy.convolution import interpolate_replace_nans
import statistics
import shutil
@@ -63,6 +65,12 @@ def remove_file(filename):
def move_files(pattern, destination):
for file in glob.glob(pattern):
shutil.move(file, os.path.join(destination,file))
def test_file(fn):
if not (os.path.isfile(fn)==True):
print(f'requested filename {fn} is not found')
cwd = os.getcwd()
print(f"Current directory: {cwd}")
sys.exit()
def get_first_file(pattern):
files = glob.glob(pattern)
@@ -541,3 +549,265 @@ def get_ds9_regions(image, src_fn, bkg_fn):
r2_bkg = regionbkg.split(",")[3].replace('\n','')
print(" r2_bkg = ", r2_bkg, "(physical, annulus)")
return x_source,y_source,r_source,x_bkg,y_bkg,r_bkg,r2_bkg
def run_evqpb(key, work_dir=None, evtfile=None, attfile=None, pimin=[200,], pimax=[12000,], label=['',], exposurefactor=10.0):
test_file(evtfile)
test_file(attfile)
# Run the SAS task evqpb over each one of _original_ event files.
evqpb_outset=work_dir+f'/EPIC_{key}_QPB.fits' # Name of the output file
inargs = [f'table={evtfile}', f'attfile={attfile}',
f'outset={evqpb_outset}',
f'exposurefactor={exposurefactor}',]
w("evqpb", inargs).run()
gti=f'EPIC_{key}_gti.fit'
test_file(gti)
expression=f'gti({gti},TIME)'
# filter the EPIC event files to create cleaned and filtered for
# flaring particle background event files for your observation
sci_clean=work_dir+f'/EPIC_{key}_filtered.fits' # Name of the output file
inargs = [f'table={evtfile}', 'withfilteredset=yes',
f'filteredset={sci_clean}',
'destruct=yes','keepfilteroutput=true',
f'expression={expression}']
w("evselect", inargs).run()
# Use the same expression to filter the FWC event files:
qpb_clean=work_dir+f'/EPIC_{key}_QPB_clean.fits' # Name of the output file
inargs = [f'table={evqpb_outset}', 'withfilteredset=yes',
f'filteredset={qpb_clean}',
'destruct=yes','keepfilteroutput=true',
f'expression={expression}']
w("evselect", inargs).run()
for idx, lab in enumerate(label):
if(key=='PN'):
expression=f'(#XMMEA_EP)&&(PATTERN<=4)&&(PI>={pimin[idx]})&&(PI<={pimax[idx]})'
else:
expression=f'(#XMMEA_EM)&&(PATTERN<=12)&&(PI>={pimin[idx]})&&(PI<={pimax[idx]})'
print(lab,expression)
sci_image=work_dir+f'/EPIC_{key}_sci_image{lab}.fits' # Name of the output file
# Extract an image for the science exposure
inargs = [f'table={sci_clean}', f'expression={expression}',
'imagebinning=binSize', f'imageset={sci_image}',
'withimageset=yes','xcolumn=X', 'ycolumn=Y',
'ximagebinsize=80', 'yimagebinsize=80',
]
w("evselect", inargs).run()
qpb_image=work_dir+f'/EPIC_{key}_qpb_image{lab}.fits' # Name of the output file
# Extract an image for the FWC exposure
inargs = [f'table={qpb_clean}', f'expression={expression}',
'imagebinning=binSize', f'imageset={qpb_image}',
'withimageset=yes','xcolumn=X', 'ycolumn=Y',
'ximagebinsize=80', 'yimagebinsize=80','zcolumn=EWEIGHT',
]
w("evselect", inargs).run()
#hdul = fits.open(qpb_image)
#primary_hdu = hdul[0]
#qpb_data = hdul[0].data
#qpb_header = hdul[0].header
#hdul.close()
#result = interpolate_replace_nans(image, kernel)
# -- Subtract background --
# skip, due to low statistics
"""
cor_image=work_dir+f'/EPIC_{key}_cor_image{lab}.fits' # Name of the output file
cmd = ['farith', f'{sci_image}', f'{qpb_image}', f'{cor_image}', 'SUB',
'copyprime=yes',
'clobber=yes']
result = subprocess.run(cmd, capture_output=True, text=True, check=True)
"""
def run_mosaic(outfile_cts=None,outfile_qpb=None,outfile_exp=None,
outfile_sub=None,outfile_flx=None,outfile_pix=None,
cts=None,qpb=None,exp=None,nn=3,devmax=5.0,cutbox=100):
ref_crd = SkyCoord(arches_ra, arches_dec, frame=FK5(), unit="deg")
nn2=nn*nn
print(outfile_cts)
print(cts[0])
print(exp[0])
# take first image as reference
ref_hdul = fits.open(cts[0])
ref_data = ref_hdul[0].data
ref_header = ref_hdul[0].header
ref_hdul.close()
ref_wcs = WCS(ref_header)
(nx,ny) = ref_data.shape
map_flx = np.zeros(ref_data.shape,dtype=np.float64)
map_cts = np.zeros(ref_data.shape,dtype=np.float64)
map_qpb = np.zeros(ref_data.shape,dtype=np.float64)
map_exp = np.zeros(ref_data.shape,dtype=np.float64)
map_pix = np.zeros(ref_data.shape,dtype=np.float64) # number of subpixels inside
# take reference exposure
exp_hdul = fits.open(exp[0])
exp_data = exp_hdul[0].data
exp_header = exp_hdul[0].header
exp_hdul.close()
exp_wcs = WCS(exp_header)
"""
for row_index in range(len(ref_data)):
for col_index in range(len(ref_data[row_index])):
refx, refy = ref_wcs.world_to_pixel(ref_crd)
if(np.absolute(refx-row_index)>cutbox):
continue
if(np.absolute(refy-col_index)>cutbox):
continue
sky = ref_wcs.pixel_to_world(row_index, col_index)
exp_xx, exp_yy = exp_wcs.world_to_pixel(sky)
exp_x = int(np.round(exp_xx))
exp_y = int(np.round(exp_yy))
if not (exp_data[exp_y, exp_x]>0):
continue
map_cts[col_index][row_index]=1.0
sep_arcmin = sky.separation(ref_crd).arcmin
if(sep_arcmin > devmax):
continue
map_cts[col_index][row_index]=2.0
ref_hdul[0].data=map_cts
ref_hdul.writeto(outfile_cts,overwrite=True)
sys.exit()
"""
for idx,in_cts in enumerate(cts):
#if not (idx==2):
# continue
hdul0 = fits.open(in_cts)
c_data = hdul0[0].data
c_header = hdul0[0].header
c_wcs = WCS(c_header)
hdul0.close()
print(f"{in_cts} filter {c_header['FILTER']}")
hdul0 = fits.open(qpb[idx])
q_data = hdul0[0].data
q_header = hdul0[0].header
q_wcs = WCS(q_header)
hdul0.close()
hdul0 = fits.open(exp[idx])
e_data = hdul0[0].data
e_header = hdul0[0].header
e_wcs = WCS(e_header)
hdul0.close()
print(exp[idx])
xx, yy = c_wcs.world_to_pixel(ref_crd)
ref_x = int(np.rint(xx))
ref_y = int(np.rint(yy))
xmin = ref_x - cutbox
xmax = ref_x + cutbox
ymin = ref_y - cutbox
ymax = ref_y + cutbox
if(xmin<0):
xmin=0
if(xmax>=nx):
xmax=nx
if(ymin<0):
ymin=0
if(ymax>=ny):
ymax=ny
#for row_index in range(len(c_data)):
# for col_index in range(len(c_data[row_index])):
for row_index in range(xmin,xmax):
for col_index in range(ymin,ymax):
c_sky = c_wcs.pixel_to_world(row_index, col_index)
exp_xx, exp_yy = e_wcs.world_to_pixel(c_sky)
exp_x = int(np.rint(np.float64(exp_xx)))
exp_y = int(np.rint(np.float64(exp_yy)))
if not (e_data[exp_y, exp_x]>0):
continue
sep_arcmin = c_sky.separation(ref_crd).arcmin
if(sep_arcmin > devmax):
continue
#xx, yy = ref_wcs.world_to_pixel(c_sky)
#refx = int(np.round(xx))
#refy = int(np.round(yy))
# do subpixeling
for ii in range(nn):
for jj in range(nn):
xx0=np.float64(row_index)-0.5+(np.float64(ii+1)-0.5)/np.float64(nn)
yy0=np.float64(col_index)-0.5+(np.float64(jj+1)-0.5)/np.float64(nn)
subpix_sky = c_wcs.pixel_to_world(xx0, yy0)
sub_xx, sub_yy = ref_wcs.world_to_pixel(subpix_sky)
refx = np.int32(np.rint(np.float64(sub_xx)))
refy = np.int32(np.rint(np.float64(sub_yy)))
exp_xx, exp_yy = e_wcs.world_to_pixel(subpix_sky)
exp_x = np.int32(np.rint(np.float64(exp_xx)))
exp_y = np.int32(np.rint(np.float64(exp_yy)))
if (refx>=0 and refx<nx and refy>=0 and refy<ny):
map_cts[refy][refx]+=float(c_data[col_index][row_index])/nn2
map_qpb[refy][refx]+=float(q_data[col_index][row_index])/nn2
map_pix[refy][refx]+=1.0/nn2
map_exp[refy][refx]+=np.float64(e_data[exp_y][exp_x])/nn2
# Convert summed sky back to normal units
for x in range(nx):
for y in range(ny):
if(map_exp[y][x]>0):
map_flx[y][x]=(map_cts[y][x]-map_qpb[y][x])/map_exp[y][x]
map_cts[y][x]=map_cts[world_to_pixely][x]/map_pix[y][x]
map_qpb[y][x]=map_qpb[y][x]/map_pix[y][x]
map_exp[y][x]=map_exp[y][x]/map_pix[y][x]
else:
map_flx[y][x]=0.0
if(outfile_pix):
ref_hdul[0].data=map_pix
ref_hdul.writeto(outfile_pix,overwrite=True)
ref_hdul[0].data=map_cts
ref_hdul.writeto(outfile_cts,overwrite=True)
ref_hdul[0].data=map_exp
ref_hdul.writeto(outfile_exp,overwrite=True)
ref_hdul[0].data=map_qpb
ref_hdul.writeto(outfile_qpb,overwrite=True)
ref_hdul[0].data=map_flx
ref_hdul.writeto(outfile_flx,overwrite=True)
#sys.exit()
"""
cmd = ['farith', f'{outfile_cts}', f'{outfile_qpb}', f'{outfile_sub}', 'SUB',
'copyprime=yes', 'clobber=yes']
result = subprocess.run(cmd, capture_output=True, text=True, check=True)
cmd = ['farith', f'{outfile_sub}', f'{outfile_exp}', f'{outfile_rat}', 'DIV',
'copyprime=yes', 'clobber=yes']
result = subprocess.run(cmd, capture_output=True, text=True, check=True)
"""