Change how DataWriter writes simulation pointing

ctapipe/conftest.py

@@ -326,6 +326,27 @@
         return output
 
 
+@pytest.fixture(scope="session")
+def dl1_divergent_file(dl1_tmp_path):
+    from ctapipe.tools.process import ProcessorTool
+
+    path = "dataset://gamma_divergent_LaPalma_baseline_20Zd_180Az_prod3_test.simtel.gz"
+    output = dl1_tmp_path / "gamma_divergent.dl1.h5"
+
+    # prevent running process multiple times in case of parallel tests
+    with FileLock(output.with_suffix(output.suffix + ".lock")):
+        if output.is_file():
+            return output
+
+        argv = [
+            f"--input={path}",
+            f"--output={output}",
+            "--EventSource.focal_length_choice=EQUIVALENT",
+        ]
+        assert run_tool(ProcessorTool(), argv=argv, cwd=dl1_tmp_path) == 0
+        return output
+
+
 @pytest.fixture(scope="session")
 def dl1_camera_frame_file(dl1_tmp_path, prod5_gamma_simtel_path):
     """

ctapipe/containers.py

@@ -229,21 +229,33 @@ def is_invalid(pixel_status):
         return (pixel_status & gain_bits) == 0
 
 
+class TelescopeConfigurationIndexContainer(Container):
+    """Index to include for per-OB telescope configuration"""
+
+    default_prefix = ""
+    obs_id = obs_id_field()
+    tel_id = tel_id_field()
+
+
 class EventIndexContainer(Container):
-    """index columns to include in event lists, common to all data levels"""
+    """Index columns to include in event lists.
 
-    default_prefix = ""  # don't want to prefix these
+    Common to all data levels
+    """
+
+    default_prefix = ""
     obs_id = obs_id_field()
     event_id = event_id_field()
 
 
 class TelEventIndexContainer(Container):
     """
-    index columns to include in telescope-wise event lists, common to all data
-    levels that have telescope-wise information
+    index columns to include in telescope-wise event lists
+
+    Common to all data levels that have telescope-wise information.
     """
 
-    default_prefix = ""  # don't want to prefix these
+    default_prefix = ""
     obs_id = obs_id_field()
     event_id = event_id_field()
     tel_id = tel_id_field()
@@ -1057,6 +1069,7 @@ class TelescopePointingContainer(Container):
     between camera and sky coordinates.
     """
 
+    default_prefix = "telescope_pointing"
     azimuth = Field(nan * u.rad, "Azimuth, measured N->E", unit=u.rad)
     altitude = Field(nan * u.rad, "Altitude", unit=u.rad)
 

ctapipe/io/datawriter.py

@@ -16,6 +16,7 @@
 from ..containers import (
     ArrayEventContainer,
     SimulatedShowerDistribution,
+    TelescopeConfigurationIndexContainer,
     TelEventIndexContainer,
 )
 from ..core import Component, Container, Field, Provenance, ToolConfigurationError
@@ -277,9 +278,8 @@ def __init__(self, event_source: EventSource, config=None, parent=None, **kwargs
         self._setup_writer()
         self._setup_outputfile()
 
-        # store last pointing to only write unique poitings
-        self._last_pointing = None
-        self._last_pointing_tel = defaultdict(lambda: (np.nan * u.deg, np.nan * u.deg))
+        # store per-ob for which telescopes we've already written the fixed pointing
+        self._constant_telescope_pointing_written = defaultdict(set)
 
     def _setup_outputfile(self):
         self._subarray.to_hdf(self._writer.h5file)
@@ -303,8 +303,10 @@ def __call__(self, event: ArrayEventContainer):
         self._at_least_one_event = True
         self.log.debug("WRITING EVENT %s", event.index)
 
-        self._write_subarray_pointing(event)
         self._write_trigger(event)
+        # write fixed pointing only for simulation, observed data will have monitoring
+        if self._is_simulation:
+            self._write_constant_pointing(event)
 
         if event.simulation is not None and event.simulation.shower is not None:
             self._writer.write(
@@ -337,6 +339,33 @@ def __call__(self, event: ArrayEventContainer):
         if self.write_muon_parameters:
             self._write_muon_telescope_events(event)
 
+    def _write_constant_pointing(self, event):
+        """
+        Write pointing configuration from event data assuming fixed pointing over the OB.
+
+        This function mainly exists due to a limitation of sim_telarray files.
+        Pointing information is only written as part of triggerred array events,
+        even though it is constant over the run. It also is not written for all
+        telescopes, only for those for which it is "known", which seem to be
+        all telescopes that participated in an array event at least once.
+
+        So we write the first pointing information for each telescope into the
+        configuration table.
+        """
+        obs_id = event.index.obs_id
+
+        for tel_id, pointing in event.pointing.tel.items():
+            if tel_id in self._constant_telescope_pointing_written[obs_id]:
+                continue
+            index = TelescopeConfigurationIndexContainer(
+                obs_id=obs_id,
+                tel_id=tel_id,
+            )
+            self._writer.write(
+                f"configuration/telescope/pointing/tel_{tel_id:03d}", (index, pointing)
+            )
+            self._constant_telescope_pointing_written[obs_id].add(tel_id)
+
     def finish(self):
         """called after all events are done"""
         self.log.info("Finishing output")
@@ -439,16 +468,7 @@ def _setup_writer(self):
             transform=tr_tel_list_to_mask,
         )
 
-        # currently the trigger info is used for the event time, but we dont'
-        # want the other bits of the trigger container in the pointing or other
-        # montitoring containers
-        writer.exclude("dl1/monitoring/subarray/pointing", "event_type")
-        writer.exclude("dl1/monitoring/subarray/pointing", "tels_with_trigger")
-        writer.exclude("dl1/monitoring/subarray/pointing", "n_trigger_pixels")
         writer.exclude("/dl1/event/telescope/trigger", "trigger_pixels")
-        writer.exclude("/dl1/monitoring/telescope/pointing/.*", "n_trigger_pixels")
-        writer.exclude("/dl1/monitoring/telescope/pointing/.*", "trigger_pixels")
-        writer.exclude("/dl1/monitoring/event/pointing/.*", "event_type")
         writer.exclude("/dl1/event/telescope/images/.*", "parameters")
         writer.exclude("/simulation/event/telescope/images/.*", "true_parameters")
 
@@ -506,15 +526,6 @@ def _setup_writer(self):
         self._writer = writer
         self.log.debug("Writer initialized: %s", self._writer)
 
-    def _write_subarray_pointing(self, event: ArrayEventContainer):
-        """store subarray pointing info in a monitoring table"""
-        pnt = event.pointing
-        current_pointing = (pnt.array_azimuth, pnt.array_altitude)
-        if current_pointing != self._last_pointing:
-            pnt.prefix = ""
-            self._writer.write("dl1/monitoring/subarray/pointing", [event.trigger, pnt])
-            self._last_pointing = current_pointing
-
     def _write_scheduling_and_observation_blocks(self):
         """write out SB and OB info"""
 
@@ -653,17 +664,6 @@ def _write_dl1_telescope_events(self, event: ArrayEventContainer):
         event
         """
 
-        # pointing info
-        for tel_id, pnt in event.pointing.tel.items():
-            current_pointing = (pnt.azimuth, pnt.altitude)
-            if current_pointing != self._last_pointing_tel[tel_id]:
-                pnt.prefix = ""
-                self._writer.write(
-                    f"dl1/monitoring/telescope/pointing/tel_{tel_id:03d}",
-                    [event.trigger.tel[tel_id], pnt],
-                )
-                self._last_pointing_tel[tel_id] = current_pointing
-
         for tel_id, dl1_camera in event.dl1.tel.items():
             tel_index = _get_tel_index(event, tel_id)
 

ctapipe/io/hdf5eventsource.py

@@ -1,12 +1,11 @@
 import logging
 from contextlib import ExitStack
-from functools import lru_cache
 from pathlib import Path
 from typing import Dict, Union
 
-import astropy.units as u
 import numpy as np
 import tables
+from astropy.table import QTable
 from astropy.utils.decorators import lazyproperty
 
 from ctapipe.atmosphere import AtmosphereDensityProfile
@@ -17,6 +16,7 @@
     CameraHillasParametersContainer,
     CameraTimingParametersContainer,
     ConcentrationContainer,
+    CoordinateFrameType,
     DispContainer,
     DL1CameraContainer,
     EventIndexContainer,
@@ -40,6 +40,7 @@
     SimulatedShowerContainer,
     SimulationConfigContainer,
     TelescopeImpactParameterContainer,
+    TelescopePointingContainer,
     TelescopeTriggerContainer,
     TelEventIndexContainer,
     TimingParametersContainer,
@@ -48,11 +49,10 @@
 from ..core import Container, Field
 from ..core.traits import UseEnum
 from ..instrument import SubarrayDescription
-from ..utils import IndexFinder
 from .astropy_helpers import read_table
 from .datalevels import DataLevel
 from .eventsource import EventSource
-from .hdf5tableio import HDF5TableReader, get_column_attrs
+from .hdf5tableio import HDF5TableReader
 from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP
 
 __all__ = ["HDF5EventSource"]
@@ -222,6 +222,14 @@
         self._obs_ids = tuple(
             self.file_.root.configuration.observation.observation_block.col("obs_id")
         )
+        pointing_key = "/configuration/telescope/pointing"
+        self._constant_telescope_pointing = {}
+        if pointing_key in self.file_.root:
+            for h5table in self.file_.root[pointing_key]._f_iter_nodes("Table"):
+                tel_id = int(h5table._v_name.partition("tel_")[-1])
+                table = QTable(read_table(self.file_, h5table._v_pathname), copy=False)
+                table.add_index("obs_id")
+                self._constant_telescope_pointing[tel_id] = table
 
     def __exit__(self, exc_type, exc_val, exc_tb):
         self.close()
@@ -574,15 +582,6 @@
             ignore_columns={"trigger_pixels"},
         )
 
-        array_pointing_finder = IndexFinder(
-            self.file_.root.dl1.monitoring.subarray.pointing.col("time")
-        )
-
-        tel_pointing_finder = {
-            table.name: IndexFinder(table.col("time"))
-            for table in self.file_.root.dl1.monitoring.telescope.pointing
-        }
-
         counter = 0
         for trigger, index in events:
             data = ArrayEventContainer(
@@ -631,8 +630,8 @@
             if len(data.trigger.tels_with_trigger) == 0:
                 continue
 
-            self._fill_array_pointing(data, array_pointing_finder)
-            self._fill_telescope_pointing(data, tel_pointing_finder)
+            self._fill_array_pointing(data)
+            self._fill_telescope_pointing(data)
 
             for tel_id in data.trigger.tel.keys():
                 key = f"tel_{tel_id:03d}"
@@ -701,67 +700,33 @@
             yield data
             counter += 1
 
-    @lazyproperty
-    def _subarray_pointing_attrs(self):
-        table = self.file_.root.dl1.monitoring.subarray.pointing
-        return get_column_attrs(table)
-
-    @lru_cache(maxsize=1000)
-    def _telescope_pointing_attrs(self, tel_id):
-        pointing_group = self.file_.root.dl1.monitoring.telescope.pointing
-        return get_column_attrs(pointing_group[f"tel_{tel_id:03d}"])
-
-    def _fill_array_pointing(self, data, array_pointing_finder):
+    def _fill_array_pointing(self, event):
         """
         Fill the array pointing information of a given event
         """
-        # Only unique pointings are stored, so reader.read() wont work as easily
-        # Thats why we match the pointings based on trigger time
-        closest_time_index = array_pointing_finder.closest(data.trigger.time.mjd)
-        table = self.file_.root.dl1.monitoring.subarray.pointing
-        array_pointing = table[closest_time_index]
-
-        data.pointing.array_azimuth = u.Quantity(
-            array_pointing["array_azimuth"],
-            self._subarray_pointing_attrs["array_azimuth"]["UNIT"],
-        )
-        data.pointing.array_altitude = u.Quantity(
-            array_pointing["array_altitude"],
-            self._subarray_pointing_attrs["array_altitude"]["UNIT"],
-        )
-        data.pointing.array_ra = u.Quantity(
-            array_pointing["array_ra"],
-            self._subarray_pointing_attrs["array_ra"]["UNIT"],
-        )
-        data.pointing.array_dec = u.Quantity(
-            array_pointing["array_dec"],
-            self._subarray_pointing_attrs["array_dec"]["UNIT"],
-        )
+        obs_id = event.index.obs_id
+        ob = self.observation_blocks[obs_id]
+        frame = ob.subarray_pointing_frame
+        if frame is CoordinateFrameType.ALTAZ:
+            event.pointing.array_azimuth = ob.subarray_pointing_lon
+            event.pointing.array_altitude = ob.subarray_pointing_lat
+        elif frame is CoordinateFrameType.ICRS:
+            event.pointing.array_ra = ob.subarray_pointing_lon
+            event.pointing.array_dec = ob.subarray_pointing_lat
+        else:
+            raise ValueError(f"Unsupported pointing frame: {frame}")
 
-    def _fill_telescope_pointing(self, data, tel_pointing_finder):
+    def _fill_telescope_pointing(self, event):
         """
         Fill the telescope pointing information of a given event
         """
-        # Same comments as to _fill_array_pointing apply
-        pointing_group = self.file_.root.dl1.monitoring.telescope.pointing
-        for tel_id in data.trigger.tel.keys():
-            key = f"tel_{tel_id:03d}"
-
-            if self.allowed_tels and tel_id not in self.allowed_tels:
+        for tel_id in event.trigger.tels_with_trigger:
+            tel_pointing = self._constant_telescope_pointing.get(tel_id)
+            if tel_pointing is None:
                 continue
 
-            tel_pointing_table = pointing_group[key]
-            closest_time_index = tel_pointing_finder[key].closest(
-                data.trigger.tel[tel_id].time.mjd
-            )
-
-            pointing_telescope = tel_pointing_table[closest_time_index]
-            attrs = self._telescope_pointing_attrs(tel_id)
-            data.pointing.tel[tel_id].azimuth = u.Quantity(
-                pointing_telescope["azimuth"],
-                attrs["azimuth"]["UNIT"],
-            )
-            data.pointing.tel[tel_id].altitude = u.Quantity(
-                pointing_telescope["altitude"],
-                attrs["altitude"]["UNIT"],
+            current = tel_pointing.loc[event.index.obs_id]
+            event.pointing.tel[tel_id] = TelescopePointingContainer(
+                altitude=current["telescope_pointing_altitude"],
+                azimuth=current["telescope_pointing_azimuth"],
             )

ctapipe/io/hdf5merger.py

@@ -32,6 +32,7 @@ class NodeType(enum.Enum):
     "/configuration/observation/scheduling_block": NodeType.TABLE,
     "/configuration/observation/observation_block": NodeType.TABLE,
     "/configuration/simulation/run": NodeType.TABLE,
+    "/configuration/telescope/pointing": NodeType.TEL_GROUP,
     "/simulation/service/shower_distribution": NodeType.TABLE,
     "/simulation/event/subarray/shower": NodeType.TABLE,
     "/simulation/event/telescope/impact": NodeType.TEL_GROUP,
@@ -280,6 +281,10 @@ def _append(self, other):
             if key in other.root:
                 self._append_table(other, other.root[key])
 
+        key = "/configuration/telescope/pointing"
+        if key in other.root:
+            self._append_table_group(other, other.root[key])
+
         # Simulation
         simulation_table_keys = [
             "/configuration/simulation/run",