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

@@ -483,6 +483,15 @@ class CoreParametersContainer(Container):
     psi = Field(nan * u.deg, "Image direction in the Tilted/Ground Frame", unit="deg")
 
 
+class TrueDispContainer(Container):
+    """
+    Container for true disp parameters
+    """
+    default_prefix = "true_disp"
+    norm = Field(nan * u.deg, "True disp norm", unit=u.deg)
+    sign = Field(nan, "True disp sign")
+
+
 class ImageParametersContainer(Container):
     """Collection of image parameters"""
 
@@ -520,6 +529,10 @@ class ImageParametersContainer(Container):
         default_factory=CoreParametersContainer,
         description="Image direction in the Tilted/Ground Frame",
     )
+    true_disp = Field(
+        default_factory=TrueDispContainer,
+        description="True disp parameters (only filled for true_parameters)",
+    )
 
 
 class DL1CameraContainer(Container):

src/ctapipe/core/traits.py

@@ -270,6 +270,8 @@ def _validate_str(self, obj, value):
             value = get_dataset_path(value.partition("dataset://")[2])
         elif url.scheme in ("", "file"):
             value = pathlib.Path(url.netloc, url.path)
+        elif os.name == "nt" and len(url.scheme) == 1:
+            value = pathlib.Path(value)
         else:
             self.error(obj, value)
 

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(

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,
+                            true_parameters.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,45 @@ def start(self):
             if self.should_compute_dl1:
                 self.process_images(event)
 
+                if event.simulation is not None:
+                    shower = event.simulation.shower
+                    for tel_id, dl1 in event.dl1.tel.items():
+                        if (
+                            tel_id not in event.simulation.tel
+                            or event.simulation.tel[tel_id].true_parameters is None
+                        ):
+                            continue
+
+                        true_param = event.simulation.tel[tel_id].true_parameters
+                        hillas = dl1.parameters.hillas
+
+                        if (
+                            hillas is not None
+                            and np.isfinite(hillas.fov_lat)
+                            and np.isfinite(hillas.fov_lon)
+                        ):
+                            pointing = event.monitoring.tel[tel_id].pointing
+                            # calculate true disp
+                            fov_lon, fov_lat = horizontal_to_telescope(
+                                alt=shower.alt,
+                                az=shower.az,
+                                pointing_alt=pointing.altitude,
+                                pointing_az=pointing.azimuth,
+                            )
+                            # numpy's trigonometric functions need radians
+                            psi = hillas.psi.to_value(u.rad)
+                            cog_lon = hillas.fov_lon
+                            cog_lat = hillas.fov_lat
+
+                            delta_lon = fov_lon - cog_lon
+                            delta_lat = fov_lat - cog_lat
+
+                            true_disp = (
+                                np.cos(psi) * delta_lon + np.sin(psi) * delta_lat
+                            )
+                            true_param.true_disp.norm = np.abs(true_disp)
+                            true_param.true_disp.sign = np.sign(true_disp.value)
+
             if self.should_compute_muon_parameters:
                 self.process_muons(event)
 

src/ctapipe/tools/tests/test_process_true_disp.py

@@ -0,0 +1,50 @@
+
+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"""
+    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_norm" in true_params.columns
+    assert "true_disp_sign" in true_params.columns
+
+    # Check that we have some valid values
+    assert np.count_nonzero(np.isfinite(true_params["true_disp_norm"])) > 0
+    assert np.count_nonzero(np.isfinite(true_params["true_disp_sign"])) > 0