diff --git a/artpipeline/.DS_Store b/artpipeline/.DS_Store new file mode 100644 index 0000000..9a874b5 Binary files /dev/null and b/artpipeline/.DS_Store differ diff --git a/artpipeline/pyproject.toml b/artpipeline/pyproject.toml new file mode 100644 index 0000000..992522d --- /dev/null +++ b/artpipeline/pyproject.toml @@ -0,0 +1,6 @@ +[build-system] +requires = [ + "setuptools>=42", +] +build-backend = "setuptools.build_meta" + diff --git a/artpipeline/setup.py b/artpipeline/setup.py new file mode 100644 index 0000000..b0bed7d --- /dev/null +++ b/artpipeline/setup.py @@ -0,0 +1,17 @@ +from setuptools import setup + +exec(open('src/artpipeline/version.py').read()) + +setup( + name="artpipeline", + version=__version__, + author="M.Pavlinsky SRG/ART-XC software team", + description="data reduction pipeline", + package_dir={"": "src"}, + packages=["artpipeline"], + entry_points={ + "console_scripts": [ + "artpipeline = artpipeline.pipeline:main", + ] + }, +) diff --git a/artpipeline/src/.DS_Store b/artpipeline/src/.DS_Store new file mode 100644 index 0000000..6975267 Binary files /dev/null and b/artpipeline/src/.DS_Store differ diff --git a/artpipeline/src/artpipeline/caldb.py b/artpipeline/src/artpipeline/caldb.py new file mode 100644 index 0000000..6947016 --- /dev/null +++ b/artpipeline/src/artpipeline/caldb.py @@ -0,0 +1,325 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2020-2025 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 +from pathlib import Path +from datetime import datetime as dt + +import numpy as np +from scipy.interpolate import interp1d + +import polars as pl + +from astropy.io import fits +# from astropy.table import Table + +from arttools.datatable import DataTable, FitsAdapter + + +# if 'ARTCALDB' in os.environ: +# ARTCALDBPATH = os.getenv('ARTCALDB') +# else: +# print('error: ARTCALDB environment variable not found') +# sys.exit(-1) + + +class CaldbInfo(object): + + def __init__(self, caldbdir): + self._caldbpath = caldbdir + self._indexfile = 'caldb.indx' + self._caldbindex = os.path.join(self._caldbpath, self._indexfile) + + with fits.open(self._caldbindex) as caldbfile: + self._version = caldbfile['CIF'].header['ARTCALDB'] + # self._index = Table(caldbfile['CIF'].data).to_pandas() + + self._index = DataTable().adapter(FitsAdapter(self._caldbindex, 'CIF')) + self._index.read() + + def get(self): + return { + 'CALDBDIR': self._caldbpath, + 'CALDBVER': self._version + } + + +class Caldb(object): + + def __init__(self, caldbdir, instname, silent=True): + + self._instname = instname + self._silent = silent + self._version = None + self._index = None + + # --- read version, index + self._caldbpath = caldbdir + self._indexfile = 'caldb.indx' + self._caldbindex = os.path.join(self._caldbpath, self._indexfile) + + with fits.open(self._caldbindex) as caldbfile: + self._version = caldbfile['CIF'].header['ARTCALDB'] + # self._index = Table(caldbfile['CIF'].data).to_pandas() + + self._index = DataTable().adapter(FitsAdapter(self._caldbindex, 'CIF')) + self._index.read() + + # print('CALDB: version {} [{}]'.format(self._version, self._caldbfilespath)) + + @property + def silent(self): + return self._silent + + def cif(self, instname, calcname): + iname = instname; cname = calcname + + if not self._silent: + print(self._silent) + print('CALDB: select: {}'.format( + 'INSTRUME==\'{}\' and CAL_CNAM==\'{}\''.format(iname, cname))) + + # iname = iname.ljust(10, ' ') + # cname = cname.ljust(20, ' ') + + # return self._index.query( + # 'INSTRUME==\'{}\' and CAL_CNAM==\'{}\''.format(iname, cname)) + + iname = iname.ljust(10, ' ') + cname = cname.ljust(20, ' ') + + return self._index.frame.filter( + pl.col('INSTRUME') == iname, + pl.col('CAL_CNAM') == cname + ) + + def file(self, calcname): + return self.instfile(self._instname, calcname) + + def instfile(self, instname, calcname): + selected = self.cif(instname, calcname) + for cdbpath, cdbfile in zip(selected['CAL_DIR'], selected['CAL_FILE']): + cdbpath = cdbpath.strip() + cdbfile = cdbfile.strip() + + if not self._silent: + print('CALDB: product: {}'.format(cdbfile)) + + p = Path(cdbpath) + p = Path(*p.parts[2:]) # chop off data/artxc_caldb + + cdbpath = p + + # return first record only, @fixme + return os.path.join(cdbpath, cdbfile) + + return None + + def path(self, folder='bcf'): + return os.path.join(self._caldbpath, folder) + + def instname(self): + return self._instname + + def version(self): + return self._version + + +class CaldbTelescope(object): + + def __init__(self, caldb): + self._caldb = caldb + self._silent = self._caldb.silent + self._oaxis_file \ + = os.path.join(self._caldb.path(), self._caldb.file('OPTAXIS')) + + with fits.open(self._oaxis_file) as hdul: + self._oaxis = { + 'X': hdul[1].data['X'][0], + 'Y': hdul[1].data['Y'][0] + } + + if not self._silent: + print('CALDB: selected OPTAXIS: X={} Y={}'.format(self._oaxis['X'], self._oaxis['Y'])) + + def oaxis(self): + return self._oaxis + + def info(self): + meta = { + 'OPTAXIS': self._oaxis + } + + +# class CaldbEnergy(object): + +# def __init__(self, instname): +# self._caldb = Caldb(instname) + +# self._thresh_file \ +# = self._caldb.file('THRESHOL') +# self._escale_file \ +# = self._caldb.file('ESCALE' ) + +# #todo: check files is not None + +# self._thresh = fits.open(os.path.join(self._caldb.path(), self._thresh_file)) +# self._escale = fits.open(os.path.join(self._caldb.path(), self._escale_file)) + +# def bot(self): +# return self._thresh['BOT'].data + +# def top(self): +# return self._thresh['TOP'].data + +# def scale(self): +# # escale_data = self._escale['ESCALE'].data +# escale_data = self._escale['E_COEFF'].data +# obt_data = escale_data['OBT'] +# c0_data = escale_data['C0' ] +# c1_data = escale_data['C1' ] + +# scale_fn = interp1d(obt_data, c0_data, bounds_error=False, fill_value=tuple(c0_data[[0, -1]])) +# shift_fn = interp1d(obt_data, c1_data, bounds_error=False, fill_value=tuple(c1_data[[0, -1]])) + +# # print last time (obt_data) +# # save it to header + +# return { +# 'C0': scale_fn, +# 'C1': shift_fn} + +# def version(self): +# return self._caldb.version() + +# def info(self): +# meta = { +# 'THRESHOL': self._thresh_file, +# 'ESCALE' : self._escale_file, +# 'CALDBVER': self.version() +# } +# return meta + +class CaldbEnergy(object): + + def __init__(self, caldbdir, instname): + self._caldb = Caldb(caldbdir, instname) + + self._chen_file \ + = self._caldb.file('CHENERGY') + + self._chen = fits.open(os.path.join(self._caldb.path(), self._chen_file)) + + def bot(self): + return self._chen['BOT'].data + + def top(self): + return self._chen['TOP'].data + + def version(self): + return self._caldb.version() + + def info(self): + meta = { + 'CHENERGY': self._chen_file, + 'CALDBVER': self.version() + } + return meta + + +class CaldbMask(object): + + def __init__(self, caldbdir, instname): + self._caldb = Caldb(caldbdir, instname) + + self._dmask_file \ + = os.path.join(self._caldb.path(), self._caldb.file('DETMASK')) + self._dmask = np.copy(fits.getdata(self._dmask_file, 1)).astype(bool).T + + # self._dmask = fits.open(os.path.join(self._caldb.path(), self._dmask_file)) + + def get(self): + return self._dmask + + def version(self): + return self._caldb.version() + + def info(self): + meta = { + 'DETMASK': self._dmask_file, + 'CALDBVER': self.version() + } + return meta + + +class CaldbOrientation(object): + + _INST_QUAT = { + 'GYRO': [ 0. , 0. , 0. , 1. ], + 'BOKZ': [ 0. , -0.707106781186548, 0. , 0.707106781186548], + 'SED1': [ 0.183012701892219, 0.683012701892219, -0.183012701892219, 0.683012701892219], + 'SED2': [-0.183012701892219, 0.683012701892219, 0.183012701892219, 0.683012701892219]} + + def __init__(self, caldbdir, instname): + self._delay_file = None + self._bs_file = None + self._caldb = Caldb(caldbdir, instname) + self._delay = self.read_time_delay() + self._instquat = self.read_instquat () + self._boresight = self.read_boresight () + + def read_time_delay(self): + delay_file = self._caldb.file('{instname}LAG'.format(instname=self._caldb.instname())) + delay_hdul = fits.open(os.path.join(self._caldb.path(), delay_file)) + delay = delay_hdul[1].data['DELAY'][0] + + self._delay_file = delay_file + + # print('CALDB: select {} DELAY={}'.format(self._caldb.instname(), delay)) + + return delay + + def read_instquat(self): + return self._INST_QUAT[self._caldb.instname()] + + def read_boresight(self, instname=None): + if instname is None: + instname = self._caldb.instname() + + bs_file = self._caldb.instfile(instname, 'BORESIGH') + bs_hdul = fits.open(os.path.join(self._caldb.path(), bs_file)) + + self._bs_file = bs_file + + return bs_hdul[1].data[0] + + def timeshift(self): + return self._delay + + def instquat(self): + return self._instquat + + def boresight(self): + return self._boresight + + def version(self): + return self._caldb.version() + + def info(self): + meta = { + 'INSTRUME': self._caldb.instname(), + 'TIMELAG' : self._delay_file, + 'BORESIGH': self._bs_file + } + return meta + + +if __name__ == '__main__': + pass diff --git a/artpipeline/src/artpipeline/energy.old.py b/artpipeline/src/artpipeline/energy.old.py new file mode 100644 index 0000000..873e3b2 --- /dev/null +++ b/artpipeline/src/artpipeline/energy.old.py @@ -0,0 +1,152 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2020-2025 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the SRG/ART-XC software project. +# ----------------------------------------------------------------------------- + +import os +import sys +import shutil + +import numpy as np +from scipy.interpolate import interp1d + +from .time import Mission + + +class Energy(object): + + def __init__(self, caldb_energy, seed=None): + self._caldb_energy = caldb_energy + self._seed = seed + + self.GRADES = np.ones(256) * (-1) + self.GRADES[18] = 0 # single central events + self.GRADES[[50, 26, 22, 19, 54, 51, 30, 27]] = [1, 2, 3, 4, 5, 6, 7, 8] # double events + self.GRADES[[58, 62, 59, 23, 55, 31, 63]] = [9, 10, 11, 12, 13, 14, 15] # triple events + + def calc(self, evtdata, hkdata): + time = evtdata['time'] + tstart = Mission(time[0]).to_datetime() + escale = self._caldb_energy.scale() + ebot = self._caldb_energy.bot() + etop = self._caldb_energy.top() + + tempdata = hkdata['td1'] + temperature = interp1d( + hkdata['time'], tempdata, bounds_error=False, kind='linear', + fill_value=(tempdata[0], tempdata[-1]))(evtdata['time']) + + if self._seed is not None: + sseq = np.random.SeedSequence(self._seed) + else: + sseq = np.random.SeedSequence() + + rng = np.random.default_rng(sseq) + seed = sseq.entropy + + energybot, sigmabot, maskbot \ + = self._calc( + evtdata, + ('raw_x' , + 'pha_bot_sub1', + 'pha_bot' , + 'pha_bot_add1'), + temperature, ebot, rng) + energytop, sigmatop, masktop \ + = self._calc( + evtdata, + ('raw_y' , + 'pha_top_sub1', + 'pha_top' , + 'pha_top_add1'), + temperature, etop, rng) + + energy = np.zeros(evtdata.size, np.double) + mask = np.zeros((8, evtdata.size), bool) + + mask[2:5,:] = masktop + mask[5:8,:] = maskbot + + emask = np.any(mask[2:,:], axis=0) + + energybot = np.sum(energybot * maskbot , axis=0) + sigmabot2 = np.sum(np.power(sigmabot, 2.) * maskbot, axis=0) + + energytop = np.sum(energytop * masktop , axis=0) + sigmatop2 = np.sum(np.power(sigmatop, 2.) * masktop, axis=0) + + energy[emask] \ + = ((energybot * sigmatop2 + energytop * sigmabot2) / (sigmabot2 + sigmatop2))[emask] + energy[emask] = self._scale(evtdata['time'][emask], energy[emask], escale) + + # energy = self._scale(energy, escale) + pi = self._en2pi(energy, emask ) + grade = self._grade(mask ) + + return energy, pi, grade, temperature, energybot, energytop, seed + + def _calc(self, evtdata, coln, temperature, energycal, rng): + strip, mask \ + = self._eventindex(evtdata[coln[0]]) + pha = np.array([evtdata[coln[1]], evtdata[coln[2]], evtdata[coln[3]]]) + + e1 = self._pha2en(pha , strip, temperature, energycal) + e2 = self._pha2en(pha + 1., strip, temperature, energycal) + + # energy = np.random.uniform(e1, e2) + energy = rng.uniform(e1, e2) + + sigma = energycal['fwhm_0'][strip] + energycal['fwhm_1'][strip] * temperature + mask = np.logical_and(mask, energy > energycal['threshold'][strip]) + + return energy, sigma, mask + + def _pha2en(self, pha, strip, temperature, energycal): + ''' Energy = c00 + c01*T + (c10 + c11*T)*PHA + (c20 + c21*T)*PHA**2 ''' + pha2 = np.power(pha, 2.) + energy = \ + energycal["c0_0"][strip] + energycal["c0_1"][strip] * temperature + \ + (energycal["c1_0"][strip] + energycal["c1_1"][strip] * temperature) * pha + \ + (energycal["c2_0"][strip] + energycal["c2_1"][strip] * temperature) * pha2 + + return energy + + def _scale(self, times, energy, escale): + ''' Еnergy = C0 + C1*Еnergy + С2*Еnergy**2 ''' + # return escale['C0'] + escale['C1'] * energy + escale['C2'] * np.power(energy, 2.) + + C0 = escale['C0'](times) + C1 = escale['C1'](times) + + return energy * C0 + C1 + + + def _en2pi(self, energy, mask): + PHA2PI = 10 # channels are in units of 100 eV + + pi = np.int16(energy * PHA2PI) - 1 + pi = np.maximum(pi, 0) + pi[~mask] = 0 + + return pi + + def _eventindex(self, strip): + coord = np.empty((3, strip.size), int) #(strip, (3, 1)) + coord[0, :] = strip - 1 + coord[1, :] = strip + coord[2, :] = strip + 1 + mask = (coord >= 0) & (coord <= 47) + coord[np.logical_not(mask)] = -1 + return coord, mask + + def _grade(self, mask): + grademask = np.packbits(mask, axis=0)[0] + return self.GRADES[grademask].astype(int) + + +if __name__ == '__main__': + pass diff --git a/artpipeline/src/artpipeline/energy.py b/artpipeline/src/artpipeline/energy.py new file mode 100755 index 0000000..df55b08 --- /dev/null +++ b/artpipeline/src/artpipeline/energy.py @@ -0,0 +1,206 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2020-2025 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the SRG/ART-XC software project. +# ----------------------------------------------------------------------------- + +import os +import sys +import shutil + +import warnings + +import numpy as np +from scipy.interpolate import interp1d + +import astropy.units as u +from astropy.time.formats import erfa +from astropy.time import Time +from astropy.time import TimezoneInfo + +# from .time import Mission + +from arttools.telescope import Teldef + +np.seterr(all="ignore") +warnings.filterwarnings('ignore') + + +class Energy(object): + + def __init__(self, caldb, seed=None): + self._caldb = caldb + self._seed = seed + + self.MJDCDB = Time(58677.6464931, format='mjd') + self.MJDREF = Time(51543.8750000, format='mjd') #corresponds to date 01.01.2000 00:00:00 (UTC+3) + self.TZ_UTC = TimezoneInfo(utc_offset=0*u.hour) #UTC time zone + self.TZ_MSK = TimezoneInfo(utc_offset=3*u.hour) #UTC+3 (Moscow) time zone + + self.GRADES = np.ones(256) * (-1) + self.GRADES[18] = 0 # single central events + self.GRADES[[50, 26, 22, 19, 54, 51, 30, 27]] = [1, 2, 3, 4, 5, 6, 7, 8] # double events + self.GRADES[[58, 62, 59, 23, 55, 31, 63]] = [9, 10, 11, 12, 13, 14, 15] # triple events + + def _calc_mjdtime(self, missiontime): + return self.MJDREF + Time(missiontime / erfa.DAYSEC, format='mjd').mjd + + def _calc_dt(self, missiontime): + mjdt = self._calc_mjdtime(missiontime) + return (mjdt - self.MJDCDB).jd + + def _eventindex(self, strip): + # coord = np.empty((3, strip.size), int) #(strip, (3, 1)) + coord = np.empty((3, strip.count()), int) #(strip, (3, 1)) + coord[0, :] = strip - 1 + coord[1, :] = strip + coord[2, :] = strip + 1 + mask = (coord >= 0) & (coord <= 47) + coord[np.logical_not(mask)] = -1 + return coord, mask + + def _calc(self, evtdata, coln, dt, dt2, random, energycal): + data_size = dt.size + strip, mask_strip \ + = self._eventindex(evtdata.column(coln[0])) + pha0 = np.array([ + evtdata.column(coln[1]), + evtdata.column(coln[2]), + evtdata.column(coln[3]) + ], dtype=np.uint32) + + C0 = energycal['C0_0'][strip] + energycal['C0_1'][strip] * dt + energycal['C0_2'][strip] * dt2 + C1 = energycal['C1_0'][strip] + energycal['C1_1'][strip] * dt + energycal['C1_2'][strip] * dt2 + C2 = energycal['C2_0'][strip] + energycal['C2_1'][strip] * dt + energycal['C2_2'][strip] * dt2 + C3 = energycal['C3_0'][strip] + energycal['C3_1'][strip] * dt + energycal['C3_2'][strip] * dt2 + C4 = energycal['C4_0'][strip] + energycal['C4_1'][strip] * dt + energycal['C4_2'][strip] * dt2 + + if random is not None: + pha_rand = random.uniform(0, 1, pha0.shape) + else: + pha_rand = np.zeros(pha0.shape) + + pha = pha0 + pha_rand + + energy3 = C0 + C1 * pha + C2 * np.power(pha, 2) + + kofA = [0.9, 1.15] + dayA = [223., 1480.] + + C1_thr = (kofA[0] - kofA[1]) / (dayA[0] - dayA[1]) + C0_thr = ((kofA[0] + kofA[1]) - (dayA[0] + dayA[1]) * C1_thr) / 2. + + kof = C0_thr + C1_thr * dt + + mask_threshold = energy3 > energycal['THRESHOLD'][strip] * kof + mask_threshold[1] = np.full_like(mask_threshold[1], True) + + mask3 = np.logical_and(mask_strip, mask_threshold) + + energy3[~mask3] = 0 + energy = np.sum(energy3, axis=0) + + # energy = np.sum(energy3 * mask3, axis=0) + + fwhm = (C3[1] + C4[1] * np.sqrt(energy)) + + return energy, fwhm, mask3 + + def _grade(self, mask): + grademask = np.packbits(mask, axis=0)[0] + + return self.GRADES[grademask].astype(int) + + def calc(self, teln, evtdata, random=None): + dt = self._calc_dt(evtdata.column('TIME')) + dt2 = np.power(dt, 2) + + energybot, fwhmbot, maskbot3 =\ + self._calc( + evtdata, + ('RAW_X' , + 'PHA_BOT_SUB1', + 'PHA_BOT' , + 'PHA_BOT_ADD1'), + dt, dt2, random, self._caldb.bot() + ) + fwhmbot_2 = np.power(fwhmbot, 2.) + + energytop, fwhmtop, masktop3 =\ + self._calc( + evtdata, + ('RAW_Y' , + 'PHA_TOP_SUB1', + 'PHA_TOP' , + 'PHA_TOP_ADD1'), + dt, dt2, random, self._caldb.top() + ) + fwhmtop_2 = np.power(fwhmtop, 2.) + + # -- calculate mask + mask = np.zeros((8, evtdata.size), bool) + mask[2:5,:] = masktop3 + mask[5:8,:] = maskbot3 + + # energy = np.zeros(dt.size, float) + energy =\ + ((energybot * fwhmtop_2 + energytop * fwhmbot_2) / (fwhmbot_2 + fwhmtop_2)) + + grades = self._grade(mask) + + return energy, energybot, energytop, grades + + +class Flag(object): + + def __init__(self, caldb): + self._caldb = caldb + + def _mask_events(self, evtdata, tstart, tstop): + times = evtdata.column('TIME') + # try: + # tstart = evtdata.tstart + # tstop = evtdata.tstop + # except KeyError: + # tstart = times[ 0] + # tstop = times[-1] + # print('WARNING: TSTART or TSTOP keywords not found') + # print('WARNING: will use event times TSTART = {}, TSTOP = {}'.format(tstart, tstop)) + + return (times > tstart) & (times < tstop) + + def _mask_shadow(self, evtdata): + evtrawx = evtdata.column('RAW_X') + evtrawy = evtdata.column('RAW_Y') + + shadowmask = self._caldb.get() + + return np.logical_not(shadowmask[evtrawx, evtrawy]) + + def _mask_edges(self, evtdata): + evtrawx = evtdata.column('RAW_X') + evtrawy = evtdata.column('RAW_Y') + + return np.any([ + evtrawx == 0, evtrawx == 47, + evtrawy == 0, evtrawy == 47], axis=0) + + def calc(self, evtdata, tstart, tstop): + flag = np.ones(evtdata.size, dtype=np.uint8) + + eventsmask = self._mask_events(evtdata, tstart, tstop) + shadowmask = self._mask_shadow(evtdata) + edgesmask = self._mask_edges (evtdata) + + flag[eventsmask] = 0 + flag[shadowmask] = 2 + flag[edgesmask ] = 3 + + return flag + + +if __name__ == '__main__': + pass diff --git a/artpipeline/src/artpipeline/events.py b/artpipeline/src/artpipeline/events.py new file mode 100755 index 0000000..1953d17 --- /dev/null +++ b/artpipeline/src/artpipeline/events.py @@ -0,0 +1,176 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2021-2025 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.gti import Gti +from arttools.datatable import DataTable, FitsAdapter + +from artpipeline.version import __version__ + + +class Events(object): + + def __init__(self, srcfile): + super().__init__() + + self._srcfile = srcfile + + self._telname = 0 + self._teln = 0 + + self._tstart = 0 + self._tstop = 0 + + self._events = DataTable() + self._hk = DataTable() + self._gti = Gti() + + @property + def size(self): + return self._events.size + + @property + def teln(self): return self._teln + + @property + def telname(self): return self._telname + + @property + def tstart(self): return self._tstart + + @property + def tstop(self): return self._tstop + + @property + def events(self): return self._events + + @property + def hk(self): return self._hk + + @property + def gti(self): return self._gti + + def _read_header(self): + with fits.open(self._srcfile) as hdul: + self._teln = hdul['events'].header['teln' ] + self._telname = hdul['events'].header['instrume'] + self._tstart = hdul['events'].header['tstart' ] + self._tstop = hdul['events'].header['tstop' ] + + def read(self): + self._read_header() + + self._events = DataTable().adapter(FitsAdapter(self._srcfile, 'events')).read() + self._hk = DataTable().adapter(FitsAdapter(self._srcfile, 'hk' )).read() + self._gti = Gti.from_hduname(self._srcfile, 'stdgti') + + return self + + def write_cl(self, dstfile, cdbinfo): + + srchdul = fits.open(self._srcfile) + + # --- primary HDU --- + hdu = fits.PrimaryHDU() + + # --- EVENTS HDU --- + a00 = np.array(self._events.column('TIME' )) + a01 = np.array(self._events.column('TIME_I' )) + a02 = np.array(self._events.column('TIME_F' )) + a03 = np.array(self._events.column('TIME_CORR' )) + a04 = np.array(self._events.column('ENERGY' )) + a05 = np.array(self._events.column('ENERGY_TOP' )) + a06 = np.array(self._events.column('ENERGY_BOT' )) + a07 = np.array(self._events.column('GRADE' )) + a08 = np.array(self._events.column('FLAG' )) + a09 = np.array(self._events.column('RA' )) + a10 = np.array(self._events.column('DEC' )) + a11 = np.array(self._events.column('RAW_X' )) + a12 = np.array(self._events.column('RAW_Y' )) + a13 = np.array(self._events.column('PHA_BOT' )) + a14 = np.array(self._events.column('PHA_BOT_ADD1')) + a15 = np.array(self._events.column('PHA_BOT_SUB1')) + a16 = np.array(self._events.column('PHA_TOP' )) + a17 = np.array(self._events.column('PHA_TOP_ADD1')) + a18 = np.array(self._events.column('PHA_TOP_SUB1')) + a19 = np.array(self._events.column('TRIGGER' )) + + c00 = fits.Column(name='TIME' , format='1D', unit='sec', array=a00 ) + c01 = fits.Column(name='TIME_I' , format='1J', unit='sec', array=a01, bzero=2147483648) + c02 = fits.Column(name='TIME_F' , format='1D', unit='sec', array=a02 ) + c03 = fits.Column(name='TIME_CORR' , format='1D', unit='sec', array=a03 ) + c04 = fits.Column(name='ENERGY' , format='1D', unit='keV', array=a04 ) + c05 = fits.Column(name='ENERGY_TOP' , format='1D', unit='keV', array=a05 ) + c06 = fits.Column(name='ENERGY_BOT' , format='1D', unit='keV', array=a06 ) + c07 = fits.Column(name='GRADE' , format='1I', unit='' , array=a07 ) + c08 = fits.Column(name='FLAG' , format='1I', unit='' , array=a08 ) + c09 = fits.Column(name='RA' , format='1D', unit='deg', array=a09 ) + c10 = fits.Column(name='DEC' , format='1D', unit='deg', array=a10 ) + c11 = fits.Column(name='RAW_X' , format='1I', unit='pix', array=a11 ) + c12 = fits.Column(name='RAW_Y' , format='1I', unit='pix', array=a12 ) + c13 = fits.Column(name='PHA_BOT' , format='1I', unit='' , array=a13 ) + c14 = fits.Column(name='PHA_BOT_ADD1', format='1I', unit='' , array=a14 ) + c15 = fits.Column(name='PHA_BOT_SUB1', format='1I', unit='' , array=a15 ) + c16 = fits.Column(name='PHA_TOP' , format='1I', unit='' , array=a16 ) + c17 = fits.Column(name='PHA_TOP_ADD1', format='1I', unit='' , array=a17 ) + c18 = fits.Column(name='PHA_TOP_SUB1', format='1I', unit='' , array=a18 ) + c19 = fits.Column(name='TRIGGER' , format='1I', unit='' , array=a19 ) + + cols = fits.ColDefs([ + c00, c01, c02, c03, c04, c05, c06, c07, c08, c09, + c10, c11, c12, c13, c14, c15, c16, c17, c18, c19 + ]) + + evthdu = fits.BinTableHDU.from_columns(cols) + evthdu.name = 'EVENTS' + + creator = 'artpipeline v.{}'.format(__version__) + + evthdu.header['creator' ] = (creator , 'program and version that created this file' ) + evthdu.header['origin' ] = ('IKI RAS' , 'institution that created this file' ) + evthdu.header['telescop'] = ('SRG/ART-XC', 'mission name' ) + + # copy header keys + COPYKEYS = [ + 'instrume', 'teln', 'detn', 'urdn', + 'mjdref', 'artday', 'tstart', 'tstop', + 'date', 'date-obs', 'time-obs', 'date-end', 'time-end', + 'maxticks', 'timedel', 'timeunit', 'timeref', 'tassign' + ] + + for key in COPYKEYS: + try: + evthdu.header[key] = srchdul['EVENTS'].header[key] + except KeyError: + pass + + + # add caldb key + for key in cdbinfo: + evthdu.header[key] = cdbinfo[key] + + # --- make file --- + # hdulist = fits.HDUList([hdu, evthdu, srchdul['']]) + hdulist = fits.HDUList([hdu, evthdu, srchdul['KVEA'], srchdul['HK'], srchdul['STDGTI']]) + + hdulist.writeto(dstfile, overwrite=True, checksum=True) + + + def setcol(self, name, arr): + self._events.setcol(name, arr) + + +if __name__ == '__main__': + pass diff --git a/artpipeline/src/artpipeline/orientation.py b/artpipeline/src/artpipeline/orientation.py new file mode 100644 index 0000000..fe5665a --- /dev/null +++ b/artpipeline/src/artpipeline/orientation.py @@ -0,0 +1,302 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2020-2025 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 numpy as np +from scipy.spatial.transform import Slerp +from scipy.spatial.transform import Rotation + +from astropy.io import fits +from astropy.nddata import bitmask + +from arttools.telescope import Teldef + + +class AttitudeStatus(object): + + INTERP = 1 + SCIREADY = 2 + PRECSTAB = 4 + SED1_ON = 8 + SED2_ON = 16 + BOKZ_ON = 32 + SED1_MEAS = 64 + SED2_MEAS = 128 + BOKZ_MEAS = 256 + + flags = np.array([ + INTERP , + SCIREADY , + PRECSTAB , + SED1_ON , + SED2_ON , + BOKZ_ON , + SED1_MEAS, + SED2_MEAS, + BOKZ_MEAS + ]) + + @classmethod + def flags_except(cls, flags): + indices = np.searchsorted(cls.flags, flags) + return np.delete(cls.flags, indices) + + @classmethod + def pack(cls, arr): + bits = [cls.flags[x.astype(bool)] for x in arr.T] + return [bitmask.interpret_bit_flags(x) for x in bits] + + @classmethod + def flagval(cls, val, flags): + return bitmask.bitfield_to_boolean_mask(val, ignore_flags=flags) + + @classmethod + def extract(cls, arr, flags): + return np.array([AttitudeStatus.flagval(x, flags) for x in arr]) + + @staticmethod + def tostr(val, fmt='{0:016b}'): + return fmt.format(val) + + +class Orientation(object): + + UNITY_QUAT = np.array([0., 0., 0., 1.]) + + VSC = { + 'X': np.array([1., 0., 0.]), + 'Y': np.array([0., 1., 0.]), + 'Z': np.array([0., 0., 1.])} + + OPTICAL_AXIS = VSC['X'] + NORTH = VSC['Z'] + + RELCORR_GTI = [ + [624390347, 624399808], + [624410643, 630954575]] + + def __init__(self, caldb): + self._caldb = caldb + + def calc_events_coords(self, oridata, evtdata, telname, random=None): + oritime = np.array(oridata['time']) + evttime = np.array(evtdata.column('TIME' )) + rawx0 = np.array(evtdata.column('RAW_X')) + rawy0 = np.array(evtdata.column('RAW_Y')) + + mask = np.logical_and(evttime > oritime[0], evttime < oritime[-1]) + + if random is not None: + rawx_rand = random.uniform(-0.5, 0.5, rawx0.shape) + rawy_rand = random.uniform(-0.5, 0.5, rawy0.shape) + else: + rawx_rand = 0 + rawy_rand = 0 + + rawx = rawx0 + rawx_rand + rawy = rawy0 + rawy_rand + + orient = Slerp(oridata['time'], oridata['quat']) + telquat = Rotation(self._caldb.read_boresight(telname)) + + evtquat = np.zeros_like(evttime) + ra = np.zeros_like(evttime) + dec = np.zeros_like(evttime) + + evtquat = orient(evttime[mask]) * telquat + + ra, dec = self.rawxy_to_radec(evtquat, rawx[mask], rawy[mask]) + + ra_deg = np.zeros_like(evttime) + dec_deg = np.zeros_like(evttime) + + ra_deg[mask] = np.rad2deg(ra) + dec_deg[mask] = np.rad2deg(dec) + + return ra_deg, dec_deg, mask + + def calc_attitudes(self, oridata, instname=None): + instquat = self.instquat(instname) + quatdata = oridata['quat_sc'] + if instname != None: + quatdata = oridata['quat'] + + orient = Slerp (oridata['time'], quatdata) + attquat = orient(oridata['time']) * Rotation(instquat) + + ra, dec, roll \ + = self.calc_radecroll(attquat) + + return np.rad2deg(ra), np.rad2deg(dec), np.rad2deg(roll) + + def calc_rawxy(self, oritime, oriquat, evttime, ra, dec, telname): + ra = np.deg2rad(ra ) + dec = np.deg2rad(dec) + + orient = Slerp(oritime, oriquat) + telquat = Rotation(self._caldb.read_boresight(telname)) + evtquat = orient(evttime) * telquat + + rawx, rawy \ + = self.radec_to_rawxy(evtquat, ra, dec) + + return rawx, rawy + + def calc_radec(self, oritime, oriquat, evttime, evtrawx, evtrawy, telname): + orient = Slerp(oritime, oriquat) + telquat = Rotation(self._caldb.read_boresight(telname)) + + evtquat = orient(evttime) * telquat + + ra, dec = self.rawxy2radec(evtquat, evtrawx, evtrawy) + + return np.rad2deg(ra), np.rad2deg(dec) + + def calc_radecroll(self, quat): + oaxis = self.vec_normalize(quat.apply(self.OPTICAL_AXIS)) + + ra, dec \ + = self.vec_to_pol(oaxis) + + YZ = self.vec_normalize(np.cross(oaxis, self.NORTH)) + roll = np.arctan2( + np.sum(YZ * quat.apply(self.VSC['Z']), axis=1), + np.sum(YZ * quat.apply(self.VSC['Y']), axis=1)) + + return ra, dec, roll + + def instquat(self, instname=None): + if instname == None: + return self.UNITY_QUAT + + if instname != 'GYRO': + return self._caldb.read_boresight(instname) + + return self.UNITY_QUAT + + def read(self, orifile): + times, quats, state = self._read_gyro(orifile) + + quat = Rotation(quats) * Rotation(self._caldb.instquat()) * Rotation(self._caldb.boresight()) + quat_sc = Rotation(quats) * Rotation(self._caldb.instquat()) + + return { + 'size' : times.size, + 'time' : times , + 'quat' : quat , + 'quat_sc': quat_sc , + 'state' : state } + + def _read_gyro(self, orifile): + ori_hdul = fits.open(orifile) + data_size = ori_hdul[1].data.size + times = ori_hdul[1].data['TIME' ] + q0 = ori_hdul[1].data['QORT_1'] + q1 = ori_hdul[1].data['QORT_2'] + q2 = ori_hdul[1].data['QORT_3'] + q3 = ori_hdul[1].data['QORT_0'] + + times -= self._caldb.timeshift() + quats = np.array([q0, q1, q2, q3]).T + + state = np.array([ + np.zeros(data_size, dtype=int), + ori_hdul[1].data['ST_SCIREADY' ], + ori_hdul[1].data['ST_PRECSTAB' ], + ori_hdul[1].data['ST_SED1_ON' ], + ori_hdul[1].data['ST_SED2_ON' ], + ori_hdul[1].data['ST_BOKZ_ON' ], + ori_hdul[1].data['ST_SED1_MEAS' ], + ori_hdul[1].data['ST_SED2_MEAS' ], + ori_hdul[1].data['ST_BOKZ_MEAS' ] + ]) + + return times, quats, state + + def vec_normalize(self, vec): + norm = np.sqrt(np.sum(vec**2, axis=-1))[..., np.newaxis] + + return vec / norm + + def rawxy_to_offset(self, rawx, rawy): + xoffset = (rawx - Teldef.centeroffset) * Teldef.d + yoffset = (rawy - Teldef.centeroffset) * Teldef.d + + return xoffset, yoffset + + def offset_to_rawxy(self, xoffset, yoffset): + x = xoffset / Teldef.d + Teldef.centeroffset + 0.5 + y = yoffset / Teldef.d + Teldef.centeroffset + 0.5 + + # return x.astype(np.int) - (x < 0), y.astype(np.int) - (y < 0) + # return x.astype(np.int), y.astype(np.int) + return x.astype(int), y.astype(int) + + def offset_to_vec(self, x, y): + vec = np.empty(shape=x.shape+(3,), dtype=np.double) + vec[..., 0] = Teldef.F + vec[..., 1] = x + vec[..., 2] = -y + + return vec + + def vec_to_offset(self, vec): + x = vec[..., 0] + y = vec[..., 1] + z = vec[..., 2] + + xoffset = y * Teldef.F / x + yoffset = -z * Teldef.F / x + + return xoffset, yoffset + + def vec_to_pol(self, vec): + x = vec[..., 0] + y = vec[..., 1] + z = vec[..., 2] + h = np.hypot(x, y) + + phi = np.arctan2(y, x) % (2. * np.pi) + theta = np.arctan2(z, h) + + return phi, theta + + def pol_to_vec(self, phi, theta): + vec = np.empty(shape=phi.shape+(3,), dtype=np.double) + + vec[..., 0] = np.cos(theta) * np.cos(phi) + vec[..., 1] = np.cos(theta) * np.sin(phi) + vec[..., 2] = np.sin(theta) + + return vec + + def rawxy_to_radec(self, quat, rawx, rawy): + xoffset, yoffset \ + = self.rawxy_to_offset(rawx, rawy) + xyvec = self.offset_to_vec(xoffset, yoffset) + ra, dec = self.vec_to_pol(quat.apply(xyvec)) + + return ra, dec + + def radec_to_rawxy(self, quat, ra, dec): + skyvec = self.pol_to_vec(ra, dec) + detvec = quat.apply(skyvec, inverse=True) + xoffset, yoffset \ + = self.vec_to_offset(detvec) + rawx, rawy \ + = self.offset_to_rawxy(xoffset, yoffset) + + return rawx, rawy + + +if __name__ == '__main__': + pass diff --git a/artpipeline/src/artpipeline/pipeline.py b/artpipeline/src/artpipeline/pipeline.py new file mode 100755 index 0000000..e82a915 --- /dev/null +++ b/artpipeline/src/artpipeline/pipeline.py @@ -0,0 +1,518 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2021-2025 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 json +import shutil +import argparse + +import datetime + +import astropy.units as u +from astropy.time.formats import erfa +from astropy.time import Time +from astropy.time import TimezoneInfo + +TZ_UTC = TimezoneInfo(utc_offset=0*u.hour) # UTC time zone +TZ_MSK = TimezoneInfo(utc_offset=3*u.hour) # UTC+3 (Moscow) time zone + +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.telescope import Teldef +from arttools.random import Random +from arttools.arttime import ArtDay, MissionTime + +from artpipeline.version import __version__ + +from artpipeline.caldb import CaldbInfo +from artpipeline.caldb import CaldbMask +from artpipeline.caldb import CaldbEnergy +from artpipeline.caldb import CaldbOrientation + +from artpipeline.events import Events +from artpipeline.energy import Flag +from artpipeline.energy import Energy +from artpipeline.skyaxis import SkyAxis +from artpipeline.skycoord import SkyCoord + + +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', + help='newfile creation stem', + action=KvAction, + default=None + ) + argparser.add_argument( + '--stem-in', + type=str, + help='input files stem', + action=KvAction, + default=None + ) + argparser.add_argument( + 'srcdir', + help='input directory', + action=KvAction + ) + argparser.add_argument( + 'tmpdir', + help='temporary directory', + action=KvAction + ) + argparser.add_argument( + '--dstdir', + help='output directory', + action=KvAction, + default=None + ) + # argparser.add_argument( + # '--user-gti', + # type=str, + # help='user gti file', + # action=KvAction, + # default=None + # ) + argparser.add_argument( + '--seed', + type=str, + help='random generator seed', + action=KvAction, + default=None + ) + argparser.add_argument( + '--no-random', + help='turn randomization off', + action='store_true', + default=False + ) + argparser.add_argument( + '--plain-energy', + help='calculate energy on a plain list of files', + action='store_true', + default=False + ) + argparser.add_argument( + '--L0', + help='process L0 files', + action='store_true', + default=False + ) + args = argparser.parse_args() + + self.stem = args.stem + self.stem_in = args.stem_in + self.srcdir = args.srcdir + self.dstdir = args.dstdir + self.tmpdir = os.path.join(args.tmpdir, str(self.pid)) + self.tmpdst = os.path.join(self.tmpdir, 'dst') + # self.usergti = args.usergti + self.seed = args.seed + self.no_random = args.no_random + self.plain_energy = args.plain_energy + self.process_l0 = args.L0 + + if self.dstdir is None: + self.dstdir = self.srcdir + + # if self.seed is not None: + # self.seed = int(self.seed) + + def print(self): + print('Input parameters:') + print(' stem = {}'.format(self.stem )) + print(' stem-in = {}'.format(self.stem_in )) + print(' srcdir = {}'.format(self.srcdir )) + print(' dstdir = {}'.format(self.dstdir )) + print(' tmpdir = {}'.format(self.tmpdir )) + # print(' user-gti = {}'.format(self.usergti )) + print(' seed = {}'.format(self.seed )) + print(' no random = {}'.format(self.no_random )) + print(' plain energy = {}'.format(self.plain_energy)) + print(' peocess L0 = {}'.format(self.process_l0 )) + print() + + caldbinfo = CaldbInfo(self.caldbdir).get() + + print('CALDB version {} [{}]'.format(caldbinfo['CALDBVER'], caldbinfo['CALDBDIR'])) + print() + + def check(self): + if not Dir(self.srcdir).exists(): + raise TaskError('error: srcdir not found') + Dir(self.dstdir).check() + Dir(self.tmpdir).check(clean=True) + Dir(self.tmpdst).check(clean=True) + + +class Pipeline(Task): + + def __init__(self): + super().__init__() + + print('*** artpipeline v.{}'.format(__version__)) + + self._name = 'artpipeline v.{}'.format(__version__) + + self._args = Args() + + self._orifile = None + + def initialize(self): + self._args.read() + self._args.print() + self._args.check() + + if not self._args.no_random: + self._random = Random(self._args.seed) + else: + self._random = None + + def finalize(self): + Dir(self._args.tmpdir).delete() + + def runmain(self): + layout = self._layout() + + self._make_orient(layout) + self._make_attitude(layout) + self._process_events(layout) + + print('*** artpipeline finished') + + def _get_stem(self, fname, sep): + # get filenme from path + fn = os.path.basename(fname) + + #find teln in fname + fi = fn.find(sep) + + return fn[:fi] + + def _list_files(self, srcdir, stem_in=None): + allfiles = [fi for fi in os.listdir(srcdir)] + allfiles.sort() + + if stem_in==None: + return allfiles + + return [fi for fi in allfiles if fi.startswith(stem_in)] + + def _plain_layout(self): + layout = {} + + layout['type'] = 'plain' + layout['post'] = '' + layout['uf' ] = {} + + # pick events + srcdir = Dir(self._args.srcdir) + eventfiles = self._list_files(srcdir.get(), self._args.stem_in) + + for TEL in Teldef.telnames: + telfiles = [] + for fname in eventfiles: + if TEL in fname: + telfiles.append((srcdir / fname).get()) + + layout['uf'][TEL] = telfiles + + return layout + + def _L0_layout(self): + layout = {} + + layout['type'] = 'L0' + layout['post'] = '_urd.L2' + layout['uf' ] = {} + + # pick L0 events + ufdir = Dir(self._args.srcdir) + eventfiles = self._list_files(ufdir.get(), self._args.stem_in) + + for TEL in Teldef.telnames: + uffiles = [] + for fname in eventfiles: + if TEL in fname: + uffiles.append((ufdir / fname).get()) + + layout['uf'][TEL] = uffiles + + # pick orientation + orientation = {} + + auxdir = Dir(self._args.srcdir) + auxfiles = self._list_files(auxdir.get(), self._args.stem_in) + + for otype in ['gyro', 'sed1', 'sed2', 'bokz']: + for auxf in auxfiles: + if auxf.endswith('_{}.L0.fits'.format(otype)): + ofile = auxdir / auxf + + if ofile.exists(): + orientation[otype] = ofile.get() + + layout['orientation'] = orientation + + return layout + + def _science_layout(self): + layout = {} + + layout['type'] = 'science' + layout['post'] = '_cl' + layout['uf' ] = {} + + # pick uf events + ufdir = Dir(self._args.srcdir) / 'event_uf' + eventfiles = self._list_files(ufdir.get(), self._args.stem_in) + + for TEL in Teldef.telnames: + uffiles = [] + for fname in eventfiles: + if TEL in fname: + uffiles.append((ufdir / fname).get()) + + layout['uf'][TEL] = uffiles + + # pick orientation + orientation = {} + + auxdir = Dir(self._args.srcdir) / 'auxiliary' + auxfiles = self._list_files(auxdir.get(), self._args.stem_in) + + for otype in ['gyro', 'sed1', 'sed2', 'bokz']: + for auxf in auxfiles: + if auxf.endswith('_{}.fits'.format(otype)): + ofile = auxdir / auxf + + if ofile.exists(): + orientation[otype] = ofile.get() + + layout['orientation'] = orientation + + return layout + + def _layout(self): + if self._args.process_l0: + return self._L0_layout() + + if self._args.plain_energy: + return self._plain_layout() + + return self._science_layout() + + def _make_orient(self, layout): + if layout['type'] == 'plain': + print('mode: plain: skipping orienation...') + return + + print('making orienation file...') + gyrofile = layout['orientation']['gyro'] + + if self._args.stem is None: + stemout = self._get_stem(gyrofile, sep='_') + else: + stemout = self._args.stem + + if self._args.process_l0: + auxdir = Dir(self._args.dstdir) + else: + auxdir = Dir(self._args.dstdir) / 'auxiliary' + + auxdir.check() + + orifile = '{}_ori.fits'.format(stemout) + orifull = (auxdir / orifile).get() + + skyaxis = SkyAxis() + skyaxis.process(gyrofile).write(orifull) + + self._orifile = orifull + + print('written: {}'.format(orifull)) + + def _make_attitude(self, layout): + if layout['type'] == 'plain': + print('mode: plain: skipping attitude...') + return + + if self._args.stem is None: + stemout = self._get_stem(self._orifile, sep='_') + else: + stemout = self._args.stem + + print('making attitude file...') + attfile = '{}_att.fits'.format(stemout) + + caldb = CaldbOrientation(self._args.caldbdir, 'GYRO') + # skycoord = SkyCoord(caldb, 'GYRO') + skycoord = SkyCoord(caldb, self._orifile, 'GYRO') + + # skycoord.attitude(self._orifile).write(attfull) + + for TEL in Teldef.telnames: + skycoord.attitude(TEL).collect() + + if self._args.process_l0: + attfull = (Dir(self._args.dstdir) / attfile).get() + else: + attfull = (Dir(self._args.dstdir) / 'auxiliary' / attfile).get() + + skycoord.write(attfull) + + print('written: {}'.format(attfull)) + + def _process_events(self, layout): + processed = datetime.datetime.now(tz=TZ_MSK).strftime("%Y-%m-%d %H:%M:%S") + + metainfo = { + 'artpipeline': __version__, + 'processed' : processed + } + + if self._args.plain_energy: + metainfo['mode'] = 'plain' + elif self._args.process_l0: + metainfo['mode'] = 'L0' + else: + metainfo['mode'] = 'science' + + caldbinfo = CaldbInfo(self._args.caldbdir) + + metainfo['caldb'] = caldbinfo.get()['CALDBVER'] + + if not self._args.no_random: + random = self._random + metainfo['seed'] = random.seed + else: + random = None + metainfo['seed'] = None + + metainfo['events'] = {} + + for TEL in Teldef.telnames: + print('processing: {}'.format(TEL)) + + for fitsname in layout['uf'].get(TEL, []): + + print('event file: {}'.format(fitsname)) + + if self._args.stem is None: + stemout = self._get_stem(fitsname, TEL) + else: + stemout = self._args.stem + + # read events + events = Events(fitsname).read() + skyflag = np.zeros_like(events.events.column('TIME')) + + if self._args.plain_energy: + times = events.events.column('TIME') + events.setcol('RA' , np.zeros_like(times)) + events.setcol('DEC', np.zeros_like(times)) + else: + # calc ra/dec + ocaldb = CaldbOrientation(self._args.caldbdir, 'GYRO') + skycoord = SkyCoord(ocaldb, self._orifile) + + skyflag = skycoord.events(self._orifile, events, random) + + # calc energies + ecaldb = CaldbEnergy(self._args.caldbdir, TEL) + energy = Energy(ecaldb) + + en, enbot, entop, grade =\ + energy.calc(events.teln, events.events, self._random) + + events.events.setcol('ENERGY' , en ) + events.events.setcol('ENERGY_BOT', enbot) + events.events.setcol('ENERGY_TOP', entop) + events.events.setcol('GRADE' , grade) + + # calc flag + mcaldb = CaldbMask(self._args.caldbdir, TEL) + flag = Flag(mcaldb) + evtflag = flag.calc(events.events, events.tstart, events.tstop) + + events.events.setcol('FLAG', skyflag + evtflag) + + # print(events.events.frame) + + if self._args.plain_energy: + stemout = self._get_stem(fitsname, '.fits') + dstfile = '{}.L2.fits'.format(stemout) + else: + dstfile = '{}{}{}.fits'.format(stemout, events.telname, layout['post']) + + if self._args.plain_energy or self._args.process_l0: + dstdir = Dir(self._args.dstdir) + else: + dstdir = Dir(self._args.dstdir) / 'event_cl' + + dstdir.check() + + dstfull = (dstdir / dstfile).get() + events.write_cl(dstfull, caldbinfo.get()) + + aday_start = MissionTime(events.tstart).to_artday() + aday_stop = MissionTime(events.tstop ).to_artday() + + metainfo['events'][dstfile] = { + 'artday' : (aday_start, aday_stop), + 'mission': (events.tstart, events.tstop), + 'events' : events.size + } + + print('written: {}'.format(dstfull)) + + # write meta + dstdir = Dir(self._args.dstdir) + metafname = dstdir.pathfor('meta.info') + with open(metafname, 'w') as metafile: + json.dump(metainfo, metafile, indent=4) + + print('written: {}'.format(metafname)) + + +def main(): + try: + Pipeline().run() + except TaskError as e: + print('error: task execution:', e) + except Error as e: + print('error: unknown:', e) + + +if __name__ == '__main__': + pass + # PipelineTask().runall() diff --git a/artpipeline/src/artpipeline/pipeline.py.new.py b/artpipeline/src/artpipeline/pipeline.py.new.py new file mode 100755 index 0000000..d9b21d3 --- /dev/null +++ b/artpipeline/src/artpipeline/pipeline.py.new.py @@ -0,0 +1,573 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2021-2025 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 json +import shutil +import argparse + +import datetime + +import astropy.units as u +from astropy.time.formats import erfa +from astropy.time import Time +from astropy.time import TimezoneInfo + +TZ_UTC = TimezoneInfo(utc_offset=0*u.hour) # UTC time zone +TZ_MSK = TimezoneInfo(utc_offset=3*u.hour) # UTC+3 (Moscow) time zone + +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.telescope import Teldef +from arttools.random import Random +from artdactools.arttime import ArtDay, MissionTime + +from artpipeline.version import __version__ + +from artpipeline.caldb import CaldbInfo +from artpipeline.caldb import CaldbMask +from artpipeline.caldb import CaldbEnergy +from artpipeline.caldb import CaldbOrientation + +from artpipeline.events import Events +from artpipeline.energy import Flag +from artpipeline.energy import Energy +from artpipeline.skyaxis import SkyAxis +from artpipeline.skycoord import SkyCoord + + +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', + help='newfile creation stem', + action=KvAction, + default=None + ) + argparser.add_argument( + '--stem-in', + type=str, + help='input files stem', + action=KvAction, + default=None + ) + argparser.add_argument( + 'srcdir', + help='input directory', + action=KvAction + ) + argparser.add_argument( + 'tmpdir', + help='temporary directory', + action=KvAction + ) + argparser.add_argument( + '--dstdir', + help='output directory', + action=KvAction, + default=None + ) + # argparser.add_argument( + # '--user-gti', + # type=str, + # help='user gti file', + # action=KvAction, + # default=None + # ) + argparser.add_argument( + '--seed', + type=str, + help='random generator seed', + action=KvAction, + default=None + ) + argparser.add_argument( + '--no-random', + help='turn randomization off', + action='store_true', + default=False + ) + argparser.add_argument( + '--plain-energy', + help='calculate energy on a plain list of files', + action='store_true', + default=False + ) + argparser.add_argument( + '--L0', + help='process L0 files', + action='store_true', + default=False + ) + args = argparser.parse_args() + + self.stem = args.stem + self.stem_in = args.stem_in + self.srcdir = args.srcdir + self.dstdir = args.dstdir + self.tmpdir = os.path.join(args.tmpdir, str(self.pid)) + self.tmpdst = os.path.join(self.tmpdir, 'dst') + # self.usergti = args.usergti + self.seed = args.seed + self.no_random = args.no_random + self.plain_energy = args.plain_energy + self.process_l0 = args.L0 + + if self.dstdir is None: + self.dstdir = self.srcdir + + # if self.seed is not None: + # self.seed = int(self.seed) + + def print(self): + print('Input parameters:') + print(' stem = {}'.format(self.stem )) + print(' stem-in = {}'.format(self.stem_in )) + print(' srcdir = {}'.format(self.srcdir )) + print(' dstdir = {}'.format(self.dstdir )) + print(' tmpdir = {}'.format(self.tmpdir )) + # print(' user-gti = {}'.format(self.usergti )) + print(' seed = {}'.format(self.seed )) + print(' no random = {}'.format(self.no_random )) + print(' plain energy = {}'.format(self.plain_energy)) + print(' peocess L0 = {}'.format(self.process_l0 )) + print() + + caldbinfo = CaldbInfo(self.caldbdir).get() + + print('CALDB version {} [{}]'.format(caldbinfo['CALDBVER'], caldbinfo['CALDBDIR'])) + print() + + def check(self): + if not Dir(self.srcdir).exists(): + raise TaskError('error: srcdir not found') + Dir(self.dstdir).check() + Dir(self.tmpdir).check(clean=True) + Dir(self.tmpdst).check(clean=True) + + +class PipeConfigL0(object): + + def __init__(self, args): + self._args = args + self._layout = self._get_layout() + + def _get_layout(self): + layout = {} + + layout['type'] = 'L0' + layout['post'] = '_urd.L2' + layout['uf' ] = {} + + # pick L0 events + ufdir = Dir(self._args.srcdir) + eventfiles = self._list_files(ufdir.get(), self._args.stem_in) + + for TEL in Teldef.telnames: + uffiles = [] + for fname in eventfiles: + if TEL in fname: + uffiles.append((ufdir / fname).get()) + + layout['uf'][TEL] = uffiles + + # pick orientation + orientation = {} + + auxdir = Dir(self._args.srcdir) + auxfiles = self._list_files(auxdir.get(), self._args.stem_in) + + for otype in ['gyro', 'sed1', 'sed2', 'bokz']: + for auxf in auxfiles: + if auxf.endswith('_{}.fits'.format(otype)): + ofile = auxdir / auxf + + if ofile.exists(): + orientation[otype] = ofile.get() + + layout['orientation'] = orientation + + return layout + + @property + def layout(self): return self._layout + + +class PipeConfigPlain(object): + + def __init__(self, args): + self._args = args + self._layout = self._get_layout() + + def _list_files(self, srcdir, stem_in=None): + allfiles = [fi for fi in os.listdir(srcdir)] + allfiles.sort() + + if stem_in==None: + return allfiles + + return [fi for fi in allfiles if fi.startswith(stem_in)] + + def _get_layout(self): + layout = {} + + layout['type'] = 'plain' + layout['post'] = '' + layout['uf' ] = {} + + # pick events + srcdir = Dir(self._args.srcdir) + eventfiles = self._list_files(srcdir.get(), self._args.stem_in) + + for TEL in Teldef.telnames: + telfiles = [] + for fname in eventfiles: + if TEL in fname: + telfiles.append((srcdir / fname).get()) + + layout['uf'][TEL] = telfiles + + return layout + + @property + def layout(self): return self._layout + + +class PipeConfigScience(object): + + def __init__(self, args): + self._args = args + self._layout = self._get_layout() + + def _list_files(self, srcdir, stem_in=None): + allfiles = [fi for fi in os.listdir(srcdir)] + allfiles.sort() + + if stem_in==None: + return allfiles + + return [fi for fi in allfiles if fi.startswith(stem_in)] + + def _get_layout(self): + layout = {} + + layout['type'] = 'science' + layout['post'] = '_cl' + layout['uf' ] = {} + + # pick uf events + ufdir = Dir(self._args.srcdir) / 'event_uf' + eventfiles = self._list_files(ufdir.get(), self._args.stem_in) + + for TEL in Teldef.telnames: + uffiles = [] + for fname in eventfiles: + if TEL in fname: + uffiles.append((ufdir / fname).get()) + + layout['uf'][TEL] = uffiles + + # pick orientation + orientation = {} + + auxdir = Dir(self._args.srcdir) / 'auxiliary' + auxfiles = self._list_files(auxdir.get(), self._args.stem_in) + + for otype in ['gyro', 'sed1', 'sed2', 'bokz']: + for auxf in auxfiles: + if auxf.endswith('_{}.fits'.format(otype)): + ofile = auxdir / auxf + + if ofile.exists(): + orientation[otype] = ofile.get() + + layout['orientation'] = orientation + + return layout + + @property + def layout(self): return self._layout + + +class EventProcContext(object): + pass + + +class Pipeline(Task): + + def __init__(self): + super().__init__() + + print('*** artpipeline v.{}'.format(__version__)) + + self._name = 'artpipeline v.{}'.format(__version__) + + self._args = Args() + + self._config = None + self._random = None + self._ofile = None + + def initialize(self): + self._args.read() + self._args.print() + self._args.check() + + if not self._args.no_random: + self._random = Random(self._args.seed) + + if self._args.plain_energy: + self._config = PipeConfigPlain(self._args) + elif self._args.process_l0: + self._config = PipeConfigL0(self._args) + else: + self._config = PipeConfigScience(self._args) + + def finalize(self): + Dir(self._args.tmpdir).delete() + + def runmain(self): + # layout = self._layout() + + layout = self._config.layout() + + self._make_orient(layout) + self._make_attitude(layout) + + self._process_events(layout) + + print('*** artpipeline finished') + + def _get_stem(self, fname, sep): + # get filenme from path + fn = os.path.basename(fname) + + #find teln in fname + fi = fn.find(sep) + + return fn[:fi] + + def _list_files(self, srcdir, stem_in=None): + allfiles = [fi for fi in os.listdir(srcdir)] + allfiles.sort() + + if stem_in==None: + return allfiles + + return [fi for fi in allfiles if fi.startswith(stem_in)] + + + def _layout(self): + if self._args.process_l0: + return self._L0_layout() + + if self._args.plain_energy: + return self._plain_layout() + + return self._science_layout() + + def _make_orient(self, layout): + if layout['type'] == 'plain': + print('mode: plain: skipping orienation...') + return + + print('making orienation file...') + gyrofile = layout['orientation']['gyro'] + + if self._args.stem is None: + stemout = self._get_stem(gyrofile, sep='_') + else: + stemout = self._args.stem + + if self._args.process_l0: + auxdir = Dir(self._args.dstdir) + else: + auxdir = Dir(self._args.dstdir) / 'auxiliary' + + auxdir.check() + + ofile = '{}_ori.fits'.format(stemout) + ofull = (auxdir / ofile).get() + + skyaxis = SkyAxis() + skyaxis.process(gyrofile).write(ofull) + + self._ofile = ofull + + print('written: {}'.format(ofull)) + + def _make_attitude(self, layout): + if layout['type'] == 'plain': + print('mode: plain: skipping attitude...') + return + + if self._args.stem is None: + stemout = self._get_stem(self._orifile, sep='_') + else: + stemout = self._args.stem + + print('making attitude file...') + attfile = '{}_att.fits'.format(stemout) + + if self._args.process_l0: + attfull = (Dir(self._args.dstdir) / attfile).get() + else: + attfull = (Dir(self._args.dstdir) / 'auxiliary' / attfile).get() + + caldb = CaldbOrientation(self._args.caldbdir, 'GYRO') + skycoord = SkyCoord(caldb, 'GYRO') + + skycoord.attitude(self._orifile).write(attfull) + + print('written: {}'.format(attfull)) + + def _process_events(self, layout): + processed = datetime.datetime.now(tz=TZ_MSK).strftime("%Y-%m-%d %H:%M:%S") + + metainfo = { + 'artpipeline': __version__, + 'processed' : processed + } + + if self._args.plain_energy: + metainfo['mode'] = 'plain' + elif self._args.process_l0: + metainfo['mode'] = 'L0' + else: + metainfo['mode'] = 'science' + + caldbinfo = CaldbInfo(self._args.caldbdir) + + metainfo['caldb'] = caldbinfo.get()['CALDBVER'] + + if not self._args.no_random: + random = self._random + metainfo['seed'] = random.seed + else: + random = None + metainfo['seed'] = None + + metainfo['events'] = {} + + for TEL in Teldef.telnames: + print('processing: {}'.format(TEL)) + + for fitsname in layout['uf'].get(TEL, []): + + print('event file: {}'.format(fitsname)) + + if self._args.stem is None: + stemout = self._get_stem(fitsname, TEL) + else: + stemout = self._args.stem + + # read events + events = Events(fitsname).read() + + if self._args.plain_energy: + times = events.events.column('TIME') + events.setcol('RA' , np.zeros_like(times)) + events.setcol('DEC', np.zeros_like(times)) + else: + # calc ra/dec + ocaldb = CaldbOrientation(self._args.caldbdir, 'GYRO') + skycoord = SkyCoord(ocaldb, 'GYRO') + + skycoord.events(self._orifile, events, random) + + # calc energies + ecaldb = CaldbEnergy(self._args.caldbdir, TEL) + energy = Energy(ecaldb) + + en, enbot, entop, grade =\ + energy.calc(events.teln, events.events, self._random) + + events.events.setcol('ENERGY' , en ) + events.events.setcol('ENERGY_BOT', enbot) + events.events.setcol('ENERGY_TOP', entop) + events.events.setcol('GRADE' , grade) + + # calc flag + mcaldb = CaldbMask(self._args.caldbdir, TEL) + flag = Flag(mcaldb) + f = flag.calc(events.events, events.tstart, events.tstop) + + events.events.setcol('FLAG', f) + + # print(events.events.frame) + + if self._args.plain_energy: + stemout = self._get_stem(fitsname, '.fits') + dstfile = '{}.L2.fits'.format(stemout) + else: + dstfile = '{}{}{}.fits'.format(stemout, events.telname, layout['post']) + + if self._args.plain_energy or self._args.process_l0: + dstdir = Dir(self._args.dstdir) + else: + dstdir = Dir(self._args.dstdir) / 'event_cl' + + dstdir.check() + + dstfull = (dstdir / dstfile).get() + events.write_cl(dstfull, caldbinfo.get()) + + aday_start = MissionTime(events.tstart).to_artday() + aday_stop = MissionTime(events.tstop ).to_artday() + + metainfo['events'][dstfile] = { + 'artday' : (aday_start, aday_stop), + 'mission': (events.tstart, events.tstop), + 'events' : events.size + } + + print('written: {}'.format(dstfull)) + + # write meta + metafname = (Dir(self._args.dstdir) / 'meta.info').get() + with open(metafname, 'w') as metafile: + json.dump(metainfo, metafile, indent=4) + + print('written: {}'.format(metafname)) + + +def main(): + try: + Pipeline().run() + except TaskError as e: + print('error: task execution:', e) + except Error as e: + print('error: unknown:', e) + + +if __name__ == '__main__': + pass + # PipelineTask().runall() diff --git a/artpipeline/src/artpipeline/skyaxis.py b/artpipeline/src/artpipeline/skyaxis.py new file mode 100755 index 0000000..fbb5c8b --- /dev/null +++ b/artpipeline/src/artpipeline/skyaxis.py @@ -0,0 +1,34 @@ +#!/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 shutil + + +class SkyAxis(object): + + def __init__(self): + super().__init__() + + self._srcfile = None + + def process(self, srcfile): + self._srcfile = srcfile + return self + + def write(self, dstfile): + + if os.path.exists(dstfile): + os.remove(dstfile) + + shutil.copy(self._srcfile, dstfile) + + +if __name__ == '__main__': + pass diff --git a/artpipeline/src/artpipeline/skycoord.py b/artpipeline/src/artpipeline/skycoord.py new file mode 100755 index 0000000..9528b16 --- /dev/null +++ b/artpipeline/src/artpipeline/skycoord.py @@ -0,0 +1,194 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2021-2025 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 numpy as np +from astropy.io import fits + +from arttools.gti import Gti +from arttools.arttime import MJDREF +from arttools.telescope import Teldef + +from artpipeline.orientation import Orientation +from artpipeline.orientation import AttitudeStatus + +from artpipeline.events import Events + +from artpipeline.version import __version__ + + +class SkyCoord(object): + + def __init__(self, caldb, orifile, instname=None): + super().__init__() + + self._caldb = caldb + self._ofile = orifile + + self._instname = instname + self._inst = None + + self._mjdref = MJDREF.mjd + + self._data = {} + self._adata = {} + + self._process_orient() + self._calc_base_attitude() + + def _process_orient(self): + self._orient = Orientation(self._caldb) + self._odata = self._orient.read(self._ofile) + self._gti = Gti.from_hduname(self._ofile, 'STDGTI') + + self._otime = self._odata['time'] + self._ostatus = AttitudeStatus.pack(self._odata['state']) + + def _calc_base_attitude(self): + if self._instname is None: + return + + att_ra, att_dec, att_roll \ + = self._orient.calc_attitudes(self._odata, self._instname) + + self._data['ATTITUDE'] = { + 'time' : self._otime, + 'ra' : att_ra , + 'dec' : att_dec , + 'roll' : att_roll , + 'status': self._ostatus + } + + x_ra, x_dec, x_roll \ + = self._orient.calc_attitudes(self._odata) + + self._data['SCX'] = { + 'time' : self._otime, + 'ra' : x_ra , + 'dec' : x_dec , + 'roll' : x_roll , + 'status': self._ostatus + } + + def attitude(self, instname=None): + + if instname is not None: + self._inst = instname + else: + self._inst = self._instname + + ra, dec, roll \ + = self._orient.calc_attitudes(self._odata, self._inst) + + self._adata[self._inst] = { + 'time' : self._otime, + 'ra' : ra , + 'dec' : dec , + 'roll' : roll , + 'status': self._ostatus + } + + return self + + def collect(self, name=None): + if self._inst is None: + return + + if name is None: + name = self._inst + + self._data[name] = self._adata.pop(self._inst, None) + self._inst = None + + def _make_hdu(self, hduname): + + if hduname not in self._data: + + hdu = fits.BinTableHDU() + hdu.name = hduname + + return hdu + + # --- HDU --- + a00 = np.array(self._data[hduname]['time' ], dtype=np.double) + a01 = np.array(self._data[hduname]['ra' ], dtype=np.double) + a02 = np.array(self._data[hduname]['dec' ], dtype=np.double) + a03 = np.array(self._data[hduname]['roll' ], dtype=np.double) + a04 = np.array(self._data[hduname]['status'], dtype=int ) + a05 = np.zeros(a00.size, dtype=int) + + c00 = fits.Column(name='TIME' , format='1D', unit='sec', array=a00) + c01 = fits.Column(name='RA' , format='1D', unit='deg', array=a01) + c02 = fits.Column(name='DEC' , format='1D', unit='deg', array=a02) + c03 = fits.Column(name='ROLL' , format='1D', unit='deg', array=a03) + c04 = fits.Column(name='STATUS', format='1I', unit='' , array=a04) + c05 = fits.Column(name='FLAG' , format='1I', unit='' , array=a05) + + cols = fits.ColDefs([ + c00, c01, c02, c03, c04, c05 + ]) + + hdu = fits.BinTableHDU.from_columns(cols) + hdu.name = hduname + + creator = 'artpipeline v.{}'.format(__version__) + + hdu.header['CREATOR' ] = (creator , 'program and version that created this file' ) + hdu.header['ORIGIN' ] = ('IKI RAS' , 'institution that created this file' ) + hdu.header['TELESCOP'] = ('SRG/ART-XC', 'mission name' ) + hdu.header['MJDREF' ] = (self._mjdref, 'MJD corresponding to SC clock start 2000.0 MSK') + hdu.header['TSTART' ] = self._gti.start + hdu.header['TSTOP' ] = self._gti.stop + hdu.header['CALDBVER'] = '{}'.format(self._caldb.version()) + + return hdu + + def write(self, attfile): + + print(self._data.keys()) + + prihdu = fits.PrimaryHDU() + + # --- make file --- + hdulist = fits.HDUList([ + prihdu, + self._make_hdu('ATTITUDE'), + self._make_hdu('T1' ), + self._make_hdu('T2' ), + self._make_hdu('T3' ), + self._make_hdu('T4' ), + self._make_hdu('T5' ), + self._make_hdu('T6' ), + self._make_hdu('T7' ), + self._make_hdu('SCX' ), + self._gti.make_hdu(self._mjdref) + ]) + + hdulist.writeto(attfile, overwrite=True, checksum=True) + + + def events(self, orifile, events, random=None): + ra, dec, mask \ + = self._orient.calc_events_coords( + self._odata, events.events, events.telname, random) + + + flag = np.zeros_like(events.events.column('TIME' )) + flag[~mask] = 100 + + events.setcol('RA' , ra ) + events.setcol('DEC' , dec ) + + return flag + + +if __name__ == '__main__': + pass diff --git a/artpipeline/src/artpipeline/version.py b/artpipeline/src/artpipeline/version.py new file mode 100755 index 0000000..a33997d --- /dev/null +++ b/artpipeline/src/artpipeline/version.py @@ -0,0 +1 @@ +__version__ = '2.1.0' diff --git a/artproducts/.DS_Store b/artproducts/.DS_Store new file mode 100644 index 0000000..e9030b3 Binary files /dev/null and b/artproducts/.DS_Store differ diff --git a/artproducts/pyproject.toml b/artproducts/pyproject.toml new file mode 100644 index 0000000..992522d --- /dev/null +++ b/artproducts/pyproject.toml @@ -0,0 +1,6 @@ +[build-system] +requires = [ + "setuptools>=42", +] +build-backend = "setuptools.build_meta" + diff --git a/artproducts/setup.py b/artproducts/setup.py new file mode 100644 index 0000000..5cb1b3c --- /dev/null +++ b/artproducts/setup.py @@ -0,0 +1,17 @@ +from setuptools import setup + +exec(open('src/artproducts/version.py').read()) + +setup( + name="artproducts", + version=__version__, + author="M.Pavlinsky SRG/ART-XC software team", + description="high-level products extraction", + package_dir={"": "src"}, + packages=["artproducts"], + entry_points={ + "console_scripts": [ + "artproducts = artproducts.products:main", + ] + }, +) diff --git a/artproducts/src/.DS_Store b/artproducts/src/.DS_Store new file mode 100644 index 0000000..1eee0ee Binary files /dev/null and b/artproducts/src/.DS_Store differ diff --git a/artproducts/src/artproducts/.DS_Store b/artproducts/src/artproducts/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/artproducts/src/artproducts/.DS_Store differ diff --git a/artproducts/src/artproducts/enspectrum.py b/artproducts/src/artproducts/enspectrum.py new file mode 100755 index 0000000..7b14ce8 --- /dev/null +++ b/artproducts/src/artproducts/enspectrum.py @@ -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 diff --git a/artproducts/src/artproducts/eventspack.py b/artproducts/src/artproducts/eventspack.py new file mode 100755 index 0000000..0c7c54e --- /dev/null +++ b/artproducts/src/artproducts/eventspack.py @@ -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 diff --git a/artproducts/src/artproducts/eventswcs.py b/artproducts/src/artproducts/eventswcs.py new file mode 100755 index 0000000..0a59136 --- /dev/null +++ b/artproducts/src/artproducts/eventswcs.py @@ -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 diff --git a/artproducts/src/artproducts/histograms.py b/artproducts/src/artproducts/histograms.py new file mode 100755 index 0000000..8b63243 --- /dev/null +++ b/artproducts/src/artproducts/histograms.py @@ -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 diff --git a/artproducts/src/artproducts/image.py b/artproducts/src/artproducts/image.py new file mode 100755 index 0000000..37bd820 --- /dev/null +++ b/artproducts/src/artproducts/image.py @@ -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 diff --git a/artproducts/src/artproducts/lightcurve.py b/artproducts/src/artproducts/lightcurve.py new file mode 100755 index 0000000..b2e6cf1 --- /dev/null +++ b/artproducts/src/artproducts/lightcurve.py @@ -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 diff --git a/artproducts/src/artproducts/products.py b/artproducts/src/artproducts/products.py new file mode 100755 index 0000000..9d86e87 --- /dev/null +++ b/artproducts/src/artproducts/products.py @@ -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(' 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(' file not set') + + if not os.path.exists(self._args.orifile): + raise TaskError(' file not found') + + if self._args.srcra is None or self._args.srcdec is None: + raise TaskError('/ 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('') + 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('') + 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('') + 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('') + 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 + diff --git a/artproducts/src/artproducts/regions.py b/artproducts/src/artproducts/regions.py new file mode 100755 index 0000000..a35d325 --- /dev/null +++ b/artproducts/src/artproducts/regions.py @@ -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 diff --git a/artproducts/src/artproducts/responces.py b/artproducts/src/artproducts/responces.py new file mode 100755 index 0000000..de54e65 --- /dev/null +++ b/artproducts/src/artproducts/responces.py @@ -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 diff --git a/artproducts/src/artproducts/version.py b/artproducts/src/artproducts/version.py new file mode 100755 index 0000000..afced14 --- /dev/null +++ b/artproducts/src/artproducts/version.py @@ -0,0 +1 @@ +__version__ = '2.0.0' diff --git a/arttools/.DS_Store b/arttools/.DS_Store new file mode 100644 index 0000000..5873551 Binary files /dev/null and b/arttools/.DS_Store differ diff --git a/arttools/pyproject.toml b/arttools/pyproject.toml new file mode 100644 index 0000000..992522d --- /dev/null +++ b/arttools/pyproject.toml @@ -0,0 +1,6 @@ +[build-system] +requires = [ + "setuptools>=42", +] +build-backend = "setuptools.build_meta" + diff --git a/arttools/setup.py b/arttools/setup.py new file mode 100644 index 0000000..e9213c5 --- /dev/null +++ b/arttools/setup.py @@ -0,0 +1,19 @@ +from setuptools import setup + +# __version__ = "1.0.0" + +exec(open('src/arttools/version.py').read()) + +setup( + name="arttools", + version=__version__, + author="M.Pavlinsky SRG/ART-XC software team", + description="SRG/ARTXC software tools", + package_dir={"": "src"}, + packages=["arttools"], + entry_points={ + "console_scripts": [ + "arttime = arttools.arttime:main" + ] + } +) diff --git a/arttools/src/.DS_Store b/arttools/src/.DS_Store new file mode 100644 index 0000000..e314a59 Binary files /dev/null and b/arttools/src/.DS_Store differ diff --git a/arttools/src/arttools/.DS_Store b/arttools/src/arttools/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/arttools/src/arttools/.DS_Store differ diff --git a/arttools/src/arttools/__init__.py b/arttools/src/arttools/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/arttools/src/arttools/argparse.py b/arttools/src/arttools/argparse.py new file mode 100755 index 0000000..8c33da1 --- /dev/null +++ b/arttools/src/arttools/argparse.py @@ -0,0 +1,27 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2021-2024 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 argparse + + +class KvAction(argparse.Action): + def __init__(self, option_strings, dest, **kwargs): + super(KvAction, self).__init__(option_strings, dest, **kwargs) + + def __call__(self, parser, namespace, values, option_string=None): + if "=" in values: + key, val = values.split("=", 1) + else: + key, val = self.dest, values + + setattr(namespace, key, val) + + +if __name__ == '__main__': + pass diff --git a/arttools/src/arttools/art.py b/arttools/src/arttools/art.py new file mode 100755 index 0000000..edd8802 --- /dev/null +++ b/arttools/src/arttools/art.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# +# ----------------------------------------------------------------------------- +# Copyright (c) 2024 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the SRG/ART-XC software project. +# ----------------------------------------------------------------------------- + + + + + +if __name__ == '__main__': + pass + + diff --git a/arttools/src/arttools/artindex.py b/arttools/src/arttools/artindex.py new file mode 100755 index 0000000..81f972f --- /dev/null +++ b/arttools/src/arttools/artindex.py @@ -0,0 +1,623 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2024 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the SRG/ART-XC software project. +# ----------------------------------------------------------------------------- + +import os +import shutil +from datetime import datetime + +import json +from dataclasses import dataclass + +import numpy as np +from astropy.io import fits + +import polars as pl + +import astropy.units as u +from astropy.time.formats import erfa +from astropy.time import Time +from astropy.time import TimezoneInfo + +from arttools.arttime import Dtime +from arttools.arttime import ArtDay + +from arttools.version import __version__ + + +MJDREFD = 51543.875 + +MJDREF = Time(51543.875, format='mjd') # 01.01.2000 00:00:00 (UTC+3) +TZ_UTC = TimezoneInfo(utc_offset=0*u.hour) # UTC time zone +TZ_MSK = TimezoneInfo(utc_offset=3*u.hour) # UTC+3 (Moscow) time zone + + +@dataclass +class SessionIndex: + def __init__(self, fname): + self.key = [] + self.npkey = np.array(self.key) + + self.stem = [] + self.station = [] + self.sequence = [] + self.data_quality = [] + self.total_quality = [] + self.processed = [] + self.start = [] + self.stop = [] + self.T1_start = [] + self.T1_stop = [] + self.T2_start = [] + self.T2_stop = [] + self.T3_start = [] + self.T3_stop = [] + self.T4_start = [] + self.T4_stop = [] + self.T5_start = [] + self.T5_stop = [] + self.T6_start = [] + self.T6_stop = [] + self.T7_start = [] + self.T7_stop = [] + self.flag = [] + + self._index_file = fname + + if os.path.isfile(self._index_file): + self.read(self._index_file) + + def read(self, fname): + with fits.open(fname) as hdul: + self.stem = hdul[1].data['STEM'].tolist() + self.station = hdul[1].data['STATION'].tolist() + self.sequence = hdul[1].data['SEQUENCE'].tolist() + self.data_quality = hdul[1].data['DATA_QUALITY'].tolist() + self.total_quality = hdul[1].data['TOTAL_QUALITY'].tolist() + self.processed = hdul[1].data['PROCESSED'].tolist() + self.start = hdul[1].data['START' ].tolist() + self.stop = hdul[1].data['STOP' ].tolist() + self.T1_start = hdul[1].data['T1_START'].tolist() + self.T1_stop = hdul[1].data['T1_STOP' ].tolist() + self.T2_start = hdul[1].data['T2_START'].tolist() + self.T2_stop = hdul[1].data['T2_STOP' ].tolist() + self.T3_start = hdul[1].data['T3_START'].tolist() + self.T3_stop = hdul[1].data['T3_STOP' ].tolist() + self.T4_start = hdul[1].data['T4_START'].tolist() + self.T4_stop = hdul[1].data['T4_STOP' ].tolist() + self.T5_start = hdul[1].data['T5_START'].tolist() + self.T5_stop = hdul[1].data['T5_STOP' ].tolist() + self.T6_start = hdul[1].data['T6_START'].tolist() + self.T6_stop = hdul[1].data['T6_STOP' ].tolist() + self.T7_start = hdul[1].data['T7_START'].tolist() + self.T7_stop = hdul[1].data['T7_STOP' ].tolist() + self.flag = hdul[1].data['FLAG'].tolist() + + for stm, sta, seq in zip(self.stem, self.station, self.sequence): + self.key.append('.'.join((stm, sta, seq))) + + self.npkey = np.array(self.key) + + def process(self, frame): + frame.make_key() + + if frame.key in self.npkey: + index = np.where(self.npkey == frame.key)[0][0] + + T1_start = frame.intervals['T1'][0] + T1_stop = frame.intervals['T1'][1] + T2_start = frame.intervals['T2'][0] + T2_stop = frame.intervals['T2'][1] + T3_start = frame.intervals['T3'][0] + T3_stop = frame.intervals['T3'][1] + T4_start = frame.intervals['T4'][0] + T4_stop = frame.intervals['T4'][1] + T5_start = frame.intervals['T5'][0] + T5_stop = frame.intervals['T5'][1] + T6_start = frame.intervals['T6'][0] + T6_stop = frame.intervals['T6'][1] + T7_start = frame.intervals['T7'][0] + T7_stop = frame.intervals['T7'][1] + + times_ok = [ + np.isclose(self.T1_start[index], T1_start), + np.isclose(self.T1_stop [index], T1_stop ), + np.isclose(self.T2_start[index], T2_start), + np.isclose(self.T2_stop [index], T2_stop ), + np.isclose(self.T3_start[index], T3_start), + np.isclose(self.T3_stop [index], T3_stop ), + np.isclose(self.T4_start[index], T4_start), + np.isclose(self.T4_stop [index], T4_stop ), + np.isclose(self.T5_start[index], T5_start), + np.isclose(self.T5_stop [index], T5_stop ), + np.isclose(self.T6_start[index], T6_start), + np.isclose(self.T6_stop [index], T6_stop ), + np.isclose(self.T7_start[index], T7_start), + np.isclose(self.T7_stop [index], T7_stop ), + ] + + quality_ok = [ + np.isclose(self.data_quality[index] , frame.data_quality ), + np.isclose(self.total_quality[index], frame.total_quality), + ] + + all_ok = np.all(times_ok) and np.all(quality_ok) + if not all_ok: + print( + "WARNING: session {}.{}.{} differs from indexed".format( + frame.stem, frame.station, frame.sequence + ) + ) + + else: + self.key.append(frame.key) + self.npkey = np.array(self.key) + + self.stem.append(frame.stem) + self.station.append(frame.station) + self.sequence.append(frame.sequence) + self.data_quality.append(frame.data_quality) + self.total_quality.append(frame.total_quality) + self.processed.append(frame.processed) + + start = np.array([ + frame.intervals["T1"][0], + frame.intervals["T2"][0], + frame.intervals["T3"][0], + frame.intervals["T4"][0], + frame.intervals["T5"][0], + frame.intervals["T6"][0], + frame.intervals["T7"][0], + ]).min() + + self.start.append(start) + + stop = np.array([ + frame.intervals["T1"][1], + frame.intervals["T2"][1], + frame.intervals["T3"][1], + frame.intervals["T4"][1], + frame.intervals["T5"][1], + frame.intervals["T6"][1], + frame.intervals["T7"][1], + ]).max() + + self.stop.append(stop ) + + self.T1_start.append(frame.intervals["T1"][0]) + self.T1_stop .append(frame.intervals["T1"][1]) + self.T2_start.append(frame.intervals["T2"][0]) + self.T2_stop .append(frame.intervals["T2"][1]) + self.T3_start.append(frame.intervals["T3"][0]) + self.T3_stop .append(frame.intervals["T3"][1]) + self.T4_start.append(frame.intervals["T4"][0]) + self.T4_stop .append(frame.intervals["T4"][1]) + self.T5_start.append(frame.intervals["T5"][0]) + self.T5_stop .append(frame.intervals["T5"][1]) + self.T6_start.append(frame.intervals["T6"][0]) + self.T6_stop .append(frame.intervals["T6"][1]) + self.T7_start.append(frame.intervals["T7"][0]) + self.T7_stop .append(frame.intervals["T7"][1]) + + # a FLAG is always 0 as we introduce session to index + self.flag.append(0) + + def index(self): + return pl.DataFrame( + { + 'key': self.key, + 'stem': self.stem, + 'station': self.station, + 'sequence': self.sequence, + 'data_quality': self.data_quality, + 'total_quality': self.total_quality, + 'processed': self.processed, + 'start' : self.start , + 'stop' : self.stop , + 'T1_start' : self.T1_start , + 'T1_stop' : self.T1_stop , + 'T2_start' : self.T2_start , + 'T2_stop' : self.T2_stop , + 'T3_start' : self.T3_start , + 'T3_stop' : self.T3_stop , + 'T4_start' : self.T4_start , + 'T4_stop' : self.T4_stop , + 'T5_start' : self.T5_start , + 'T5_stop' : self.T5_stop , + 'T6_start' : self.T6_start , + 'T6_stop' : self.T6_stop , + 'T7_start' : self.T7_start , + 'T7_stop' : self.T7_stop , + 'flag' : self.flag , + } + ).sort('key') + + def write_index_fits(self): + df = self.index() + + hdu = fits.PrimaryHDU() + + c00 = fits.Column(name='STEM', format='20A', unit='', array=np.array(df['stem'])) + c01 = fits.Column( + name='STATION', format='2A', unit='', array=np.array(df['station']) + ) + c02 = fits.Column( + name='SEQUENCE', format='3A', unit='', array=np.array(df['sequence']) + ) + c03 = fits.Column(name='FLAG' , format='1I' , unit='' , array=np.array(df['flag' ])) + c04 = fits.Column( + name='DATA_QUALITY', + format='1D', + unit='%', + array=np.array(df['data_quality']), + ) + c05 = fits.Column( + name='TOTAL_QUALITY', + format='1D', + unit='%', + array=np.array(df['total_quality'])) + c06 = fits.Column(name='PROCESSED', format='19A', unit='date', array=np.array(df['processed'])) + c07 = fits.Column(name='START' , format='1D' , unit='sec' , array=np.array(df['start' ])) + c08 = fits.Column(name='STOP' , format='1D' , unit='sec' , array=np.array(df['stop' ])) + c09 = fits.Column(name='T1_START' , format='1D' , unit='sec' , array=np.array(df['T1_start' ])) + c10 = fits.Column(name='T1_STOP' , format='1D' , unit='sec' , array=np.array(df['T1_stop' ])) + c11 = fits.Column(name='T2_START' , format='1D' , unit='sec' , array=np.array(df['T2_start' ])) + c12 = fits.Column(name='T2_STOP' , format='1D' , unit='sec' , array=np.array(df['T2_stop' ])) + c13 = fits.Column(name='T3_START' , format='1D' , unit='sec' , array=np.array(df['T3_start' ])) + c14 = fits.Column(name='T3_STOP' , format='1D' , unit='sec' , array=np.array(df['T3_stop' ])) + c15 = fits.Column(name='T4_START' , format='1D' , unit='sec' , array=np.array(df['T4_start' ])) + c16 = fits.Column(name='T4_STOP' , format='1D' , unit='sec' , array=np.array(df['T4_stop' ])) + c17 = fits.Column(name='T5_START' , format='1D' , unit='sec' , array=np.array(df['T5_start' ])) + c18 = fits.Column(name='T5_STOP' , format='1D' , unit='sec' , array=np.array(df['T5_stop' ])) + c19 = fits.Column(name='T6_START' , format='1D' , unit='sec' , array=np.array(df['T6_start' ])) + c20 = fits.Column(name='T6_STOP' , format='1D' , unit='sec' , array=np.array(df['T6_stop' ])) + c21 = fits.Column(name='T7_START' , format='1D' , unit='sec' , array=np.array(df['T7_start' ])) + c22 = fits.Column(name='T7_STOP' , format='1D' , unit='sec' , array=np.array(df['T7_stop' ])) + + cols = fits.ColDefs( + [ + c00, c01, c02, c03, c04, c05, c06, c07, c08, c09, + c10, c11, c12, c13, c14, c15, c16, c17, c18, c19, + c20, c21, c22 + ] + ) + + indexhdu = fits.BinTableHDU.from_columns(cols) + indexhdu.name = "INDEX" + + creator = 'artindex v.{}'.format(__version__) + indexhdu.header['creator' ] = (creator, 'program and version that created this file') + + dtnow = datetime.now(tz=TZ_MSK) + + indexhdu.header['mjdref' ] = (MJDREFD , 'MJD corresponding to SC clock start 2000.0 MSK' ) + indexhdu.header['date' ] = (dtnow.strftime('%Y/%m/%d'), 'date that FITS file was created' ) + + hdulist = fits.HDUList([hdu, indexhdu]) + hdulist.writeto(self._index_file, overwrite=True) + + print("written: {}".format(self._index_file)) + + +@dataclass +class SessionIndexFrame: + def __init__(self): + self.key = "" + + self.stem = "" + self.station = "00" + self.sequence = "000" + self.data_quality = 0 + self.total_quality = 0 + + self.processed = "" + + self.intervals = { + "T1": (0, 0), + "T2": (0, 0), + "T3": (0, 0), + "T4": (0, 0), + "T5": (0, 0), + "T6": (0, 0), + "T7": (0, 0), + } + + self.flag = 0 + + def make_key(self): + self.key = ".".join((self.stem, self.station, self.sequence)) + + + +@dataclass +class ArtdayIndex: + def __init__(self, fname): + self.artday = [] + self.msktime = [] + self.complete = [] + self.status = [] + self.errors = [] + self.processed = [] + self.toolver = [] + self.sciproc = [] + self.scitool = [] + self.scicaldb = [] + self.start = [] + self.stop = [] + self.daylen = [] + self.t1 = [] + self.t2 = [] + self.t3 = [] + self.t4 = [] + self.t5 = [] + self.t6 = [] + self.t7 = [] + self.gyro = [] + self.bokz = [] + self.sed1 = [] + self.sed2 = [] + + self._index = fname + + if os.path.isfile(self._index): + self.read(self._index) + + def read(self, fname): + with fits.open(fname) as hdul: + self.artday = hdul[1].data['ARTDAY' ].tolist() + self.msktime = hdul[1].data['MSKTIME' ].tolist() + self.complete = hdul[1].data['COMPLETE' ].tolist() + self.status = hdul[1].data['STATUS' ].tolist() + self.errors = hdul[1].data['ERRORS' ].tolist() + self.processed = hdul[1].data['PROCESSED' ].tolist() + self.toolver = hdul[1].data['TOOLVER' ].tolist() + self.sciproc = hdul[1].data['SCI_PROCESSED'].tolist() + self.scitool = hdul[1].data['SCI_TOOLVER' ].tolist() + self.scicaldb = hdul[1].data['SCI_CALDB' ].tolist() + self.start = hdul[1].data['START' ].tolist() + self.stop = hdul[1].data['STOP' ].tolist() + self.daylen = hdul[1].data['DAYLEN' ].tolist() + self.t1 = hdul[1].data['T1' ].tolist() + self.t2 = hdul[1].data['T2' ].tolist() + self.t3 = hdul[1].data['T3' ].tolist() + self.t4 = hdul[1].data['T4' ].tolist() + self.t5 = hdul[1].data['T5' ].tolist() + self.t6 = hdul[1].data['T6' ].tolist() + self.t7 = hdul[1].data['T7' ].tolist() + self.gyro = hdul[1].data['GYRO' ].tolist() + self.bokz = hdul[1].data['BOKZ' ].tolist() + self.sed1 = hdul[1].data['SED1' ].tolist() + self.sed2 = hdul[1].data['SED2' ].tolist() + + def process(self, frame): + + if frame['artday'] not in self.artday: + self.artday.append(frame['artday']) + self.msktime.append(frame['msktime']) + self.complete.append(frame['complete']) + self.status.append(-100) + self.errors.append(frame['errors']) + self.processed.append(frame['processed']) + self.toolver.append(frame['artfitstool']) + self.sciproc.append(' ') + self.scitool.append(' ') + self.scicaldb.append(' ') + self.start.append(frame['start']) + self.stop.append(frame['stop']) + self.daylen.append(frame['daylen']) + self.t1.append(frame['T1']) + self.t2.append(frame['T2']) + self.t3.append(frame['T3']) + self.t4.append(frame['T4']) + self.t5.append(frame['T5']) + self.t6.append(frame['T6']) + self.t7.append(frame['T7']) + self.gyro.append(frame['GYRO']) + self.bokz.append(frame['BOKZ']) + self.sed1.append(frame['SED1']) + self.sed2.append(frame['SED2']) + else: + index = np.where(np.array(self.artday) == frame['artday']) + if not index: + print('warning: failed to update index for artday={}'.format(frame['artday'])) + return + + index = index[0][0] + + self.complete[index] = frame['complete'] + self.errors[index] = frame['errors'] + self.processed[index] = frame['processed'] + self.toolver[index] = frame['artfitstool'] + self.start[index] = frame['start'] + self.stop[index] = frame['stop'] + self.daylen[index] = frame['daylen'] + self.t1[index] = frame['T1'] + self.t2[index] = frame['T2'] + self.t3[index] = frame['T3'] + self.t4[index] = frame['T4'] + self.t5[index] = frame['T5'] + self.t6[index] = frame['T6'] + self.t7[index] = frame['T7'] + self.gyro[index] = frame['GYRO'] + self.bokz[index] = frame['BOKZ'] + self.sed1[index] = frame['SED1'] + self.sed2[index] = frame['SED2'] + + + def index(self): + return pl.DataFrame( + { + 'artday' : self.artday , + 'msktime' : self.msktime , + 'complete' : self.complete , + 'status' : self.status , + 'errors' : self.errors , + 'processed' : self.processed, + 'toolver' : self.toolver , + 'sciproc' : self.sciproc , + 'scitool' : self.scitool , + 'scicaldb' : self.scicaldb , + 'start' : self.start , + 'stop' : self.stop , + 'daylen' : self.daylen , + 'T1' : self.t1 , + 'T2' : self.t2 , + 'T3' : self.t3 , + 'T4' : self.t4 , + 'T5' : self.t5 , + 'T6' : self.t6 , + 'T7' : self.t7 , + 'GYRO' : self.gyro , + 'BOKZ' : self.bokz , + 'SED1' : self.sed1 , + 'SED2' : self.sed2 , + } + , schema={ + 'artday' : pl.String , + 'msktime' : pl.String , + 'complete' : pl.Int16 , + 'status' : pl.Int16 , + 'errors' : pl.Int16 , + 'processed' : pl.String , + 'toolver' : pl.String , + 'sciproc' : pl.String , + 'scitool' : pl.String , + 'scicaldb' : pl.String , + 'start' : pl.Float64, + 'stop' : pl.Float64, + 'daylen' : pl.Float64, + 'T1' : pl.Int16 , + 'T2' : pl.Int16 , + 'T3' : pl.Int16 , + 'T4' : pl.Int16 , + 'T5' : pl.Int16 , + 'T6' : pl.Int16 , + 'T7' : pl.Int16 , + 'GYRO' : pl.Int16 , + 'BOKZ' : pl.Int16 , + 'SED1' : pl.Int16 , + 'SED2' : pl.Int16 + } + ).sort('artday') + + def _calculate_status(self, df): + status = np.array(df['status']) + errors = np.array(df['errors']) + artdays = np.array(df['artday'], dtype=int) + + dtime = Dtime(datetime.utcnow().isoformat(), tzone=TZ_UTC) + nowday = ArtDay(dtime.to_artday()).get() + nowdays = np.full_like(artdays, nowday) + + deltadays = nowdays - artdays + + mask_maychange = status != 1 + mask_haserrors = errors != 0 + mask_freshdays = deltadays < 4 + + newstatus = np.zeros_like(artdays) + newstatus[mask_maychange & mask_freshdays & ~mask_haserrors] = -2 + newstatus[mask_maychange & mask_freshdays & mask_haserrors] = -5 + newstatus[mask_maychange & ~mask_freshdays & mask_haserrors] = -1 + newstatus[~mask_maychange] = 1 + + return newstatus + + def write_index_fits(self): + df = self.index() + status = self._calculate_status(df) + + # print(df) + + hdu = fits.PrimaryHDU() + + # status = np.full_like(df['complete' ], -1) + + c00 = fits.Column(name='ARTDAY' , format='6A' , unit='' , array=np.array(df['artday' ])) + c01 = fits.Column(name='MSKTIME' , format='19A', unit='date', array=np.array(df['msktime' ])) + c02 = fits.Column(name='COMPLETE' , format='1I' , unit='' , array=np.array(df['complete' ])) + c03 = fits.Column(name='STATUS' , format='1I' , unit='' , array=np.array(status )) + c04 = fits.Column(name='ERRORS' , format='1I' , unit='' , array=np.array(df['errors' ])) + c05 = fits.Column(name='PROCESSED' , format='19A', unit='date', array=np.array(df['processed'])) + c06 = fits.Column(name='TOOLVER' , format='10A', unit='date', array=np.array(df['toolver' ])) + c07 = fits.Column(name='SCI_PROCESSED', format='19A', unit='' , array=np.array(df['sciproc' ])) + c08 = fits.Column(name='SCI_TOOLVER' , format='10A', unit='' , array=np.array(df['scitool' ])) + c09 = fits.Column(name='SCI_CALDB' , format='10A', unit='' , array=np.array(df['scicaldb' ])) + c10 = fits.Column(name='START' , format='1D' , unit='sec' , array=np.array(df['start' ])) + c11 = fits.Column(name='STOP' , format='1D' , unit='sec' , array=np.array(df['stop' ])) + c12 = fits.Column(name='DAYLEN' , format='1D' , unit='sec' , array=np.array(df['daylen' ])) + c13 = fits.Column(name='T1' , format='1I' , unit='' , array=np.array(df['T1' ])) + c14 = fits.Column(name='T2' , format='1I' , unit='' , array=np.array(df['T2' ])) + c15 = fits.Column(name='T3' , format='1I' , unit='' , array=np.array(df['T3' ])) + c16 = fits.Column(name='T4' , format='1I' , unit='' , array=np.array(df['T4' ])) + c17 = fits.Column(name='T5' , format='1I' , unit='' , array=np.array(df['T5' ])) + c18 = fits.Column(name='T6' , format='1I' , unit='' , array=np.array(df['T6' ])) + c19 = fits.Column(name='T7' , format='1I' , unit='' , array=np.array(df['T7' ])) + c20 = fits.Column(name='GYRO' , format='1I' , unit='' , array=np.array(df['GYRO' ])) + c21 = fits.Column(name='BOKZ' , format='1I' , unit='' , array=np.array(df['BOKZ' ])) + c22 = fits.Column(name='SED1' , format='1I' , unit='' , array=np.array(df['SED1' ])) + c23 = fits.Column(name='SED2' , format='1I' , unit='' , array=np.array(df['SED2' ])) + + cols = fits.ColDefs( + [ + c00, c01, c02, c03, c04, c05, c06, c07, c08, c09, + c10, c11, c12, c13, c14, c15, c16, c17, c18, c19, + c20, c21, c22, c23 + ] + ) + + indexhdu = fits.BinTableHDU.from_columns(cols) + indexhdu.name = "INDEX" + + creator = 'artindex v.{}'.format(__version__) + indexhdu.header['creator' ] = (creator, 'program and version that created this file') + + + dtnow = datetime.now(tz=TZ_MSK) + + indexhdu.header['mjdref' ] = (MJDREFD , 'MJD corresponding to SC clock start 2000.0 MSK' ) + indexhdu.header['date' ] = (dtnow.strftime('%Y/%m/%d'), 'date that FITS file was created' ) + + + hdulist = fits.HDUList([hdu, indexhdu]) + hdulist.writeto(self._index, overwrite=True) + + print("written: {}".format(self._index)) + + +@dataclass +class SessionIndexFrame: + def __init__(self): + self.key = "" + + self.stem = "" + self.station = "00" + self.sequence = "000" + self.data_quality = 0 + self.total_quality = 0 + + self.processed = "" + + self.intervals = { + "T1": (0, 0), + "T2": (0, 0), + "T3": (0, 0), + "T4": (0, 0), + "T5": (0, 0), + "T6": (0, 0), + "T7": (0, 0), + } + + self.flag = 0 + + def make_key(self): + self.key = ".".join((self.stem, self.station, self.sequence)) + + +if __name__ == "__main__": + pass diff --git a/arttools/src/arttools/arttime.py b/arttools/src/arttools/arttime.py new file mode 100755 index 0000000..8e58f74 --- /dev/null +++ b/arttools/src/arttools/arttime.py @@ -0,0 +1,388 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2021-2024 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the M. Pavlinksy SRG/ART-XC telescope software. +# ----------------------------------------------------------------------------- + +import os +import sys +import shutil +import argparse + +import numpy as np + +from datetime import datetime + +import astropy.units as u +from astropy.io import fits +from astropy.time.formats import erfa +from astropy.time import Time +from astropy.time import TimezoneInfo + +from arttools.version import __version__ + +from arttools.argparse import KvAction +from arttools.toolargs import ToolArgs + + +MJDREF = Time(51543.875, format='mjd') #corresponds to date 01.01.2000 00:00:00 (UTC+3) +TZ_UTC = TimezoneInfo(utc_offset=0*u.hour) #UTC time zone +TZ_MSK = TimezoneInfo(utc_offset=3*u.hour) #UTC+3 (Moscow) time zone + + +class ArtTime(object): + + def __init__(self, artday=None, missiontime=None, mjd=None, datetime_utc=None, datetime_msk=None): + + self._has_value = 0 + self._try_aday(artday) + self._try_mtime(missiontime) + self._try_mjd(mjd) + self._try_dtime_utc(datetime_utc) + self._try_dtime_msk(datetime_msk) + + if self._has_value > 1: + print('Error: ArtTime needs only one argument') + exit() + + if self._has_value == 0: + datetime_utc = datetime.utcnow() + self._try_dtime_utc(datetime_utc.isoformat()) + + def _try_mtime(self, missiontime): + if missiontime is None: + return + + self._has_value = self._has_value+1 + + self._missiontime = MissionTime(missiontime) + self._artday = ArtDay(self._missiontime.to_artday()) + self._mjd = MjdDay(self._missiontime.to_mjd()) + self._datetime_utc = Dtime(dtime=str(self._missiontime.to_datetime(tzone=TZ_UTC))) + self._datetime_msk = Dtime(dtime=str(self._missiontime.to_datetime(tzone=TZ_MSK))) + return + + def _try_aday(self, artday): + if artday is None: + return + + self._has_value = self._has_value+1 + + self._artday = ArtDay(artday) + self._missiontime = MissionTime(self._artday.to_mission()) + self._mjd = MjdDay(self._artday.to_mjd()) + self._datetime_utc = Dtime(dtime=str(self._artday.to_datetime(tzone=TZ_UTC))) + self._datetime_msk = Dtime(dtime=str(self._artday.to_datetime(tzone=TZ_MSK))) + + def _try_mjd(self, mjd): + if mjd is None: + return + + self._has_value = self._has_value+1 + + self._mjd = MjdDay(mjd) + self._artday = ArtDay(self._mjd.to_artday()) + self._missiontime = MissionTime(self._mjd.to_mission()) + self._datetime_utc = Dtime(dtime=str(self._mjd.to_datetime(tzone=TZ_UTC))) + self._datetime_msk = Dtime(dtime=str(self._mjd.to_datetime(tzone=TZ_MSK))) + + def _try_dtime_utc(self, datetime_utc): + if datetime_utc is None: + return + + self._has_value = self._has_value+1 + + self._datetime_utc = Dtime(datetime_utc, tzone=TZ_UTC) + self._datetime_msk = Dtime(str(self._datetime_utc.to_msk())) + self._missiontime = MissionTime(self._datetime_utc.to_mission()) + self._artday = ArtDay(self._datetime_utc.to_artday()) + self._mjd = MjdDay(self._datetime_utc.to_mjd()) + + def _try_dtime_msk(self, datetime_msk): + if datetime_msk is None: + return + + self._has_value = self._has_value+1 + + self._datetime_msk = Dtime(datetime_msk, tzone=TZ_MSK) + self._datetime_utc = Dtime(str(self._datetime_msk.to_utc())) + self._missiontime = MissionTime(self._datetime_msk.to_mission()) + self._artday = ArtDay(self._datetime_msk.to_artday()) + self._mjd = MjdDay(self._datetime_msk.to_mjd()) + + @property + def missiontime(self): + return self._missiontime + + @property + def artday(self): + return self._artday + + @property + def mjd(self): + return self._mjd + + @property + def datetime_utc(self): + return self._datetime_utc + + @property + def datetime_msk(self): + return self._datetime_msk + + def __str__(self): + aday_str = 'Artday --> ' + str(self._artday ) + mtime_str = 'Mission time --> ' + str(self._missiontime ) + mjd_str = 'MJD --> ' + str(self._mjd ) + date_msk_str = 'Date MSK --> ' + str(self._datetime_msk) + date_utc_str = 'Date UTC --> ' + str(self._datetime_utc) + + return '\n'.join((aday_str, mtime_str, mjd_str, date_msk_str, date_utc_str)) + + +class MissionTime(object): + + def __init__(self, missiontime): + self._missiontime = missiontime + + def get(self): + return self._missiontime + + @property + def time(self): + return self._missiontime + + def __str__(self): + return str(self._missiontime) + + def to_datetime(self, tzone=TZ_UTC): + days = self.to_mjd() + return days.to_datetime(timezone=tzone) + + def to_artday(self): + ARTDAYLEN = erfa.DAYSEC + return self._missiontime / ARTDAYLEN + + def to_mjd(self): + mjd_value = MJDREF.mjd + Time(self._missiontime / erfa.DAYSEC, format='mjd').mjd + return Time(mjd_value, format='mjd') + + +class ArtDay(object): + + def __init__(self, artday): + self._artday = artday + + def get(self): + return self._artday + + @property + def time(self): + return self._artday + + def __str__(self): + return str(self._artday) + + def to_mission(self): + ARTDAYLEN = erfa.DAYSEC + return self._artday * ARTDAYLEN + + def to_datetime(self, tzone=TZ_UTC): + days = self.to_mjd() + return days.to_datetime(timezone=tzone) + + def to_mjd(self): + mjd_value = MJDREF.mjd + Time(self._artday, format='mjd').mjd + return Time(mjd_value, format='mjd') + + def to_str(self): + return '{:0>6}'.format(self._artday) + + +class MjdDay(object): + + def __init__(self, mjd): + self._mjd = Time(mjd, format='mjd') + + def get(self): + return self._mjd + + @property + def time(self): + return self._mjd + + def __str__(self): + return str(self._mjd) + + def to_mission(self): + ARTDAYLEN = erfa.DAYSEC + return (self._mjd.mjd - MJDREF.mjd) * ARTDAYLEN + + def to_artday(self): + return self._mjd.mjd - MJDREF.mjd + + def to_datetime(self, tzone=TZ_UTC): + return Time(self._mjd, format='mjd').to_datetime(timezone=tzone) + + +class Dtime(object): + + def __init__(self, dtime, tzone=None): + dtime = datetime.fromisoformat(dtime) + if tzone is None: + self._dtime = dtime + else: + self._dtime = dtime.replace(tzinfo=tzone) + self._tzone = tzone + + def get(self): + return self._dtime + + @property + def time(self): + return self._dtime + + def __str__(self): + return str(self._dtime) + + def to_mjd(self): + days = Time(self._dtime) + return days + + def to_artday(self): + return self.to_mjd().mjd - MJDREF.mjd + + def to_mission(self): + ARTDAYLEN = erfa.DAYSEC + return self.to_artday() * ARTDAYLEN + + def to_msk(self): + dtime_utc = Time(self._dtime) + return dtime_utc.to_datetime(timezone=TZ_MSK) + + def to_utc(self): + dtime_msk = Time(self._dtime) + return dtime_msk.to_datetime(timezone=TZ_UTC) + + +class Args(ToolArgs): + + def __init__(self): + super().__init__() + + self.artday = None + self.missiontime = None + self.mjd = None + self.datetime_utc = None + self.datetime_msk = None + + def read(self): + argparser = argparse.ArgumentParser( + description=('The arttime tool converts artdays, mission time, MSK and UTC date time' + ' and MJD to each other. Calling arttime without parameters will use the current ' + 'local time as input.')) + argparser.add_argument( + '--artday', '--aday', + help=('ART-XC days.' + 'Example: arttime --artday=8000'), + action=KvAction) + argparser.add_argument( + '--missiontime', '--mtime', + help=('Onboard ART-XC seconds. ' + 'Example: arttime --missiontime=691200000'), + action=KvAction) + argparser.add_argument( + '--mjd', + help=('MJD. ' + 'Example: arttime --mjd=59543.875'), + action=KvAction) + argparser.add_argument( + '--datetime_utc', '--dtime_utc', + help=('Date time based on UTC time zone. ' + 'Example: arttime --datetime_utc=2021-11-25T21:00:00'), + action=KvAction) + argparser.add_argument( + '--datetime_msk', '--dtime_msk', + help=('Date time based on Moscow time zone (UTC+3). ' + 'Example: arttime --datetime_msk=2021-11-26T00:00:00'), + action=KvAction) + args = argparser.parse_args() + + if args.artday is not None: + self.artday = float(args.artday) + if args.missiontime is not None: + self.missiontime = float(args.missiontime) + if args.mjd is not None: + self.mjd = float(args.mjd) + self.datetime_utc = args.datetime_utc + self.datetime_msk = args.datetime_msk + + def print(self): + logstr = '\n'.join(( + 'Input parameters:', + ' artday = {}'.format(self.artday ), + ' missiontime = {}'.format(self.missiontime ), + ' mjd = {}'.format(self.mjd ), + ' datetime_utc = {}'.format(self.datetime_utc), + ' datetime_msk = {}'.format(self.datetime_msk) + )) + print(logstr) + + +class Task(object): + + def __init__(self): + pass + + def __enter__(self): + self.initialize() + return self + + def __exit__(self, ex_type, ex_value, ex_traceback): + self.finalize() + + def initialize(self): + print('running: arttime v.{}'.format(__version__) ) + print('================================================================================') + + self._args = Args() + self._args.read() + + self._args.print() + + def finalize(self): + print('--------------------------------------------------------------------------------') + + def run(self): + if self._args.artday is not None: + print(ArtTime(artday=self._args.artday)) + + if self._args.missiontime is not None: + print(ArtTime(missiontime=self._args.missiontime)) + + if self._args.mjd is not None: + print(ArtTime(mjd=self._args.mjd)) + + if self._args.datetime_utc is not None: + print(ArtTime(datetime_utc=self._args.datetime_utc)) + + if self._args.datetime_msk is not None: + print(ArtTime(datetime_msk=self._args.datetime_msk)) + + if len(sys.argv) <= 1: + datetime_utc = datetime.utcnow() + print('*** Convertation for the current local time:') + print(ArtTime(datetime_utc=datetime_utc.isoformat())) + + +def main(): + with Task() as task: + task.run() + + +if __name__ == '__main__': + main() + diff --git a/arttools/src/arttools/cmd.py b/arttools/src/arttools/cmd.py new file mode 100755 index 0000000..71734c8 --- /dev/null +++ b/arttools/src/arttools/cmd.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2025 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the SRG/ART-XC software project. +# ----------------------------------------------------------------------------- + +import sys +import asyncio as aio + +from arttools.errors import Error + + +class InternalError(Error): + pass + + +class UncleanExecutionError(Error): + pass + + +class ExecutionResult(object): + + def __init__(self, cmd, retcode, stdout_data, stderr_data): + self.command = cmd + self.return_code = retcode + self.stdout_data = stdout_data + self.stderr_data = stderr_data + + +class Command(object): + + def __init__(self, cmd, stdout=sys.stdout, stderr=sys.stderr): + self._cmd = cmd + self._stdout = stdout + self._stderr = stderr + + self._result = None + + async def _run(self): + process = await aio.create_subprocess_shell( + cmd=self._cmd, + stdout=self._stdout, + stderr=self._stderr + ) + + stdout_data, stderr_data = await process.communicate() + + if stdout_data is not None: + stdout_data = stdout_data.decode() + + if stderr_data is not None: + stderr_data = stderr_data.decode() + + self._result = \ + ExecutionResult( + self._cmd, + process.returncode, + stdout_data, + stderr_data + ) + + def run(self): + aio.run(self._run()) + + return self._result + +if __name__ == "__main__": + pass + diff --git a/arttools/src/arttools/datatable.py b/arttools/src/arttools/datatable.py new file mode 100755 index 0000000..8ec6fb1 --- /dev/null +++ b/arttools/src/arttools/datatable.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2025 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 polars as pl + +import numpy as np +from astropy.io import fits + + +class DataTable(object): + + def __init__(self): + self._dataframe = None + self._adapter = None + + @property + def frame(self): return self._dataframe + + @property + def size(self): return self._dataframe.height + + def adapter(self, adapter): + self._adapter = adapter + + return self + + def read(self): + if self._adapter is not None: + self._dataframe = self._adapter.read() + + return self + + def column(self, name): + return self._dataframe.get_column(name) + + def setcol(self, name, arr): + self._dataframe =\ + self._dataframe.with_columns( + pl.Series(name, arr) + ) + + +class FitsAdapter(object): + + def __init__(self, filename, tablename): + self._FORMATS = { + 'D': pl.Float64, + 'J': pl.UInt32 , + 'I': pl.UInt16 , + 'B': pl.UInt8 , + 'A': pl.String + } + + self._filename = filename + self._tablename = tablename + + def _schema_for(self, columns): + schema = [] + for col in columns: + colformat = ''.join([ch for ch in col.format if ch.isalpha()]) + schema.append((col.name, self._FORMATS[colformat])) + + return schema + + def read(self): + with fits.open(self._filename) as hdul: + hdu = hdul[self._tablename] + + schema = self._schema_for(hdu.data.columns) + + data = {} + for col in hdu.data.columns: + colformat = ''.join([ch for ch in col.format if ch.isalpha()]) + + if colformat == 'D': data[col.name] = np.array(hdu.data[col.name], dtype=np.float64) + elif colformat == 'I': data[col.name] = np.array(hdu.data[col.name], dtype=np.uint16 ) + else: + data[col.name] = hdu.data[col.name] + + return pl.DataFrame(data, schema=schema) + + def write(self): + pass + + +class Adapter(object): + pass + + +if __name__ == '__main__': + pass + + diff --git a/arttools/src/arttools/errors.py b/arttools/src/arttools/errors.py new file mode 100755 index 0000000..49060d0 --- /dev/null +++ b/arttools/src/arttools/errors.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2025 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the MVN software project. +# ----------------------------------------------------------------------------- + +class Error(Exception): + pass + +class ConfigNotFoundError(Error): + pass + +class RegionsError(Error): + pass + +if __name__ == "__main__": + pass + diff --git a/arttools/src/arttools/events.py b/arttools/src/arttools/events.py new file mode 100755 index 0000000..83e584b --- /dev/null +++ b/arttools/src/arttools/events.py @@ -0,0 +1,404 @@ +#!/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 shutil + +import datetime + +import numpy as np +from scipy.interpolate import interp1d + +from astropy.io import fits + +import astropy.units as u +from astropy.time.formats import erfa +from astropy.time import Time +from astropy.time import TimezoneInfo + +TZ_UTC = TimezoneInfo(utc_offset=0*u.hour) # UTC time zone +TZ_MSK = TimezoneInfo(utc_offset=3*u.hour) # UTC+3 (Moscow) time zone +from arttools.arttime import MissionTime + +from arttools.gti import Gti +from arttools.datatable import DataTable, FitsAdapter + +from arttools.version import __version__ + + +class Events(object): + _deadtimes = { + 1: 770e-6, + 2: 770e-6, + 3: 770e-6, + 4: 770e-6, + 5: 770e-6, + 6: 1060e-6, + 7: 770e-6, + } + + _eventlist_keys = [ + 'TIME', + 'TIME_I', + 'TIME_F', + 'TIME_CORR', + 'ENERGY', + 'ENERGY_TOP', + 'ENERGY_BOT', + 'GRADE', + 'FLAG', + 'RA', + 'DEC', + 'RAW_X', + 'RAW_Y', + 'PHA_BOT', + 'PHA_BOT_ADD1', + 'PHA_BOT_SUB1', + 'PHA_TOP', + 'PHA_TOP_ADD1', + 'PHA_TOP_SUB1', + 'TRIGGER', + ] + + _eventlist_layout = { + 'TIME' : {'format': '1D', 'unit': 'sec' }, + 'TIME_I' : {'format': '1J', 'unit': 'sec', 'bzero': 2147483648}, + 'TIME_F' : {'format': '1D', 'unit': 'sec' }, + 'TIME_CORR' : {'format': '1D', 'unit': 'sec' }, + 'ENERGY' : {'format': '1D', 'unit': 'keV' }, + 'ENERGY_TOP' : {'format': '1D', 'unit': 'keV' }, + 'ENERGY_BOT' : {'format': '1D', 'unit': 'keV' }, + 'GRADE' : {'format': '1I', 'unit': '' }, + 'FLAG' : {'format': '1I', 'unit': '' }, + 'RA' : {'format': '1D', 'unit': 'deg' }, + 'DEC' : {'format': '1D', 'unit': 'deg' }, + 'RAW_X' : {'format': '1I', 'unit': 'pix' }, + 'RAW_Y' : {'format': '1I', 'unit': 'pix' }, + 'PHA_BOT' : {'format': '1I', 'unit': '' }, + 'PHA_BOT_ADD1': {'format': '1I', 'unit': '' }, + 'PHA_BOT_SUB1': {'format': '1I', 'unit': '' }, + 'PHA_TOP' : {'format': '1I', 'unit': '' }, + 'PHA_TOP_ADD1': {'format': '1I', 'unit': '' }, + 'PHA_TOP_SUB1': {'format': '1I', 'unit': '' }, + 'TRIGGER' : {'format': '1I', 'unit': '' }, + } + + def __init__(self, srcfile): + super().__init__() + + self._srcfile = srcfile + + self._telname = 0 + self._teln = 0 + + self._tstart = 0 + self._tstop = 0 + + self._events = DataTable() + self._hk = DataTable() + self._gti = Gti() + + @property + def empty(self): + return self.size == 0 + + @property + def size(self): + return self._events.size + + @property + def teln(self): + return self._teln + + @property + def telname(self): + return self._telname + + @property + def mjdref(self): + return self._mjdref + + @property + def tstart(self): + return self._tstart + + @property + def tstop(self): + return self._tstop + + # @property + # def deadtime(self): + # return self._deadtime + + @property + def events(self): + return self._events + + @property + def hk(self): + return self._hk + + @property + def gti(self): + return self._gti + + @property + def filename(self): + return self._srcfile + + def _telescope_deadtime(self): + return self._deadtimes[self._teln] + + def deadtime(self, start, stop): + interval = stop - start + + hktimes = np.array(self.hk.column('TIME' )) + hkevt = np.array(self.hk.column('EVENTS')).astype(np.int64) + + X = hktimes[1:] + Y = np.diff(hkevt) / np.diff(hktimes) + hkmask = Y > 0 + + X = X[hkmask] + Y = Y[hkmask] + + datamask = (X >= start) & (X <= stop) + interp = interp1d(X, Y, kind='linear', bounds_error=False, fill_value='extrapolate') + + begin = interp(start).item() + end = interp(stop).item() + middle = Y[datamask] + + cnttimes = np.array([start, *X[datamask], stop]) + countrate = np.array([begin, *middle, end]) + + return np.average(((countrate[1:] + countrate[:-1]) / 2.), weights=np.diff(cnttimes)) * interval * self._telescope_deadtime() + + def bindeadtime(self, edges, bins): + hktimes = np.array(self.hk.column('TIME' )) + hkevt = np.array(self.hk.column('EVENTS')).astype(np.int64) + + binrate = np.ones_like(bins) + + evtcountfn = interp1d(hktimes, hkevt, + kind='linear', bounds_error=False, fill_value='extrapolate') + + edge_events = evtcountfn(edges) + + # 2) make mask for event counter reversals (negative diffs) + edgediff = np.diff(edge_events) + edgemask = edgediff > 0 + + # # 3) interpolate filtered (2) points + edgecountfn = interp1d( + edges[1:][edgemask], edgediff[edgemask], + kind='linear', bounds_error=False, fill_value='extrapolate') + + # 4) calculate deadtimes and fractional exposure for bins + bindeadt = edgecountfn(edges) * self._telescope_deadtime() + binrate -= bindeadt[1:] / np.diff(edges) + + return bindeadt[1:], binrate + + def _read_header(self): + with fits.open(self._srcfile) as hdul: + self._teln = hdul['events'].header['teln' ] + self._telname = hdul['events'].header['instrume'] + self._mjdref = hdul['events'].header['mjdref' ] + self._tstart = hdul['events'].header['tstart' ] + self._tstop = hdul['events'].header['tstop' ] + + def read(self): + self._read_header() + + self._events = DataTable().adapter(FitsAdapter(self._srcfile, 'events')).read() + self._hk = DataTable().adapter(FitsAdapter(self._srcfile, 'hk' )).read() + self._gti = Gti.from_hduname(self._srcfile, 'stdgti') + + return self + + def filter(self, gti): + pass + + def write_cl(self, dstfile, cdbinfo=[]): + + srchdul = fits.open(self._srcfile) + + # --- primary HDU --- + hdu = fits.PrimaryHDU() + + # --- EVENTS HDU --- + a00 = np.array(self._events.column('TIME' )) + a01 = np.array(self._events.column('TIME_I' )) + a02 = np.array(self._events.column('TIME_F' )) + a03 = np.array(self._events.column('TIME_CORR' )) + a04 = np.array(self._events.column('ENERGY' )) + a05 = np.array(self._events.column('ENERGY_TOP' )) + a06 = np.array(self._events.column('ENERGY_BOT' )) + a07 = np.array(self._events.column('GRADE' )) + a08 = np.array(self._events.column('FLAG' )) + a09 = np.array(self._events.column('RA' )) + a10 = np.array(self._events.column('DEC' )) + a11 = np.array(self._events.column('RAW_X' )) + a12 = np.array(self._events.column('RAW_Y' )) + a13 = np.array(self._events.column('PHA_BOT' )) + a14 = np.array(self._events.column('PHA_BOT_ADD1')) + a15 = np.array(self._events.column('PHA_BOT_SUB1')) + a16 = np.array(self._events.column('PHA_TOP' )) + a17 = np.array(self._events.column('PHA_TOP_ADD1')) + a18 = np.array(self._events.column('PHA_TOP_SUB1')) + a19 = np.array(self._events.column('TRIGGER' )) + + c00 = fits.Column(name='TIME' , format='1D', unit='sec', array=a00 ) + c01 = fits.Column(name='TIME_I' , format='1J', unit='sec', array=a01, bzero=2147483648) + c02 = fits.Column(name='TIME_F' , format='1D', unit='sec', array=a02 ) + c03 = fits.Column(name='TIME_CORR' , format='1D', unit='sec', array=a03 ) + c04 = fits.Column(name='ENERGY' , format='1D', unit='keV', array=a04 ) + c05 = fits.Column(name='ENERGY_TOP' , format='1D', unit='keV', array=a05 ) + c06 = fits.Column(name='ENERGY_BOT' , format='1D', unit='keV', array=a06 ) + c07 = fits.Column(name='GRADE' , format='1I', unit='' , array=a07 ) + c08 = fits.Column(name='FLAG' , format='1I', unit='' , array=a08 ) + c09 = fits.Column(name='RA' , format='1D', unit='deg', array=a09 ) + c10 = fits.Column(name='DEC' , format='1D', unit='deg', array=a10 ) + c11 = fits.Column(name='RAW_X' , format='1I', unit='pix', array=a11 ) + c12 = fits.Column(name='RAW_Y' , format='1I', unit='pix', array=a12 ) + c13 = fits.Column(name='PHA_BOT' , format='1I', unit='' , array=a13 ) + c14 = fits.Column(name='PHA_BOT_ADD1', format='1I', unit='' , array=a14 ) + c15 = fits.Column(name='PHA_BOT_SUB1', format='1I', unit='' , array=a15 ) + c16 = fits.Column(name='PHA_TOP' , format='1I', unit='' , array=a16 ) + c17 = fits.Column(name='PHA_TOP_ADD1', format='1I', unit='' , array=a17 ) + c18 = fits.Column(name='PHA_TOP_SUB1', format='1I', unit='' , array=a18 ) + c19 = fits.Column(name='TRIGGER' , format='1I', unit='' , array=a19 ) + + cols = fits.ColDefs([ + c00, c01, c02, c03, c04, c05, c06, c07, c08, c09, + c10, c11, c12, c13, c14, c15, c16, c17, c18, c19 + ]) + + evthdu = fits.BinTableHDU.from_columns(cols) + evthdu.name = 'EVENTS' + + creator = 'artpipeline v.{}'.format(__version__) + + evthdu.header['creator' ] = (creator , 'program and version that created this file' ) + evthdu.header['origin' ] = ('IKI RAS' , 'institution that created this file' ) + evthdu.header['telescop'] = ('SRG/ART-XC', 'mission name' ) + + # copy header keys + COPYKEYS = [ + 'instrume', 'teln', 'detn', 'urdn', + 'mjdref', 'artday', 'tstart', 'tstop', + 'date', 'date-obs', 'time-obs', 'date-end', 'time-end', + 'maxticks', 'timedel', 'timeunit', 'timeref', 'tassign' + ] + + for key in COPYKEYS: + try: + evthdu.header[key] = srchdul['EVENTS'].header[key] + except KeyError: + pass + + + # add caldb key + for key in cdbinfo: + evthdu.header[key] = cdbinfo[key] + + # --- make file --- + hdulist = fits.HDUList([hdu, evthdu, srchdul['KVEA'], srchdul['HK'], srchdul['STDGTI']]) + + hdulist.writeto(dstfile, overwrite=True, checksum=True) + + def modify_and_write_cl(self, dstfile, creator, keys, layout, data, position, mask, gti, cdbinfo=[]): + + srchdul = fits.open(self._srcfile) + mjdref = srchdul['EVENTS'].header['MJDREF'] + + # -- prepare keys + newkeys = self._eventlist_keys[:position] + keys + self._eventlist_keys[position:] + + # -- prepare layout + newlayout = self._eventlist_layout + newlayout.update(layout) + + # -- prepare data + newdata = { } + for key in self._eventlist_keys: + newdata[key] = np.array(self._events.column(key)) + + for key in keys: + newdata[key] = np.array(data[key]) + + # -- prepare fits structure + hdu = fits.PrimaryHDU() + + cols = [] + for key in newkeys: + lay = newlayout[key] + + if 'bzero' in lay: + col = fits.Column(name=key, format=lay['format'], unit=lay['unit'], array=newdata[key][mask], bzero=lay['bzero']) + else: + col = fits.Column(name=key, format=lay['format'], unit=lay['unit'], array=newdata[key][mask]) + + cols.append(col) + + evthdu = fits.BinTableHDU.from_columns(cols) + evthdu.name = 'EVENTS' + + evthdu.header['creator' ] = (creator , 'program and version that created this file' ) + evthdu.header['origin' ] = ('IKI RAS' , 'institution that created this file' ) + evthdu.header['telescop'] = ('SRG/ART-XC', 'mission name' ) + + # copy header keys + COPYKEYS = [ + 'instrume', 'teln', 'detn', 'urdn', + 'maxticks', 'timedel', 'timeunit', 'timeref', 'tassign' + ] + + for key in COPYKEYS: + try: + evthdu.header[key] = srchdul['EVENTS'].header[key] + except KeyError: + pass + + dtnow = datetime.datetime.now(tz=TZ_MSK) + dtobs = MissionTime(self.gti.start).to_datetime(TZ_MSK) + dtend = MissionTime(self.gti.stop ).to_datetime(TZ_MSK) + + deadtime = self.deadtime(self.gti.start, self.gti.stop) + + evthdu.header['MJDREF' ] = (mjdref , 'MJD corresponding to SC clock start 2000.0 MSK') + + evthdu.header['tstart' ] = (self.gti.start , 'start time of the observation [sec]' ) + evthdu.header['tstop' ] = (self.gti.stop , 'end time of the observation [sec]' ) + + evthdu.header['date' ] = (dtnow.strftime('%Y/%m/%d') , 'date that FITS file was created' ) + evthdu.header['date-obs'] = (dtobs.strftime('%Y/%m/%d') , 'start date of the observation MSK') + evthdu.header['date-end'] = (dtend.strftime('%Y/%m/%d') , 'end date of the observation MSK' ) + evthdu.header['time-obs'] = (dtobs.strftime('%H:%M:%S') , 'start time of the observation MSK') + evthdu.header['time-end'] = (dtend.strftime('%H:%M:%S') , 'end time of the observation MSK' ) + + evthdu.header['onsource'] = (self.gti.onsource , '[sec]' ) + evthdu.header['exposure'] = (self.gti.exposure - deadtime, '[sec]' ) + + # add caldb key + # for key in cdbinfo: + # evthdu.header[key] = cdbinfo[key] + + # --- make file --- + hdulist = fits.HDUList([hdu, evthdu, srchdul['KVEA'], srchdul['HK'], gti.make_hdu(mjdref)]) + + hdulist.writeto(dstfile, overwrite=True, checksum=True) + + def setcol(self, name, arr): + self._events.setcol(name, arr) + + +if __name__ == '__main__': + pass diff --git a/arttools/src/arttools/fstools.py b/arttools/src/arttools/fstools.py new file mode 100755 index 0000000..2c76cb3 --- /dev/null +++ b/arttools/src/arttools/fstools.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2025 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 + +from enum import Enum + + +class Dir(object): + + class Type(Enum): + path = 1 + file = 2 + + def __init__(self, path): + self._path = path + + @property + def value(self): + return self._path + + def get(self): + return self._path + + def pathfor(self, name): + return os.path.join(self._path, name) + + def __truediv__(self, path): + return Dir(os.path.join(self._path, path)) + + def ls(self, type=None, filt=None): + flist = [f for f in os.listdir(self._path)] + flist.sort() + + if filt is not None: + flist = [f for f in flist if filt(f)] + + typefn = { + self.Type.path: os.path.isdir , + self.Type.file: os.path.isfile, + }.get(type, None) + + if typefn is not None: + flist = [f for f in flist if typefn(os.path.join(self._path, f))] + + return flist + + def lsdirs(self, filt=None): + return self.ls(type=self.Type.path, filt=filt) + + def lsfiles(self, filt=None): + return self.ls(type=self.Type.file, filt=filt) + + def empty(self): + if not os.path.isdir(self._path): + return False + + return len(os.listdir(self._path)) == 0 + + def delete(self): + if self.exists(): + shutil.rmtree(self._path) + + def exists(self): + if self._path is None: + return False + + return os.path.exists(self._path) + + def check(self, clean=False): + if clean: + self.delete() + + if not self.exists(): + os.makedirs(self._path) + + +if __name__ == "__main__": + pass + diff --git a/arttools/src/arttools/gti.py b/arttools/src/arttools/gti.py new file mode 100644 index 0000000..74cb876 --- /dev/null +++ b/arttools/src/arttools/gti.py @@ -0,0 +1,254 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2020-2025 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the SRG/ART-XC software project. +# ----------------------------------------------------------------------------- + +import numpy as np +from astropy.io import fits + +class Gti(object): + + def __init__(self, array=np.zeros((0,2), dtype=np.double), start=None, stop=None): + self._array = np.reshape(np.asarray(array), (-1,2)) + + self._start = self._ntime(start, default=lambda: self._array.min(axis=0)[0]) + self._stop = self._ntime(stop , default=lambda: self._array.max(axis=0)[1]) + + @property + def array(self): + return self._array + + @property + def empty(self): + return self._array.size == 0 + + @property + def size(self) : + return self._array.shape[0] + + @property + def start(self): + return self._start + + @property + def stop(self) : + return self._stop + + @property + def exposure(self): + return -np.sum(np.subtract.reduce(self._array, axis=1)) + + @property + def onsource(self): + return self._stop - self._start + + def _ntime(self, val, default): + if val is not None: + return val + if self._array.size == 0: + return None + return default() + + def mask(self, times): + masks = [] + for istart, istop in self._array: + masks.append(np.logical_and(times >= istart, times < istop)) + + return np.any(masks, axis=0) + + @staticmethod + def from_mask(times, mask): + slices = np.ma.clump_masked(np.ma.masked_where(mask, mask)) + gtis = [[times[s.start], times[s.stop - 1]] for s in slices] + return Gti(gtis) + + @staticmethod + def from_hduname(fitsfile, hduname): + with fits.open(fitsfile) as hdul: + hdu = hdul[hduname] + + array = np.array([hdu.data['START'], hdu.data['STOP']]).T + try: + start = hdu.header['TSTART'] + stop = hdu.header['TSTOP' ] + except KeyError: + start = hdu.data['START'][ 0] + stop = hdu.data['STOP' ][-1] + + return Gti(array, start, stop) + + @staticmethod + def from_hdu(hdu): + array = np.array([hdu.data['START'], hdu.data['STOP']]).T + try: + start = hdu.header['TSTART'] + stop = hdu.header['TSTOP' ] + except KeyError: + start = hdu.data['START'][ 0] + stop = hdu.data['STOP' ][-1] + + return Gti(array, start, stop) + + def clone(self): + return Gti(self._array, self._start, self._stop) + + def normalize(self, edge_delta=0.5): + if self._array.size == 0: + return + + gti_edges_dt = self._array[1:, 0] - self._array[:-1, 1] + gti_edges_mask = gti_edges_dt > edge_delta / 2. + gti_indices = np.ravel(np.where(gti_edges_mask)) + + new_gtis = [] + current_index = 0 + + for index in gti_indices: + start = self._array[current_index, 0] + stop = self._array[index, 1] + new_gtis.append([start, stop]) + + current_index = index + 1 + + start = self._array[current_index, 0] + new_gtis.append([start, self._stop]) + + return Gti(new_gtis) + + def __eq__(self, other): + if self is other: return True + if other is None : return False + + return np.array_equal(self._array, other.array) \ + and (self._start == other.start) and (self._stop == other.stop) + + def __add__(self, other): + if self is other: return self + if other is None : return self + if other.empty : return self + if self.empty : return other + + array = self._union(other.array) + start = min(self._start, other.start) + stop = max(self._stop , other.stop ) + + return Gti(array, start, stop) + + def __sub__(self, other): + if self is other: return Gti() + if self.empty : return Gti() + if other is None : return self + if other.empty : return self + + array = self._diff(other.array) + + return Gti(array, self._start, self._stop) + + def __mul__(self, other): + if self is other: return self + if other is None : return self + if self.empty : return Gti() + if other.empty : return Gti() + + array = self._intersect(other.array) + start = min(self._start, other.start) + stop = max(self._stop , other.stop ) + + return Gti(array, start, stop) + + def _union(self, array): + array = np.concatenate((self._array, array), axis=0) + array = np.sort(array, axis=0) + + starts = array[:,0] + stops = array[:,1] + + mask = np.ones(array.shape[0] + 1, dtype=bool) + mask[ 1:-1] = starts[1:] >= stops[:-1] + + return np.vstack((starts[:][mask[:-1]], stops[:][mask[1:]])).T + + def _diff(self, array): + starts = self._array[:,0] + stops = self._array[:,1] + + for item in array: + istart = item[0] + istop = item[1] + + mask_lower = np.logical_and(istart >= starts, istart <= stops) + mask_upper = np.logical_and(istop >= starts, istop <= stops) + + starts[mask_lower] = istart + stops [mask_upper] = istop + + return np.dstack((starts, stops)) + + def _intersect_helper(self, arr1, arr2): + newarray = [] + + for start1, stop1 in arr1: + for start2, stop2 in arr2: + + start = max(start1, start2) + stop = min(stop1 , stop2 ) + + has_intersections = stop > start + if has_intersections: + newarray.append([start, stop]) + + return newarray + + def _intersect(self, array): + newarray = [] + + newarray1 = self._intersect_helper(self._array, array) + newarray2 = self._intersect_helper(array, self._array) + + newarray.extend(newarray1) + newarray.extend(newarray2) + + return np.unique(np.array(newarray), axis=0) + + def make_hdu(self, mjdref): + a1 = np.array(self._array[:, 0]) + a2 = np.array(self._array[:, 1]) + c1 = fits.Column(name='START', format='1D', unit='sec', array=a1) + c2 = fits.Column(name='STOP' , format='1D', unit='sec', array=a2) + cols = fits.ColDefs([c1, c2]) + + gtihdu = fits.BinTableHDU.from_columns(cols) + gtihdu.name = 'STDGTI' + gtihdu.header['TSTART'] = self._start + gtihdu.header['TSTOP' ] = self._stop + gtihdu.header['MJDREF'] = (mjdref, 'MJD corresponding to SC clock start 2000.0 MSK') + + return gtihdu + + def write_fits(self, gtifile, mjdref): + # --- GTI file structure + hdu = fits.PrimaryHDU() + + a1 = np.array(self._array[:, 0]) + a2 = np.array(self._array[:, 1]) + c1 = fits.Column(name='START', format='1D', unit='sec', array=a1) + c2 = fits.Column(name='STOP' , format='1D', unit='sec', array=a2) + cols = fits.ColDefs([c1, c2]) + + gtihdu = fits.BinTableHDU.from_columns(cols) + gtihdu.name = 'STDGTI' + gtihdu.header['TSTART'] = self._start + gtihdu.header['TSTOP' ] = self._stop + gtihdu.header['MJDREF'] = (mjdref, 'MJD corresponding to SC clock start 2000.0 MSK') + + # --- wirte GTI + hdulist = fits.HDUList([hdu, gtihdu]) + hdulist.writeto(gtifile) + + +if __name__ == '__main__': + pass diff --git a/arttools/src/arttools/heatools.py b/arttools/src/arttools/heatools.py new file mode 100755 index 0000000..2ae5afa --- /dev/null +++ b/arttools/src/arttools/heatools.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2024 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 + +from arttools.cmd import Command, ExecutionResult +from arttools.errors import Error + + +class HeatoolsError(Error): + pass + + +class Heatools(object): + _ftcopy = 'ftcopy {src}[{ext}] {dst} copyall=no history=no clobber=yes' + _ftappend = 'ftappend {srcfile}[{srcext}]{filter} {dstfile} history=yes' + _maketime = 'maketime {srcfile}[{srcext}] {gtifile} "{filter_expr}" time=TIME prefr=0 postfr=0 compact=no histkw=no' + _ftmgtime = 'ftmgtime {gtifilelist} {dstgti} merge=and emptygti=apply history=yes clobber=yes' + _ftselect = 'ftselect "{srcfile}[{srcext}]" {dstfile} \'gtifilter("{gtifile}")\' copyall=no history=yes clobber=yes' + _ftchksum = 'ftchecksum {srcfile} update=yes chatter=0' + + @classmethod + def check(cls, envvar='HEADAS'): + if not envvar in os.environ: + raise HeatoolsError('error: Heasoft environment not found (is it initialzed?)') + + @classmethod + def ftcopy(cls, src, dst, ext): + Command( + cls._ftcopy.format( + src=src, + dst=dst, + ext=ext + ) + ).run() + + @classmethod + def ftappend(cls, srcfile, srcext, dstfile, filt=''): + Command( + cls._ftappend.format( + srcfile=srcfile, + srcext=srcext, + filter=filt, + dstfile=dstfile + ) + ).run() + + def maketime(self, srcfile, srcext, gtifile, filter_expr): + Command( + self._maketime.format( + srcfile=srcfile, + srcext=srcext, + gtifile=gtifile, + filter_expr=filter_expr + ) + ).run() + + def ftmgtime(self, gtifilelist, dstgti): + Command( + self._ftmgtime.format( + gtifilelist=gtifilelist, + dstgti=dstgti + ) + ).run() + + def ftselect(self, srcfile, srcext, dstfile, gtifile): + Command( + self._ftselect.format( + srcfile=srcfile, + srcext=srcext, + dstfile=dstfile, + gtifile=gtifile + ) + ).run() + + def ftchecksum(self, srcfile): + Command( + self._ftchksum.format( + srcfile=srcfile + ) + ).run() + + +if __name__ == '__main__': + pass + + diff --git a/arttools/src/arttools/printer.py b/arttools/src/arttools/printer.py new file mode 100755 index 0000000..c4bf474 --- /dev/null +++ b/arttools/src/arttools/printer.py @@ -0,0 +1,137 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# +# ----------------------------------------------------------------------------- +# Copyright (c) 2024 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the SRG/ART-XC software project. +# ----------------------------------------------------------------------------- + +import sys +import logging +from enum import Enum + +from datetime import datetime + + +# class Log(object): + +# @staticmethod +# def configure(who, logfile): +# import logging + +# logging.basicConfig( +# level=logging.INFO, +# format='{}: %(asctime)s [%(levelname)-7s] %(message)s'.format(who), +# handlers=[ +# logging.FileHandler(logfile), +# logging.StreamHandler() +# ] +# ) + + +class PrintType(Enum): + LOG = 'log' + INFO = 'info' + WARNING = 'warning' + ERROR = 'error' + + +class PrintFormatter(object): + + def __init__(self): + self._newline = '\n' + self._line_width = 80 + self._line_type_1 = '-' * self._line_width + self._line_type_2 = '=' * self._line_width + self._line_type_3 = '*' * 5 + self._offset_w = 0 + self._offset_str = ' ' * self._offset_w + self._separator = self._newline + self._offset_str + + def offset(self, offset): + self._offset_w = offset + self._offset_str = ' ' * self._offset_w + self._separator = self._newline + self._offset_str + + def header(self, toolname, version): + return self.multiline([ + '[' + self._timestr() + '] running: {} v.{}'.format(toolname, version), + self._line_type_2 + ]) + + def footer(self, toolname, version, message): + return self.multiline([ + '{} v.{}: exited with status: {}'.format(toolname, version, message), + self._line_type_1 + ]) + + def log(self, message): + return self._offset_str + message + + def message(self, level, message): + return self._format(level, message) + + def multiline(self, vmessage): + return self._separator.join(vmessage) + + def _format(self, level, message): + return '[{}] {:>7}: {}'.format(self._timestr(), level, message) + + def _timestr(self): + return str(datetime.now()) + + +class Printer(object): + + def __init__(self, toolname, version, logfile=None): + super().__init__() + self._toolname = toolname + self._version = version + self._formatter = PrintFormatter() + + # -- setup logger + logfmt = logging.Formatter('%(message)s') + self._logger = logging.getLogger(self._toolname) + self._logger.setLevel(logging.INFO) + self._logger.handlers = [] # reset log handlers + + logsh = logging.StreamHandler(sys.stdout) + logsh.setFormatter(logfmt) + self._logger.addHandler(logsh) + + if logfile is not None: + logfh = logging.FileHandler(logfile) + logfh.setFormatter(logfmt) + self._logger.addHandler(logfh) + + def info(self, message): + self._logger.info( + self._formatter.message('INFO' , message)) + + def warning(self, message): + self._logger.warning( + self._formatter.message('WARNING', message)) + + def error(self, message): + self._logger.error( + self._formatter.message('ERROR' , message)) + + def log(self, message): + self._logger.info( + self._formatter.log(message)) + + def header(self): + self._logger.info( + self._formatter.header(self._toolname, self._version)) + + def footer(self, message): + self._logger.info( + self._formatter.footer(self._toolname, self._version, message)) + + +if __name__ == '__main__': + pass + + diff --git a/arttools/src/arttools/random.py b/arttools/src/arttools/random.py new file mode 100755 index 0000000..63844c7 --- /dev/null +++ b/arttools/src/arttools/random.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2025 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 + + +class Random(object): + + def __init__(self, seed=None): + super().__init__() + + self._seed = seed + self._rng = None + + self._initialize() + + def _initialize(self): + if self._seed is None: + sseq = np.random.SeedSequence() + self._seed = sseq.entropy + else: + sseq = np.random.SeedSequence(self._seed) + + self._rng = np.random.default_rng(sseq) + + @property + def seed(self): return self._seed + + def uniform(self, low, high, size=None): + return self._rng.uniform(low, high, size) + + +if __name__ == '__main__': + pass + + diff --git a/arttools/src/arttools/task.py b/arttools/src/arttools/task.py new file mode 100755 index 0000000..ad4e03f --- /dev/null +++ b/arttools/src/arttools/task.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2025 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 + +from abc import ABC, abstractmethod + +from arttools.cmd import Command +from arttools.errors import Error + + +class TaskError(Error): + pass + + +class TaskArgs(ABC): + + def __init__(self, cfgfile=None): + super().__init__() + + self.pid = os.getpid() + self.env = {} + self.cfg = {} + + self.read_env() + self.read_cfg(cfgfile) + + @abstractmethod + def read_env(self): + pass + + @abstractmethod + def read_cfg(self, cfgfile): + pass + + @abstractmethod + def read_args(self): + pass + + def read(self): + self.read_args() + + @abstractmethod + def check(self): + pass + + +class TaskCmds(ABC): + + def __init__(self): + super().__init__() + + def run(self): + pass + + +class Task(ABC): + + def __init__(self): + super().__init__() + + @abstractmethod + def initialize(self): + pass + + @abstractmethod + def finalize(self): + pass + + @abstractmethod + def runmain(self): + pass + + def run(self): + try: + self.initialize() + self.runmain() + self.finalize() + except TaskError as e: + print('error: task execution:', e) + except Error as e: + print('error: unknown:', e) + +if __name__ == "__main__": + pass + diff --git a/arttools/src/arttools/telescope.py b/arttools/src/arttools/telescope.py new file mode 100644 index 0000000..6f765d4 --- /dev/null +++ b/arttools/src/arttools/telescope.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2016-2025 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 Teldef(object): + F = 2693. + d = .595 + centeroffset = 23.5 + assemblies = [ + {'detn': 33, 'urdn': 28, 'teln': 1}, + {'detn': 34, 'urdn': 22, 'teln': 2}, + {'detn': 24, 'urdn': 23, 'teln': 3}, + {'detn': 25, 'urdn': 24, 'teln': 4}, + {'detn': 22, 'urdn': 25, 'teln': 5}, + {'detn': 26, 'urdn': 26, 'teln': 6}, + {'detn': 23, 'urdn': 30, 'teln': 7} + ] + + urdalist = ['02', '04', '08', '10', '20', '40', '80'] + urdnlist = [ 28 , 22 , 23 , 24 , 25 , 26 , 30 ] + detnlist = [ 33 , 34 , 24 , 24 , 22 , 26 , 23 ] + telnlist = [ 1 , 2 , 3 , 4 , 5 , 6 , 7 ] + telnames = [''.join(('T', str(teln))) for teln in telnlist] + + @classmethod + def telname(cls, teln): + if teln not in cls.telnlist: + return 'T0' + return ''.join(('T', str(teln))) + + @classmethod + def teln(cls, detn, urdn): + for item in cls.assemblies: + if item['detn'] == detn and item['urdn'] == urdn: + return item['teln'] + + return 0 + +if __name__ == '__main__': + pass diff --git a/arttools/src/arttools/toolargs.py b/arttools/src/arttools/toolargs.py new file mode 100755 index 0000000..47bf68f --- /dev/null +++ b/arttools/src/arttools/toolargs.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2016-2024 IKI RAS (http://iki.rssi.ru/) +# Space Research Institute of the Russian Academy of Science +# +# This file is part of the M. Pavlinksy SRG/ART-XC telescope software. +# ----------------------------------------------------------------------------- + +import os +import shutil + + +class ToolArgs(object): + + def __init__(self, configfile=None): + self._configfile = configfile + + self.pid = os.getpid() + self.config = {} + + self._read_env() + + def _read_env(self): + if self._configfile is None: + return + + def make_tmpdir(self, path): + return os.path.join(path, str(self.pid)) + + def clean_path(self, path): + if os.path.exists(path): + shutil.rmtree(path) + + def check_path(self, path, clean=False): + if clean and os.path.exists(path): + shutil.rmtree(path) + if not os.path.exists(path): + os.makedirs(path) + + +if __name__ == '__main__': + pass + + diff --git a/arttools/src/arttools/version.py b/arttools/src/arttools/version.py new file mode 100755 index 0000000..2940e26 --- /dev/null +++ b/arttools/src/arttools/version.py @@ -0,0 +1 @@ +__version__ = '2.0.8' diff --git a/arttools/src/arttools/wcs.py b/arttools/src/arttools/wcs.py new file mode 100644 index 0000000..05919b2 --- /dev/null +++ b/arttools/src/arttools/wcs.py @@ -0,0 +1,143 @@ +#!/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 json + +import numpy as np + +from astropy.io import fits +from astropy.wcs import WCS + +import astropy.units as u + +from astropy.coordinates import Angle +from astropy.coordinates import SkyCoord + + +class Fov(object): + + def __init__(self, ra_cen, dec_cen): + self._ra_cen = ra_cen + self._dec_cen = dec_cen + + fov_size = 0.66 * u.deg + pix_size = 45.5728027182 * u.arcsec / u.pix + + fov_wcs, fov_shape = \ + self._make_fov_wcs( + self._ra_cen, + self._dec_cen, + fov_size, + pix_size + ) + + self._wcs = Wcs() + self._wcs.init(fov_wcs) + + self._shape = fov_shape + + @property + def wcs(self): + return self._wcs + + @property + def ra(self): + return self._ra_cen + + @property + def dec(self): + return self._dec_cen + + def _make_fov_wcs(self, + ra_cen, dec_cen, + fov_size, pix_size, + ref_pix=None, + projection='TAN' ): + + center = SkyCoord(ra_cen, dec_cen, unit=u.deg) + + fov_width = Angle(fov_size, unit=u.deg) + fov_height = Angle(fov_size, unit=u.deg) + + pix_width = Angle(pix_size) + pix_height = Angle(pix_size) + + nx = int(np.ceil((fov_width / pix_width ).decompose().value)) + ny = int(np.ceil((fov_height / pix_height).decompose().value)) + + if ref_pix is None: + ref_pix = [(nx + 1) / 2., (ny + 1) / 2.] + + wcs = Wcs(naxis=2) + + wcs.wcs.radesys = 'FK5' + wcs.wcs.equinox = 2000. + + wcs.wcs.crval = [center.ra.deg, center.dec.deg] + wcs.wcs.crpix = ref_pix + + wcs.wcs.cdelt = [-pix_width.value, pix_height.value] + + return wcs, (ny, nx) + + +class Wcs(object): + + def __init__(self): + self._strip_size = 45.5728027182 + + def init(self, wcs): + self._wcs = wcs + self._projection = self._make_projection(self._wcs) + + return self + + def from_other(self, projection): + self._wcs = None + self._projection = projection + + return self + + @staticmethod + def from_file(fname): + with open(fname) as wcsfile: + return Wcs().init(json.load(wcsfile)) + + @staticmethod + def from_eventfile(fname): + with fits.open(fname) as hdul: + projection = WCS(hdul['EVENTS'].header) + return Wcs().from_other(projection) + + def get(self): + return self._projection + + def world2pix(self, ra, dec): + y, x = ( + self._projection.all_world2pix( + np.array([ra, dec]).T, 1) - 0.5 + ).astype(int).T + + return x, y + + def _make_projection(self, wcs, subpixels=5): + pixsize = self._strip_size / 3600. / subpixels + + projection = WCS(naxis=2) + projection.wcs.crpix = [wcs['nx' ], wcs['ny' ]] + projection.wcs.crval = [wcs['crval1'], wcs['crval2']] + projection.wcs.cdelt = [pixsize, pixsize] + projection.wcs.ctype = ["RA---TAN", "DEC--TAN"] + projection.wcs.radesys = "FK5" + + return projection + + +if __name__ == '__main__': + pass diff --git a/arttools/tests/test_artdatatable.py b/arttools/tests/test_artdatatable.py new file mode 100644 index 0000000..6e0564f --- /dev/null +++ b/arttools/tests/test_artdatatable.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2020 Space Research Institute (http://iki.rssi.ru/) +# +# This file is part of ARTXC software project. +# ----------------------------------------------------------------------------- + +import os, sys +import unittest +import unittest.mock as mock + +import numpy as np + +sys.path.append('../src') + +from pyartxc.gti import Gti + + +class Test_(unittest.TestCase): + + def setUp(self): + pass + + def test_create_empty(self): + pass + + +if __name__ == '__main__': + unittest.main() diff --git a/arttools/tests/test_arttools_gti.py b/arttools/tests/test_arttools_gti.py new file mode 100755 index 0000000..4776657 --- /dev/null +++ b/arttools/tests/test_arttools_gti.py @@ -0,0 +1,208 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# ----------------------------------------------------------------------------- +# Copyright (c) 2020-2025 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, sys +import unittest +import unittest.mock as mock + +import numpy as np + +sys.path.append(os.path.join(os.getcwd(), '../src/')) + + +import arttools +from arttools.gti import Gti + + +class TestGti(unittest.TestCase): + + def setUp(self): + self.gti_0 = Gti() + self.gti_1 = Gti([[100., 200.],]) + self.gti_2 = Gti([[100., 200.], [300, 305]]) + self.gti_3 = Gti([100., 200., 300., 400.]) + + def test_create_empty(self): + self.assertTrue (self.gti_0.empty) + self.assertEqual(self.gti_0.size, 0) + + def test_create_from_intervals_1(self): + self.assertFalse(self.gti_1.empty) + self.assertEqual(self.gti_1.size, 1) + self.assertEqual(self.gti_1.exposure(), 100.) + + def test_create_from_intervals_2(self): + self.assertFalse(self.gti_2.empty) + self.assertEqual(self.gti_2.size, 2) + self.assertEqual(self.gti_2.exposure(), 105.) + + def test_create_from_intervals_3(self): + with self.assertRaises(ValueError): + Gti([100., 200., 300.]) + + def test_startstop(self): + self.assertEqual(self.gti_0.start, None) + self.assertEqual(self.gti_0.stop , None) + self.assertEqual(self.gti_2.start, 100.) + self.assertEqual(self.gti_2.stop , 305.) + + def test_create_from_timemask(self): + times = np.array([100, 200, 300, 400, 500, 600]) + param = np.array([-99, -98, -10, 0, -98, -99]) + + gti = Gti.from_timemask(times, param < -95) + + self.assertFalse(gti.empty) + self.assertEqual(gti.size , 2) + self.assertEqual(gti.exposure(), 200) + self.assertEqual(gti.start , 100) + self.assertEqual(gti.stop , 600) + + def test_clone(self): + gti_clone = self.gti_2.clone() + + self.assertFalse(gti_clone.empty) + self.assertEqual(gti_clone.size, 2) + self.assertEqual(gti_clone.exposure(), 105.) + + def test_equal(self): + gti = Gti([100., 200.]) + + self.assertTrue (gti == gti ) + self.assertTrue (gti == self.gti_1) + self.assertFalse(gti == self.gti_2) + self.assertFalse(gti == None ) + + def test_union_1(self): + gti_1 = Gti([100., 200.]) + gti_2 = Gti() + + res_gti_1 = gti_1 + None + res_gti_2 = gti_1 + gti_1 + res_gti_3 = gti_1 + gti_2 + res_gti_4 = gti_2 + gti_1 + + self.assertTrue(res_gti_1 == gti_1 ) + self.assertTrue(res_gti_2 == gti_1 ) + self.assertTrue(res_gti_3 == gti_1 ) + self.assertTrue(res_gti_4 == gti_1 ) + self.assertTrue(res_gti_3 == res_gti_4) + + def test_union_2(self): + gti_1 = Gti([100., 200.]) + gti_2 = Gti([150., 250.]) + + res_gti_1 = gti_1 + gti_2 + res_gti_2 = gti_2 + gti_1 + + self.assertEqual(res_gti_1.size , 1) + self.assertEqual(res_gti_1.exposure(), 150) + self.assertEqual(res_gti_1.start , 100) + self.assertEqual(res_gti_1.stop , 250) + self.assertTrue (res_gti_1 == res_gti_2) + + def test_union_3(self): + gti_1 = Gti([100., 200., 300., 400.]) + gti_2 = Gti([150., 250., 350., 450., 370., 470.]) + gti_3 = Gti([170., 270.]) + + res_gti = gti_1 + gti_2 + gti_3 + + self.assertEqual(res_gti.size , 2) + self.assertEqual(res_gti.exposure(), 340) + self.assertEqual(res_gti.start , 100) + self.assertEqual(res_gti.stop , 470) + + def test_subtraction_1(self): + gti_1 = Gti([100., 200.]) + gti_2 = Gti() + + res_gti_1 = gti_1 - None + res_gti_2 = gti_1 - gti_1 + res_gti_3 = gti_1 - gti_2 + res_gti_4 = gti_2 - gti_1 + + self.assertTrue(res_gti_1 == gti_1 ) + self.assertTrue(res_gti_2 == self.gti_0) + self.assertTrue(res_gti_3 == gti_1 ) + self.assertTrue(res_gti_4 == self.gti_0) + + def test_subtraction_2(self): + gti_1 = Gti([100., 200.]) + gti_2 = Gti([150., 250.]) + + res_gti = gti_1 - gti_2 + + self.assertEqual(res_gti.size , 1) + self.assertEqual(res_gti.start , 100) + self.assertEqual(res_gti.stop , 200) + self.assertEqual(res_gti.exposure(), 50) + self.assertTrue( + np.array_equal(res_gti.array, np.array([[150., 200.]]))) + + + def test_subtraction_3(self): + gti_2 = Gti([50., 120., 170., 250.]) + gti_1 = Gti([100., 200.]) + + res_gti = gti_2 - gti_1 + + self.assertEqual(res_gti.size , 2) + self.assertEqual(res_gti.start , 50) + self.assertEqual(res_gti.stop , 250) + self.assertEqual(res_gti.exposure(), 50) + self.assertTrue( + np.array_equal(res_gti.array, np.array([[100., 120.], [170., 200.]]))) + + def test_intersection_1(self): + gti_1 = Gti([100., 200.]) + gti_2 = Gti() + + res_gti_1 = gti_1 * None + res_gti_2 = gti_1 * gti_1 + res_gti_3 = gti_1 * gti_2 + res_gti_4 = gti_2 * gti_1 + + self.assertTrue(res_gti_1 == self.gti_0) + self.assertTrue(res_gti_2 == gti_1 ) + self.assertTrue(res_gti_3 == self.gti_0) + self.assertTrue(res_gti_4 == self.gti_0) + + def test_intersection_2(self): + gti_1 = Gti([100., 200.]) + gti_2 = Gti([50., 120., 170., 250.]) + + res_gti = gti_2 * gti_1 + + self.assertEqual(res_gti.size , 2) + self.assertEqual(res_gti.start , 50) + self.assertEqual(res_gti.stop , 250) + self.assertEqual(res_gti.exposure(), 50) + self.assertTrue( + np.array_equal(res_gti.array, np.array([[100., 120.], [170., 200.]]))) + + def test_intersection_3(self): + gti_1 = Gti([100., 200.]) + gti_2 = Gti([50., 120., 170., 250.]) + + res_gti = gti_1 * gti_2 + + self.assertEqual(res_gti.size , 2) + self.assertEqual(res_gti.start , 50) + self.assertEqual(res_gti.stop , 250) + self.assertEqual(res_gti.exposure(), 50) + self.assertTrue( + np.array_equal(res_gti.array, np.array([[100., 120.], [170., 200.]]))) + + # def test_inversion(self): + # pass + + +if __name__ == '__main__': + unittest.main() diff --git a/config/.DS_Store b/config/.DS_Store new file mode 100644 index 0000000..c9cbc73 Binary files /dev/null and b/config/.DS_Store differ diff --git a/config/artinit-dev-ashty.sh b/config/artinit-dev-ashty.sh new file mode 100755 index 0000000..72d4b9a --- /dev/null +++ b/config/artinit-dev-ashty.sh @@ -0,0 +1,3 @@ +export ARTCONFIG="/home/ashty/devel/projects/iki/artxc/software/artpipeline/config/artpipeline.cfg" + +conda activate /home/ashty/devel/projects/iki/artxc/env/newartenv diff --git a/config/artinit.sh b/config/artinit.sh new file mode 100755 index 0000000..72d4b9a --- /dev/null +++ b/config/artinit.sh @@ -0,0 +1,3 @@ +export ARTCONFIG="/home/ashty/devel/projects/iki/artxc/software/artpipeline/config/artpipeline.cfg" + +conda activate /home/ashty/devel/projects/iki/artxc/env/newartenv diff --git a/config/artpipeline.cfg b/config/artpipeline.cfg new file mode 100644 index 0000000..5124dda --- /dev/null +++ b/config/artpipeline.cfg @@ -0,0 +1,10 @@ +{ + "inbox": "/home/ashty/data/iki/mvn_sandbox/mvn/inbox", + "newbox": "/home/ashty/data/iki/mvn_sandbox/new", + "tmpdir": "/home/ashty/data/iki/mvn_sandbox/temp", + "mvndir": "/home/ashty/data/iki/mvn_sandbox/mvn", + "datadir": "/home/ashty/data/iki/mvn_sandbox/mvn/data", + "indexdir": "/home/ashty/data/iki/mvn_sandbox/mvn/index", + "dnlogsdir": "/home/ashty/data/iki/mvn_sandbox/mvn/dnlogs", + "vvlbasedir": "/home/ashty/data/iki/mvn_sandbox/mvn/vvlbase" +} diff --git a/config/devel/activate-conda.sh b/config/devel/activate-conda.sh new file mode 100755 index 0000000..4ddb08e --- /dev/null +++ b/config/devel/activate-conda.sh @@ -0,0 +1 @@ +conda activate ${PWD}/../mvnenv/ \ No newline at end of file diff --git a/config/devel/activate.sh b/config/devel/activate.sh new file mode 100755 index 0000000..4ddb08e --- /dev/null +++ b/config/devel/activate.sh @@ -0,0 +1 @@ +conda activate ${PWD}/../mvnenv/ \ No newline at end of file diff --git a/config/devel/install.sh b/config/devel/install.sh new file mode 100755 index 0000000..b7df081 --- /dev/null +++ b/config/devel/install.sh @@ -0,0 +1,5 @@ +#!/usr/bin/env bash + + +pip install mvntools/ +pip install mvndac/ diff --git a/config/devel/uninstall.sh b/config/devel/uninstall.sh new file mode 100755 index 0000000..b84ace5 --- /dev/null +++ b/config/devel/uninstall.sh @@ -0,0 +1,8 @@ +#!/usr/bin/env bash + + +pip uninstall mvntools --yes +rm -rf mvntools/build + +pip uninstall mvndac --yes +rm -rf mvndac/build diff --git a/config/setup/conda/configure.sh b/config/setup/conda/configure.sh new file mode 100755 index 0000000..beea544 --- /dev/null +++ b/config/setup/conda/configure.sh @@ -0,0 +1,11 @@ +#!/usr/bin/env bash + +echo "cleaning previous environment..." +rm -rf .venv + +echo "creating new environment..." +conda create --yes --prefix $(pwd)/env python=3.11 pip + +echo "finished" + + diff --git a/config/setup/conda/environment.txt b/config/setup/conda/environment.txt new file mode 100755 index 0000000..ee9464f --- /dev/null +++ b/config/setup/conda/environment.txt @@ -0,0 +1,7 @@ + +rm -rf .venv + +conda create --yes --prefix \$PWD/.venv/ python=3.11 pip +conda activate \$PWD/.venv/ + +pip install -r requirements.txt diff --git a/config/setup/configure-dev.sh b/config/setup/configure-dev.sh new file mode 100755 index 0000000..175deef --- /dev/null +++ b/config/setup/configure-dev.sh @@ -0,0 +1,17 @@ +#!/usr/bin/env bash + + +echo "cleaning previous environment..." +rm -rf .venv + +echo "creating new environment..." +python3 -m venv .venv +source .venv/bin/activate + +echo "install pip..." +python3 -m ensurepip --default-pip + +echo "install requirements..." +pip install -r requirements-dev.txt + +echo "finished" diff --git a/config/setup/configure.sh b/config/setup/configure.sh new file mode 100755 index 0000000..019c86f --- /dev/null +++ b/config/setup/configure.sh @@ -0,0 +1,18 @@ +#!/usr/bin/env bash + +echo "cleaning previous environment..." +rm -rf .venv + +echo "creating new environment..." +python3 -m venv .venv +source .venv/bin/activate + +echo "install pip..." +python3 -m ensurepip --default-pip + +echo "install requirements..." +pip install -r requirements.txt + +echo "finished" + + diff --git a/requirements-dev.txt b/requirements-dev.txt new file mode 100644 index 0000000..6f1b5db --- /dev/null +++ b/requirements-dev.txt @@ -0,0 +1,10 @@ +pydantic +refurb +ruff +typer[all] +rich +numpy +scipy +polars +astropy +regions \ No newline at end of file diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..99a66ef --- /dev/null +++ b/requirements.txt @@ -0,0 +1,8 @@ +pydantic +typer[all] +rich +numpy +scipy +polars +astropy +regions