Add a component to reduce waveform readout window

docs/changes/2780.feature.rst

@@ -0,0 +1,3 @@
+Add a new component ``ReadoutWindowReducer`` and options
+to enable it in ``ctapipe-process`` to study performance
+with shorter readout windows.

src/ctapipe/image/reducer.py

@@ -3,6 +3,7 @@
 """
 
 from abc import abstractmethod
+from copy import deepcopy
 
 import numpy as np
 
@@ -214,3 +215,39 @@ def select_pixels(self, waveforms, tel_id=None, selected_gain_channel=None):
             mask = dilate(camera_geom, mask)
 
         return mask
+
+
+class ReadoutWindowReducer(TelescopeComponent):
+    """
+    Reduce the readout window size of telescope using a fixed window.
+    """
+
+    window_start = IntTelescopeParameter(
+        default_value=None,
+        allow_none=True,
+        help="Start sample of readout window",
+    ).tag(config=True)
+
+    window_end = IntTelescopeParameter(
+        default_value=None,
+        allow_none=True,
+        help="Last sample of readout window (non-inclusive)",
+    ).tag(config=True)
+
+    def __init__(self, subarray, **kwargs):
+        # we mutate the subarray, make sure we do not modify it for someone else
+        super().__init__(deepcopy(subarray), **kwargs)
+
+        for tel_id, tel in self.subarray.tel.items():
+            start = self.window_start.tel[tel_id]
+            end = self.window_end.tel[tel_id]
+            n_samples = len(np.arange(tel.camera.readout.n_samples)[start:end])
+            tel.camera.readout.n_samples = n_samples
+
+    def __call__(self, event):
+        for container in (event.r0, event.r1, event.dl0):
+            for tel_id, tel_container in container.tel.items():
+                if tel_container.waveform is not None:
+                    start = self.window_start.tel[tel_id]
+                    end = self.window_end.tel[tel_id]
+                    tel_container.waveform = tel_container.waveform[..., start:end]

src/ctapipe/image/tests/test_reducer.py

@@ -1,3 +1,5 @@
+from copy import deepcopy
+
 import astropy.units as u
 import numpy as np
 import pytest
@@ -6,6 +8,7 @@
 
 from ctapipe.image.reducer import NullDataVolumeReducer, TailCutsDataVolumeReducer
 from ctapipe.instrument import SubarrayDescription
+from ctapipe.io import EventSource
 
 
 @pytest.fixture(scope="module")
@@ -92,3 +95,72 @@ def test_tailcuts_data_volume_reducer(subarray_lst):
 
     assert (reduced_waveforms != 0).sum() == (1 + 4 + 14) * n_samples
     assert_array_equal(expected_waveforms, reduced_waveforms)
+
+
+def test_readout_window_reducer(prod5_gamma_simtel_path):
+    from ctapipe.image.reducer import ReadoutWindowReducer
+
+    def check(subarray, event, event_reduced, event_no_op):
+        checked = 0
+        for dl in ("r0", "r1", "dl0"):
+            for tel_id in getattr(event, dl).tel.keys():
+                if str(subarray.tel[tel_id]).startswith("LST"):
+                    start = 10
+                    end = 20
+                    n_samples = end - start
+                elif "Nectar" in str(subarray.tel[tel_id]):
+                    start = 15
+                    end = 40
+                    n_samples = end - start
+                else:
+                    start = None
+                    end = None
+                    n_samples = subarray.tel[tel_id].camera.readout.n_samples
+
+                original_waveform = getattr(event, dl).tel[tel_id].waveform
+                reduced_waveform = getattr(event_reduced, dl).tel[tel_id].waveform
+                no_op_waveform = getattr(event_no_op, dl).tel[tel_id].waveform
+
+                assert reduced_waveform.ndim == 3
+                assert reduced_waveform.shape[-1] == n_samples
+                np.testing.assert_array_equal(original_waveform, no_op_waveform)
+                np.testing.assert_array_equal(
+                    original_waveform[..., start:end], reduced_waveform
+                )
+
+                checked += 1
+        return checked
+
+    with EventSource(prod5_gamma_simtel_path) as source:
+        reducer_no_op = ReadoutWindowReducer(source.subarray)
+        reducer = ReadoutWindowReducer(
+            source.subarray,
+            window_start=[
+                ("type", "*", None),
+                ("type", "LST*", 10),
+                ("type", "*Nectar*", 15),
+            ],
+            window_end=[
+                ("type", "*", None),
+                ("type", "LST*", 20),
+                ("type", "*Nectar*", 40),
+            ],
+        )
+
+        # make sure we didn't modify the original subarray
+        assert source.subarray.tel[1].camera.readout.n_samples == 40
+        assert reducer_no_op.subarray.tel[1].camera.readout.n_samples == 40
+        # new subarray should have reduced window
+        assert reducer.subarray.tel[1].camera.readout.n_samples == 10
+
+        n_checked = 0
+        for event in source:
+            event_no_op = deepcopy(event)
+            event_reduced = deepcopy(event)
+
+            reducer(event_reduced)
+            reducer_no_op(event_no_op)
+
+            n_checked += check(source.subarray, event, event_reduced, event_no_op)
+
+        assert n_checked > 0

src/ctapipe/tools/process.py

@@ -12,6 +12,7 @@
 from ..image import ImageCleaner, ImageModifier, ImageProcessor
 from ..image.extractor import ImageExtractor
 from ..image.muon import MuonProcessor
+from ..image.reducer import ReadoutWindowReducer
 from ..instrument import SoftwareTrigger
 from ..io import (
     DataLevel,
@@ -75,6 +76,11 @@ class ProcessorTool(Tool):
         default_value=False,
     ).tag(config=True)
 
+    reduce_readout_window = Bool(
+        default_value=False,
+        help="If True, use ReadoutWindowReducer on waveforms.",
+    ).tag(config=True)
+
     aliases = {
         ("i", "input"): "EventSource.input_url",
         ("o", "output"): "DataWriter.output_path",
@@ -152,6 +158,7 @@ class ProcessorTool(Tool):
             metadata.Instrument,
             metadata.Contact,
             SoftwareTrigger,
+            ReadoutWindowReducer,
         ]
         + classes_with_traits(EventSource)
         + classes_with_traits(ImageCleaner)
@@ -179,6 +186,15 @@ def setup(self):
             sys.exit(1)
 
         subarray = self.event_source.subarray
+
+        if self.reduce_readout_window:
+            self.readout_window_reducer = ReadoutWindowReducer(
+                subarray=subarray, parent=self
+            )
+            subarray = self.readout_window_reducer.subarray
+        else:
+            self.readout_window_reducer = None
+
         self.software_trigger = SoftwareTrigger(parent=self, subarray=subarray)
         self.calibrate = CameraCalibrator(parent=self, subarray=subarray)
         self.process_images = ImageProcessor(subarray=subarray, parent=self)
@@ -300,6 +316,9 @@ def start(self):
                 )
                 continue
 
+            if self.readout_window_reducer is not None:
+                self.readout_window_reducer(event)
+
             if self.should_calibrate:
                 self.calibrate(event)
 

src/ctapipe/tools/tests/test_process.py

@@ -595,3 +595,46 @@ def test_prod6_issues(tmp_path):
         images = loader.read_telescope_events([32], true_images=True)
         images.add_index("event_id")
         np.testing.assert_array_equal(images.loc[1664106]["true_image"], -1)
+
+
+def test_readout_window_reducer(tmp_path, provenance):
+    """Test ReadoutWindowReducer in process tool."""
+    output = tmp_path / "test_reduced_window.dl1.h5"
+
+    config = {
+        "ProcessorTool": {
+            "reduce_readout_window": True,
+            "ReadoutWindowReducer": {
+                "window_start": [
+                    ("type", "*", None),
+                    ("type", "LST*", 10),
+                ],
+                "window_end": [
+                    ("type", "*", None),
+                    ("type", "LST*", 20),
+                ],
+            },
+            "DataWriter": {
+                "write_r1_waveforms": True,
+                "write_dl1_images": True,
+                "write_dl1_parameters": False,
+            },
+        }
+    }
+    config_path = tmp_path / "config.json"
+    config_path.write_text(json.dumps(config))
+
+    provenance_log = tmp_path / "provenance.log"
+    input_path = get_dataset_path("gamma_prod5.simtel.zst")
+    run_tool(
+        ProcessorTool(),
+        argv=[
+            f"--config={config_path}",
+            f"--input={input_path}",
+            f"--output={output}",
+            f"--provenance-log={provenance_log}",
+            "--overwrite",
+        ],
+        cwd=tmp_path,
+        raises=True,
+    )