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 9a633efa..c64f9652 100644 --- a/setup.cfg +++ b/setup.cfg @@ -18,20 +18,20 @@ 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 - ctapipe~=0.19.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 1f814c90..763e805a 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -2,11 +2,12 @@ """ 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 ( ReflectorShape, @@ -25,7 +26,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 @@ -50,11 +52,20 @@ PixelStatus, TriggerBits, ) +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 @@ -93,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 @@ -118,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 @@ -133,6 +141,16 @@ 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): + if set_dvr_bits: + not_broken = get_channel_info(pixel_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 + + class LSTEventSource(EventSource): """ EventSource for LST R0 data. @@ -300,6 +318,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 @@ -309,6 +329,10 @@ 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) self._subarray = self.create_subarray(self.tel_id, reference_location) self.r0_r1_calibrator = LSTR0Corrections( @@ -381,6 +405,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, ) @@ -437,24 +467,37 @@ 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 + + # reorder to nominal pixel order + 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 + n_samples = zfits_event.num_samples - # TODO: decide from applied preprocessing options, not gains + readout_shape = (n_channels, n_pixels, n_samples) + 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 + waveform = reordered_waveform + + # FIXME, check using evb_preprocessing and make ctapipe support 2 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), - ) - array_event.r1.tel[self.tel_id] = R1CameraContainer() + 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) - - shape = (zfits_event.num_pixels, zfits_event.num_samples) - array_event.r1.tel[self.tel_id] = R1CameraContainer( - waveform=zfits_event.waveform.reshape(shape), - selected_gain_channel=selected_gain_channel, - ) - array_event.r0.tel[self.tel_id] = R0CameraContainer() + waveform = waveform[0] array_event.lst.tel[self.tel_id] = self.fill_lst_from_ctar1(zfits_event) @@ -462,15 +505,38 @@ 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) + + r1 = R1CameraContainer( + waveform=waveform, + selected_gain_channel=selected_gain_channel, + ) + + 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 + + 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, + 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 @@ -530,28 +596,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: @@ -986,6 +1055,16 @@ 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 = _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, + ) + return r0, r1 def fill_r0r1_container(self, array_event, zfits_event): @@ -997,15 +1076,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( @@ -1013,7 +1089,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/calibration.py b/src/ctapipe_io_lst/calibration.py index 845b00fa..2f69fab7 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. @@ -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': @@ -257,7 +259,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: @@ -480,6 +482,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 +1010,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, diff --git a/src/ctapipe_io_lst/evb_preprocessing.py b/src/ctapipe_io_lst/evb_preprocessing.py new file mode 100644 index 00000000..da1f0f8b --- /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 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 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) diff --git a/src/ctapipe_io_lst/tests/test_cta_r1.py b/src/ctapipe_io_lst/tests/test_cta_r1.py index 7ebc4d29..ae378553 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 @@ -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 @@ -141,9 +142,14 @@ 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)), ) ) + # 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 +157,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") @@ -278,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