diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index 61aa8c49..597cf1cb 100644 --- a/src/ctapipe_io_lst/__init__.py +++ b/src/ctapipe_io_lst/__init__.py @@ -52,7 +52,7 @@ PixelStatus, TriggerBits, ) -from .evb_preprocessing import get_processings_for_trigger_bits, EVBPreprocessing +from .evb_preprocessing import get_processings_for_trigger_bits, EVBPreprocessingFlag __all__ = [ @@ -416,10 +416,8 @@ def scheduling_blocks(self): @property def datalevels(self): if self.cta_r1: - if EVBPreprocessing.CALIBRATION in self.evb_preprocessing[TriggerBits.MONO]: + if EVBPreprocessingFlag.PE_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) @@ -555,6 +553,7 @@ def fill_lst_from_ctar1(self, zfits_event): evt.chips_flags = debug.chips_flags evt.charges_hg = debug.charges_gain1 evt.charges_lg = debug.charges_gain2 + evt.tdp_action = debug.tdp_action # unpack Dragon counters counters = debug.counters.view(DRAGON_COUNTERS_DTYPE) @@ -618,6 +617,8 @@ def _generator(self): self.log.warning('Event with event_id=0 found, skipping') continue + + # container for LST data array_event = LSTArrayEventContainer( count=count, @@ -652,11 +653,19 @@ def _generator(self): self.fill_pointing_info(array_event) # apply low level corrections - if self.apply_drs4_corrections: + self.r0_r1_calibrator.update_first_capacitors(array_event) + tdp_action = array_event.lst.tel[self.tel_id].evt.tdp_action + is_calibrated = False + if tdp_action is not None: + tdp_action = EVBPreprocessingFlag(int(tdp_action)) + is_calibrated = EVBPreprocessingFlag.PE_CALIBRATION in tdp_action + + if self.apply_drs4_corrections and not is_calibrated: self.r0_r1_calibrator.apply_drs4_corrections(array_event) - # flat field tagging is performed on r1 data, so can only # be done after the drs4 corrections are applied + # it also assumes uncalibrated data, so cannot be done if EVB + # already calibrated the data if self.use_flatfield_heuristic: self.tag_flatfield_events(array_event) @@ -664,7 +673,7 @@ def _generator(self): self.check_interleaved_pedestal(array_event) # gain select and calibrate to pe - if self.r0_r1_calibrator.calibration_path is not None: + if not is_calibrated and self.r0_r1_calibrator.calibration_path is not None: # skip flatfield and pedestal events if asked if ( self.calibrate_flatfields_and_pedestals @@ -672,6 +681,9 @@ def _generator(self): ): self.r0_r1_calibrator.calibrate(array_event) + # dl1 and drs4 timeshift needs to be filled always + self.r0_r1_calibrator.fill_time_correction(array_event) + yield array_event @staticmethod diff --git a/src/ctapipe_io_lst/calibration.py b/src/ctapipe_io_lst/calibration.py index 2f69fab7..c456bbb7 100644 --- a/src/ctapipe_io_lst/calibration.py +++ b/src/ctapipe_io_lst/calibration.py @@ -222,7 +222,6 @@ def __init__(self, subarray, config=None, parent=None, **kwargs): self.mon_data = self._read_calibration_file(self.calibration_path) def apply_drs4_corrections(self, event: LSTArrayEventContainer): - self.update_first_capacitors(event) for tel_id in event.trigger.tels_with_trigger: r1 = event.r1.tel[tel_id] @@ -277,7 +276,7 @@ def update_first_capacitors(self, event: LSTArrayEventContainer): ) def calibrate(self, event: LSTArrayEventContainer): - for tel_id in event.r0.tel: + for tel_id in event.trigger.tels_with_trigger: r1 = event.r1.tel[tel_id] # if `apply_drs4_corrections` is False, we did not fill in the # waveform yet. @@ -311,6 +310,21 @@ def calibrate(self, event: LSTArrayEventContainer): else: r1.waveform[invalid_pixels[r1.selected_gain_channel, PIXEL_INDEX]] = 0.0 + # needed for charge scaling in ctapipe dl1 calib + if r1.selected_gain_channel is not None: + relative_factor = np.empty(N_PIXELS) + relative_factor[r1.selected_gain_channel == HIGH_GAIN] = self.calib_scale_high_gain.tel[tel_id] + relative_factor[r1.selected_gain_channel == LOW_GAIN] = self.calib_scale_low_gain.tel[tel_id] + else: + relative_factor = np.empty((N_GAINS, N_PIXELS)) + relative_factor[HIGH_GAIN] = self.calib_scale_high_gain.tel[tel_id] + relative_factor[LOW_GAIN] = self.calib_scale_low_gain.tel[tel_id] + + 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 # first drs4 time shift (zeros if no calib file was given) time_shift = self.get_drs4_time_correction( @@ -330,17 +344,6 @@ def calibrate(self, event: LSTArrayEventContainer): event.calibration.tel[tel_id].dl1.time_shift = time_shift - # needed for charge scaling in ctpaipe dl1 calib - if r1.selected_gain_channel is not None: - relative_factor = np.empty(N_PIXELS) - relative_factor[r1.selected_gain_channel == HIGH_GAIN] = self.calib_scale_high_gain.tel[tel_id] - relative_factor[r1.selected_gain_channel == LOW_GAIN] = self.calib_scale_low_gain.tel[tel_id] - else: - relative_factor = np.empty((N_GAINS, N_PIXELS)) - relative_factor[HIGH_GAIN] = self.calib_scale_high_gain.tel[tel_id] - relative_factor[LOW_GAIN] = self.calib_scale_low_gain.tel[tel_id] - - event.calibration.tel[tel_id].dl1.relative_factor = relative_factor @staticmethod def _read_calibration_file(path): diff --git a/src/ctapipe_io_lst/evb_preprocessing.py b/src/ctapipe_io_lst/evb_preprocessing.py index da1f0f8b..8f47a673 100644 --- a/src/ctapipe_io_lst/evb_preprocessing.py +++ b/src/ctapipe_io_lst/evb_preprocessing.py @@ -1,4 +1,4 @@ -from enum import IntEnum +from enum import IntEnum, IntFlag from .constants import TriggerBits from collections import defaultdict @@ -6,20 +6,40 @@ 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. + The values of this Enum is the bit index of this step in the tdp_action value. - See EVB ICD section in + Notes + ----- + This was supposed to be documented in the EVB ICD: https://edms.cern.ch/ui/file/2411710/2.6/LSTMST-ICD-20191206.pdf + But that document doest match the current EVB code. + + As of 2024-01-26, there is a bug in EVB that the ttype_pattern and + tdp_action arrays are actually mixed up in the camera_configuration + object. + """ + # pre-processing flags + GAIN_SELECTION = 0 # PPF0 + BASELINE_SUBTRACTION = 1 # PPF1 + DELTA_T_CORRECTION = 2 # PPF2 + SPIKE_REMOVAL = 3 # PPF3 + + # processing flags + PEDESTAL_SUBTRACTION = 5 # PF0 + PE_CALIBRATION = 6 # PF0 + + +class EVBPreprocessingFlag(IntFlag): + """ + IntFlag version of the EVBPreprocessing, as stored in Event.debug.tdp_action """ - 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 + GAIN_SELECTION = 1 << EVBPreprocessing.GAIN_SELECTION + BASELINE_SUBTRACTION = 1 << EVBPreprocessing.BASELINE_SUBTRACTION + DELTA_T_CORRECTION = 1 << EVBPreprocessing.DELTA_T_CORRECTION + SPIKE_REMOVAL = 1 << EVBPreprocessing.SPIKE_REMOVAL + + PEDESTAL_SUBTRACTION = 1 << EVBPreprocessing.PEDESTAL_SUBTRACTION + PE_CALIBRATION = 1 << EVBPreprocessing.PE_CALIBRATION def get_processings_for_trigger_bits(camera_configuration): @@ -27,21 +47,23 @@ 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} + # EVB has a bug, it stores the tdp_action in the wrong field + tdp_action = camera_configuration.debug.ttype_pattern + + # first entry is default handling + # only look at the first byte for now (maximum 6 bits defied above) + default = EVBPreprocessingFlag(int(tdp_action[0]) & 0xff) actions = defaultdict(lambda: default) - # the following bits refer to the entries in tdp_type + # the following entries refer to the entries in tdp_type + # but with offset of 1, because 0 is the default 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) - } + # only look at the first byte for now (maximum 6 bits defied above) + actions[TriggerBits(int(trigger_bits))] = EVBPreprocessingFlag(int(tdp_action[i]) & 0xff) return actions diff --git a/src/ctapipe_io_lst/tests/test_calib.py b/src/ctapipe_io_lst/tests/test_calib.py index eae32368..0b61d9a3 100644 --- a/src/ctapipe_io_lst/tests/test_calib.py +++ b/src/ctapipe_io_lst/tests/test_calib.py @@ -1,12 +1,14 @@ -import pytest -import pickle -from ctapipe_io_lst.constants import HIGH_GAIN import os +from importlib_resources import files, as_file from pathlib import Path -from traitlets.config import Config +import pickle + import numpy as np +import pytest import tables -from importlib_resources import files, as_file +from traitlets.config import Config + +from ctapipe_io_lst.constants import HIGH_GAIN resource_dir = files('ctapipe_io_lst') / 'tests/resources' @@ -299,3 +301,43 @@ def test_spike_positions(): for key, pos in positions.items(): assert sorted(pos) == sorted(expected_positions[key]) + + +def test_calibrate_precalibrated(): + from ctapipe_io_lst import LSTEventSource + + # file with pe calibration applied by EVB + test_file = "20231219/LST-1.1.Run16255.0000_first50.fits.fz" + + path = test_data / "real/R0" / test_file + config = Config({ + 'LSTEventSource': { + 'pointing_information': False, + 'apply_drs4_corrections': True, + 'LSTR0Corrections': { + 'drs4_pedestal_path': test_drs4_pedestal_path, + 'drs4_time_calibration_path': test_time_calib_path, + 'calibration_path': test_calib_path, + }, + }, + }) + + + previous = None + with LSTEventSource(path, config=config) as source: + n_read = 0 + for event in source: + n_read += 1 + + time_shift = event.calibration.tel[1].dl1.time_shift + # test we filled the timeshift although the data is pre-calibrated + assert time_shift is not None + assert np.any(time_shift != 0) + + # check that time_shift varies from event to event due to drs4 time shift + # regression test for the bug of not updating first_capacitors + if previous is not None and time_shift.shape == previous.shape: + assert np.any(time_shift != previous) + previous = time_shift + + assert n_read == 200 diff --git a/src/ctapipe_io_lst/tests/test_cta_r1.py b/src/ctapipe_io_lst/tests/test_cta_r1.py index 04a53ca7..870fa7ed 100644 --- a/src/ctapipe_io_lst/tests/test_cta_r1.py +++ b/src/ctapipe_io_lst/tests/test_cta_r1.py @@ -1,25 +1,28 @@ import os +import socket from pathlib import Path from contextlib import ExitStack + import pytest import numpy as np -from ctapipe.io import EventSource from astropy.time import Time import astropy.units as u +from traitlets.config import Config + +from ctapipe.io import EventSource +from ctapipe.image.toymodel import WaveformModel, Gaussian +from ctapipe.containers import EventType import protozfits from protozfits.R1v1_pb2 import CameraConfiguration, Event, TelescopeDataStream from protozfits.R1v1_debug_pb2 import DebugEvent, DebugCameraConfiguration from protozfits.CoreMessages_pb2 import AnyArray -from traitlets.config import Config + from ctapipe_io_lst import LSTEventSource 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.constants import CLOCK_FREQUENCY_KHZ, TriggerBits from ctapipe_io_lst.event_time import time_to_cta_high +from ctapipe_io_lst.evb_preprocessing import EVBPreprocessingFlag test_data = Path(os.getenv('LSTCHAIN_TEST_DATA', 'test_data')) @@ -126,6 +129,30 @@ def dummy_cta_r1(dummy_cta_r1_dir): ucts_address = np.frombuffer(socket.inet_pton(socket.AF_INET, "10.0.100.123"), dtype=np.uint32)[0] run_start = Time("2023-05-16T16:06:31.123") + + # trigger dependent processing definition + # tdp_type is a list of TriggerBits + # tdp_action is the bit pattern of pre-processing steps to apply + # to the corresponding tdp_type. + tdp_type = np.zeros(15, dtype=np.uint16) + tdp_action = np.zeros(16, dtype=np.uint16) + + # tdp_action[0] is the default for unknown trigger bits + tdp_action[0] = EVBPreprocessingFlag.BASELINE_SUBTRACTION | EVBPreprocessingFlag.GAIN_SELECTION + # tdp_action[i] corresponds to tdp_type[i - 1] + # physics (mono for now) + tdp_type[0] = TriggerBits.MONO + tdp_action[1] = EVBPreprocessingFlag.BASELINE_SUBTRACTION | EVBPreprocessingFlag.GAIN_SELECTION + # flat field events + tdp_type[1] = TriggerBits.CALIBRATION + tdp_action[2] = EVBPreprocessingFlag.BASELINE_SUBTRACTION + tdp_type[2] = TriggerBits.CALIBRATION | TriggerBits.MONO + tdp_action[3] = EVBPreprocessingFlag.BASELINE_SUBTRACTION + # pedestal events + tdp_type[3] = TriggerBits.PEDESTAL + tdp_action[4] = EVBPreprocessingFlag.BASELINE_SUBTRACTION + + camera_config = CameraConfiguration( tel_id=1, local_run_id=10000, @@ -138,12 +165,15 @@ def dummy_cta_r1(dummy_cta_r1_dir): num_samples_nominal=num_samples, pixel_id_map=to_anyarray(np.arange(num_pixels).astype(np.uint16)), module_id_map=to_anyarray(np.arange(num_modules).astype(np.uint16)), + debug=DebugCameraConfiguration( cs_serial="???", evb_version="evb-dummy", cdhs_version="evb-dummy", - tdp_type=to_anyarray(np.zeros(15, dtype=np.uint16)), + tdp_type=to_anyarray(tdp_type), tdp_action=to_anyarray(np.zeros(16, dtype=np.uint16)), + # at the moment, the tdp_action and ttype_pattern fields are mixed up in EVB + ttype_pattern=to_anyarray(tdp_action), ) ) @@ -304,3 +334,27 @@ def test_drs4_calibration(dummy_cta_r1): n_events += 1 assert n_events == 100 + + +test_files = [ + "20231214/LST-1.1.Run16102.0000_first50.fits.fz", # has only baseline enabled + "20231218/LST-1.1.Run16231.0000_first50.fits.fz", # baseline + gain selection for physics + "20231219/LST-1.1.Run16255.0000_first50.fits.fz", # all corrections + gain selection for physics +] + + +@pytest.mark.parametrize("test_file", test_files) +def test_read_real_files(test_file): + config = Config({ + 'LSTEventSource': { + 'pointing_information': False, + 'apply_drs4_corrections': False, + }, + }) + + with EventSource(test_data / "real/R0" / test_file, config=config) as source: + n_read = 0 + for _ in source: + n_read += 1 + + assert n_read == 200 diff --git a/src/ctapipe_io_lst/tests/test_evb_preprocessing.py b/src/ctapipe_io_lst/tests/test_evb_preprocessing.py new file mode 100644 index 00000000..ea97f1af --- /dev/null +++ b/src/ctapipe_io_lst/tests/test_evb_preprocessing.py @@ -0,0 +1,56 @@ +import os +from pathlib import Path +from protozfits import File +from ctapipe_io_lst import TriggerBits +from ctapipe_io_lst.evb_preprocessing import EVBPreprocessingFlag + +import pytest + + +test_data = Path(os.getenv('LSTCHAIN_TEST_DATA', 'test_data')).absolute() + +test_files = [ + "20231214/LST-1.1.Run16102.0000_first50.fits.fz", # has only baseline enabled + "20231218/LST-1.1.Run16231.0000_first50.fits.fz", # baseline + gain selection for physics + "20231219/LST-1.1.Run16255.0000_first50.fits.fz", # all corrections + gain selection for physics +] + +all_corrections = EVBPreprocessingFlag.BASELINE_SUBTRACTION | EVBPreprocessingFlag.DELTA_T_CORRECTION | EVBPreprocessingFlag.PE_CALIBRATION | EVBPreprocessingFlag.PEDESTAL_SUBTRACTION + +expected = [ + { + TriggerBits.MONO: EVBPreprocessingFlag.BASELINE_SUBTRACTION, + TriggerBits.STEREO: EVBPreprocessingFlag.BASELINE_SUBTRACTION, + TriggerBits.MONO|TriggerBits.PEDESTAL: EVBPreprocessingFlag.BASELINE_SUBTRACTION, + TriggerBits.MONO|TriggerBits.CALIBRATION: EVBPreprocessingFlag.BASELINE_SUBTRACTION, + TriggerBits.PEDESTAL: EVBPreprocessingFlag.BASELINE_SUBTRACTION, + TriggerBits.CALIBRATION: EVBPreprocessingFlag.BASELINE_SUBTRACTION, + }, + { + TriggerBits.MONO: EVBPreprocessingFlag.BASELINE_SUBTRACTION | EVBPreprocessingFlag.GAIN_SELECTION, + TriggerBits.STEREO: EVBPreprocessingFlag.BASELINE_SUBTRACTION | EVBPreprocessingFlag.GAIN_SELECTION, + TriggerBits.MONO|TriggerBits.PEDESTAL: EVBPreprocessingFlag.BASELINE_SUBTRACTION | EVBPreprocessingFlag.GAIN_SELECTION, + TriggerBits.MONO|TriggerBits.CALIBRATION: EVBPreprocessingFlag.BASELINE_SUBTRACTION, + TriggerBits.PEDESTAL: EVBPreprocessingFlag.BASELINE_SUBTRACTION, + TriggerBits.CALIBRATION: EVBPreprocessingFlag.BASELINE_SUBTRACTION, + }, + { + TriggerBits.MONO: all_corrections | EVBPreprocessingFlag.GAIN_SELECTION, + TriggerBits.STEREO: all_corrections| EVBPreprocessingFlag.GAIN_SELECTION, + TriggerBits.MONO|TriggerBits.PEDESTAL: all_corrections | EVBPreprocessingFlag.GAIN_SELECTION, + TriggerBits.MONO|TriggerBits.CALIBRATION: all_corrections, + TriggerBits.PEDESTAL: all_corrections, + TriggerBits.CALIBRATION: all_corrections, + }, +] +@pytest.mark.parametrize(("test_file", "expected"), zip(test_files, expected)) +def test_get_processings_for_trigger_bits(test_file, expected): + from ctapipe_io_lst.evb_preprocessing import get_processings_for_trigger_bits + + path = test_data / "real/R0/" / test_file + + with File(str(path)) as f: + camera_config = f.CameraConfiguration[0] + + result = get_processings_for_trigger_bits(camera_config) + assert result == expected diff --git a/src/ctapipe_io_lst/tests/test_lsteventsource.py b/src/ctapipe_io_lst/tests/test_lsteventsource.py index 85f5bc03..424077f4 100644 --- a/src/ctapipe_io_lst/tests/test_lsteventsource.py +++ b/src/ctapipe_io_lst/tests/test_lsteventsource.py @@ -24,6 +24,11 @@ test_drive_report = test_data / 'real/monitoring/DrivePositioning/DrivePosition_log_20200218.txt' +calib_version = "ctapipe-v0.17" +calib_path = test_data / 'real/monitoring/PixelCalibration/Cat-A/' +test_calib_path = calib_path / f'calibration/20200218/{calib_version}/calibration_filters_52.Run02006.0000.h5' +test_time_calib_path = calib_path / f'drs4_time_sampling_from_FF/20191124/{calib_version}/time_calibration.Run01625.0000.h5' + # ADC_SAMPLES_SHAPE = (2, 14, 40) @@ -358,8 +363,6 @@ def test_reference_position(): assert u.isclose(position.height, LST1_LOCATION.height) - - @pytest.mark.parametrize("timeshift", (-5, 70)) def test_time_correction(timeshift): from ctapipe_io_lst import LSTEventSource @@ -381,3 +384,26 @@ def test_time_correction(timeshift): dt_tel = event_shifted.trigger.tel[1].time - event.trigger.tel[1].time assert u.isclose(dt.to_value(u.s), timeshift) assert u.isclose(dt_tel.to_value(u.s), timeshift) + + +def test_evb_calibrated_data(): + from ctapipe_io_lst import LSTEventSource + input_url = test_data / 'real/R0/20231219/LST-1.1.Run16255.0000_first50.fits.fz' + + config = { + 'LSTEventSource': { + "pointing_information": False, + 'LSTR0Corrections': { + 'drs4_time_calibration_path': str(test_time_calib_path), + 'calibration_path': str(test_calib_path), + }, + }, + } + + with LSTEventSource(input_url, config=Config(config)) as source: + read_events = 0 + for e in source: + read_events += 1 + assert np.all(e.calibration.tel[1].dl1.time_shift != 0) + + assert read_events == 200