Read pixel time shift if present and supported by ctapipe

src/ctapipe_io_zfits/conftest.py

@@ -81,6 +81,27 @@ def get_module_and_pixel_id_map(n_modules, n_pixels_module, missing_modules=None
         },
         id="standard",
     ),
+    pytest.param(
+        {
+            "obs_start": Time("2025-02-04T20:45:31"),
+            "sb_creator_id": 2,
+            "sb_id": 124,
+            "obs_id": 790,
+            "pixel_time_shift": True,
+        },
+        id="pixel_time_shift",
+    ),
+    pytest.param(
+        {
+            "obs_start": Time("2025-02-04T20:45:31"),
+            "sb_creator_id": 2,
+            "sb_id": 124,
+            "obs_id": 791,
+            "dvr": True,
+            "pixel_time_shift": True,
+        },
+        id="dvred_pixel_time_shift",
+    ),
     pytest.param(
         {
             "missing_modules": [50, 200],
@@ -156,6 +177,12 @@ def dummy_dl0(dl0_base, request):
         n_modules=265, n_pixels_module=7, missing_modules=missing_modules
     )
     n_pixels = len(pixel_id_map)
+    pixel_stored = np.ones(n_pixels, dtype=bool)
+    if config.get("dvr", False):
+        # Keep every second pixel, so that the event contains fewer pixels than
+        # the camera configuration and exercises DVR reordering.
+        pixel_stored[::2] = False
+    n_pixels_stored = pixel_stored.sum()
 
     camera_configuration = DL0_Telescope.CameraConfiguration(
         tel_id=1,
@@ -245,7 +272,18 @@ def convert_waveform(waveform):
             # TODO: randomize event to test actually parsing it
 
             # TODO: fill actual signal into waveform, not just 0
-            waveform = rng.normal(0.0, 1.0, size=(1, n_pixels, 40)).astype(np.float32)
+            waveform = rng.normal(0.0, 1.0, size=(1, n_pixels_stored, 40)).astype(
+                np.float32
+            )
+
+            additional_fields = {}
+            if config.get("pixel_time_shift", False):
+                if config.get("dvr", False):
+                    time_shift = np.arange(1, n_pixels_stored + 1, dtype=np.int16)
+                else:
+                    time_shift = rng.normal(0, 0.5, size=(n_pixels_stored,))
+                    time_shift = np.round(100 * time_shift).astype(np.int16)
+                additional_fields["pixel_time_shift"] = numpy_to_any_array(time_shift)
 
             lst_event_files[sdh_id].write_message(
                 DL0_Telescope.Event(
@@ -256,12 +294,13 @@ def convert_waveform(waveform):
                     event_time_qns=int(time_qns),
                     # identified as signal, low gain stored, high gain stored
                     pixel_status=numpy_to_any_array(
-                        np.full(n_pixels, 0b00001101, dtype=np.uint8)
+                        np.where(pixel_stored, 0b00001101, 0).astype(np.uint8)
                     ),
                     waveform=numpy_to_any_array(convert_waveform(waveform)),
                     num_channels=1,
                     num_samples=40,
-                    num_pixels_survived=n_pixels,
+                    num_pixels_survived=n_pixels_stored,
+                    **additional_fields,
                 )
             )
             events_written[sdh_id] += 1

src/ctapipe_io_zfits/dl0.py

@@ -4,6 +4,7 @@
 from contextlib import ExitStack
 
 import numpy as np
+from ctapipe import __version__ as ctapipe_version
 from ctapipe.containers import (
     ArrayEventContainer,
     CameraCalibrationContainer,
@@ -20,6 +21,7 @@
 from ctapipe.instrument import SubarrayDescription
 from ctapipe.io import DataLevel, EventSource
 from ctapipe.io.simteleventsource import GainChannel
+from packaging.version import Version
 from protozfits import File
 
 from .instrument import build_subarray_description, get_array_elements_by_id
@@ -33,7 +35,11 @@
 
 log = logging.getLogger(__name__)
 
+CTAPIPE_VERSION = Version(ctapipe_version)
+CTAPIPE_GE_0_31 = CTAPIPE_VERSION >= Version("0.31.0a0")
+
 ARRAY_ELEMENTS = get_array_elements_by_id()
+TEN_PS_TO_NS = np.float32(0.01)
 
 
 def _is_compatible(input_url, extname, allowed_protos):
@@ -96,15 +102,14 @@ def _fill_dl0_container(
 
     pixel_status = tel_event.pixel_status
     # FIXME: seems ACADA doesn't set pixels to "stored" when no DVR is applied
-    if n_pixels_stored == camera_config.num_pixels and np.all(
-        PixelStatus.get_dvr_status(pixel_status) == 0
-    ):
+    all_pixels_stored = n_pixels_stored == camera_config.num_pixels
+    if all_pixels_stored and np.all(PixelStatus.get_dvr_status(pixel_status) == 0):
         pixel_status = pixel_status | PixelStatus.DVR_1
 
     pixel_stored = PixelStatus.get_dvr_status(pixel_status) != 0
-    n_pixels_nominal = camera_geometry.n_pixels
 
     # fill not readout pixels with 0, reorder pixels
+    n_pixels_nominal = camera_geometry.n_pixels
     waveform = np.full(
         (n_channels, n_pixels_nominal, n_samples), dvr_fill_value, dtype=np.float32
     )
@@ -131,6 +136,17 @@ def _fill_dl0_container(
     else:
         selected_gain_channel = None
 
+    extra_fields = {}
+    if CTAPIPE_GE_0_31 and tel_event.pixel_time_shift is not None:
+        # pixel_time_shift is stored as int16, in 10 ps increments. Convert to ns.
+        pixel_time_shift = tel_event.pixel_time_shift.astype(np.float32) * TEN_PS_TO_NS
+        pixel_time_shift = pixel_time_shift.reshape((n_channels, n_pixels_stored))
+        pixel_time_shift_reordered = np.zeros((n_channels, n_pixels_nominal))
+        pixel_time_shift_reordered[..., camera_config.pixel_id_map[pixel_stored]] = (
+            pixel_time_shift
+        )
+        extra_fields["pixel_time_shift"] = pixel_time_shift_reordered
+
     return DL0CameraContainer(
         pixel_status=pixel_status_reordered,
         event_type=EventType(int(tel_event.event_type)),
@@ -141,6 +157,7 @@ def _fill_dl0_container(
         ),
         waveform=waveform,
         first_cell_id=tel_event.first_cell_id,
+        **extra_fields,
     )
 
 

src/ctapipe_io_zfits/tests/test_dl0.py

@@ -6,6 +6,8 @@
 from ctapipe.io import EventSource, TableLoader
 from ctapipe.tools.process import ProcessorTool
 
+from ctapipe_io_zfits.dl0 import CTAPIPE_GE_0_31
+
 
 def test_is_compatible(dummy_dl0):
     from ctapipe_io_zfits import ProtozfitsDL0EventSource
@@ -46,6 +48,17 @@ def test_subarray_events(dummy_dl0):
             n_read += 1
             time = time + 0.001 * u.s
 
+            if CTAPIPE_GE_0_31 and dummy_dl0.get("pixel_time_shift"):
+                pixel_time_shift = array_event.dl0.tel[1].pixel_time_shift
+                assert pixel_time_shift is not None
+
+                if dummy_dl0.get("dvr"):
+                    pixel_stored = array_event.dl0.tel[1].pixel_status != 0
+                    np.testing.assert_array_equal(
+                        pixel_time_shift[0, ~pixel_stored], 0.0
+                    )
+                    assert np.all(pixel_time_shift[0, pixel_stored] > 0.0)
+
         assert n_read == 100