initial project code commit

This commit is contained in:
2026-04-22 13:15:55 +03:00
parent e624299842
commit 9dcc2f9473
69 changed files with 7387 additions and 0 deletions

BIN
artproducts/src/.DS_Store vendored Normal file

Binary file not shown.

BIN
artproducts/src/artproducts/.DS_Store vendored Normal file

Binary file not shown.

View File

@@ -0,0 +1,284 @@
#!/usr/bin/env python3
# -*- coding: utf8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2021-2026 IKI RAS (http://iki.rssi.ru/)
# Space Research Institute of the Russian Academy of Science
#
# This file is part of the M. Pavlinsky SRG/ART-XC software project.
# -----------------------------------------------------------------------------
import os
import numpy as np
from regions import Regions
from astropy.io import fits
from arttools.telescope import Teldef
from artproducts.version import __version__
class EnergySpectrum(object):
def __init__(self, events, make_individual=False):
self._events = events
self._start = events.tstart
self._stop = events.tstop
self._phadata = {}
self._phameta = {}
self._make_individual = make_individual
def create(self, region):
all_phadata = {
'pha' : [],
'bins' : [],
'edges': [],
'gti' : []
}
self._phameta['all'] = {
'exposure': 0.0,
'deadtime': 0.0,
'livetime': 0.0,
'events' : 0,
'selected': 0,
'filter' : None
}
for telname, evtdata in self._events:
phadata, phameta = self._make_pha(evtdata, region)
self._phadata[telname] = phadata
self._phameta[telname] = phameta
for key in ['pha', 'bins', 'edges', 'gti']:
all_phadata[key].append(phadata[key])
for key in ['exposure', 'deadtime', 'livetime', 'events', 'selected']:
self._phameta['all'][key] += phameta[key]
self._phameta['all']['projection'] = phameta['projection']
self._phameta['all']['filter' ] = phameta['filter' ]
all_bins = all_phadata['bins' ][0]
all_edges = all_phadata['edges'][0]
all_pha = np.sum(all_phadata['pha'], axis=0)
all_gti = np.prod(all_phadata['gti'])
self._phadata['all'] = {
'pha' : all_pha ,
'bins' : all_bins ,
'edges': all_edges,
'gti' : all_gti
}
def _make_pha(self, evtdata, region):
data = evtdata['events' ]
projection = evtdata['projection']
gti = data.gti
exp, deadt, livet \
= self._make_timekeys(data, gti)
mask, filter_expr \
= self._make_mask(data, region, projection)
times = np.array(data.events.column('TIME'))
bins = np.linspace(0.1, 65., 650)
pha, edges = np.histogram(times[mask], bins=bins)
phadata = {
'pha' : pha ,
'bins' : bins ,
'edges': edges,
'gti' : gti
}
phameta = {
'projection': projection,
'exposure' : exp ,
'deadtime' : deadt,
'livetime' : livet,
'events' : len(times),
'selected' : len(times[mask]),
'filter' : filter_expr
}
return phadata, phameta
def _make_timekeys(self, evtdata, gti):
exposure = gti.exposure
deadtime = 0
for gtistart, gtistop in gti.array:
deadtime += evtdata.deadtime(gtistart, gtistop)
livetime = exposure - deadtime
return exposure, deadtime, livetime
def _make_goodevt_mask(self, evtdata, grade=9, flag=0):
maskexpr = '(GRADE < {grade}) && (FLAG == {flag})'
maskexpr = maskexpr.format(grade=grade, flag=flag)
grademask = evtdata.events.column('GRADE') < grade
flagmask = evtdata.events.column('FLAG' ) == flag
return np.all([grademask, flagmask], axis=0), maskexpr
def _make_mask(self, evtdata, region, projection):
goodmask, goodexpr \
= self._make_goodevt_mask(evtdata)
regionmask = region.mask(projection,
np.array(evtdata.events.column('RA' )),
np.array(evtdata.events.column('DEC')))
expr = '{}'.format(goodexpr)
return np.all([goodmask, regionmask], axis=0), expr
def _write_pha(self, name, dstfname, arffile=None, rmffile=None, backfile=None, backscale=None):
if not name in self._phadata:
return False
mjdref = 51543.875
phadata = self._phadata[name]
phameta = self._phameta[name]
hdu = fits.PrimaryHDU()
a1 = np.arange(649)
a2 = np.array(phadata['pha'])
a3 = np.zeros(a1.size)
a4 = np.ones (a1.size)
c1 = fits.Column(name='CHANNEL' , format='1J', unit='' , array=a1)
c2 = fits.Column(name='COUNTS' , format='1J', unit='count', array=a2)
c3 = fits.Column(name='QUALITY' , format='1I', unit='' , array=a3)
c4 = fits.Column(name='GROUPING', format='1I', unit='' , array=a4)
cols = fits.ColDefs([c1, c2, c3, c4])
hdu1 = fits.BinTableHDU.from_columns(cols)
hdu1.name = 'SPECTRUM'
hdu1.header['MISSION' ] = 'SRG'
hdu1.header['TELESCOP'] = 'ART-XC'
hdu1.header['INSTRUME'] = name
hdu1.header['EXPOSURE'] = phameta['exposure'] - phameta['deadtime']
hdu1.header['AREASCAL'] = 1.
if backfile is not None:
hdu1.header['BACKFILE'] = backfile
if backscale is not None:
hdu1.header['BACKSCAL'] = backscale
else:
hdu1.header['BACKSCAL'] = 1.
hdu1.header['CORRFILE'] = ''
hdu1.header['CORRSCAL'] = -1
if arffile is not None:
hdu1.header['ANCRFILE'] = arffile
if rmffile is not None:
hdu1.header['RESPFILE'] = rmffile
hdu1.header['HDUCLASS'] = 'OGIP'
hdu1.header['HDUCLAS1'] = 'SPECTRUM'
hdu1.header['HDUVERS' ] = '1.2.1'
hdu1.header['POISSERR'] = True
hdu1.header['CHANTYPE'] = 'PI'
hdu1.header['DETCHANS'] = 649
hdu1.header['TLMIN1' ] = 0
hdu1.header['TLMAX1' ] = 648
hdu2 = phadata['gti'].make_hdu(mjdref)
hdulist = fits.HDUList([hdu, hdu1, hdu2])
hdulist.writeto(dstfname, overwrite=True, checksum=True)
return True
def write(self, dstfname, respfiles=None, backfiles=None, backscale=None):
written = {
'files': []
}
allfname = dstfname.format('_all')
backname = None
if backfiles is not None:
if 'all' in backfiles:
backname = backfiles['all']
arffile = None
rmffile = None
if respfiles is not None:
if 'all' in respfiles['arf']:
arffile = respfiles['arf']['all']
if 'all' in respfiles['rmf']:
rmffile = respfiles['rmf']['all']
# def _write_pha(self, name, dstfname, arffile=None, respfile=None, backfile=None, backscale=None):
if self._write_pha(
'all',
allfname,
arffile=arffile,
rmffile=rmffile,
backfile=backname,
backscale=backscale):
written['all' ] = os.path.basename(allfname)
written['files'].append(allfname)
if not self._make_individual:
return written
for telname in Teldef.telnames:
print(telname)
telfname = dstfname.format(telname)
backname = None
if backfiles is not None:
if telname in backfiles:
backname = backfiles[telname]
arffile = None
rmffile = None
if respfiles is not None:
if telname in respfiles['arf']:
arffile = respfiles['arf'][telname]
if telname in respfiles['rmf']:
rmffile = respfiles['rmf'][telname]
if self._write_pha(
telname,
telfname,
arffile=arffile,
rmffile=rmffile,
backfile=backname,
backscale=backscale):
written[telname] = os.path.basename(telfname)
written['files'].append(telfname)
return written
if __name__ == '__main__':
pass

View File

@@ -0,0 +1,128 @@
#!/usr/bin/env python3
# -*- coding: utf8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2026 IKI RAS (http://iki.rssi.ru/)
# Space Research Institute of the Russian Academy of Science
#
# This file is part of the M. Pavlinsky SRG/ART-XC software project.
# -----------------------------------------------------------------------------
import numpy as np
from arttools.wcs import Wcs
from arttools.events import Events
from arttools.telescope import Teldef
class Calculate(object):
def __init__(self, method, parameters, datafields=None):
self._method = method
self._parameters = parameters
self._datafields = datafields
self._result = {
'all': 0
}
@property
def result(self):
return self._result
def apply(self, telname, data):
params = [data[p] for p in self._parameters]
value = self._method(*params)
self._result[telname] = value
self._result['all' ] += value
class EventsPack(object):
def __init__(self, evtdir):
self._evtdir = evtdir
self._epoch0 = 0
self._tstart = 0
self._tstop = 0
self._count = 0
self._data = {}
def __iter__(self):
return iter(self._data.items())
@property
def empty(self):
return not self._data
@property
def count(self):
return self._count
@property
def epoch0(self):
return self._epoch0
@property
def tstart(self):
return self._tstart
@property
def tstop(self):
return self._tstop
def _event_files(self, stem):
for telname in Teldef.telnames:
evtfile = '{}{}_cl.fits'.format(stem, telname)
yield self._evtdir.pathfor(evtfile)
def read(self, stem, epoch0):
starts = []
stops = []
count = 0
for eventfile in self._event_files(stem):
events = Events(eventfile).read()
if events.empty:
continue
starts.append(events.tstart)
stops.append(events.tstop)
count += events.size
projection = Wcs.from_eventfile(eventfile)
self._data[events.telname] = {
'events' : events,
'projection': projection.get()
}
self._tstart = np.max(starts)
self._tstop = np.min(stops )
if epoch0 is not None:
self._epoch0 = epoch0
else:
self._epoch0 = self._tstart
self._count = count
return self
def call(self, callee):
for telname in Teldef.telnames:
if telname not in self._data:
continue
callee.apply(telname, self._data[telname])
return callee.result
if __name__ == '__main__':
pass

View File

@@ -0,0 +1,117 @@
#!/usr/bin/env python3
# -*- coding: utf8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2021-2026 IKI RAS (http://iki.rssi.ru/)
# Space Research Institute of the Russian Academy of Science
#
# This file is part of the M. Pavlinsky SRG/ART-XC software project.
# -----------------------------------------------------------------------------
import warnings
warnings.simplefilter('ignore')
import numpy as np
from astropy.io import fits
from arttools.wcs import Wcs
from arttools.fstools import Dir
from arttools.gti import Gti
from arttools.events import Events
from artproducts.version import __version__
class EventsWcs(object):
def __init__(self, srcfile, wcs):
self._wcs = wcs
self._projection = Wcs.from_file(wcs)
self._eventfile = Events(srcfile)
@property
def file(self):
return self._eventfile
@property
def tstart(self):
return self._eventfile.tstart
@property
def tstop(self):
return self._eventfile.tstop
@property
def gti(self):
return self._eventfile.gti
def deadtime(self, start, stop):
return self._eventfile.deadtime(start, stop)
def _wcs_mask(self, x, y):
xsize = int(self._projection.get().wcs.crpix[0]*2 + 1)
ysize = int(self._projection.get().wcs.crpix[1]*2 + 1)
asize = x * ysize + xsize
wcs_mask = np.all([x > 0, y > 0, asize < xsize * ysize], axis=0)
return wcs_mask
def prepare(self, fitwcs=False, usergti=None):
self._eventfile.read()
ra = np.array(self._eventfile.events.column('RA' ))
dec = np.array(self._eventfile.events.column('DEC' ))
times = np.array(self._eventfile.events.column('TIME'))
evt_gti = self._eventfile.gti
self._x = np.zeros_like(ra )
self._y = np.zeros_like(dec)
self._x, self._y \
= self._projection.world2pix(ra, dec)
if fitwcs:
wcs_mask = self._wcs_mask(self._x, self._y)
else:
wcs_mask = np.full_like(times, True)
wcs_gti = Gti.from_mask(times, wcs_mask).normalize(edge_delta=10.)
gti = evt_gti * wcs_gti
if usergti is not None:
self._gti = gti * Gti.from_hduname(usergti, 'stdgti')
else:
self._gti = gti
self._mask = self._gti.mask(times)
def read(self):
self._eventfile.read()
def write(self, dstfile, creator=''):
pos = 13
keys = ['X', 'Y']
layout = {
'X': {'format': '1D', 'unit': ''},
'Y': {'format': '1D', 'unit': ''},
}
data = {
'X': self._x,
'Y': self._y
}
self._eventfile.modify_and_write_cl(dstfile, creator, keys, layout, data, pos, self._mask, self._gti)
with fits.open(dstfile, mode='update') as dsthdul:
dsthdul['EVENTS'].header.update(self._projection.get().to_header())
dsthdul['EVENTS'].header['MJDREF'] = (self._eventfile.mjdref, 'MJD corresponding to SC clock start 2000.0 MSK')
if __name__ == '__main__':
pass

View File

@@ -0,0 +1,26 @@
#!/usr/bin/env python3
# -*- coding: utf8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2021-2026 IKI RAS (http://iki.rssi.ru/)
# Space Research Institute of the Russian Academy of Science
#
# This file is part of the M. Pavlinsky SRG/ART-XC software project.
# -----------------------------------------------------------------------------
class Hist(object):
def __init__(self, vstart, vstop):
pass
class TimeHist(object):
def __init__(self, epoch0, start, stop, binsize):
pass
if __name__ == '__main__':
pass

View File

@@ -0,0 +1,149 @@
#!/usr/bin/env python3
# -*- coding: utf8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2021-2026 IKI RAS (http://iki.rssi.ru/)
# Space Research Institute of the Russian Academy of Science
#
# This file is part of the M. Pavlinsky SRG/ART-XC software project.
# -----------------------------------------------------------------------------
# import warnings
# warnings.simplefilter('ignore')
import numpy as np
from astropy.io import fits
from arttools.wcs import Wcs
from arttools.fstools import Dir
from arttools.events import Events
from artproducts.version import __version__
from artproducts.eventwcs import EventWcs
class Image(object):
def __init__(self, eventfiles, wcs):
self._wcs = wcs
self._projection = Wcs.from_file(wcs)
self._eventfiles = eventfiles
self._elow = 0
self._ehigh = 0
self._mjdref = 0
self._exposure = 0
self._deadtime = 0
self._tstart = 0
self._tstop = 0
self._onsource = 0
def make(self, elow, ehigh):
self._image = self._make_image()
self._elow = elow
self._ehigh = ehigh
image_exp = []
image_deadt = []
image_start = []
image_stop = []
for evtfile in self._eventfiles:
events = EventWcs(evtfile, self._wcs)
events.read()
self._mjdref = events.file.mjdref
self._image = self._events_to_image(events, elow, ehigh, self._image)
start = events.tstart
stop = events.tstop
image_exp.append(events.gti.exposure)
image_deadt.append(events.deadtime(start, stop))
image_start.append(start)
image_stop.append(stop)
self._exposure = np.sum(image_exp )
self._deadtime = np.sum(image_deadt)
self._tstart = np.min(image_start)
self._tstop = np.max(image_stop )
self._onsource = self._tstop - self._tstart
print('onsource [sec] = {}'.format(self._onsource))
print('exposure [sec] = {} (average per head {})'.format(self._exposure, self._exposure / 7.))
print('deadtime [sec] = {} (average per head {})'.format(self._deadtime, self._deadtime / 7.))
print('tstart [sec] = {}'.format(self._tstart ))
print('tstop [sec] = {}'.format(self._tstop ))
def write(self, dstfile):
fits.PrimaryHDU(
data=self._image,
header=self._projection.get().to_header()
).writeto(dstfile, overwrite=True, checksum=True)
gtistart = self._tstart
gtistop = self._tstop
deadtime = self._deadtime
exposure = self._exposure - self._deadtime
onsource = self._onsource
with fits.open(dstfile, mode='update') as hdul:
hdul[0].header['MJDREF' ] = (self._mjdref, '[d] MJDREF')
hdul[0].header['TSTART' ] = (gtistart , '[sec]' )
hdul[0].header['TSTOP' ] = (gtistop , '[sec]' )
hdul[0].header['ONSOURCE'] = (onsource , '[sec]' )
hdul[0].header['EXPOSURE'] = (exposure , '[sec]' )
hdul[0].header['E_LOW' ] = (self._elow , '[keV]' )
hdul[0].header['E_HIGH' ] = (self._ehigh , '[keV]' )
print('written: {}'.format(dstfile))
def _make_mask(self, events, elow, ehigh):
en = events.file.events.column('ENERGY')
grade = events.file.events.column('GRADE' )
flag = events.file.events.column('FLAG' )
mask = np.all([
en >= elow ,
en <= ehigh,
grade < 9 ,
flag == 0
], axis=0)
return mask
def _make_image(self):
xsize = int(self._projection.get().wcs.crpix[0]*2 + 1)
ysize = int(self._projection.get().wcs.crpix[1]*2 + 1)
image = np.zeros((ysize, xsize), dtype=np.double)
return image
def _events_to_image(self, events, elow, ehigh, image):
times = np.array(events.file.events.column('TIME'))
ra = np.array(events.file.events.column('RA' ))
dec = np.array(events.file.events.column('DEC' ))
xsize = int(self._projection.get().wcs.crpix[0]*2 + 1)
ysize = int(self._projection.get().wcs.crpix[1]*2 + 1)
x, y = self._projection.world2pix(ra, dec)
asize = x * ysize + xsize
xy_mask = np.all([x > 0, y > 0, asize < xsize * ysize], axis=0)
evt_mask = events.gti.mask(times)
img_mask = self._make_mask(events, elow, ehigh)
masklist = [evt_mask, img_mask, xy_mask]
mask = np.all(masklist, axis=0)
np.add.at(image, (x[mask], y[mask]), 1)
return image
if __name__ == '__main__':
pass

View File

@@ -0,0 +1,327 @@
#!/usr/bin/env python3
# -*- coding: utf8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2021-2026 IKI RAS (http://iki.rssi.ru/)
# Space Research Institute of the Russian Academy of Science
#
# This file is part of the M. Pavlinsky SRG/ART-XC software project.
# -----------------------------------------------------------------------------
import os
import numpy as np
from astropy.io import fits
from regions import Regions
from arttools.gti import Gti
from arttools.wcs import Wcs
from arttools.telescope import Teldef
from artproducts.version import __version__
class Lightcurve(object):
def __init__(self, events, binsize, make_individual=False):
self._events = events
self._epoch0 = events.epoch0
self._start = events.tstart
self._stop = events.tstop
self._binsize = binsize
self._lcdata = {}
self._lcmeta = {}
self._make_inividual = make_individual
@property
def meta(self):
return self._lcmeta
@property
def get(self):
return self._lcdata
def create(self, region, elow, ehigh):
all_lcdata = {
'time' : [],
'count' : [],
'error' : [],
'exposure': [],
'flag' : [],
'gti' : []
}
self._lcmeta['all'] = {
'exposure': 0.0,
'deadtime': 0.0,
'livetime': 0.0,
'elow' : elow ,
'ehigh' : ehigh,
'events' : 0,
'selected': 0,
'filter' : None
}
for telname, evtdata in self._events:
lcdata, lcmeta = self._make_lcurve(evtdata, region, elow, ehigh)
self._lcdata[telname] = lcdata
self._lcmeta[telname] = lcmeta
for key in ['time', 'count', 'error', 'exposure', 'flag', 'gti']:
all_lcdata[key].append(lcdata[key])
for key in ['exposure', 'deadtime', 'livetime', 'events', 'selected']:
self._lcmeta['all'][key] += lcmeta[key]
self._lcmeta['all']['filter'] = lcmeta['filter']
all_time = all_lcdata['time'][0]
all_count = np.sum(all_lcdata['count'], axis=0)
all_error = np.sqrt(np.sum(np.power(all_lcdata['error'], 2.), axis=0))
all_exposure = np.sum(all_lcdata['exposure'], axis=0)
all_flag = np.zeros_like(all_lcdata['flag'][0], dtype=int)
one_mask = np.any([ff == 1 for ff in all_lcdata['flag']], axis=0)
two_mask = np.any([ff == 2 for ff in all_lcdata['flag']], axis=0)
all_flag[one_mask] = 1
all_flag[two_mask] = 2
all_gti = np.prod(all_lcdata['gti'])
self._lcdata['all'] = {
'time' : all_time,
'count' : all_count,
'error' : all_error,
'exposure': all_exposure,
'flag' : all_flag,
'gti' : all_gti
}
def _make_bins(self):
dt = self._start - self._epoch0
steps = np.ceil(dt / self._binsize)
t0 = self._epoch0 + self._binsize * steps
edges = np.arange(t0, self._stop, self._binsize)
bins = (edges[ :-1] + edges[ 1: ]) / 2.
return edges, bins
def _make_flag(self, edges, bins, gti):
flag = np.ones_like(bins, dtype=int)
edge_intervals = zip(edges[:-1], edges[1:])
edge_intervals = Gti([*edge_intervals])
index = 0
for interval in edge_intervals.array:
interval_gti = Gti(interval) * gti
flag[index] = 2 if np.diff(interval) > interval_gti.exposure else 0
index += 1
return flag
def _make_lcurve(self, evtdata, region, elow, ehigh):
data = evtdata['events' ]
proj = evtdata['projection']
edges, bins = self._make_bins()
gti = data.gti * Gti([self._start, self._stop])
exp, deadt, livet \
= self._make_timekeys(data, gti)
mask, filter_expr \
= self._make_mask(data, region, proj, elow, ehigh)
times = np.array(data.events.column('TIME'))
counts, _ = np.histogram(times[mask], bins=edges)
errors = np.sqrt(counts)
binexp = np.diff(edges)
bindeadt, _ = data.bindeadtime(edges, bins)
flag = self._make_flag(edges, bins, gti)
exposure = binexp - bindeadt
lcdata = {
'time' : bins,
'count' : counts,
'error' : errors,
'exposure': exposure,
'flag' : flag,
'gti' : gti
}
lcmeta = {
'exposure': exp ,
'deadtime': deadt,
'livetime': livet,
'elow' : elow ,
'ehigh' : ehigh,
'events' : len(times),
'selected': len(times[mask]),
'filter' : filter_expr
}
return lcdata, lcmeta
def _make_timekeys(self, evtdata, gti):
exposure = gti.exposure
deadtime = 0
for gtistart, gtistop in gti.array:
deadtime += evtdata.deadtime(gtistart, gtistop)
livetime = exposure - deadtime
return exposure, deadtime, livetime
def _make_goodevt_mask(self, evtdata, grade=9, flag=0):
maskexpr = '(GRADE < {grade}) && (FLAG == {flag})'
maskexpr = maskexpr.format(grade=grade, flag=flag)
grademask = evtdata.events.column('GRADE') < grade
flagmask = evtdata.events.column('FLAG' ) == flag
return np.all([grademask, flagmask], axis=0), maskexpr
def _make_energy_mask(self, evtdata, elow, ehigh):
maskexpr = '(ENERGY >= {elow} && ENERGY <= {ehigh})'
maskexpr = maskexpr.format(elow=elow, ehigh=ehigh)
return np.all([
evtdata.events.column('ENERGY') >= elow ,
evtdata.events.column('ENERGY') <= ehigh
], axis = 0), maskexpr
def _make_mask(self, evtdata, region, projection, elow, ehigh):
goodmask, goodexpr \
= self._make_goodevt_mask(evtdata)
enmask, enexpr \
= self._make_energy_mask(evtdata, elow, ehigh)
regionmask = region.mask(projection,
np.array(evtdata.events.column('RA' )),
np.array(evtdata.events.column('DEC')))
expr = '{} && {}'.format(goodexpr, enexpr)
return np.all([goodmask, enmask, regionmask], axis=0), expr
def _write_curve(self, name, dstfname, backfile=None):
if not name in self._lcdata:
return False
mjdref = 51543.875
lcdata = self._lcdata[name]
lcmeta = self._lcmeta[name]
time = lcdata['time' ]
rate = lcdata['count'] / lcdata['exposure']
error = lcdata['error'] / lcdata['exposure']
flag = lcdata['flag' ]
gti = lcdata['gti' ]
# # --- DEBUG ---
# import matplotlib.pyplot as plt
# plt.title(name)
# # plt.plot(lcdata['time'], rate)
# # plt.errorbar(lcdata['time'], rate, yerr=error)
# # plt.plot(lcdata['time'], lcdata['count'])
# plt.errorbar(lcdata['time'], lcdata['count'], yerr=lcdata['error'])
# # plt.plot(hkdata['times'][1:], countdiff)
# # plt.plot(edges[1:], fracexp)
# plt.show()
# # ---
hdu = fits.PrimaryHDU()
# --- STDGTI HDU ---
a1 = np.array(time )
a2 = np.array(rate )
a3 = np.array(error)
a4 = np.array(flag )
c1 = fits.Column(name='TIME' , format='1D', unit='sec' , array=a1)
c2 = fits.Column(name='RATE' , format='1E', unit='counts/sec', array=a2)
c3 = fits.Column(name='ERROR', format='1E', unit='' , array=a3)
c4 = fits.Column(name='FLAG' , format='1I', unit='' , array=a4)
cols = fits.ColDefs([c1, c2, c3, c4])
hdu1 = fits.BinTableHDU.from_columns(cols)
hdu1.name = 'RATE'
hdu1.header['ORIGIN' ] = 'IKI RAS'
hdu1.header['CREATOR' ] = 'artproducts v.{}'.format(__version__)
hdu1.header['TELESCOP'] = 'SRG/ART-XC'
hdu1.header['INSTRUME'] = name
hdu1.header['MJDREF' ] = mjdref
hdu1.header['TSTART' ] = (self._start , '[sec]')
hdu1.header['TSTOP' ] = (self._stop , '[sec]')
hdu1.header['EXPOSURE'] = (lcmeta['exposure'], '[sec]')
hdu1.header['DEADTIME'] = (lcmeta['deadtime'], '[sec]')
hdu1.header['LIVETIME'] = (lcmeta['livetime'], '[sec]')
hdu1.header['E_LOW' ] = (lcmeta['elow' ], '[keV]')
hdu1.header['E_HIGH' ] = (lcmeta['ehigh' ], '[keV]')
hdu1.header['BIN' ] = (self._binsize , '[sec]')
hdu1.header['TIMEDEL' ] = (self._binsize , '[sec]')
if backfile is not None:
hdu1.header['BACKFILE'] = backfile
hdu2 = gti.make_hdu(mjdref)
# --- make FITS file ---
hdulist = fits.HDUList([hdu, hdu1, hdu2])
hdulist.writeto(dstfname, overwrite=True, checksum=True)
return True
def write(self, dstfname, backfiles=None):
written = {
'files': []
}
allfname = dstfname.format('_all')
backname = None
if backfiles is not None:
if 'all' in backfiles:
backname = backfiles['all']
if self._write_curve('all', allfname, backname):
written['all' ] = os.path.basename(allfname)
written['files'].append(allfname)
if not self._make_inividual:
return written
for telname in Teldef.telnames:
telfname = dstfname.format(telname)
backname = None
if backfiles is not None:
if telname in backfiles:
backname = backfiles[telname]
if self._write_curve(telname, telfname, backname):
written[telname] = os.path.basename(telfname)
written['files'].append(telfname)
return written
if __name__ == '__main__':
pass

View File

@@ -0,0 +1,486 @@
#!/usr/bin/env python3
# -*- coding: utf8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2021-2026 IKI RAS (http://iki.rssi.ru/)
# Space Research Institute of the Russian Academy of Science
#
# This file is part of the M. Pavlinsky SRG/ART-XC software project.
# -----------------------------------------------------------------------------
import os
import sys
import shutil
import argparse
import numpy as np
from arttools.cmd import Command, ExecutionResult
from arttools.task import (
Task, TaskError, TaskArgs, TaskCmds
)
from arttools.errors import Error
from arttools.fstools import Dir
from arttools.argparse import KvAction
from arttools.errors import RegionsError
from artpipeline.caldb import CaldbInfo
from artpipeline.caldb import CaldbOrientation
from artpipeline.orientation import Orientation
from artproducts.version import __version__
from artproducts.eventswcs import EventsWcs
from artproducts.regions import Region
from artproducts.eventspack import EventsPack, Calculate
from artproducts.responces import Arf
from artproducts.responces import Rmf
from artproducts.image import Image
from artproducts.lightcurve import Lightcurve
from artproducts.enspectrum import EnergySpectrum
class Args(TaskArgs):
def __init__(self):
super().__init__()
def read_env(self):
if not 'ARTCALDB' in os.environ:
raise TaskError('ARTCALDB not found (environment variable not set?)')
self.caldbdir = os.environ['ARTCALDB']
def read_cfg(self, cfgfile):
pass
def read_args(self):
argparser = argparse.ArgumentParser()
argparser.add_argument(
'stem',
action=KvAction,
help='newfile creation stem')
argparser.add_argument(
'srcdir',
action=KvAction,
help='input directory')
argparser.add_argument(
'dstdir',
action=KvAction,
help='output directory')
argparser.add_argument(
'tmpdir',
action=KvAction,
help='temporary directory')
argparser.add_argument(
'--prepare-skytile',
default=False,
action='store_true',
help='input files stem')
argparser.add_argument(
'--fit-to-wcs',
default=False,
action='store_true',
help='process only wcs-bound events')
argparser.add_argument(
'--wcs',
default=None,
action=KvAction,
help='WCS region')
argparser.add_argument(
'--srcra',
default=None,
action=KvAction,
help='source position RA')
argparser.add_argument(
'--srcdec',
default=None,
action=KvAction,
help='source position Dec')
argparser.add_argument(
'--srcreg',
default=None,
action=KvAction,
help='source region')
argparser.add_argument(
'--bkgreg',
default=None,
action=KvAction,
help='background region')
argparser.add_argument(
'--elow',
default=4.,
action=KvAction,
help='low energy bound')
argparser.add_argument(
'--ehigh',
default=12.,
action=KvAction,
help='high energy bound')
argparser.add_argument(
'--binsize',
default=1.,
action=KvAction,
help='time series bin size')
argparser.add_argument(
'--epoch0',
default=None,
action=KvAction,
help='zero epoch for light curves')
argparser.add_argument(
'--orient',
default=None,
action=KvAction,
help='orientation file')
argparser.add_argument(
'--usergti',
type=str,
default=None,
help='user gti file',
action=KvAction)
argparser.add_argument(
'--make-individual',
default=False,
action='store_true',
help='make individual products for each telescope')
args = argparser.parse_args()
self.stem = args.stem
self.srcdir = Dir(args.srcdir)
self.dstdir = Dir(args.dstdir)
self.tmpdir = Dir(args.tmpdir) / str(self.pid)
self.tmpdst = self.tmpdir / 'dst'
self.prepare = args.prepare_skytile
self.fitwcs = args.fit_to_wcs
self.wcs = args.wcs
self.srcreg = args.srcreg
self.bkgreg = args.bkgreg
self.srcra = args.srcra
self.srcdec = args.srcdec
self.elow = float(args.elow )
self.ehigh = float(args.ehigh )
self.binsize = float(args.binsize)
self.epoch0 = None
if args.epoch0 is not None:
self.epoch0 = float(args.epoch0)
self.orifile = args.orient
self.usergti = args.usergti
self.make_individual = args.make_individual
def print(self):
print('Input parameters:')
print(' stem = {}'.format(self.stem ))
print(' srcdir = {}'.format(self.srcdir.get() ))
print(' dstdir = {}'.format(self.dstdir.get() ))
print(' tmpdir = {}'.format(self.tmpdir.get() ))
print(' prepare = {}'.format(self.prepare ))
print(' fitwcs = {}'.format(self.fitwcs ))
print(' wcs = {}'.format(self.wcs ))
print(' srcra = {}'.format(self.srcra ))
print(' srcdec = {}'.format(self.srcdec ))
print(' srcreg = {}'.format(self.srcreg ))
print(' bkgreg = {}'.format(self.bkgreg ))
print(' elow = {}'.format(self.elow ))
print(' ehigh = {}'.format(self.ehigh ))
print(' binsize = {}'.format(self.binsize ))
print(' epoch0 = {}'.format(self.epoch0 ))
print(' orient = {}'.format(self.orifile ))
print(' usergti = {}'.format(self.usergti ))
print(' individual = {}'.format(self.make_individual))
print()
def check(self):
if not self.srcdir.exists():
raise TaskError('<srcdir> not found')
self.dstdir.check()
self.tmpdir.check(clean=True)
self.tmpdst.check(clean=True)
class Products(Task):
def __init__(self):
super().__init__()
self._name = 'artproducts v.{}'.format(__version__)
print('*** {}'.format(self._name))
self._args = Args()
def initialize(self):
self._args.read()
self._args.print()
self._args.check()
self._imgdir = self._args.dstdir / 'img'
self._phadir = self._args.dstdir / 'pha'
self._lcdir = self._args.dstdir / 'lc'
def finalize(self):
self._args.tmpdir.delete()
def runmain(self):
if self._args.prepare:
self._evtdir = self._args.dstdir / 'events'
self._prepare()
else:
if self._args.orifile is None:
raise TaskError('<orient> file not set')
if not os.path.exists(self._args.orifile):
raise TaskError('<orient> file not found')
if self._args.srcra is None or self._args.srcdec is None:
raise TaskError('<srcra>/<srcdec> not set')
self._evtdir = self._args.srcdir
self._products()
print('*** artproducts v.{} finished'.format(__version__))
def _prepare(self):
self._evtdir.check(clean=True)
self._imgdir.check(clean=True)
self._phadir.check(clean=True)
self._lcdir.check(clean=True)
evtfiles = []
print('Preprocessing event files:')
for telname in Teldef.telnames:
evtfile = ''.join((self._args.stem, telname, '_cl.fits'))
evtdst = self._evtdir.pathfor(evtfile)
evtlist = EventsWcs(self._args.srcdir.pathfor(evtfile), self._args.wcs)
evtlist.prepare(self._args.fitwcs, self._args.usergti)
evtlist.write(evtdst, creator=self._name)
print('written: {}'.format(evtdst))
evtfiles.append(evtdst)
print()
print('Make image:')
imgfile = ''.join((self._args.stem, '_all.img'))
imgdst = self._imgdir.pathfor(imgfile)
image = Image(evtfiles, self._args.wcs)
image.make(self._args.elow, self._args.ehigh)
image.write(imgdst)
def _products(self):
cdbinfo = CaldbInfo(self._args.caldbdir).get()
print('CALDB info: version={}, path={}'.format(cdbinfo['CALDBVER'], cdbinfo['CALDBDIR']))
self._orient \
= Orientation(
CaldbOrientation(self._args.caldbdir, 'GYRO'))
oridata = self._orient.read(self._args.orifile)
print()
print('Reading events...')
events = EventsPack(self._evtdir).read(self._args.stem, self._args.epoch0)
if events.empty:
raise TaskError('event files empty')
print(' GTI = [{}, {}]'.format(events.tstart, events.tstop))
print(' epoch0 = {}'.format(events.epoch0))
print(' total events = {}'.format(events.count))
print()
print('Reading regions...')
srcreg = Region(self._args.srcreg)
bkgreg = Region(self._args.bkgreg)
srcarea = events.call(Calculate(srcreg.area, ['projection']))
bkgarea = events.call(Calculate(bkgreg.area, ['projection']))
areascale = 1.
if bkgarea['all'] > 0:
backscale = srcarea['all'] / bkgarea['all']
else:
backscale = srcarea['all']
print(' source area = {}'.format(srcarea['all']))
print(' background area = {}'.format(bkgarea['all']))
print(' area scale = {}'.format(areascale))
print(' back scale = {}'.format(backscale))
print()
print('Making lightcurves...')
self._make_lcurves(events, self._args.binsize, srcreg, bkgreg)
print()
print('Making energy spectra...')
self._make_enspectra(events, oridata, srcreg, bkgreg, areascale, backscale)
def _make_lcurves(self, events, binsize, srcreg, bkgreg):
dstdir = self._args.dstdir / 'lc'
dstdir.check()
print('background:')
lcurvebkgfile = ''.join((self._args.stem, '{}', '_bkg.lc'))
lcurvebkgfile = dstdir.pathfor(lcurvebkgfile)
lcurvebkg = Lightcurve(events, binsize, make_individual=self._args.make_individual)
lcurvebkg.create(bkgreg, self._args.elow, self._args.ehigh)
meta = lcurvebkg.meta['all']
total = meta['events']
selected = meta['selected']
per = np.round((selected / total) * 100, 2)
print(' filter = {}'.format(meta['filter']))
print(' selected events = {} of {} ({}%)'.format(selected, total, per))
writtenbkg = lcurvebkg.write(lcurvebkgfile)
print('written:')
for item in writtenbkg['files']:
print(' {}'.format(item))
print('source:')
lcurvesrcfile = ''.join((self._args.stem, '{}', '_src.lc'))
lcurvesrcfile = dstdir.pathfor(lcurvesrcfile)
lcurvesrc = Lightcurve(events, binsize, make_individual=self._args.make_individual)
lcurvesrc.create(srcreg, self._args.elow, self._args.ehigh)
meta = lcurvesrc.meta['all']
total = meta['events']
selected = meta['selected']
per = np.round((selected / total) * 100, 2)
print(' filter = {}'.format(meta['filter']))
print(' selected events = {} of {} ({}%)'.format(selected, total, per))
writtensrc = lcurvesrc.write(lcurvesrcfile, writtenbkg)
print('written:')
for item in writtensrc['files']:
print(' {}'.format(item))
def _make_responses(self, events, oridata, srcreg):
dstdir = self._args.dstdir / 'pha'
dstdir.check()
dstarffile = ''.join((self._args.stem, '{}', '.arf'))
dstarffile = dstdir.pathfor(dstarffile)
arf = Arf(self._args.caldbdir, self._orient, make_individual=self._args.make_individual)
arf.create(events, oridata, srcreg, self._args.srcra, self._args.srcdec)
print('Warning: average ARF responce will *NOT* be generated (use external procedure)')
arfwritten = arf.write(dstarffile)
print('written arf responces:')
if not arfwritten['files']:
print('<none>')
else:
for arfitem in arfwritten['files']:
print(arfitem)
dstdir = self._args.dstdir / 'pha'
dstdir.check()
dstrmffile = ''.join((self._args.stem, '{}', '.rmf'))
dstrmffile = dstdir.pathfor(dstrmffile)
rmf = Rmf(self._args.caldbdir, make_individual=self._args.make_individual)
rmf.create()
print('Warning: average RMF responce will *NOT* be generated (use external procedure)')
rmfwritten = rmf.write(dstrmffile)
print('written rmf responces:')
if not rmfwritten['files']:
print('<none>')
else:
for rmfitem in rmfwritten['files']:
print(rmfitem)
return arfwritten, rmfwritten
def _make_enspectra(self, events, oridata, srcreg, bkgreg, areascale, backscale):
# responces
arffiles, rmffiles = self._make_responses(events, oridata, srcreg)
responces = {
'arf': arffiles,
'rmf': rmffiles
}
# background spectra
enspecbkgfile = ''.join((self._args.stem, '{}', '_bkg.pha'))
enspecbkgfile = (self._args.dstdir / 'pha').pathfor(enspecbkgfile)
spectrumbkg = EnergySpectrum(events, make_individual=self._args.make_individual)
spectrumbkg.create(bkgreg)
writtenbkg = spectrumbkg.write(enspecbkgfile)
print('written background pha:')
if not writtenbkg['files']:
print('<none>')
else:
for phaitem in writtenbkg['files']:
print(phaitem)
# source spectra
enspecsrcfile = ''.join((self._args.stem, '{}', '_src.pha'))
enspecsrcfile = (self._args.dstdir / 'pha').pathfor(enspecsrcfile)
spectrumsrc = EnergySpectrum(events, make_individual=self._args.make_individual)
spectrumsrc.create(srcreg)
writtensrc = spectrumsrc.write(enspecsrcfile, responces, writtenbkg, backscale)
print('written source pha:')
if not writtensrc['files']:
print('<none>')
else:
for phaitem in writtensrc['files']:
print(phaitem)
def main():
try:
Products().run()
except RegionsError as e:
print('error: regions:', e)
except TaskError as e:
print('error: task execution:', e)
except Exception as e:
print('error: unknown:', e)
# Products().run()
if __name__ == '__main__':
pass

View File

@@ -0,0 +1,69 @@
#!/usr/bin/env python3
# -*- coding: utf8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2021-2026 IKI RAS (http://iki.rssi.ru/)
# Space Research Institute of the Russian Academy of Science
#
# This file is part of the M. Pavlinsky SRG/ART-XC software project.
# -----------------------------------------------------------------------------
from regions import Regions
from regions import CompoundSkyRegion
from regions.shapes.circle import CircleSkyRegion
from regions.shapes.annulus import CircleAnnulusSkyRegion
from astropy.coordinates import SkyCoord
from arttools.errors import RegionsError
class Region(object):
def __init__(self, regfile):
self._regfile = regfile
# -- read regions
self._reglist = Regions.read(self._regfile)
self._regunit = self._reglist[0].center.info.unit
# --- check
if not self._reglist:
raise RegionsError('no regions found')
for region in self._reglist:
region_allowed = \
isinstance(region, (CircleSkyRegion, CircleAnnulusSkyRegion))
if not region_allowed:
raise RegionsError('unsupported region: {}'.format(type(region)))
# -- make union region
self._region = self._reglist[0]
for region in self._reglist[1:]:
self._region = self._region.union(region)
def area(self, projection):
area = 0
for region in self._reglist:
area += region.to_pixel(projection).area
return area
def mask(self, projection, ra, dec):
radec = SkyCoord(ra, dec, unit=self._regunit, frame='fk5')
return self._region.contains(radec, projection)
# masks = []
# for r in self._regions:
# pass
# rmask = r.contains(radec, projection)
# masks.append(rmask)
# return np.any(masks, axis=0)
if __name__ == '__main__':
pass

View File

@@ -0,0 +1,316 @@
#!/usr/bin/env python3
# -*- coding: utf8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2026 IKI RAS (http://iki.rssi.ru/)
# Space Research Institute of the Russian Academy of Science
#
# This file is part of the M. Pavlinsky SRG/ART-XC software project.
# -----------------------------------------------------------------------------
import os
import shutil
import numpy as np
from astropy.io import fits
from arttools.fstools import Dir
from arttools.telescope import Teldef
from artpipeline.caldb import Caldb
from artpipeline.caldb import CaldbTelescope
class Arf(object):
EBINS = np.array([4 , 5 , 6, 8, 10, 14, 18 , 25])
EECF = np.array([4.5, 5.5, 7, 9, 12, 16, 21.5 ])
PCORR = {
'T1': [-5.53812028e-06, 5.44531758e-04, -1.85116554e-02, 2.49542143e-01, -4.78704144e-02],
'T2': [-9.23161584e-06, 8.50470939e-04, -2.71210316e-02, 3.46807627e-01, -5.10534726e-01],
'T3': [-6.39175857e-06, 6.37572616e-04, -2.21113466e-02, 3.02875927e-01, -3.00549874e-01],
'T4': [-3.79975811e-06, 3.78351940e-04, -1.38558825e-02, 2.11893093e-01, 5.14311717e-02],
'T5': [-8.51754314e-06, 8.25307741e-04, -2.78721554e-02, 3.75886273e-01, -6.00306320e-01],
'T6': [-5.44416755e-06, 5.33468306e-04, -1.81215482e-02, 2.45438278e-01, -1.13582748e-01],
'T7': [-7.31913644e-06, 6.80886499e-04, -2.23323962e-02, 2.98891459e-01, -3.07578086e-01]
}
ECF_R3 = {
'T1': [0.8405893291928408, 0.8578133072455523, 0.8518679621578527, 0.8435798995261736, 0.8289480999358520, 0.7258092846885180, 0.6826893958459460],
'T2': [0.7616948204797905, 0.7551934086121032, 0.7265351319204454, 0.7101925259375811, 0.6654367283015024, 0.5911152050970271, 0.6112236775368495],
'T3': [0.8291727938992725, 0.8257931679150124, 0.8188083629048590, 0.8057277406342708, 0.7642634099905733, 0.6574702032925820, 0.6553824446937162],
'T4': [0.8587518633143223, 0.8629811713390994, 0.8603473794365152, 0.8509428156196533, 0.8255339443817531, 0.7794067227678055, 0.7313918593446697],
'T5': [0.8710124678275055, 0.8724032722913310, 0.8690978752281667, 0.8561072162925373, 0.8341782730242377, 0.7799308511134618, 0.6982241619541473],
'T6': [0.8621286963648814, 0.8566493477571664, 0.8479740523247319, 0.8344740169388141, 0.8119985053434914, 0.7497446386004684, 0.6533235819021327],
'T7': [0.8361233953579142, 0.8314687457264499, 0.8167286887863530, 0.8073078896029924, 0.7795029443860121, 0.6759143436928101, 0.6415284037194291]
}
ECF_R5 ={
'T1': [0.8674891551004330, 0.8822134383381507, 0.8788946213372089, 0.8744636104265867, 0.8597222695685300, 0.7557929438315966, 0.7189331472867164],
'T2': [0.7936437182568867, 0.7891579186682454, 0.7674259998254548, 0.7589714843682488, 0.7158317895510184, 0.6411584008631043, 0.6751573402559861],
'T3': [0.8582713402059252, 0.8542245360657326, 0.8490572253823582, 0.8385523961057851, 0.7992921226502082, 0.7082381667579745, 0.7148249035011437],
'T4': [0.8858430082172500, 0.8867274026307359, 0.8870164916860136, 0.8793272870044280, 0.8561885232880430, 0.8040219587253561, 0.7604390487299408],
'T5': [0.8930065538441408, 0.8950301833952440, 0.8942787134906093, 0.8854633646996730, 0.8671659219235113, 0.8201161803780476, 0.7338377592101047],
'T6': [0.8876123076470380, 0.8835789100530228, 0.8807700953185790, 0.8736529044505901, 0.8565833527159418, 0.7983704564521831, 0.7009986224989362],
'T7': [0.8694210354890501, 0.8620577130319129, 0.8498935555726601, 0.8429003979100774, 0.8204142620537092, 0.7239238390896476, 0.6748621998062087]
}
def __init__(self, caldbpath, orient, make_individual=False):
self._caldbdir = Dir(caldbpath)
self._orient = orient
self._make_individual = make_individual
srcvigndir = self._caldbdir / 'ext' / 'vign'
self._srcvignfile = srcvigndir.pathfor('art-xc_vignea_full_q200.fits')
srcarfdir = self._caldbdir / 'ext' / 'arf'
self._srcarffile = srcarfdir.pathfor('artxc_arf_20191030_v004.arf')
self._arcmin2deg = 1. / 60.
self._pix2arcmin = 45. / 60.
self._respdata = {}
self._respmeta = {}
def _read_src_vign(self, vignfile):
hdul = fits.open(vignfile)
mapshape = hdul['All events' ].data[0][1].shape
vignmaps = hdul['All events' ].data
vignxoff = hdul['Offset angles'].data['X']
vignyoff = hdul['Offset angles'].data['Y']
return {
'vignmaps': vignmaps,
'mapshape': mapshape,
'xoff' : vignxoff,
'yoff' : vignyoff
}
def _read_src_arf(self, arffile):
hdul = fits.open(arffile)
elow = hdul[1].data['ENERG_LO']
ehigh = hdul[1].data['ENERG_HI']
resp = hdul[1].data['SPECRESP']
return {
'elow' : elow ,
'ehigh' : ehigh,
'specresp': resp
}
def _calc_detxy(self, telname, evtdata, oridata, ra, dec):
evttimes = np.array(evtdata.events.column('TIME'), dtype=np.double)
evtra = np.full(shape=evttimes.shape, fill_value=ra , dtype=np.double)
evtdec = np.full(shape=evttimes.shape, fill_value=dec, dtype=np.double)
detx, dety \
= self._orient.calc_rawxy(
oridata['time'],
oridata['quat'],
evttimes ,
evtra ,
evtdec ,
telname
)
return np.median(detx).astype(int), np.median(dety).astype(int)
def _calc_offset(self, telname, x, y):
telcaldb = \
CaldbTelescope(
Caldb(self._caldbdir.get(), telname))
oapixdist = np.hypot(
float(x) - telcaldb.oaxis()['X'],
float(y) - telcaldb.oaxis()['Y'])
oaoffset = oapixdist * self._pix2arcmin
return oaoffset
def _make_masks(self, mapshape, xoff, yoff, offset):
tolerance = 0.1
xoff = Teldef.F * np.tan(np.deg2rad(xoff * self._arcmin2deg))
yoff = Teldef.F * np.tan(np.deg2rad(yoff * self._arcmin2deg))
cxy = np.array(np.meshgrid(xoff, yoff)).T.reshape(-1, 2)
cx = cxy.T[0]
cy = cxy.T[1]
offsetmap = np.rad2deg(np.arctan(
np.hypot(cx, cy) / Teldef.F)) / self._arcmin2deg
cenmask = np.abs(offsetmap - 0. ) <= tolerance
offmask = np.abs(offsetmap - offset) <= tolerance
return cenmask.reshape(mapshape), offmask.reshape(mapshape)
def _make_vign_resp(self, vigndata, offset):
cenmask, offmask \
= self._make_masks(
vigndata['mapshape'],
vigndata['xoff' ],
vigndata['yoff' ],
offset )
energies = vigndata['vignmaps']['E' ]
effareas = vigndata['vignmaps']['EFFAREA']
scales = []
for ea in effareas:
scales.append(np.mean(ea[offmask]) / np.mean(ea[cenmask]))
return energies, scales
def _make_arf_resp(self, telname, arfdata, vignscale):
arfen = (arfdata['elow'] + arfdata['ehigh']) / 2.
arfresp = arfdata['specresp'] / 7.
vignen = vignscale[0]
vignresp = vignscale[1]
scaled = np.interp(arfen, vignen, vignresp)
newresp = arfresp * scaled
# --- polycorr ---
pcorr = self.PCORR[telname]
newresp = newresp * np.polyval(pcorr, arfen)
# ---
# --- ecfcorr ---
ecf = np.array(self.ECF_R3[telname])
newresp = newresp * np.interp(arfen, self.EECF, ecf)
# ---
return newresp
def create(self, events, oridata, srcreg, srcra, srcdec):
srcarfdata = self._read_src_arf(self._srcarffile)
srcvigndata = self._read_src_vign(self._srcvignfile)
for telname, evtdata in events:
events = evtdata['events' ]
projection = evtdata['projection']
detx, dety = self._calc_detxy(
telname,
events ,
oridata,
srcra ,
srcdec
)
offset = self._calc_offset(telname, detx, dety)
es, vs = self._make_vign_resp(srcvigndata, offset)
arfresp = self._make_arf_resp(telname, srcarfdata, (es, vs))
self._respdata[telname] = arfresp
self._respmeta[telname] = {
'ra' : srcra ,
'dec' : srcdec,
'detx' : detx ,
'dety' : dety ,
'offset': offset
}
logstr = '{}: ra = {}, dec = {}, x = {}, y = {}, offset = {}'.format(
telname, srcra, srcdec, detx, dety, offset
)
print(logstr)
def _write_resp(self, telname, srcfile, dstfile):
if not os.path.exists(srcfile):
return False
if telname not in self._respdata:
return False
if os.path.exists(dstfile):
os.remove(dstfile)
shutil.copy(srcfile, dstfile)
with fits.open(dstfile, mode='update') as hdul:
hdul[1].data["SPECRESP"] = self._respdata[telname]
hdul.flush()
return True
def write(self, dstarf):
written = {
'files': []
}
if not self._make_individual:
return written
for telname in Teldef.telnames:
if telname not in self._respdata:
continue
dstarfname = dstarf.format(telname)
if not self._write_resp(telname, self._srcarffile, dstarfname):
continue
written[telname] = os.path.basename(dstarfname)
written['files'].append(dstarfname)
return written
class Rmf(object):
def __init__(self, caldbpath, make_individual=False):
self._caldbdir = Dir(caldbpath)
self._make_individual = make_individual
srcrmfdir = self._caldbdir / 'ext' / 'rmf'
self._srcrmffile = srcrmfdir.pathfor('artxc_rmf_{}_20191030_v003.rmf')
def create(self):
pass
def _write_resp(self, telname, srcfile, dstfile):
if not os.path.exists(srcfile):
return False
if os.path.exists(dstfile):
os.remove(dstfile)
shutil.copy(srcfile, dstfile)
return True
def write(self, dstrmf):
written = {
'files': []
}
if not self._make_individual:
return written
for telname in Teldef.telnames:
srcrmfname = self._srcrmffile.format(telname)
dstrmfname = dstrmf.format(telname)
if not self._write_resp(telname, srcrmfname, dstrmfname):
continue
written[telname] = os.path.basename(dstrmfname)
written['files'].append(dstrmfname)
return written
if __name__ == '__main__':
pass

View File

@@ -0,0 +1 @@
__version__ = '2.0.0'