Merge pull request #225 from cta-observatory/ctapipe_0.21.2

Add compatibility with ctapipe 0.21

.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:

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

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/*

src/ctapipe_io_lst/__init__.py

@@ -7,7 +7,6 @@
 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 .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 @@
 ]
 
 
-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 @@ def fill_from_cta_r1(self, array_event, zfits_event):
             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 @@ def _generator(self):
             # 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 @@ def fill_trigger_info(self, array_event):
         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 @@ def fill_r0r1_camera_container(self, zfits_event):
             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 @@ def check_interleaved_pedestal(self, array_event):
             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

src/ctapipe_io_lst/calibration.py

@@ -16,6 +16,7 @@
 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 @@ def calibrate(self, event: LSTArrayEventContainer):
             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 @@ def fill_time_correction(self, event):
                 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 @@ def get_drs4_time_correction(self, tel_id, first_capacitors, selected_gain_chann
         """
 
         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 @@ def get_drs4_time_correction(self, tel_id, first_capacitors, selected_gain_chann
         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],
             )

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)

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 @@
 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

src/ctapipe_io_lst/tests/test_cta_r1.py

@@ -18,7 +18,7 @@
 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 @@
 
 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)
 

src/ctapipe_io_lst/tests/test_lsteventsource.py

@@ -12,7 +12,7 @@
 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,15 +261,22 @@ 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)
             pixel_idx = np.arange(N_PIXELS)
             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