feat: Implement true_disp calculation in ctapipe-process + test

docs/changes/2950.feature.rst

@@ -0,0 +1 @@
+Implement calculation of ``true_disp`` parameters (norm and sign) in ``ctapipe-process`` for simulated events.

src/ctapipe/containers.py

@@ -788,6 +788,12 @@ class SimulatedCameraContainer(Container):
         description="true impact parameter",
     )
 
+    true_disp = Field(
+        nan * u.deg,
+        description="True disp parameter",
+        unit=u.deg,
+    )
+
 
 class SimulatedEventContainer(Container):
     shower = Field(

src/ctapipe/image/image_processor.py

@@ -251,8 +251,10 @@ def _process_telescope_event(self, event):
                     peak_time=None,  # true image from simulation has no peak time
                     default=DEFAULT_TRUE_IMAGE_PARAMETERS,
                 )
+                from ctapipe.core import Container
+
                 for container in sim_camera.true_parameters.values():
-                    if not container.prefix.startswith("true_"):
+                    if isinstance(container, Container) and not container.prefix.startswith("true_"):
                         container.prefix = f"true_{container.prefix}"
 
                 self.log.debug(
@@ -261,3 +263,34 @@ def _process_telescope_event(self, event):
                         recursive=True
                     ),
                 )
+
+                if (
+                    dl1_camera.parameters.hillas is not None
+                    and np.isfinite(dl1_camera.parameters.hillas.fov_lat)
+                    and np.isfinite(dl1_camera.parameters.hillas.fov_lon)
+                ):
+                    pointing = event.monitoring.tel[tel_id].pointing
+                    shower = event.simulation.shower
+
+                    from ctapipe.reco.preprocessing import (
+                        horizontal_to_telescope,
+                        calculate_true_disp,
+                    )
+
+                    fov_lon, fov_lat = horizontal_to_telescope(
+                        alt=shower.alt,
+                        az=shower.az,
+                        pointing_alt=pointing.altitude,
+                        pointing_az=pointing.azimuth,
+                    )
+
+                    hillas = dl1_camera.parameters.hillas
+
+                    sim_camera.true_disp = calculate_true_disp(
+                        fov_lon=fov_lon,
+                        fov_lat=fov_lat,
+                        hillas_psi=hillas.psi,
+                        hillas_lon=hillas.fov_lon,
+                        hillas_lat=hillas.fov_lat,
+                    )
+

src/ctapipe/io/datawriter.py

@@ -660,6 +660,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/reco/preprocessing.py

@@ -29,6 +29,7 @@
     "table_to_X",
     "horizontal_to_telescope",
     "telescope_to_horizontal",
+    "calculate_true_disp",
 ]
 
 
@@ -154,3 +155,35 @@ def telescope_to_horizontal(lon, lat, pointing_alt, pointing_az):
         horizontal_coord = tel_coord.transform_to(AltAz())
 
     return horizontal_coord.alt.to(u.deg), horizontal_coord.az.to(u.deg)
+
+
+@u.quantity_input(
+    fov_lon=u.deg, fov_lat=u.deg, hillas_psi=u.rad, hillas_lon=u.deg, hillas_lat=u.deg
+)
+def calculate_true_disp(fov_lon, fov_lat, hillas_psi, hillas_lon, hillas_lat):
+    """
+    Calculate the true disp parameter (distance between hillas cog and source position)
+
+    Parameters
+    ----------
+    fov_lon : u.Quantity
+        Source longitude in telescope frame
+    fov_lat : u.Quantity
+        Source latitude in telescope frame
+    hillas_psi : u.Quantity
+        Hillas psi parameter (angle between major axis and x-axis)
+    hillas_lon : u.Quantity
+        Hillas center longitude
+    hillas_lat : u.Quantity
+        Hillas center latitude
+
+    Returns
+    -------
+    true_disp : u.Quantity
+        The true disp parameter (signed distance)
+    """
+    delta_lon = fov_lon - hillas_lon
+    delta_lat = fov_lat - hillas_lat
+
+    true_disp = np.cos(hillas_psi) * delta_lon + np.sin(hillas_psi) * delta_lat
+    return true_disp

src/ctapipe/tools/process.py

@@ -5,6 +5,8 @@
 # pylint: disable=W0201
 import sys
 
+import astropy.units as u
+import numpy as np
 from tqdm.auto import tqdm
 
 from ..calib import CameraCalibrator, GainSelector
@@ -30,6 +32,7 @@
     DL2_EVENT_STATISTICS_GROUP,
 )
 from ..reco import Reconstructor, ShowerProcessor
+from ..reco.preprocessing import horizontal_to_telescope
 from ..utils import EventTypeFilter
 
 COMPATIBLE_DATALEVELS = [
@@ -362,6 +365,8 @@ def start(self):
             if self.should_compute_dl1:
                 self.process_images(event)
 
+
+
             if self.should_compute_muon_parameters:
                 self.process_muons(event)
 

src/ctapipe/tools/tests/test_process_true_disp.py

@@ -0,0 +1,49 @@
+
+import numpy as np
+import pandas as pd
+import pytest
+import tables
+
+from ctapipe.core import run_tool
+from ctapipe.tools.process import ProcessorTool
+from ctapipe.utils import get_dataset_path, resource_file
+from ctapipe.io.hdf5dataformat import SIMULATION_PARAMETERS_GROUP
+
+def test_true_disp_calculation(tmp_path, dl1_image_file):
+    """check true_disp calculation in ctapipe-process"""
+    print("DEBUG: Starting test_true_disp_calculation", flush=True)
+    config = resource_file("stage1_config.json")
+
+    output_file = tmp_path / "true_disp_test.dl1.h5"
+
+    run_tool(
+        ProcessorTool(),
+        argv=[
+            f"--config={config}",
+            f"--input={dl1_image_file}",
+            f"--output={output_file}",
+            "--write-parameters",
+            "--overwrite",
+        ],
+        cwd=tmp_path,
+        raises=True,
+    )
+
+    # check if fields exist and are not all NaN
+    # We need to find a telescope that has events
+
+    with tables.open_file(output_file, mode="r") as testfile:
+        # Check if true parameters group exists for at least one telescope
+        sim_params_group = testfile.root.simulation.event.telescope.parameters
+        assert len(sim_params_group._v_children) > 0
+
+        first_tel_group = list(sim_params_group._v_children.keys())[0]
+
+    true_params = pd.read_hdf(
+        output_file, f"{SIMULATION_PARAMETERS_GROUP}/{first_tel_group}"
+    )
+
+    assert "true_disp" in true_params.columns
+
+    # Check that we have some valid values
+    assert np.count_nonzero(np.isfinite(true_params["true_disp"])) > 0