Repository: cta-observatory/ctapipe_io_lst

Issue:
ctapipe_io_lst not compatible with ctapipe v0.21.1
The current ctapipe_io_lst is not compatible with the latest version of ctapipe (to write v0.6.0 dl1 files) and leads to the following error.

Traceback (most recent call last):
  File "/home/bastien.lacave/miniforge3/envs/gammapy-1.2/lib/python3.10/site-packages/ctapipe/core/tool.py", line 426, in run
    self.start()
  File "/home/bastien.lacave/miniforge3/envs/gammapy-1.2/lib/python3.10/site-packages/ctapipe/tools/process.py", line 314, in start
    self.calibrate(event)
  File "/home/bastien.lacave/miniforge3/envs/gammapy-1.2/lib/python3.10/site-packages/ctapipe/calib/camera/calibrator.py", line 328, in __call__
    self._calibrate_dl0(event, tel_id)
  File "/home/bastien.lacave/miniforge3/envs/gammapy-1.2/lib/python3.10/site-packages/ctapipe/calib/camera/calibrator.py", line 202, in _calibrate_dl0
    dl0_waveform[:, ~signal_pixels] = 0
IndexError: boolean index did not match indexed array along dimension 1; dimension is 36 but corresponding boolean dimension is 1855
2024-06-17 12:12:40,147 [1;33mWARNING[0m [ctapipe.ctapipe-process.DataWriter] (datawriter.finish): No events have been written to the output file

Discussion:

Patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 4396dff..91142d3 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -20,6 +20,10 @@ jobs:
       matrix:
         python-version: ["3.9", "3.10", "3.11"]
         ctapipe-version: ["0.19.3", "0.20.0"]
+        include:
+          # ctapipe 0.21 requires >= 3.10
+          - python-version: "3.12"
+            ctapipe-version: "0.21.2"
 
     defaults:
       run:
diff --git a/environment.yml b/environment.yml
index 495ec70..f9df547 100644
--- a/environment.yml
+++ b/environment.yml
@@ -1,11 +1,10 @@
 name: lstio
 channels:
   - conda-forge
-  - default
 dependencies:
-  - astropy=5.2
+  - astropy >=5.2,<7
   - python=3.11 # nail the python version, so conda does not try upgrading / dowgrading
-  - ctapipe=0.19
+  - ctapipe
   - protozfits=2.4
   - eventio
   - corsikaio
diff --git a/setup.cfg b/setup.cfg
index 21b9e90..66d3406 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -30,11 +30,11 @@ package_dir =
 python_requires = >=3.9
 zip_safe = False
 install_requires=
-    astropy~=5.2
-    ctapipe >=0.19.0,<0.21.0a0
-    protozfits~=2.4
-    numpy>=1.20
-    scipy<1.14
+    astropy >=5.2,<7.0.0a0
+    ctapipe >=0.19.0,<0.22.0a0
+    protozfits ~=2.4
+    numpy >=1.20
+    scipy <1.14
 
 [options.package_data]
 * = resources/*
diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py
index 6c33ab3..a60128d 100644
--- a/src/ctapipe_io_lst/__init__.py
+++ b/src/ctapipe_io_lst/__init__.py
@@ -7,7 +7,6 @@ from ctapipe.instrument.subarray import EarthLocation
 import logging
 import numpy as np
 from astropy import units as u
-from ctapipe import __version__ as ctapipe_version
 from ctapipe.core import Provenance
 from ctapipe.instrument import (
     ReflectorShape,
@@ -53,6 +52,7 @@ from .constants import (
 )
 
 from .evb_preprocessing import get_processings_for_trigger_bits, EVBPreprocessingFlag
+from .compat import CTAPIPE_GE_0_20, CTAPIPE_GE_0_21
 
 
 __all__ = [
@@ -61,10 +61,6 @@ __all__ = [
 ]
 
 
-CTAPIPE_VERSION = tuple(int(v) for v in ctapipe_version.split(".")[:3])
-CTAPIPE_0_20 = CTAPIPE_VERSION >= (0, 20)
-
-
 log = logging.getLogger(__name__)
 
 
@@ -523,7 +519,7 @@ class LSTEventSource(EventSource):
             selected_gain_channel=selected_gain_channel,
         )
 
-        if CTAPIPE_0_20:
+        if CTAPIPE_GE_0_20:
             r1.pixel_status = pixel_status
             r1.event_type = EventType(zfits_event.event_type)
             r1.event_time = trigger.time
@@ -688,6 +684,16 @@ class LSTEventSource(EventSource):
             # dl1 and drs4 timeshift needs to be filled always
             self.r0_r1_calibrator.fill_time_correction(array_event)
 
+            # since ctapipe 0.21, waveform is always 3d, also for gain selected data
+            # FIXME: this is the easiest solution to keep compatibility for ctapipe < 0.21
+            # once we drop all version < 0.21, the proper solution would be to directly fill
+            # the correct shape
+            if CTAPIPE_GE_0_21:
+                for c in (array_event.r0, array_event.r1):
+                    for tel_c in c.tel.values():
+                        if tel_c.waveform is not None and tel_c.waveform.ndim == 2:
+                            tel_c.waveform = tel_c.waveform[np.newaxis, ...]
+
             yield array_event
 
     @staticmethod
@@ -959,7 +965,7 @@ class LSTEventSource(EventSource):
         if trigger.event_type == EventType.UNKNOWN:
             self.log.warning(f'Event {array_event.index.event_id} has unknown event type, trigger: {trigger_bits:08b}')
 
-        if CTAPIPE_0_20:
+        if CTAPIPE_GE_0_20:
             array_event.r1.tel[tel_id].event_type = trigger.event_type
 
     def tag_flatfield_events(self, array_event):
@@ -1092,7 +1098,7 @@ class LSTEventSource(EventSource):
             r0 = R0CameraContainer(waveform=reordered_waveform)
             r1 = R1CameraContainer()
 
-        if CTAPIPE_0_20:
+        if CTAPIPE_GE_0_20:
             # reorder to nominal pixel order
             pixel_status = _reorder_pixel_status(
                 zfits_event.pixel_status, pixel_id_map, set_dvr_bits=not self.dvr_applied
@@ -1186,5 +1192,5 @@ class LSTEventSource(EventSource):
             return
 
         array_event.trigger.event_type = event_type
-        if CTAPIPE_0_20:
+        if CTAPIPE_GE_0_20:
             array_event.r1.tel[self.tel_id].event_type = event_type
diff --git a/src/ctapipe_io_lst/calibration.py b/src/ctapipe_io_lst/calibration.py
index c456bbb..8f5a9fa 100644
--- a/src/ctapipe_io_lst/calibration.py
+++ b/src/ctapipe_io_lst/calibration.py
@@ -16,6 +16,7 @@ from ctapipe.containers import FlatFieldContainer, MonitoringCameraContainer, Mo
 from ctapipe.io import HDF5TableReader, read_table
 from traitlets import Enum
 
+from .compat import CTAPIPE_GE_0_21
 from .containers import LSTArrayEventContainer
 
 
@@ -323,6 +324,7 @@ class LSTR0Corrections(TelescopeComponent):
             event.calibration.tel[tel_id].dl1.relative_factor = relative_factor
 
     def fill_time_correction(self, event):
+
         for tel_id in event.trigger.tels_with_trigger:
             r1 = event.r1.tel[tel_id]
             # store calibration data needed for dl1 calibration in ctapipe
@@ -337,10 +339,10 @@ class LSTR0Corrections(TelescopeComponent):
                 time_corr = self.mon_data.tel[tel_id].calibration.time_correction
                 # time_shift is subtracted in ctapipe,
                 # but time_correction should be added
-                if r1.selected_gain_channel is not None:
-                    time_shift -= time_corr[r1.selected_gain_channel, PIXEL_INDEX].to_value(u.ns)
-                else:
+                if CTAPIPE_GE_0_21 or r1.selected_gain_channel is None:
                     time_shift -= time_corr.to_value(u.ns)
+                else:
+                    time_shift -= time_corr[r1.selected_gain_channel, PIXEL_INDEX].to_value(u.ns)
 
             event.calibration.tel[tel_id].dl1.time_shift = time_shift
 
@@ -393,7 +395,7 @@ class LSTR0Corrections(TelescopeComponent):
         """
 
         if self.drs4_time_calibration_path.tel[tel_id] is None:
-            if selected_gain_channel is None:
+            if CTAPIPE_GE_0_21 or selected_gain_channel is None:
                 return np.zeros((N_GAINS, N_PIXELS))
             else:
                 return np.zeros(N_PIXELS)
@@ -402,16 +404,16 @@ class LSTR0Corrections(TelescopeComponent):
         if tel_id not in self.fan:
             self.load_drs4_time_calibration_file_for_tel(tel_id)
 
-        if selected_gain_channel is not None:
-            return calc_drs4_time_correction_gain_selected(
+        if CTAPIPE_GE_0_21 or selected_gain_channel is None:
+            return calc_drs4_time_correction_both_gains(
                 first_capacitors,
-                selected_gain_channel,
                 self.fan[tel_id],
                 self.fbn[tel_id],
             )
         else:
-            return calc_drs4_time_correction_both_gains(
+            return calc_drs4_time_correction_gain_selected(
                 first_capacitors,
+                selected_gain_channel,
                 self.fan[tel_id],
                 self.fbn[tel_id],
             )
diff --git a/src/ctapipe_io_lst/compat.py b/src/ctapipe_io_lst/compat.py
new file mode 100644
index 0000000..e5c199f
--- /dev/null
+++ b/src/ctapipe_io_lst/compat.py
@@ -0,0 +1,5 @@
+from ctapipe import __version__ as ctapipe_version
+
+CTAPIPE_VERSION = tuple(int(v) for v in ctapipe_version.split(".")[:3])
+CTAPIPE_GE_0_20 = CTAPIPE_VERSION >= (0, 20)
+CTAPIPE_GE_0_21 = CTAPIPE_VERSION >= (0, 21)


Test patch:
diff --git a/src/ctapipe_io_lst/tests/test_calib.py b/src/ctapipe_io_lst/tests/test_calib.py
index 0b61d9a..3dd55fd 100644
--- a/src/ctapipe_io_lst/tests/test_calib.py
+++ b/src/ctapipe_io_lst/tests/test_calib.py
@@ -1,5 +1,5 @@
 import os
-from importlib_resources import files, as_file
+from importlib.resources import files, as_file
 from pathlib import Path
 import pickle
 
@@ -9,6 +9,7 @@ import tables
 from traitlets.config import Config
 
 from ctapipe_io_lst.constants import HIGH_GAIN
+from ctapipe_io_lst.compat import CTAPIPE_GE_0_21
 
 
 resource_dir = files('ctapipe_io_lst') / 'tests/resources'
@@ -203,7 +204,10 @@ def test_missing_module():
             assert np.count_nonzero(waveform == 0) >= N_PIXELS_MODULE * (N_SAMPLES - 4)
 
             # waveforms in failing pixels must be all 0
-            assert np.all(waveform[failing_pixels[HIGH_GAIN]] == 0)
+            if CTAPIPE_GE_0_21:
+                np.testing.assert_equal(waveform[:, failing_pixels[HIGH_GAIN]], 0)
+            else:
+                np.testing.assert_equal(waveform[failing_pixels[HIGH_GAIN]], 0)
 
 def test_no_gain_selection():
     from ctapipe_io_lst import LSTEventSource
diff --git a/src/ctapipe_io_lst/tests/test_cta_r1.py b/src/ctapipe_io_lst/tests/test_cta_r1.py
index 870fa7e..794fea0 100644
--- a/src/ctapipe_io_lst/tests/test_cta_r1.py
+++ b/src/ctapipe_io_lst/tests/test_cta_r1.py
@@ -18,7 +18,7 @@ from protozfits.R1v1_pb2 import CameraConfiguration, Event, TelescopeDataStream
 from protozfits.R1v1_debug_pb2 import DebugEvent, DebugCameraConfiguration
 from protozfits.CoreMessages_pb2 import AnyArray
 
-from ctapipe_io_lst import LSTEventSource
+from ctapipe_io_lst import LSTEventSource, CTAPIPE_GE_0_21
 from ctapipe_io_lst.anyarray_dtypes import CDTS_AFTER_37201_DTYPE, TIB_DTYPE
 from ctapipe_io_lst.constants import CLOCK_FREQUENCY_KHZ, TriggerBits
 from ctapipe_io_lst.event_time import time_to_cta_high
@@ -47,8 +47,12 @@ DTYPE_TO_ANYARRAY_TYPE = {v: k for k, v in ANY_ARRAY_TYPE_TO_NUMPY_TYPE.items()}
 
 subarray = LSTEventSource.create_subarray(tel_id=1)
 GEOMETRY = subarray.tel[1].camera.geometry
+pulse_shape = subarray.tel[1].camera.readout.reference_pulse_shape[0]
+if CTAPIPE_GE_0_21:
+    pulse_shape = pulse_shape[np.newaxis, ...]
+
 waveform_model = WaveformModel(
-    reference_pulse=subarray.tel[1].camera.readout.reference_pulse_shape[0],
+    reference_pulse=pulse_shape,
     reference_pulse_sample_width=subarray.tel[1].camera.readout.reference_pulse_sample_width,
     sample_width=1 * u.ns,
 )
@@ -87,6 +91,8 @@ def create_shower(rng):
 
 def create_waveform(image, peak_time, num_samples=40, gains=(86, 5), offset=400):
     r1 = waveform_model.get_waveform(image, peak_time, num_samples)
+    if CTAPIPE_GE_0_21:
+        r1 = r1[0]
     return np.array([r1 * gain + offset for gain in gains]).astype(np.uint16)
 
 
@@ -328,7 +334,10 @@ def test_drs4_calibration(dummy_cta_r1):
 
             assert e.r1.tel[1].waveform.dtype == np.float32
             if e.trigger.event_type is EventType.SUBARRAY:
-                assert e.r1.tel[1].waveform.shape == (1855, 36)
+                if CTAPIPE_GE_0_21:
+                    assert e.r1.tel[1].waveform.shape == (1, 1855, 36)
+                else:
+                    assert e.r1.tel[1].waveform.shape == (1855, 36)
             else:
                 assert e.r1.tel[1].waveform.shape == (2, 1855, 36)
 
diff --git a/src/ctapipe_io_lst/tests/test_lsteventsource.py b/src/ctapipe_io_lst/tests/test_lsteventsource.py
index 1c77617..ed39657 100644
--- a/src/ctapipe_io_lst/tests/test_lsteventsource.py
+++ b/src/ctapipe_io_lst/tests/test_lsteventsource.py
@@ -12,7 +12,7 @@ from ctapipe.containers import CoordinateFrameType, EventType, PointingMode
 from ctapipe.calib.camera.gainselection import ThresholdGainSelector
 
 from ctapipe_io_lst.constants import LST1_LOCATION, N_GAINS, N_PIXELS_MODULE, N_SAMPLES, N_PIXELS
-from ctapipe_io_lst import CTAPIPE_0_20, TriggerBits, PixelStatus
+from ctapipe_io_lst import CTAPIPE_GE_0_20, CTAPIPE_GE_0_21, TriggerBits, PixelStatus
 
 test_data = Path(os.getenv('LSTCHAIN_TEST_DATA', 'test_data')).absolute()
 test_r0_dir = test_data / 'real/R0/20200218'
@@ -166,7 +166,7 @@ def test_missing_modules_r1v1():
         for gain, pixel in zip(missing_gain, missing_pixel):
             np.testing.assert_equal(event.r0.tel[1].waveform[gain, pixel], 0.0)
 
-        if CTAPIPE_0_20:
+        if CTAPIPE_GE_0_20:
             np.testing.assert_equal(event.lst.tel[1].evt.pixel_status, event.r1.tel[1].pixel_status)
 
     assert n_events == 40
@@ -209,7 +209,10 @@ def test_gain_selected():
             if event.r0.tel[1].waveform is not None:
                 assert event.r0.tel[1].waveform.shape == (N_GAINS, N_PIXELS, N_SAMPLES)
 
-            assert event.r1.tel[1].waveform.shape == (N_PIXELS, N_SAMPLES - 4)
+            if CTAPIPE_GE_0_21:
+                assert event.r1.tel[1].waveform.shape == (1, N_PIXELS, N_SAMPLES - 4)
+            else:
+                assert event.r1.tel[1].waveform.shape == (N_PIXELS, N_SAMPLES - 4)
 
             # compare to original file
             selected_gain = gain_selector(original_event.r1.tel[1].waveform)
@@ -258,7 +261,10 @@ def test_dvr():
             if dvr_event.r0.tel[1].waveform is not None:
                 assert dvr_event.r0.tel[1].waveform.shape == (N_GAINS, N_PIXELS, N_SAMPLES)
 
-            assert dvr_event.r1.tel[1].waveform.shape == (N_PIXELS, N_SAMPLES - 4)
+            if CTAPIPE_GE_0_21:
+                assert dvr_event.r1.tel[1].waveform.shape == (1, N_PIXELS, N_SAMPLES - 4)
+            else:
+                assert dvr_event.r1.tel[1].waveform.shape == (N_PIXELS, N_SAMPLES - 4)
 
             # compare to original file
             selected_gain = gain_selector(original_event.r1.tel[1].waveform)
@@ -266,7 +272,11 @@ def test_dvr():
             waveform = original_event.r1.tel[1].waveform[selected_gain, pixel_idx]
 
             readout_pixels = (dvr_event.lst.tel[1].evt.pixel_status & PixelStatus.DVR_STATUS) > 0
-            assert np.allclose(dvr_event.r1.tel[1].waveform[readout_pixels], waveform[readout_pixels])
+
+            if CTAPIPE_GE_0_21:
+                assert np.allclose(dvr_event.r1.tel[1].waveform[:, readout_pixels], waveform[readout_pixels])
+            else:
+                assert np.allclose(dvr_event.r1.tel[1].waveform[readout_pixels], waveform[readout_pixels])
 
     assert dvr_event.count == 199
 
@@ -354,7 +364,7 @@ def test_pedestal_events(tmp_path):
             else:
                 assert event.trigger.event_type != EventType.SKY_PEDESTAL
 
-            if CTAPIPE_0_20:
+            if CTAPIPE_GE_0_20:
                 assert event.r1.tel[1].event_type == event.trigger.event_type