Calculate and store true_disp per-telescope when running ctapipe-process

docs/changes/3003.feature.rst

@@ -0,0 +1,5 @@
+Added calculation and storage of ``true_disp_norm`` and ``true_disp_sign`` per-telescope
+when running ``ctapipe-process``. These values are computed in ``ImageProcessor``
+(telescope frame only) and written to ``simulation/event/telescope/parameters``.
+A new ``TrueDispContainer`` is added to ``SimulatedCameraContainer``, and
+``HDF5EventSource`` reads these columns back with NaN defaults for older files.

docs/conf.py

@@ -131,6 +131,8 @@ def add_reference_type(prefix, objs):
         "traitlets.traitlets.ObserveHandler",
         "traitlets.traitlets.T",
         "traitlets.traitlets.G",
+        "traitlets.traitlets.K",
+        "traitlets.traitlets.V",
         "Sentinel",
         "ObserveHandler",
         "dict[K, V]",

docs/user-guide/data_format/index.rst

@@ -128,8 +128,10 @@ Simulation Data Format
       - simulated camera images
       - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.SimulatedCameraContainer`
     * - ``/simulation/event/telescope/parameters/tel_{TEL_ID:03d}``
-      - Parameters derived form the simulated camera images
-      - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ImageParametersContainer`
+      - Parameters derived from the simulated camera images, including true Hillas parameters
+        and true disp values (:py:class:`~ctapipe.containers.TrueDispContainer`) for training
+        and evaluating angular reconstruction algorithms
+      - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ImageParametersContainer`, :py:class:`~ctapipe.containers.TrueDispContainer`
     * - ``/simulation/service/shower_distribution``
       - simulated shower distribution histograms
       - :py:class:`~ctapipe.containers.SimulatedShowerDistribution`

src/ctapipe/containers.py

@@ -50,6 +50,7 @@
     "ReconstructedEnergyContainer",
     "ReconstructedGeometryContainer",
     "DispContainer",
+    "TrueDispContainer",
     "SimulatedCameraContainer",
     "SimulatedShowerContainer",
     "SimulatedShowerDistribution",
@@ -764,6 +765,26 @@ class TelescopeImpactParameterContainer(Container):
 
 
 # Simulation containers
+class TrueDispContainer(Container):
+    """
+    Container for the true (simulated) disp values, used for debugging angular
+    reconstruction and for training the DispReconstructor.
+    """
+
+    default_prefix = "true_disp"
+
+    norm = Field(
+        nan * u.deg,
+        "True norm of the disp parameter (distance from shower COG to true source position)",
+        unit=u.deg,
+    )
+    sign = Field(
+        np.float32(nan),
+        "True sign of the disp parameter (+1, -1, or 0 in degenerate case where "
+        "true source position coincides exactly with the shower CoG)",
+    )
+
+
 class SimulatedShowerContainer(Container):
     default_prefix = "true"
     energy = Field(nan * u.TeV, "Simulated Energy", unit=u.TeV)
@@ -811,6 +832,12 @@ class SimulatedCameraContainer(Container):
         type=ImageParametersContainer,
     )
 
+    true_disp = Field(
+        default_factory=TrueDispContainer,
+        description="True disp parameters (norm and sign) derived from the true shower direction",
+        type=TrueDispContainer,
+    )
+
     impact = Field(
         default_factory=TelescopeImpactParameterContainer,
         description="true impact parameter",

src/ctapipe/image/image_processor.py

@@ -2,11 +2,13 @@
 High level image processing  (ImageProcessor Component)
 """
 
+import warnings
 from copy import deepcopy
 
+import astropy.units as u
 import numpy as np
 
-from ctapipe.coordinates import TelescopeFrame
+from ctapipe.coordinates import MissingFrameAttributeWarning, TelescopeFrame
 
 from ..containers import (
     ArrayEventContainer,
@@ -16,6 +18,7 @@
     IntensityStatisticsContainer,
     PeakTimeStatisticsContainer,
     TimingParametersContainer,
+    TrueDispContainer,
 )
 from ..core import QualityQuery, TelescopeComponent
 from ..core.traits import Bool, BoolTelescopeParameter, ComponentName, List
@@ -213,6 +216,72 @@ def _parameterize_image(
         # parameterization
         return default
 
+    def _calculate_true_disp(self, event, sim_camera):
+        """
+        Calculate the true disp parameters (norm and sign) for a simulated telescope event.
+
+        This method should only be called when ``use_telescope_frame=True``, as it
+        requires the hillas parameters in the telescope frame (fov_lon, fov_lat, psi).
+
+        Parameters
+        ----------
+        event : ArrayEventContainer
+            The event container with simulation and monitoring data.
+        sim_camera : SimulatedCameraContainer
+            The simulated camera container with true_parameters already set.
+
+        Returns
+        -------
+        TrueDispContainer
+            Container with true_disp_norm and true_disp_sign values, or default
+            NaN-filled container if the calculation is not possible (missing shower
+            info or pointing data).
+        """
+        from astropy.coordinates import AltAz
+
+        shower = event.simulation.shower
+        if shower is None:
+            return TrueDispContainer()
+
+        pointing_alt = event.monitoring.pointing.array_altitude
+        pointing_az = event.monitoring.pointing.array_azimuth
+
+        if np.isnan(pointing_alt.value) or np.isnan(pointing_az.value):
+            return TrueDispContainer()
+
+        # Clamp altitude to [-90, 90] deg to avoid floating-point rounding errors
+        # that can push values slightly outside the valid range for astropy AltAz
+        pointing_alt = np.clip(pointing_alt.to_value(u.deg), -90.0, 90.0) * u.deg
+        shower_alt = np.clip(shower.alt.to_value(u.deg), -90.0, 90.0) * u.deg
+
+        hillas = sim_camera.true_parameters.hillas
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", MissingFrameAttributeWarning)
+            horizontal_coord = AltAz(alt=shower_alt, az=shower.az)
+            pointing = AltAz(alt=pointing_alt, az=pointing_az)
+            tel_frame = TelescopeFrame(telescope_pointing=pointing)
+            tel_coord = horizontal_coord.transform_to(tel_frame)
+
+        fov_lon = tel_coord.fov_lon.to(u.deg)
+        fov_lat = tel_coord.fov_lat.to(u.deg)
+
+        psi = hillas.psi.to_value(u.rad)
+        cog_lon = hillas.fov_lon.to(u.deg)
+        cog_lat = hillas.fov_lat.to(u.deg)
+
+        delta_lon = fov_lon - cog_lon
+        delta_lat = fov_lat - cog_lat
+
+        # Project displacement vector onto the shower major axis (defined by psi)
+        # to determine on which end of the shower axis the true source lies
+        true_disp_projected = np.cos(psi) * delta_lon + np.sin(psi) * delta_lat
+        true_sign = np.float32(np.sign(true_disp_projected.value))
+        # Norm is the full Euclidean distance from COG to true source position
+        true_norm = np.sqrt(delta_lon**2 + delta_lat**2)
+
+        return TrueDispContainer(norm=true_norm, sign=true_sign)
+
     def _process_telescope_event(self, event):
         """
         Loop over telescopes and process the calibrated images into parameters
@@ -255,6 +324,9 @@ def _process_telescope_event(self, event):
                     if not container.prefix.startswith("true_"):
                         container.prefix = f"true_{container.prefix}"
 
+                if self.use_telescope_frame:
+                    sim_camera.true_disp = self._calculate_true_disp(event, sim_camera)
+
                 self.log.debug(
                     "sim params: %s",
                     event.simulation.tel[tel_id].true_parameters.as_dict(

src/ctapipe/image/tests/test_image_processor.py

@@ -7,14 +7,99 @@
 import astropy.units as u
 import numpy as np
 import pytest
+from astropy.coordinates import EarthLocation
 from numpy import isfinite
 
 from ctapipe.calib import CameraCalibrator
 from ctapipe.containers import (
+    ArrayEventContainer,
     CameraHillasParametersContainer,
     CameraTimingParametersContainer,
+    DL1CameraContainer,
+    SimulatedCameraContainer,
+    SimulatedEventContainer,
+    SimulatedShowerContainer,
 )
+from ctapipe.coordinates import CameraFrame
 from ctapipe.image import ImageCleaner, ImageProcessor
+from ctapipe.instrument import CameraGeometry, SubarrayDescription, TelescopeDescription
+from ctapipe.instrument.camera import CameraDescription, CameraReadout
+from ctapipe.instrument.optics import OpticsDescription, ReflectorShape, SizeType
+
+
+@pytest.fixture(scope="module")
+def synthetic_subarray():
+    """
+    Create a fully synthetic subarray (no network access needed).
+    Uses a rectangular camera with a 1m focal length.
+    """
+    focal_length = 1.0 * u.m
+    geom = CameraGeometry.make_rectangular(10, 10)
+    geom.frame = CameraFrame(focal_length=focal_length)
+    n_pix = geom.n_pixels
+
+    readout = CameraReadout(
+        name="TestCam",
+        sampling_rate=1 * u.GHz,
+        reference_pulse_shape=np.ones((1, 3)),
+        reference_pulse_sample_width=1 * u.ns,
+        n_channels=1,
+        n_pixels=n_pix,
+        n_samples=50,
+    )
+    cam_desc = CameraDescription("TestCam", geom, readout)
+    optics = OpticsDescription(
+        name="testoptics",
+        size_type=SizeType.MST,
+        n_mirrors=1,
+        n_mirror_tiles=1,
+        mirror_area=1.0 * u.m**2,
+        equivalent_focal_length=focal_length,
+        effective_focal_length=focal_length,
+        reflector_shape=ReflectorShape.PARABOLIC,
+    )
+    tel_desc = TelescopeDescription(name="TestTel", optics=optics, camera=cam_desc)
+    reference_location = EarthLocation(
+        lon=-17 * u.deg, lat=28 * u.deg, height=2200 * u.m
+    )
+    return SubarrayDescription(
+        name="test",
+        tel_positions={1: np.array([0.0, 0.0, 0.0]) * u.m},
+        tel_descriptions={1: tel_desc},
+        reference_location=reference_location,
+    )
+
+
+@pytest.fixture
+def synthetic_event_with_simulation(synthetic_subarray):
+    """
+    Create a synthetic ArrayEventContainer with simulation data for tel_id=1.
+    The image has a simple signal in the central pixels so that Hillas parameters
+    can be computed.
+    """
+    n_pix = synthetic_subarray.tel[1].camera.geometry.n_pixels
+    image = np.zeros(n_pix)
+    # Put signal in the central pixels (rows 3-7, cols 3-7 of the 10x10 grid)
+    # Signal falls off in one direction to give a clear major axis
+    for i in range(5):
+        image[30 + i * 10 + 3 : 30 + i * 10 + 7] = 100.0 * np.exp(-((i - 2) ** 2) / 2)
+
+    true_image = (image * 0.5).astype(np.int32)
+
+    event = ArrayEventContainer()
+    event.monitoring.pointing.array_altitude = 70.0 * u.deg
+    event.monitoring.pointing.array_azimuth = 0.0 * u.deg
+
+    event.simulation = SimulatedEventContainer()
+    event.simulation.shower = SimulatedShowerContainer(
+        alt=70.5 * u.deg,
+        az=0.5 * u.deg,
+    )
+
+    event.dl1.tel[1] = DL1CameraContainer(image=image, peak_time=None)
+    event.monitoring.tel[1]  # initialize monitoring
+    event.simulation.tel[1] = SimulatedCameraContainer(true_image=true_image)
+    return event
 
 
 @pytest.mark.parametrize("cleaner", ImageCleaner.non_abstract_subclasses().values())
@@ -94,3 +179,48 @@ def test_image_processor_camera_frame(cleaner, example_event, example_subarray):
         assert isinstance(dl1.parameters.timing, CameraTimingParametersContainer)
         assert np.isnan(dl1.parameters.hillas.length.value)
         assert dl1.parameters.hillas.length.unit == u.m
+
+
+def test_true_disp_calculation(synthetic_subarray, synthetic_event_with_simulation):
+    """Test that true_disp_norm and true_disp_sign are calculated for simulation events."""
+    process_images = ImageProcessor(subarray=synthetic_subarray)
+    process_images(synthetic_event_with_simulation)
+
+    sim_camera = synthetic_event_with_simulation.simulation.tel[1]
+    true_disp = sim_camera.true_disp
+
+    # Hillas parameters should have been computed for the non-trivial image
+    if isfinite(sim_camera.true_parameters.hillas.fov_lon.value):
+        assert isfinite(true_disp.norm.value)
+        assert true_disp.norm >= 0 * u.deg
+        assert true_disp.norm.unit == u.deg
+        assert true_disp.sign in {-1.0, 0.0, 1.0}
+
+
+def test_true_disp_requires_telescope_frame(
+    synthetic_subarray, synthetic_event_with_simulation
+):
+    """Test that true_disp remains NaN when use_telescope_frame=False."""
+    process_images = ImageProcessor(
+        subarray=synthetic_subarray, use_telescope_frame=False
+    )
+    process_images(synthetic_event_with_simulation)
+
+    sim_camera = synthetic_event_with_simulation.simulation.tel[1]
+    # When not in telescope frame, true_disp should remain at defaults (NaN)
+    assert np.isnan(sim_camera.true_disp.norm.value)
+    assert np.isnan(sim_camera.true_disp.sign)
+
+
+def test_true_disp_no_simulation(synthetic_subarray, synthetic_event_with_simulation):
+    """Test that _calculate_true_disp handles missing simulation shower gracefully."""
+    event = deepcopy(synthetic_event_with_simulation)
+    event.simulation.shower = None
+
+    process_images = ImageProcessor(subarray=synthetic_subarray)
+    process_images(event)
+
+    sim_camera = event.simulation.tel[1]
+    # Without shower, true_disp should remain at defaults (NaN)
+    assert np.isnan(sim_camera.true_disp.norm.value)
+    assert np.isnan(sim_camera.true_disp.sign)

src/ctapipe/io/datawriter.py

@@ -662,6 +662,7 @@ def _write_dl1_telescope_events(self, event: ArrayEventContainer):
                             true_parameters.concentration,
                             true_parameters.morphology,
                             true_parameters.intensity_statistics,
+                            event.simulation.tel[tel_id].true_disp,
                         ],
                     )
 

src/ctapipe/io/hdf5eventsource.py

@@ -44,6 +44,7 @@
     TelEventIndexContainer,
     TimingParametersContainer,
     TriggerContainer,
+    TrueDispContainer,
 )
 from ..core import Container, Field, Provenance
 from ..core.traits import UseEnum
@@ -663,13 +664,15 @@ def _init_dl1_true_parameter_readers(self):
                     ConcentrationContainer,
                     MorphologyContainer,
                     IntensityStatisticsContainer,
+                    TrueDispContainer,
                 ],
                 prefixes=[
                     f"true_{hillas_prefix}",
                     "true_leakage",
                     "true_concentration",
                     "true_morphology",
                     "true_intensity",
+                    "true_disp",
                 ],
             )
             for table in self.file_.root.dl1.event.telescope.parameters
@@ -1002,6 +1005,7 @@ def _fill_parameters_for_tel(
             morphology=simulated_params[3],
             intensity_statistics=simulated_params[4],
         )
+        simulated.true_disp = simulated_params[5]
 
     def _fill_muons_for_tel(self, data, tel_id, key, muon_readers):
         if not self.has_muon_parameters: