Change how DataWriter writes simulation pointing

Instead of writing the unique pointings with their timestamp,
which has no physical meaning in case of simulations,
we now write the first pointing per observation block into
the "/configuration/telescope/pointing/tel_nnn" tables.

ctapipe/conftest.py

@@ -326,6 +326,27 @@ def dl1_file(dl1_tmp_path, prod5_gamma_simtel_path):
         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,11 +1,10 @@
 import logging
-from functools import lru_cache
 from pathlib import Path
 from typing import Dict
 
-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
@@ -16,6 +15,7 @@
     CameraHillasParametersContainer,
     CameraTimingParametersContainer,
     ConcentrationContainer,
+    CoordinateFrameType,
     DispContainer,
     DL1CameraContainer,
     EventIndexContainer,
@@ -39,6 +39,7 @@
     SimulatedShowerContainer,
     SimulationConfigContainer,
     TelescopeImpactParameterContainer,
+    TelescopePointingContainer,
     TelescopeTriggerContainer,
     TelEventIndexContainer,
     TimingParametersContainer,
@@ -47,11 +48,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"]
@@ -217,6 +217,14 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs):
         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()
@@ -569,15 +577,6 @@ def _generator(self):
             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(
@@ -626,8 +625,8 @@ def _generator(self):
             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}"
@@ -696,67 +695,33 @@ def _generator(self):
             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",