diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 4396dff1..91142d39 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 495ec70f..f9df547e 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 21b9e900..66d3406c 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 6c33ab34..a60128da 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/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 diff --git a/src/ctapipe_io_lst/calibration.py b/src/ctapipe_io_lst/calibration.py index c456bbb7..8f5a9fae 100644 --- a/src/ctapipe_io_lst/calibration.py +++ b/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], ) diff --git a/src/ctapipe_io_lst/compat.py b/src/ctapipe_io_lst/compat.py new file mode 100644 index 00000000..e5c199f1 --- /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) diff --git a/src/ctapipe_io_lst/tests/test_calib.py b/src/ctapipe_io_lst/tests/test_calib.py index 0b61d9a3..3dd55fdf 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 @@ 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 870fa7ed..794fea0f 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_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) diff --git a/src/ctapipe_io_lst/tests/test_lsteventsource.py b/src/ctapipe_io_lst/tests/test_lsteventsource.py index c26bf0e6..61be5c38 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.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