-
Notifications
You must be signed in to change notification settings - Fork 1
/
L0toL15.py
120 lines (97 loc) · 4.78 KB
/
L0toL15.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
from astropy.io import fits
from astropy.table import Table
import argparse
import sys
import os
import numpy as np
from math import pi
import pandas
import copy
import arttools
from arttools._det_spatial import get_shadowed_pix_mask_for_urddata
from arttools.energy import get_events_energy
from arttools.caldb import get_shadowmask, get_energycal
from arttools.time import get_gti
from astropy.table import Table
import socket
import smtplib
from email.message import EmailMessage
parser = argparse.ArgumentParser(description="process L0 data to L1 format")
parser.add_argument("stem", help="part of the L0 files name, which are euqal to them")
ARTCALDBPATH = os.environ["ARTCALDB"]
indexfname = "artxc_index.fits"
#caldbindex = pandas.DataFrame(fits.getdata(indexfname, "CIF"))
URDTOTEL = {28: "T1",
22: "T2",
23: "T3",
24: "T4",
25: "T5",
26: "T6",
30: "T7"}
if __name__ == "__main__":
try:
sock = socket.socket()
sock.connect(("10.5.2.24", 8081))
sock.send(str.encode(sys.argv[1]))
except ConnectionRefusedError as conerr:
msg = EmailMessage()
msg["to"] = "san@iki.rssi.ru"
msg["from"] = "art@cosmos.ru"
msg["subbj"] = "makeL1b fail to start daemon"
msg.set_content(conerr)
smtplib.SMTP("localhost").send_message(msg)
finally:
if len(sys.argv) != 4 or "-h" in sys.argv:
print("description run like that 'python3 L1toL15.py stem outdir'"\
", where stem is srg_20190727_214739_000")
raise ValueError("wrong arguments")
fname = sys.argv[1]
stem = fname.rsplit(".")[0]
outdir = sys.argv[2]
attfname = sys.argv[3]
if os.path.abspath(outdir) == os.path.abspath(os.path.dirname(stem)):
raise ValueError("The L0 files will be overwriten")
if not os.path.exists(outdir):
os.mkdir(outdir)
attdata = arttools.plot.get_attdata(attfname)
urdfname = fname
urdfile = fits.open(urdfname)
urddata = Table(urdfile["EVENTS"].data).as_array()
flag = np.ones(urddata.size, np.uint8)
locgti = get_gti(urdfile, gtiextname="KVEA_GTI")
locgti.merge_joint()
caldbfile = get_energycal(urdfile)
bkigti = arttools.time.make_bki_gti(urdfile)
flag[(locgti & attdata.gti).mask_outofgti_times(urddata["TIME"])] = 0
RA, DEC = np.empty(urddata.size, np.double), np.empty(urddata.size, np.double)
attmask = attdata.gti.mask_outofgti_times(urddata["TIME"])
RA[attmask], DEC[attmask] = arttools.orientation.get_photons_sky_coord(urddata[attmask], urdfile["EVENTS"].header["URDN"], attdata)
ENERGY, xc, yc, grades = get_events_energy(urddata, urdfile["HK"].data, caldbfile)
shadow = get_shadowmask(urdfile)
maskshadow = get_shadowed_pix_mask_for_urddata(urddata, shadow)
flag[np.logical_not(maskshadow)] = 2
flag[np.any([urddata["RAW_X"] == 0, urddata["RAW_X"] == 47, urddata["RAW_Y"] == 0, urddata["RAW_Y"] == 47], axis=0)] = 3
h = copy.copy(urdfile["EVENTS"].header)
h.pop("NAXIS2")
if bkigti.exposure > 0:
bkimask = bkigti.mask_outofgti_times(urddata["TIME"])
flag[bkimask] += 100
cols = urdfile["EVENTS"].data.columns
cols.add_col(fits.Column(name="ENERGY", array=ENERGY, format="1D", unit="keV"))
cols.add_col(fits.Column(name="RA", array=np.copy(RA*180./pi), format="1D", unit="deg"))
cols.add_col(fits.Column(name="DEC", array=np.copy(DEC*180./pi), format="1D", unit="deg"))
cols.add_col(fits.Column(name="GRADE", array=grades, format="I"))
cols.add_col(fits.Column(name="FLAG", array=flag, format="I"))
cols["TIME_I"].array = urddata["TIME_I"]
newurdtable = fits.BinTableHDU.from_columns(cols, header=h)
newurdtable.name = "EVENTS"
gtitable = fits.BinTableHDU(Table(locgti.arr, names=("START", "STOP")), header=urdfile["GTI"].header)
newfile = fits.HDUList([urdfile[0], newurdtable, urdfile["HK"], gtitable])
newfile.writeto(os.path.join(outdir, os.path.basename(urdfname)), overwrite=True)
if not os.path.exists(os.path.join(outdir, os.path.basename(attfname))):
attfile = fits.open(attfname)
hdus = [attfile[0], ]
for telescope in arttools.telescope.TELESCOPES:
ra, dec, roll = arttools.orientation.quat_to_pol_and_roll(attdata(attdata.times)*arttools.caldb.ARTQUATS[telescope])
hdus.append(fits.BinTableHDU(Table(np.array([attdata.times, ra*180/pi, dec*180/pi, roll*180/pi]).T, names=("TIME", "RA", "DEC", "ROLL")), header={"TELESCP":telescope}, name=telescope))
fits.HDUList(hdus).writeto(os.path.join(outdir, os.path.basename(attfname)))