From 4be86da8be6f54bd752e04144d1306c2f4e05afe Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 10 Jul 2023 10:09:04 +0000 Subject: [PATCH 01/19] Add public method to only update the last readout times without applying dt correction --- src/ctapipe_io_lst/calibration.py | 34 +++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/src/ctapipe_io_lst/calibration.py b/src/ctapipe_io_lst/calibration.py index 845b00fa..137b5f24 100644 --- a/src/ctapipe_io_lst/calibration.py +++ b/src/ctapipe_io_lst/calibration.py @@ -480,6 +480,14 @@ def subtract_pedestal(self, event, tel_id): event.r1.tel[tel_id].selected_gain_channel, ) + def update_last_readout_times(self, event, tel_id): + lst = event.lst.tel[tel_id] + update_last_readout_times( + local_clock_counter=lst.evt.local_clock_counter, + first_capacitors=self.first_cap[tel_id], + last_readout_time=self.last_readout_time[tel_id], + expected_pixels_id=lst.svc.pixel_ids, + ) def time_lapse_corr(self, event, tel_id): """ @@ -1000,6 +1008,32 @@ def apply_timelapse_correction( ) +@njit(cache=True) +def update_last_readout_times( + local_clock_counter, + first_capacitors, + last_readout_time, + expected_pixels_id, +): + """ + Update the last readout time for all pixels / capcacitors + """ + n_modules = len(expected_pixels_id) // N_PIXELS_MODULE + for gain in range(N_GAINS): + for module in range(n_modules): + time_now = local_clock_counter[module] + for pixel_in_module in range(N_PIXELS_MODULE): + pixel_index = module * N_PIXELS_MODULE + pixel_in_module + pixel_id = expected_pixels_id[pixel_index] + + update_last_readout_time( + pixel_in_module=pixel_in_module, + first_capacitor=first_capacitors[gain, pixel_id], + time_now=time_now, + last_readout_time=last_readout_time[gain, pixel_id], + ) + + @njit(cache=True) def apply_timelapse_correction_gain_selected( waveform, From cbaa57b0d3bc08290168418e00ab4ec5df611cb8 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Fri, 15 Sep 2023 12:08:57 +0200 Subject: [PATCH 02/19] Set dvr bit for input r1 --- src/ctapipe_io_lst/__init__.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index 1f814c90..c8bcb706 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -466,11 +466,15 @@ def fill_from_cta_r1(self, array_event, zfits_event): def fill_lst_from_ctar1(self, zfits_event): evt = LSTEventContainer( - pixel_status=zfits_event.pixel_status, + pixel_status=zfits_event.pixel_status.copy(), first_capacitor_id=zfits_event.first_cell_id, calibration_monitoring_id=zfits_event.calibration_monitoring_id, local_clock_counter=zfits_event.module_hires_local_clock_counter, ) + # set bits for dvr if not already set + if not self.dvr_applied: + not_broken = get_channel_info(evt.pixel_status) != 0 + evt.pixel_status[not_broken] |= PixelStatus.DVR_STATUS_0 if zfits_event.debug is not None: debug = zfits_event.debug From 021ebd161dd121e4b94392bc0469febcaf07f458 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Fri, 15 Sep 2023 10:12:36 +0000 Subject: [PATCH 03/19] add datastream to MultiFiles --- src/ctapipe_io_lst/multifiles.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/ctapipe_io_lst/multifiles.py b/src/ctapipe_io_lst/multifiles.py index 91bb11b0..c94f215b 100644 --- a/src/ctapipe_io_lst/multifiles.py +++ b/src/ctapipe_io_lst/multifiles.py @@ -106,6 +106,7 @@ def __init__(self, path, *args, **kwargs): self._events_tables = {} self._headers = {} self.camera_config = None + self.data_stream = None self.dvr_applied = None if self.all_streams and file_info is not None: @@ -182,7 +183,10 @@ def _load_next_subrun(self, stream): # new files use CameraConfiguration self.cta_r1 = True config = next(file_.CameraConfiguration) - + + if self.data_stream is None: + self.data_stream = next(file_.DataStream) + if self.camera_config is None: self.camera_config = config From b35752bcd004a48e0a71aa946a436e55c83fbbce Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Fri, 15 Sep 2023 10:15:35 +0000 Subject: [PATCH 04/19] ctapipe 0.20 --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 9a633efa..4bb9d909 100644 --- a/setup.cfg +++ b/setup.cfg @@ -31,7 +31,7 @@ python_requires = >=3.8 zip_safe = False install_requires= astropy~=5.2 - ctapipe~=0.19.0 + ctapipe~=0.20.0 protozfits~=2.2 numpy>=1.20 From fec629211470f2cc642fa5a82d2ae510b49c96b3 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Fri, 15 Sep 2023 14:29:47 +0200 Subject: [PATCH 05/19] Reorder pixels by id --- src/ctapipe_io_lst/__init__.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index c8bcb706..50a92b7d 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -310,6 +310,8 @@ def __init__(self, input_url=None, **kwargs): self.n_samples = self.camera_config.num_samples self.lst_service = self.fill_lst_service_container(self.tel_id, self.camera_config) + self.reverse_pixel_id_map = np.argsort(self.pixel_id_map) + self._subarray = self.create_subarray(self.tel_id, reference_location) self.r0_r1_calibrator = LSTR0Corrections( subarray=self._subarray, parent=self @@ -438,20 +440,24 @@ def create_subarray(tel_id=1, reference_location=None): def fill_from_cta_r1(self, array_event, zfits_event): tel_id = self.tel_id + # FIXME: missing modules / pixels + # FIXME: DVR? should not happen in r1 but dl0, but our own converter uses the old R1 + # TODO: decide from applied preprocessing options, not gains if zfits_event.num_channels == 2: shape = (zfits_event.num_channels, zfits_event.num_pixels, zfits_event.num_samples) array_event.r0.tel[self.tel_id] = R0CameraContainer( - waveform=zfits_event.waveform.reshape(shape), + waveform=zfits_event.waveform.reshape(shape)[:, self.reverse_pixel_id_map], ) array_event.r1.tel[self.tel_id] = R1CameraContainer() + else: has_high_gain = (zfits_event.pixel_status & PixelStatus.HIGH_GAIN_STORED).astype(bool) selected_gain_channel = np.where(has_high_gain, 0, 1) shape = (zfits_event.num_pixels, zfits_event.num_samples) array_event.r1.tel[self.tel_id] = R1CameraContainer( - waveform=zfits_event.waveform.reshape(shape), + waveform=zfits_event.waveform.reshape(shape)[self.reverse_pixel_id_map], selected_gain_channel=selected_gain_channel, ) array_event.r0.tel[self.tel_id] = R0CameraContainer() From 107f179dac761c187bcbccb159fc7b54e3de899d Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 20 Sep 2023 15:49:07 +0200 Subject: [PATCH 06/19] Use a fresh container for each new event --- src/ctapipe_io_lst/__init__.py | 42 ++++++++++++++----------- src/ctapipe_io_lst/tests/test_cta_r1.py | 7 ++++- 2 files changed, 30 insertions(+), 19 deletions(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index 50a92b7d..b240132d 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -25,7 +25,8 @@ from ctapipe.core.traits import Bool, Float, Enum, Path from ctapipe.containers import ( CoordinateFrameType, PixelStatusContainer, EventType, PointingMode, R0CameraContainer, R1CameraContainer, - SchedulingBlockContainer, ObservationBlockContainer, + SchedulingBlockContainer, ObservationBlockContainer, MonitoringContainer, MonitoringCameraContainer, + EventIndexContainer, ) from ctapipe.coordinates import CameraFrame @@ -443,7 +444,6 @@ def fill_from_cta_r1(self, array_event, zfits_event): # FIXME: missing modules / pixels # FIXME: DVR? should not happen in r1 but dl0, but our own converter uses the old R1 - # TODO: decide from applied preprocessing options, not gains if zfits_event.num_channels == 2: shape = (zfits_event.num_channels, zfits_event.num_pixels, zfits_event.num_samples) array_event.r0.tel[self.tel_id] = R0CameraContainer( @@ -540,28 +540,31 @@ def fill_lst_from_ctar1(self, zfits_event): def _generator(self): - # container for LST data - array_event = LSTArrayEventContainer() - array_event.meta['input_url'] = self.input_url - array_event.meta['max_events'] = self.max_events - array_event.meta['origin'] = 'LSTCAM' - # also add service container to the event section # initialize general monitoring container - self.initialize_mon_container(array_event) + mon = self.initialize_mon_container() # loop on events for count, zfits_event in enumerate(self.multi_file): - array_event.count = count - array_event.index.event_id = zfits_event.event_id - array_event.index.obs_id = self.local_run_id - # Skip "empty" events that occur at the end of some runs if zfits_event.event_id == 0: self.log.warning('Event with event_id=0 found, skipping') continue + # container for LST data + array_event = LSTArrayEventContainer( + count=count, + index=EventIndexContainer( + obs_id=self.local_run_id, + event_id=zfits_event.event_id, + ), + mon=mon, + ) + array_event.meta['input_url'] = self.input_url + array_event.meta['max_events'] = self.max_events + array_event.meta['origin'] = 'LSTCAM' + array_event.lst.tel[self.tel_id].svc = self.lst_service if self.cta_r1: @@ -1007,15 +1010,12 @@ def fill_r0r1_container(self, array_event, zfits_event): array_event.r0.tel[self.tel_id] = r0 array_event.r1.tel[self.tel_id] = r1 - def initialize_mon_container(self, array_event): + def initialize_mon_container(self): """ Fill with MonitoringContainer. For the moment, initialize only the PixelStatusContainer """ - container = array_event.mon - mon_camera_container = container.tel[self.tel_id] - shape = (N_GAINS, N_PIXELS) # all pixels broken by default status_container = PixelStatusContainer( @@ -1023,7 +1023,13 @@ def initialize_mon_container(self, array_event): pedestal_failing_pixels=np.zeros(shape, dtype=bool), flatfield_failing_pixels=np.zeros(shape, dtype=bool), ) - mon_camera_container.pixel_status = status_container + + camera_container = MonitoringCameraContainer( + pixel_status = status_container, + ) + container = MonitoringContainer() + container.tel[self.tel_id] = camera_container + return container def fill_mon_container(self, array_event, zfits_event): """ diff --git a/src/ctapipe_io_lst/tests/test_cta_r1.py b/src/ctapipe_io_lst/tests/test_cta_r1.py index 7ebc4d29..0e8d0f5b 100644 --- a/src/ctapipe_io_lst/tests/test_cta_r1.py +++ b/src/ctapipe_io_lst/tests/test_cta_r1.py @@ -8,7 +8,7 @@ import astropy.units as u import protozfits -from protozfits.CTA_R1_pb2 import CameraConfiguration, Event +from protozfits.CTA_R1_pb2 import CameraConfiguration, Event, TelescopeDataStream from protozfits.Debug_R1_pb2 import DebugEvent, DebugCameraConfiguration from protozfits.CoreMessages_pb2 import AnyArray from traitlets.config import Config @@ -144,6 +144,9 @@ def dummy_cta_r1(dummy_cta_r1_dir): ) ) + # assume evb did no pe calibration + data_stream = TelescopeDataStream(sb_id=10000, obs_id=10000, waveform_offset=400, waveform_scale=1) + rng = np.random.default_rng() streams = [] @@ -151,6 +154,8 @@ def dummy_cta_r1(dummy_cta_r1_dir): for stream_path in stream_paths: stream = stack.enter_context(protozfits.ProtobufZOFits(n_tiles=5, rows_per_tile=20, compression_block_size_kb=64 * 1024)) stream.open(str(stream_path)) + stream.move_to_new_table("DataStream") + stream.write_message(data_stream) stream.move_to_new_table("CameraConfiguration") stream.write_message(camera_config) stream.move_to_new_table("Events") From 303ee29604e3fbfdd17b848e2b570270dfff251d Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 20 Sep 2023 15:29:51 +0000 Subject: [PATCH 07/19] Add code to parse evb pre-processing of event types --- src/ctapipe_io_lst/evb_preprocessing.py | 47 +++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 src/ctapipe_io_lst/evb_preprocessing.py diff --git a/src/ctapipe_io_lst/evb_preprocessing.py b/src/ctapipe_io_lst/evb_preprocessing.py new file mode 100644 index 00000000..ae97ef8c --- /dev/null +++ b/src/ctapipe_io_lst/evb_preprocessing.py @@ -0,0 +1,47 @@ +from enum import IntEnum +from .constants import TriggerBits +from collections import defaultdict + +class EVBPreprocessing(IntEnum): + """ + The preprocessing steps that can be applied by EVB. + + The values of this Enum is the index of this step in the tdp_action array. + + See EVB ICD section in + https://edms.cern.ch/ui/file/2411710/2.6/LSTMST-ICD-20191206.pdf + """ + GAIN_SELECTION = 0 + BASELINE_SUBTRACTION = 1 + DELTA_T_CORRECTION = 2 + KEEP_ALL_EVENTS = 3 + MUON_SEARCH = 4 + PEDESTAL_SUBTRACTION = 5 + CALIBRATION = 6 + PEDESTAL_SUM = 7 + CALIBRATION_SUM = 8 + + +def get_processings_for_trigger_bits(camera_configuration): + """ + Parse the tdp_action/type information into a dict mapping + """ + tdp_type = camera_configuration.debug.tdp_type + tdp_action = camera_configuration.debug.tdp_action + + # first bit (no shift) is default handling + default = {step for step in EVBPreprocessing if tdp_action[step] & 1} + actions = defaultdict(lambda: default) + + # the following bits refer to the entries in tdp_type + for i, trigger_bits in enumerate(tdp_type, start=1): + # all-zero trigger bits can be ignored + if trigger_bits == 0: + continue + + actions[TriggerBits(int(trigger_bits))] = { + step for step in EVBPreprocessing + if tdp_action[step] & (1 << i) + } + + return actions From ac2cef8dba9e019d46e0315437143f9b99c10883 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 26 Sep 2023 18:24:07 +0200 Subject: [PATCH 08/19] Fix event type for CTAr1 case --- src/ctapipe_io_lst/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index b240132d..cdeb9cd7 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -468,7 +468,7 @@ def fill_from_cta_r1(self, array_event, zfits_event): trigger.time = cta_high_res_to_time(zfits_event.event_time_s, zfits_event.event_time_qns) trigger.tels_with_trigger = [tel_id] trigger.tel[tel_id].time = trigger.time - trigger.event_type = self._event_type_from_trigger_bits(zfits_event.event_type) + trigger.event_type = EventType(zfits_event.event_type) def fill_lst_from_ctar1(self, zfits_event): evt = LSTEventContainer( From 63150286c2830cf29ebb46399df0be3815fc0cc9 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 28 Sep 2023 14:58:28 +0000 Subject: [PATCH 09/19] Make it work with ctapipe-process --- src/ctapipe_io_lst/__init__.py | 59 +++++++++++++++++++------ src/ctapipe_io_lst/calibration.py | 2 +- src/ctapipe_io_lst/evb_preprocessing.py | 2 +- 3 files changed, 48 insertions(+), 15 deletions(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index cdeb9cd7..c3e3cd4e 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -51,6 +51,8 @@ PixelStatus, TriggerBits, ) +from .evb_preprocessing import get_processings_for_trigger_bits, EVBPreprocessing + log = logging.getLogger(__name__) @@ -301,6 +303,8 @@ def __init__(self, input_url=None, **kwargs): self.n_pixels = self.camera_config.num_pixels self.n_samples = self.camera_config.num_samples_nominal self.lst_service = self.fill_lst_service_container_ctar1(self.camera_config) + self.evb_preprocessing = get_processings_for_trigger_bits(self.camera_config) + self.data_stream = self.multi_file.data_stream else: self.tel_id = self.camera_config.telescope_id self.local_run_id = self.camera_config.configuration_id @@ -310,6 +314,8 @@ def __init__(self, input_url=None, **kwargs): self.n_pixels = self.camera_config.num_pixels self.n_samples = self.camera_config.num_samples self.lst_service = self.fill_lst_service_container(self.tel_id, self.camera_config) + self.evb_preprocessing = None + self.data_stream = None self.reverse_pixel_id_map = np.argsort(self.pixel_id_map) @@ -384,6 +390,12 @@ def scheduling_blocks(self): @property def datalevels(self): + if self.cta_r1: + if EVBPreprocessing.CALIBRATION in self.evb_preprocessing[TriggerBits.MONO]: + return (DataLevel.R1, ) + else: + return (DataLevel.R0, ) + if self.r0_r1_calibrator.calibration_path is not None: return (DataLevel.R0, DataLevel.R1) return (DataLevel.R0, ) @@ -440,27 +452,40 @@ def create_subarray(tel_id=1, reference_location=None): def fill_from_cta_r1(self, array_event, zfits_event): tel_id = self.tel_id + scale = self.data_stream.waveform_scale + offset = self.data_stream.waveform_offset + pixel_id_map = self.camera_config.pixel_id_map # FIXME: missing modules / pixels # FIXME: DVR? should not happen in r1 but dl0, but our own converter uses the old R1 - if zfits_event.num_channels == 2: - shape = (zfits_event.num_channels, zfits_event.num_pixels, zfits_event.num_samples) - array_event.r0.tel[self.tel_id] = R0CameraContainer( - waveform=zfits_event.waveform.reshape(shape)[:, self.reverse_pixel_id_map], - ) - array_event.r1.tel[self.tel_id] = R1CameraContainer() + # reorder to nominal pixel order + pixel_status = np.zeros(N_PIXELS, dtype=zfits_event.pixel_status.dtype) + pixel_status[pixel_id_map] = zfits_event.pixel_status + + # set dvr bits, so that calibration code doesn't reset them... + not_broken = get_channel_info(pixel_status) != 0 + pixel_status[not_broken] |= PixelStatus.DVR_STATUS_0 + n_channels = zfits_event.num_channels + n_pixels = zfits_event.num_pixels + n_samples = zfits_event.num_samples + + readout_shape = (n_channels, n_pixels, n_samples) + waveform = zfits_event.waveform.reshape(readout_shape) + waveform = (waveform.astype(np.float32) - offset) / scale + + reordered_waveform = np.full((n_channels, N_PIXELS, n_samples), 0.0, dtype=np.float32) + reordered_waveform[:, pixel_id_map] = waveform + + # FIXME, check using evb_preprocessing and make ctapipe support 2 gains + if zfits_event.num_channels == 2: + selected_gain_channel = None else: - has_high_gain = (zfits_event.pixel_status & PixelStatus.HIGH_GAIN_STORED).astype(bool) + has_high_gain = (pixel_status & PixelStatus.HIGH_GAIN_STORED).astype(bool) selected_gain_channel = np.where(has_high_gain, 0, 1) + waveform = waveform[0] - shape = (zfits_event.num_pixels, zfits_event.num_samples) - array_event.r1.tel[self.tel_id] = R1CameraContainer( - waveform=zfits_event.waveform.reshape(shape)[self.reverse_pixel_id_map], - selected_gain_channel=selected_gain_channel, - ) - array_event.r0.tel[self.tel_id] = R0CameraContainer() array_event.lst.tel[self.tel_id] = self.fill_lst_from_ctar1(zfits_event) @@ -470,6 +495,14 @@ def fill_from_cta_r1(self, array_event, zfits_event): trigger.tel[tel_id].time = trigger.time trigger.event_type = EventType(zfits_event.event_type) + array_event.r1.tel[self.tel_id] = R1CameraContainer( + waveform=waveform, + selected_gain_channel=selected_gain_channel, + pixel_status=pixel_status, + event_type=EventType(zfits_event.event_type), + event_time=trigger.time, + ) + def fill_lst_from_ctar1(self, zfits_event): evt = LSTEventContainer( pixel_status=zfits_event.pixel_status.copy(), diff --git a/src/ctapipe_io_lst/calibration.py b/src/ctapipe_io_lst/calibration.py index 845b00fa..1a1d7e87 100644 --- a/src/ctapipe_io_lst/calibration.py +++ b/src/ctapipe_io_lst/calibration.py @@ -257,7 +257,7 @@ def apply_drs4_corrections(self, event: LSTArrayEventContainer): r1.waveform -= self.offset.tel[tel_id] broken_pixels = event.mon.tel[tel_id].pixel_status.hardware_failing_pixels - dvred_pixels = (event.lst.tel[tel_id].evt.pixel_status & PixelStatus.DVR_STATUS ) == 0 + dvred_pixels = (event.lst.tel[tel_id].evt.pixel_status & PixelStatus.DVR_STATUS) == 0 invalid_pixels = broken_pixels | dvred_pixels if r1.selected_gain_channel is None: diff --git a/src/ctapipe_io_lst/evb_preprocessing.py b/src/ctapipe_io_lst/evb_preprocessing.py index ae97ef8c..da1f0f8b 100644 --- a/src/ctapipe_io_lst/evb_preprocessing.py +++ b/src/ctapipe_io_lst/evb_preprocessing.py @@ -38,7 +38,7 @@ def get_processings_for_trigger_bits(camera_configuration): # all-zero trigger bits can be ignored if trigger_bits == 0: continue - + actions[TriggerBits(int(trigger_bits))] = { step for step in EVBPreprocessing if tdp_action[step] & (1 << i) From a23b1f0b45c7e505a62563e7660186cfa7b02746 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 28 Sep 2023 15:21:36 +0000 Subject: [PATCH 10/19] Fix application of scale/offset --- src/ctapipe_io_lst/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index c3e3cd4e..09b73e58 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -473,7 +473,7 @@ def fill_from_cta_r1(self, array_event, zfits_event): readout_shape = (n_channels, n_pixels, n_samples) waveform = zfits_event.waveform.reshape(readout_shape) - waveform = (waveform.astype(np.float32) - offset) / scale + waveform = waveform.astype(np.float32) / scale - offset reordered_waveform = np.full((n_channels, N_PIXELS, n_samples), 0.0, dtype=np.float32) reordered_waveform[:, pixel_id_map] = waveform From 01c5fa4319017e60d7c16eca07faa4e6f5e70632 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 28 Sep 2023 16:10:38 +0000 Subject: [PATCH 11/19] Fix looping over r0 instead of triggered telescope ids, assign reordered waveform --- src/ctapipe_io_lst/__init__.py | 2 +- src/ctapipe_io_lst/calibration.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index 09b73e58..99a29cc3 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -477,6 +477,7 @@ def fill_from_cta_r1(self, array_event, zfits_event): reordered_waveform = np.full((n_channels, N_PIXELS, n_samples), 0.0, dtype=np.float32) reordered_waveform[:, pixel_id_map] = waveform + waveform = reordered_waveform # FIXME, check using evb_preprocessing and make ctapipe support 2 gains if zfits_event.num_channels == 2: @@ -486,7 +487,6 @@ def fill_from_cta_r1(self, array_event, zfits_event): selected_gain_channel = np.where(has_high_gain, 0, 1) waveform = waveform[0] - array_event.lst.tel[self.tel_id] = self.fill_lst_from_ctar1(zfits_event) trigger = array_event.trigger diff --git a/src/ctapipe_io_lst/calibration.py b/src/ctapipe_io_lst/calibration.py index 1a1d7e87..0f87fc4b 100644 --- a/src/ctapipe_io_lst/calibration.py +++ b/src/ctapipe_io_lst/calibration.py @@ -224,11 +224,11 @@ def __init__(self, subarray, config=None, parent=None, **kwargs): def apply_drs4_corrections(self, event: LSTArrayEventContainer): self.update_first_capacitors(event) - for tel_id, r0 in event.r0.tel.items(): + for tel_id in event.trigger.tels_with_trigger: r1 = event.r1.tel[tel_id] # If r1 was not yet filled, copy of r0 converted if r1.waveform is None: - r1.waveform = r0.waveform + r1.waveform = event.r0.tel[tel_id].waveform # float32 can represent all values of uint16 exactly, # so this does not loose precision. From 8cc56827993c1be302b2038f66b779dd68f1cd2e Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 23 Oct 2023 14:51:03 +0000 Subject: [PATCH 12/19] Also update last_readout_times when not applying delta t correction --- src/ctapipe_io_lst/calibration.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/ctapipe_io_lst/calibration.py b/src/ctapipe_io_lst/calibration.py index ad01186b..2f69fab7 100644 --- a/src/ctapipe_io_lst/calibration.py +++ b/src/ctapipe_io_lst/calibration.py @@ -240,6 +240,8 @@ def apply_drs4_corrections(self, event: LSTArrayEventContainer): if self.apply_timelapse_correction: self.time_lapse_corr(event, tel_id) + else: + self.update_last_readout_times(event, tel_id) if self.apply_spike_correction: if self.spike_correction_method == 'subtraction': From 6bc361b2772ff97c9b0ca8d73967f9fb420ca3ec Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 24 Oct 2023 16:07:55 +0200 Subject: [PATCH 13/19] Support both ctapipe 0.19 and 0.20 --- .github/workflows/ci.yml | 4 ++-- setup.cfg | 2 +- src/ctapipe_io_lst/__init__.py | 23 +++++++++++++++++------ src/ctapipe_io_lst/tests/test_cta_r1.py | 2 ++ 4 files changed, 22 insertions(+), 9 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 34630b96..8cfe8420 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -18,8 +18,8 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ["3.8", "3.9", "3.10", "3.11"] - ctapipe-version: ["0.19.2"] + python-version: ["3.9", "3.10", "3.11"] + ctapipe-version: ["0.19.3", "0.20.0"] defaults: run: diff --git a/setup.cfg b/setup.cfg index 4bb9d909..3650958b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -31,7 +31,7 @@ python_requires = >=3.8 zip_safe = False install_requires= astropy~=5.2 - ctapipe~=0.20.0 + ctapipe >=0.19.0,<0.21.0a0 protozfits~=2.2 numpy>=1.20 diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index 99a29cc3..9feeb066 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -7,6 +7,7 @@ import numpy as np from astropy import units as u from pkg_resources import resource_filename +from ctapipe import __version__ as ctapipe_version from ctapipe.core import Provenance from ctapipe.instrument import ( ReflectorShape, @@ -53,11 +54,17 @@ from .evb_preprocessing import get_processings_for_trigger_bits, EVBPreprocessing +__all__ = [ + 'LSTEventSource', + '__version__', +] -log = logging.getLogger(__name__) + +CTAPIPE_VERSION = tuple(int(v) for v in ctapipe_version.split(".")[:3]) +CTAPIPE_0_20 = CTAPIPE_VERSION >= (0, 20) -__all__ = ['LSTEventSource', '__version__'] +log = logging.getLogger(__name__) # Date from which the flatfield heuristic will be switch off by default @@ -495,14 +502,18 @@ def fill_from_cta_r1(self, array_event, zfits_event): trigger.tel[tel_id].time = trigger.time trigger.event_type = EventType(zfits_event.event_type) - array_event.r1.tel[self.tel_id] = R1CameraContainer( + r1 = R1CameraContainer( waveform=waveform, selected_gain_channel=selected_gain_channel, - pixel_status=pixel_status, - event_type=EventType(zfits_event.event_type), - event_time=trigger.time, ) + if CTAPIPE_0_20: + r1.pixel_status = pixel_status + r1.event_type = EventType(zfits_event.event_type) + r1.event_time = trigger.time + + array_event.r1.tel[self.tel_id] = r1 + def fill_lst_from_ctar1(self, zfits_event): evt = LSTEventContainer( pixel_status=zfits_event.pixel_status.copy(), diff --git a/src/ctapipe_io_lst/tests/test_cta_r1.py b/src/ctapipe_io_lst/tests/test_cta_r1.py index 0e8d0f5b..4379c769 100644 --- a/src/ctapipe_io_lst/tests/test_cta_r1.py +++ b/src/ctapipe_io_lst/tests/test_cta_r1.py @@ -141,6 +141,8 @@ def dummy_cta_r1(dummy_cta_r1_dir): cs_serial="???", evb_version="evb-dummy", cdhs_version="evb-dummy", + tdp_type=to_anyarray(np.zeros(15, dtype=np.uint16)), + tdp_action=to_anyarray(np.zeros(16, dtype=np.uint16)), ) ) From c20ce92fff3d6a6424098d59d0c840e48a18471d Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 25 Oct 2023 10:07:02 +0200 Subject: [PATCH 14/19] Use importlib resources instead of deprecated pkg_resources --- src/ctapipe_io_lst/__init__.py | 22 ++++++++++------------ src/ctapipe_io_lst/tests/test_calib.py | 11 +++++------ 2 files changed, 15 insertions(+), 18 deletions(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index 9feeb066..00b2e1c5 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -2,11 +2,11 @@ """ EventSource for LSTCam protobuf-fits.fz-files. """ +from importlib.resources import files, as_file from ctapipe.instrument.subarray import EarthLocation import logging import numpy as np from astropy import units as u -from pkg_resources import resource_filename from ctapipe import __version__ as ctapipe_version from ctapipe.core import Provenance from ctapipe.instrument import ( @@ -54,6 +54,7 @@ from .evb_preprocessing import get_processings_for_trigger_bits, EVBPreprocessing + __all__ = [ 'LSTEventSource', '__version__', @@ -103,11 +104,9 @@ def get_channel_info(pixel_status): def load_camera_geometry(): ''' Load camera geometry from bundled resources of this repo ''' - f = resource_filename( - 'ctapipe_io_lst', 'resources/LSTCam.camgeom.fits.gz' - ) - Provenance().add_input_file(f, role="CameraGeometry") - cam = CameraGeometry.from_table(f) + with as_file(files("ctapipe_io_lst") / "resources/LSTCam.camgeom.fits.gz") as path: + Provenance().add_input_file(path, role="CameraGeometry") + cam = CameraGeometry.from_table(path) cam.frame = CameraFrame(focal_length=OPTICS.effective_focal_length) return cam @@ -128,13 +127,12 @@ def read_pulse_shapes(): # temporary replace the reference pulse shape # ("oversampled_pulse_LST_8dynode_pix6_20200204.dat") # with a dummy one in order to disable the charge corrections in the charge extractor - infilename = resource_filename( - 'ctapipe_io_lst', - 'resources/oversampled_pulse_LST_8dynode_pix6_20200204.dat' - ) - data = np.genfromtxt(infilename, dtype='float', comments='#') - Provenance().add_input_file(infilename, role="PulseShapes") + pulse_shape_path = 'resources/oversampled_pulse_LST_8dynode_pix6_20200204.dat' + with as_file(files("ctapipe_io_lst") / pulse_shape_path) as path: + data = np.genfromtxt(path, dtype='float', comments='#') + Provenance().add_input_file(path, role="PulseShapes") + daq_time_per_sample = data[0, 0] * u.ns pulse_shape_time_step = data[0, 1] * u.ns diff --git a/src/ctapipe_io_lst/tests/test_calib.py b/src/ctapipe_io_lst/tests/test_calib.py index 4ad59a00..80a58e5e 100644 --- a/src/ctapipe_io_lst/tests/test_calib.py +++ b/src/ctapipe_io_lst/tests/test_calib.py @@ -5,12 +5,10 @@ from traitlets.config import Config import numpy as np import tables -import pkg_resources +from importlib_resources import files, as_file -resource_dir = Path(pkg_resources.resource_filename( - 'ctapipe_io_lst', 'tests/resources' -)) +resource_dir = files('ctapipe_io_lst') / 'tests/resources' test_data = Path(os.getenv('LSTCHAIN_TEST_DATA', 'test_data')).absolute() test_r0_path = test_data / 'real/R0/20200218/LST-1.1.Run02008.0000_first50.fits.fz' test_r0_calib_path = test_data / 'real/R0/20200218/LST-1.1.Run02006.0004.fits.fz' @@ -41,8 +39,9 @@ def test_get_first_capacitor(): first_capacitor_id = event.lst.tel[tel_id].evt.first_capacitor_id - with tables.open_file(resource_dir / 'first_caps.hdf5', 'r') as f: - expected = f.root.first_capacitor_for_modules[:] + with as_file(resource_dir / 'first_caps.hdf5') as path: + with tables.open_file(path, 'r') as f: + expected = f.root.first_capacitor_for_modules[:] first_caps = get_first_capacitors_for_pixels(first_capacitor_id) From 5730ee7e4aa2ba2d869fe7dcadd592a6031e1cb2 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 25 Oct 2023 10:07:41 +0200 Subject: [PATCH 15/19] Set pixel status also for old data format --- src/ctapipe_io_lst/__init__.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index 00b2e1c5..d9971027 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -1041,6 +1041,20 @@ def fill_r0r1_camera_container(self, zfits_event): r0 = R0CameraContainer(waveform=reordered_waveform) r1 = R1CameraContainer() + if CTAPIPE_0_20: + # reorder to nominal pixel order + pixel_status = np.zeros(N_PIXELS, dtype=zfits_event.pixel_status.dtype) + pixel_status[pixel_id_map] = zfits_event.pixel_status + + # set dvr bits, so that calibration code doesn't reset them... + not_broken = get_channel_info(pixel_status) != 0 + pixel_status[not_broken] |= PixelStatus.DVR_STATUS_0 + + r1.pixel_status = pixel_status + r1.event_time = cta_high_res_to_time( + zfits_event.trigger_time_s, zfits_event.trigger_time_qns, + ) + return r0, r1 def fill_r0r1_container(self, array_event, zfits_event): From 407598e61d9995d6fab3d8e9278633d34cc0ef23 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 25 Oct 2023 10:11:56 +0200 Subject: [PATCH 16/19] Refactor common code into function --- src/ctapipe_io_lst/__init__.py | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index d9971027..bf23dfc3 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -141,6 +141,18 @@ def read_pulse_shapes(): return daq_time_per_sample, pulse_shape_time_step, data[1:].T +def _reorder_pixel_status(pixel_status, pixel_id_map, set_dvr_bits=True): + reordered_pixel_status = np.zeros(N_PIXELS, dtype=pixel_status.dtype) + reordered_pixel_status[pixel_id_map] = pixel_status + + # set dvr bits, so that calibration code doesn't reset them... + if set_dvr_bits: + not_broken = get_channel_info(pixel_status) != 0 + reordered_pixel_status[not_broken] |= PixelStatus.DVR_STATUS_0 + + return reordered_pixel_status + + class LSTEventSource(EventSource): """ EventSource for LST R0 data. @@ -465,12 +477,9 @@ def fill_from_cta_r1(self, array_event, zfits_event): # FIXME: DVR? should not happen in r1 but dl0, but our own converter uses the old R1 # reorder to nominal pixel order - pixel_status = np.zeros(N_PIXELS, dtype=zfits_event.pixel_status.dtype) - pixel_status[pixel_id_map] = zfits_event.pixel_status - - # set dvr bits, so that calibration code doesn't reset them... - not_broken = get_channel_info(pixel_status) != 0 - pixel_status[not_broken] |= PixelStatus.DVR_STATUS_0 + pixel_status = _reorder_pixel_status( + zfits_event.pixel_status, pixel_id_map, set_dvr_bits=True + ) n_channels = zfits_event.num_channels n_pixels = zfits_event.num_pixels @@ -1043,13 +1052,9 @@ def fill_r0r1_camera_container(self, zfits_event): if CTAPIPE_0_20: # reorder to nominal pixel order - pixel_status = np.zeros(N_PIXELS, dtype=zfits_event.pixel_status.dtype) - pixel_status[pixel_id_map] = zfits_event.pixel_status - - # set dvr bits, so that calibration code doesn't reset them... - not_broken = get_channel_info(pixel_status) != 0 - pixel_status[not_broken] |= PixelStatus.DVR_STATUS_0 - + pixel_status = _reorder_pixel_status( + zfits_event.pixel_status, pixel_id_map, set_dvr_bits=True + ) r1.pixel_status = pixel_status r1.event_time = cta_high_res_to_time( zfits_event.trigger_time_s, zfits_event.trigger_time_qns, From 7a3d1f06f68ce75861c0ba597809970fed8f5f40 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 25 Oct 2023 10:15:43 +0200 Subject: [PATCH 17/19] Update required python versions --- setup.cfg | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/setup.cfg b/setup.cfg index 3650958b..c64f9652 100644 --- a/setup.cfg +++ b/setup.cfg @@ -18,16 +18,16 @@ classifiers = Topic :: Scientific/Engineering :: Astronomy Topic :: Scientific/Engineering :: Physics Programming Language :: Python :: 3 :: Only - Programming Language :: Python :: 3.8 Programming Language :: Python :: 3.9 Programming Language :: Python :: 3.10 + Programming Language :: Python :: 3.11 [options] packages = find: package_dir = = src -python_requires = >=3.8 +python_requires = >=3.9 zip_safe = False install_requires= astropy~=5.2 From 24e61381ceb9de3885636c2e45f9905d10d0926b Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 25 Oct 2023 10:31:01 +0200 Subject: [PATCH 18/19] Fix _reorder_pixel_status for case of missing modules --- src/ctapipe_io_lst/__init__.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index bf23dfc3..59981de1 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -142,14 +142,12 @@ def read_pulse_shapes(): def _reorder_pixel_status(pixel_status, pixel_id_map, set_dvr_bits=True): - reordered_pixel_status = np.zeros(N_PIXELS, dtype=pixel_status.dtype) - reordered_pixel_status[pixel_id_map] = pixel_status - - # set dvr bits, so that calibration code doesn't reset them... if set_dvr_bits: not_broken = get_channel_info(pixel_status) != 0 - reordered_pixel_status[not_broken] |= PixelStatus.DVR_STATUS_0 + pixel_status = pixel_status[not_broken] | PixelStatus.DVR_STATUS_0 + reordered_pixel_status = np.zeros(N_PIXELS, dtype=pixel_status.dtype) + reordered_pixel_status[pixel_id_map] = pixel_status return reordered_pixel_status From e868560dca3ee8f30994905378b1e268ec5667fe Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 25 Oct 2023 11:31:32 +0200 Subject: [PATCH 19/19] Fill R0 data in case EVB has not applied pe calibration --- src/ctapipe_io_lst/__init__.py | 11 +++++++++-- src/ctapipe_io_lst/tests/test_cta_r1.py | 16 +++++++++++++++- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index 59981de1..763e805a 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -484,8 +484,8 @@ def fill_from_cta_r1(self, array_event, zfits_event): n_samples = zfits_event.num_samples readout_shape = (n_channels, n_pixels, n_samples) - waveform = zfits_event.waveform.reshape(readout_shape) - waveform = waveform.astype(np.float32) / scale - offset + raw_waveform = zfits_event.waveform.reshape(readout_shape) + waveform = raw_waveform.astype(np.float32) / scale - offset reordered_waveform = np.full((n_channels, N_PIXELS, n_samples), 0.0, dtype=np.float32) reordered_waveform[:, pixel_id_map] = waveform @@ -519,6 +519,13 @@ def fill_from_cta_r1(self, array_event, zfits_event): array_event.r1.tel[self.tel_id] = r1 + if DataLevel.R0 in self.datalevels: + reordered_raw_waveform = np.full((n_channels, N_PIXELS, n_samples), 0, dtype=np.uint16) + reordered_raw_waveform[:, pixel_id_map] = raw_waveform + array_event.r0.tel[self.tel_id] = R0CameraContainer( + waveform=reordered_raw_waveform, + ) + def fill_lst_from_ctar1(self, zfits_event): evt = LSTEventContainer( pixel_status=zfits_event.pixel_status.copy(), diff --git a/src/ctapipe_io_lst/tests/test_cta_r1.py b/src/ctapipe_io_lst/tests/test_cta_r1.py index 4379c769..ae378553 100644 --- a/src/ctapipe_io_lst/tests/test_cta_r1.py +++ b/src/ctapipe_io_lst/tests/test_cta_r1.py @@ -16,6 +16,7 @@ from ctapipe_io_lst.anyarray_dtypes import CDTS_AFTER_37201_DTYPE, TIB_DTYPE from ctapipe_io_lst.constants import CLOCK_FREQUENCY_KHZ from ctapipe.image.toymodel import WaveformModel, Gaussian +from ctapipe.containers import EventType import socket from ctapipe_io_lst.event_time import time_to_cta_high @@ -285,8 +286,21 @@ def test_drs4_calibration(dummy_cta_r1): }) with EventSource(dummy_cta_r1, config=config) as source: - n_events = 0 for e in source: + if e.trigger.event_type is EventType.SUBARRAY: + n_channels = 1 + else: + n_channels = 2 + + assert e.r0.tel[1].waveform.dtype == np.uint16 + assert e.r0.tel[1].waveform.shape == (n_channels, 1855, 40) + + 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) + else: + assert e.r1.tel[1].waveform.shape == (2, 1855, 36) + n_events += 1 assert n_events == 100