From 617570407051fee13070a89c2da7c9e631fa5f06 Mon Sep 17 00:00:00 2001 From: TjarkMiener Date: Wed, 3 Jul 2024 15:19:43 +0200 Subject: [PATCH] use astropy to get the sky offset and separation from tel poiting and source --- dl1_data_handler/reader.py | 26 ++++++--- dl1_data_handler/transforms.py | 97 ++++++++++++++++------------------ 2 files changed, 64 insertions(+), 59 deletions(-) diff --git a/dl1_data_handler/reader.py b/dl1_data_handler/reader.py index 9e80bb75..db313b90 100644 --- a/dl1_data_handler/reader.py +++ b/dl1_data_handler/reader.py @@ -9,6 +9,7 @@ from dl1_data_handler.processor import DL1DataProcessor import astropy.units as u +from astropy.coordinates import SkyCoord from astropy.table import ( Table, join, # let us merge tables horizontally @@ -129,7 +130,7 @@ def __init__( # Telescope pointings self.telescope_pointings = {} - self.fix_pointing_alt, self.fix_pointing_az = None, None + self.fix_pointing = None if self.process_type == "Observation": for tel_id in self.tel_ids: with lock: @@ -147,17 +148,26 @@ def __init__( # Only fix telescope pointings valid for MCs! # No divergent pointing implemented! - self.fix_pointing_alt = self.telescope_pointings[f"tel_{tel_id:03d}"][ + fix_pointing_alt = self.telescope_pointings[f"tel_{tel_id:03d}"][ "telescope_pointing_altitude" - ][0] - self.fix_pointing_az = self.telescope_pointings[f"tel_{tel_id:03d}"][ + ] + fix_pointing_az = self.telescope_pointings[f"tel_{tel_id:03d}"][ "telescope_pointing_azimuth" - ][0] - # Set the telescope pointing to the delta Alt/Az tranform + ] + # Reading the pointing for the first obs_id assuming fix tel pointing + fix_pointing_az = fix_pointing_az[0] * fix_pointing_az.unit + fix_pointing_alt = fix_pointing_alt[0] * fix_pointing_alt.unit + self.fix_pointing = SkyCoord( + fix_pointing_az.to_value(u.deg), + fix_pointing_alt.to_value(u.deg), + frame="altaz", + unit="deg", + ) + # Set the telescope pointing of the SkyOffsetSeparation tranform to the fix pointing if transforms is not None: for transform in transforms: - if transform.name == "deltaAltAz": - transform.set_pointing(self.fix_pointing_alt, self.fix_pointing_az) + if transform.name == "SkyOffsetSeparation": + transform.set_pointing(self.fix_pointing) # AI-based trigger system self.trigger_settings = trigger_settings diff --git a/dl1_data_handler/transforms.py b/dl1_data_handler/transforms.py index 83eeddb2..48b15cad 100644 --- a/dl1_data_handler/transforms.py +++ b/dl1_data_handler/transforms.py @@ -1,14 +1,14 @@ +import astropy.units as u +from astropy.coordinates import SkyCoord import numpy as np import itertools from .processor import Transform class ShowerPrimaryID(Transform): - def __init__( - self, name="true_shower_primary_id", particle_id_col_name="true_shower_primary_id" - ): + def __init__(self): super().__init__() - self.particle_id_col_name = particle_id_col_name + self.particle_id_col_name = "true_shower_primary_id" self.shower_primary_id_to_class = { 0: 1, 101: 0, @@ -22,7 +22,7 @@ def __init__( 255: "hadron", } - self.name = name + self.name = "true_shower_primary_id" self.dtype = np.dtype("int8") def describe(self, description): @@ -57,19 +57,18 @@ def transform(self, example): return example -class MCEnergy(Transform): - def __init__(self, name="energy", energy_col_name="true_energy", unit="log(TeV)"): +class LogEnergy(Transform): + def __init__(self): super().__init__() - self.name = name - self.energy_col_name = energy_col_name + self.name = "energy" self.shape = 1 self.dtype = np.dtype("float32") - self.unit = unit + self.unit = "log(TeV)" def describe(self, description): self.description = [ {**des, "name": self.name, "dtype": self.dtype, "unit": self.unit} - if des["name"] == self.energy_col_name + if des["name"] == "true_energy" else des for des in description ] @@ -77,36 +76,21 @@ def describe(self, description): def __call__(self, example): for i, (val, des) in enumerate(zip(example, self.description)): - if des["base_name"] == self.energy_col_name: - example[i] = ( - np.array([np.log10(val)]) - if self.unit == "log(TeV)" - else np.array([val]) - ) + if des["base_name"] == "true_energy": + example[i] = np.log10(val) return example -class DeltaAltAz(Transform): - def __init__( - self, - base_name="direction", - alt_col_name="true_alt", - az_col_name="true_az", - deg2rad=True, - north_pointing_correction=True, - ): +class SkyOffsetSeparation(Transform): + def __init__(self, transform_to_rad=False): super().__init__() - self.name = "deltaAltAz" - self.base_name = base_name - self.alt_col_name = alt_col_name - self.az_col_name = az_col_name - self.deg2rad = deg2rad - self.north_pointing_correction = north_pointing_correction - self.shape = 2 + self.name = "SkyOffsetSeparation" + self.base_name = "direction" + self.shape = 3 self.dtype = np.dtype("float32") - self.unit = "rad" - self.pointing_alt, self.pointing_az = None, None + self.unit = u.rad if transform_to_rad else u.deg + self.fix_pointing = None def describe(self, description): self.description = description @@ -117,36 +101,47 @@ def describe(self, description): "base_name": self.base_name, "shape": self.shape, "dtype": self.dtype, - "unit": self.unit, + "unit": str(self.unit), } ) return self.description - def set_pointing(self, pointing_alt, pointing_az): - self.pointing_alt = pointing_alt - self.pointing_az = pointing_az - return + def set_pointing(self, fix_pointing): + self.fix_pointing = fix_pointing def __call__(self, example): for i, (val, des) in enumerate( itertools.zip_longest(example, self.description) ): - if des["base_name"] == self.alt_col_name: - alt = np.radians(example[i]) if self.deg2rad else example[i] - alt -= self.pointing_alt - elif des["base_name"] == self.az_col_name: - az = np.radians(example[i]) if self.deg2rad else example[i] - if self.north_pointing_correction and az > 3*np.pi/2: - az -= 2 * np.pi - az -= self.pointing_az + if des["base_name"] == "true_alt": + alt = example[i] + elif des["base_name"] == "true_az": + az = example[i] elif des["base_name"] == self.base_name: - example.append(np.array([alt, az])) + true_direction = SkyCoord( + az * u.deg, + alt * u.deg, + frame="altaz", + unit="deg", + ) + sky_offset = self.fix_pointing.spherical_offsets_to(true_direction) + angular_separation = self.fix_pointing.separation(true_direction) + example.append( + np.array( + [ + sky_offset[0].to_value(self.unit), + sky_offset[1].to_value(self.unit), + angular_separation.to_value(self.unit), + ] + ) + ) return example + class CoreXY(Transform): - def __init__(self, name="impact"): + def __init__(self): super().__init__() - self.name = name + self.name = "impact" self.shape = 2 self.dtype = np.dtype("float32") self.unit = "km"