diff --git a/environment.yml b/environment.yml index ebce1d7f..d97b4c72 100644 --- a/environment.yml +++ b/environment.yml @@ -4,7 +4,7 @@ channels: - default dependencies: - astropy=5.2 - - python=3.10 # nail the python version, so conda does not try upgrading / dowgrading + - python=3.11 # nail the python version, so conda does not try upgrading / dowgrading - ctapipe=0.19 - eventio - corsikaio diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index f4e1349a..e99d1e61 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -17,7 +17,6 @@ OpticsDescription, SizeType, ) -from enum import IntFlag, auto from astropy.time import Time from ctapipe.io import EventSource, read_table @@ -46,6 +45,7 @@ ) from .constants import ( HIGH_GAIN, LST_LOCATIONS, N_GAINS, N_PIXELS, N_SAMPLES, LST1_LOCATION, REFERENCE_LOCATION, + PixelStatus, TriggerBits, ) @@ -56,40 +56,6 @@ NO_FF_HEURISTIC_DATE = Time("2022-01-01T00:00:00") -class TriggerBits(IntFlag): - ''' - See TIB User manual - ''' - UNKNOWN = 0 - MONO = auto() - STEREO = auto() - CALIBRATION = auto() - SINGLE_PE = auto() - SOFTWARE = auto() - PEDESTAL = auto() - SLOW_CONTROL = auto() - - PHYSICS = MONO | STEREO - OTHER = CALIBRATION | SINGLE_PE | SOFTWARE | PEDESTAL | SLOW_CONTROL - - -class PixelStatus(IntFlag): - ''' - Pixel status information - - See Section A.5 of the CTA R1 Data Model: - https://forge.in2p3.fr/dmsf/files/8627 - ''' - RESERVED_0 = auto() - RESERVED_1 = auto() - HIGH_GAIN_STORED = auto() - LOW_GAIN_STORED = auto() - SATURATED = auto() - PIXEL_TRIGGER_1 = auto() - PIXEL_TRIGGER_2 = auto() - PIXEL_TRIGGER_3 = auto() - - BOTH_GAINS_STORED = HIGH_GAIN_STORED | LOW_GAIN_STORED #: LST Optics Description OPTICS = OpticsDescription( @@ -117,7 +83,7 @@ def get_channel_info(pixel_status): 2: low-gain read out 3: both gains read out ''' - return (pixel_status & 0b1100) >> 2 + return (pixel_status & PixelStatus.BOTH_GAINS_STORED) >> 2 def load_camera_geometry(): @@ -132,7 +98,6 @@ def load_camera_geometry(): def read_pulse_shapes(): - ''' Reads in the data on the pulse shapes and readout speed, from an external file @@ -315,6 +280,7 @@ def __init__(self, input_url=None, **kwargs): self.run_id = self.camera_config.configuration_id self.tel_id = self.camera_config.telescope_id self.run_start = Time(self.camera_config.date, format='unix') + self.dvr_applied = self.multi_file.dvr_applied reference_location = EarthLocation( lon=self.reference_position_lon * u.deg, @@ -572,7 +538,15 @@ def fill_lst_event_container(self, array_event, zfits_event): lst_evt.configuration_id = zfits_event.configuration_id lst_evt.event_id = zfits_event.event_id lst_evt.tel_event_id = zfits_event.tel_event_id - lst_evt.pixel_status = zfits_event.pixel_status + + lst_evt.pixel_status = np.zeros(N_PIXELS, dtype=zfits_event.pixel_status.dtype) + lst_evt.pixel_status[self.camera_config.expected_pixels_id] = zfits_event.pixel_status + + # set bits for dvr if not already set + if not self.dvr_applied: + not_broken = get_channel_info(lst_evt.pixel_status) != 0 + lst_evt.pixel_status[not_broken] |= PixelStatus.DVR_STATUS_0 + lst_evt.ped_id = zfits_event.ped_id lst_evt.module_status = zfits_event.lstcam.module_status lst_evt.extdevices_presence = zfits_event.lstcam.extdevices_presence @@ -793,6 +767,11 @@ def fill_r0r1_camera_container(self, zfits_event): dtype = zfits_event.waveform.dtype fill = np.iinfo(dtype).max + if self.dvr_applied: + stored_pixels = (zfits_event.pixel_status & PixelStatus.DVR_STATUS) > 0 + else: + stored_pixels = slice(None) # all pixels stored + # we assume that either all pixels are gain selected or none # only broken pixels are allowed to be missing completely if gain_selected: @@ -813,7 +792,7 @@ def fill_r0r1_camera_container(self, zfits_event): ) reordered_waveform = np.full((N_PIXELS, N_SAMPLES), fill, dtype=dtype) - reordered_waveform[expected_pixels] = waveform + reordered_waveform[expected_pixels[stored_pixels]] = waveform reordered_selected_gain = np.full(N_PIXELS, -1, dtype=np.int8) reordered_selected_gain[expected_pixels] = selected_gain @@ -824,11 +803,10 @@ def fill_r0r1_camera_container(self, zfits_event): selected_gain_channel=reordered_selected_gain, ) else: - reshaped_waveform = zfits_event.waveform.reshape(N_GAINS, n_pixels, n_samples) + reshaped_waveform = zfits_event.waveform.reshape(N_GAINS, -1, n_samples) # re-order the waveform following the expected_pixels_id values - # could also just do waveform = reshaped_waveform[np.argsort(expected_ids)] reordered_waveform = np.full((N_GAINS, N_PIXELS, N_SAMPLES), fill, dtype=dtype) - reordered_waveform[:, expected_pixels, :] = reshaped_waveform + reordered_waveform[:, expected_pixels[stored_pixels], :] = reshaped_waveform r0 = R0CameraContainer(waveform=reordered_waveform) r1 = R1CameraContainer() diff --git a/src/ctapipe_io_lst/calibration.py b/src/ctapipe_io_lst/calibration.py index 4c0b897b..845b00fa 100644 --- a/src/ctapipe_io_lst/calibration.py +++ b/src/ctapipe_io_lst/calibration.py @@ -14,16 +14,17 @@ from ctapipe.calib.camera.gainselection import ThresholdGainSelector from ctapipe.containers import FlatFieldContainer, MonitoringCameraContainer, MonitoringContainer, PedestalContainer, PixelStatusContainer, WaveformCalibrationContainer from ctapipe.io import HDF5TableReader, read_table -from .containers import LSTArrayEventContainer from traitlets import Enum +from .containers import LSTArrayEventContainer + from .constants import ( N_GAINS, N_PIXELS, N_MODULES, N_SAMPLES, LOW_GAIN, HIGH_GAIN, N_PIXELS_MODULE, N_CAPACITORS_PIXEL, N_CAPACITORS_CHANNEL, CLOCK_FREQUENCY_KHZ, CHANNEL_ORDER_LOW_GAIN, CHANNEL_ORDER_HIGH_GAIN, N_CHANNELS_MODULE, - PIXEL_INDEX, + PIXEL_INDEX, PixelStatus ) __all__ = [ @@ -255,12 +256,14 @@ def apply_drs4_corrections(self, event: LSTArrayEventContainer): if self.offset.tel[tel_id] != 0: r1.waveform -= self.offset.tel[tel_id] - mon = event.mon.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 + invalid_pixels = broken_pixels | dvred_pixels + if r1.selected_gain_channel is None: - r1.waveform[mon.pixel_status.hardware_failing_pixels] = 0.0 + r1.waveform[invalid_pixels] = 0.0 else: - broken = mon.pixel_status.hardware_failing_pixels[r1.selected_gain_channel, PIXEL_INDEX] - r1.waveform[broken] = 0.0 + r1.waveform[invalid_pixels[r1.selected_gain_channel, PIXEL_INDEX]] = 0.0 def update_first_capacitors(self, event: LSTArrayEventContainer): @@ -298,10 +301,13 @@ def calibrate(self, event: LSTArrayEventContainer): ) 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 + invalid_pixels = broken_pixels | dvred_pixels + if r1.selected_gain_channel is None: - r1.waveform[broken_pixels] = 0.0 + r1.waveform[invalid_pixels] = 0.0 else: - r1.waveform[broken_pixels[r1.selected_gain_channel, PIXEL_INDEX]] = 0.0 + r1.waveform[invalid_pixels[r1.selected_gain_channel, PIXEL_INDEX]] = 0.0 # store calibration data needed for dl1 calibration in ctapipe # first drs4 time shift (zeros if no calib file was given) diff --git a/src/ctapipe_io_lst/constants.py b/src/ctapipe_io_lst/constants.py index fb7d4838..4d0ff749 100644 --- a/src/ctapipe_io_lst/constants.py +++ b/src/ctapipe_io_lst/constants.py @@ -1,6 +1,7 @@ import numpy as np import astropy.units as u from astropy.coordinates import EarthLocation +from enum import IntFlag, auto N_GAINS = 2 N_MODULES = 265 @@ -43,3 +44,40 @@ LST_LOCATIONS = { 1: LST1_LOCATION, } + + +class TriggerBits(IntFlag): + ''' + See TIB User manual + ''' + UNKNOWN = 0 + MONO = auto() + STEREO = auto() + CALIBRATION = auto() + SINGLE_PE = auto() + SOFTWARE = auto() + PEDESTAL = auto() + SLOW_CONTROL = auto() + + PHYSICS = MONO | STEREO + OTHER = CALIBRATION | SINGLE_PE | SOFTWARE | PEDESTAL | SLOW_CONTROL + + +class PixelStatus(IntFlag): + ''' + Pixel status information + + See Section A.5 of the CTA R1 Data Model: + https://forge.in2p3.fr/dmsf/files/8627 + ''' + DVR_STATUS_0 = auto() + DVR_STATUS_1 = auto() + HIGH_GAIN_STORED = auto() + LOW_GAIN_STORED = auto() + SATURATED = auto() + PIXEL_TRIGGER_1 = auto() + PIXEL_TRIGGER_2 = auto() + PIXEL_TRIGGER_3 = auto() + + BOTH_GAINS_STORED = HIGH_GAIN_STORED | LOW_GAIN_STORED + DVR_STATUS = DVR_STATUS_0 | DVR_STATUS_1 diff --git a/src/ctapipe_io_lst/containers.py b/src/ctapipe_io_lst/containers.py index 08fb9bd5..87e4c16f 100644 --- a/src/ctapipe_io_lst/containers.py +++ b/src/ctapipe_io_lst/containers.py @@ -1,10 +1,9 @@ """ Container structures for data that should be read or written to disk """ -from astropy import units as u -from numpy import nan from ctapipe.core import Container, Field, Map from ctapipe.containers import ArrayEventContainer +from functools import partial __all__ = [ @@ -99,8 +98,8 @@ class LSTCameraContainer(Container): """ Container for Fields that are specific to each LST camera """ - evt = Field(LSTEventContainer(), "LST specific event Information") - svc = Field(LSTServiceContainer(), "LST specific camera_config Information") + evt = Field(default_factory=LSTEventContainer, description="LST specific event Information") + svc = Field(default_factory=LSTServiceContainer, description="LST specific camera_config Information") class LSTContainer(Container): @@ -110,12 +109,13 @@ class LSTContainer(Container): # create the camera container tel = Field( - Map(LSTCameraContainer), - "map of tel_id to LSTTelContainer") + default_factory=partial(Map, LSTCameraContainer), + description="map of tel_id to LSTTelContainer" + ) class LSTArrayEventContainer(ArrayEventContainer): """ Data container including LST and monitoring information """ - lst = Field(LSTContainer(), "LST specific Information") + lst = Field(default_factory=LSTContainer, description="LST specific Information") diff --git a/src/ctapipe_io_lst/multifiles.py b/src/ctapipe_io_lst/multifiles.py index 8d2f6c7d..69003be3 100644 --- a/src/ctapipe_io_lst/multifiles.py +++ b/src/ctapipe_io_lst/multifiles.py @@ -105,7 +105,9 @@ def __init__(self, path, *args, **kwargs): self._files = {} self._events = PriorityQueue() self._events_tables = {} + self._headers = {} self.camera_config = None + self.dvr_applied = None if self.all_streams and file_info is not None: for stream in count(1): @@ -161,6 +163,12 @@ def _load_next_subrun(self, stream): self._files[stream] = File(str(path)) self.log.info("Opened file %s", path) self._events_tables[stream] = self._files[stream].Events + self._headers[stream] = self._events_tables[stream].header + dvr_applied = self._headers[stream].get("LSTDVR", False) + if self.dvr_applied is None: + self.dvr_applied = dvr_applied + elif dvr_applied != self.dvr_applied: + raise IOError("Mixing subruns / streams with and without DVR applied is not supported") # load first event from each stream event = next(self._events_tables[stream]) diff --git a/src/ctapipe_io_lst/tests/test_lsteventsource.py b/src/ctapipe_io_lst/tests/test_lsteventsource.py index 34a01c4c..814f9c27 100644 --- a/src/ctapipe_io_lst/tests/test_lsteventsource.py +++ b/src/ctapipe_io_lst/tests/test_lsteventsource.py @@ -12,10 +12,11 @@ from ctapipe.calib.camera.gainselection import ThresholdGainSelector from ctapipe_io_lst.constants import N_GAINS, N_PIXELS_MODULE, N_SAMPLES, N_PIXELS -from ctapipe_io_lst import TriggerBits +from ctapipe_io_lst import TriggerBits, PixelStatus test_data = Path(os.getenv('LSTCHAIN_TEST_DATA', 'test_data')).absolute() test_r0_dir = test_data / 'real/R0/20200218' +test_r0_dvr_dir = test_data / 'real/R0DVR' test_r0_path = test_r0_dir / 'LST-1.1.Run02006.0004.fits.fz' test_r0_path_all_streams = test_r0_dir / 'LST-1.1.Run02008.0000_first50.fits.fz' @@ -184,6 +185,57 @@ def test_gain_selected(): assert event.count == 199 + +def test_dvr(): + from ctapipe_io_lst import LSTEventSource + + config = Config(dict( + LSTEventSource=dict( + default_trigger_type='tib', # ucts unreliable in this run + apply_drs4_corrections=True, + pointing_information=False, + use_flatfield_heuristic=True, + LSTR0Corrections=dict( + apply_drs4_pedestal_correction=False, + apply_spike_correction=False, + apply_timelapse_correction=False, + offset=400, + ) + ) + )) + + dvr_source = LSTEventSource( + test_r0_dvr_dir / 'LST-1.1.Run02008.0100_first50.fits.fz', + config=config, + ) + original_source = LSTEventSource( + test_r0_dir / 'LST-1.1.Run02008.0100_first50.fits.fz', + config=config, + ) + gain_selector = ThresholdGainSelector(threshold=3500) + for dvr_event, original_event in zip(dvr_source, original_source): + if dvr_event.trigger.event_type in {EventType.FLATFIELD, EventType.SKY_PEDESTAL}: + assert 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 is not None + assert dvr_event.r1.tel[1].waveform.shape == (N_GAINS, N_PIXELS, N_SAMPLES - 4) + else: + 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) + + # 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]) + + assert dvr_event.count == 199 + + def test_pointing_info(): from ctapipe_io_lst import LSTEventSource diff --git a/src/ctapipe_io_lst/tests/test_multifile.py b/src/ctapipe_io_lst/tests/test_multifile.py index 80e47a1e..f10bbf0d 100644 --- a/src/ctapipe_io_lst/tests/test_multifile.py +++ b/src/ctapipe_io_lst/tests/test_multifile.py @@ -13,6 +13,7 @@ def test_multifile_streams(): with MultiFiles(path) as multi_files: assert multi_files.n_open_files == 4 + assert multi_files.dvr_applied is False event_count = 0 for event in multi_files: