#!/usr/bin/env python
# S. Caroff
"""
ONLINE analysis:
Produce online alt az from RA DEC of the pointing position
Usage:
$> online_alt_az_dl1_maker
--input-file "./dl1_5202_0_10000evt.h5"
--RA 83.254
--DEC 102.35
"""
import argparse
import time
from pathlib import Path
import astropy.units as u
import numpy as np
import tables
from astropy.coordinates import AltAz, EarthLocation, SkyCoord
from astropy.time import Time
from astropy.utils import iers
from rta_to_lstchain_formator.utils.rta_argparse import EnumAction
from rta_to_lstchain_formator.utils.rta_logging import LOGGING_LEVELS, init_logging
iers.conf.auto_download = False
[docs]
def online_alt_az_dl1_maker_arg_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
description="Add alt-az values in-place in DL1 params in a file, from RA-DEC values already in the file",
formatter_class=argparse.ArgumentDefaultsHelpFormatter,
)
parser.add_argument(
"--DEC", type=float, dest="DEC", help="Declination of the telescope", default=None, required=True
)
parser.add_argument(
"--input-file",
"-f",
type=str,
dest="input_file",
help='DL1 file, example : "dl1_5202_0_10000evt.h5"',
default=None,
required=True,
)
parser.add_argument(
"--log_file",
"-l",
type=Path,
dest="log_file_path",
help="Path to a file where this process will log its execution",
)
parser.add_argument(
"--log_level",
type=LOGGING_LEVELS,
action=EnumAction,
dest="logging_level",
help="Logging level: all messages with a level severity above or equal the specified level will be logged.",
default=LOGGING_LEVELS.WARNING,
)
parser.add_argument(
"--nbins",
type=int,
dest="nbins",
help="number of bins for the interpolation of alt az",
default=10,
required=False,
)
parser.add_argument(
"--nb_tries",
type=int,
dest="nb_tries",
help="Number of times the software will attempt to open the input file",
default=20,
required=False,
)
parser.add_argument(
"--RA", type=float, dest="RA", help="Right Ascension of the telescope", default=None, required=True
)
parser.add_argument(
"--timeout",
type=float,
dest="open_timeout",
help="Time to wait between input file opening attempts",
default=0.5,
required=False,
)
return parser
[docs]
def main():
parser = online_alt_az_dl1_maker_arg_parser()
args = parser.parse_args()
init_logging(
log_header_str="RTA_online_altaz_maker", log_filename=args.log_file_path, log_level=args.logging_level
)
input_file = args.input_file # "dl3_LST-1.Run04190.*.fits"
obs_location = EarthLocation.from_geodetic(-17.891485 * u.deg, 28.761518 * u.deg, 2187 * u.m)
# Failing to open a file on fefs doesn't necessarily mean we won't succeed next time !
for _ in range(args.nb_tries - 1):
try:
hfile = tables.open_file(input_file, mode="r+")
break
except Exception:
print("Failed to open input file {}, re-trying in {}".format(args.input_file, args.open_timeout))
time.sleep(args.open_timeout)
else:
# try 1 last time, let tables error raise if failing again
hfile = tables.open_file(input_file, mode="r+")
Trigger_time = hfile.root.dl1.event.telescope.parameters.tel_001.col("trigger_time")
is_good_event_table = hfile.root.dl1.event.telescope.parameters.tel_001.col("is_good_event")
selection_mask = np.logical_and(Trigger_time < 4000000000.0, is_good_event_table > 0.5)
is_good_event_table = np.multiply(selection_mask, 1)
time_table = np.linspace(Trigger_time[selection_mask][0], Trigger_time[selection_mask][-1], args.nbins)
pointing_coord = SkyCoord(ra=args.RA, dec=args.DEC, unit=u.deg)
event_altaz = pointing_coord.transform_to(
AltAz(
obstime=Time(time_table, format="unix"),
location=obs_location,
obswl=0.35 * u.micron,
relative_humidity=0.5,
temperature=10 * u.deg_C,
pressure=790 * u.hPa,
)
)
alt_table = np.interp(Trigger_time, time_table, event_altaz.alt.rad)
az_table = np.interp(Trigger_time, time_table, event_altaz.az.rad)
table = hfile.root.dl1.event.telescope.parameters.tel_001
table.modify_column(start=0, stop=len(alt_table), step=1, column=alt_table, colname="alt_tel")
table.modify_column(start=0, stop=len(az_table), step=1, column=az_table, colname="az_tel")
table.modify_column(
start=0, stop=len(is_good_event_table), step=1, column=is_good_event_table, colname="is_good_event"
)
table.flush()
hfile.close()
if __name__ == "__main__":
main()