diff --git a/src/ctapipe_io_lst/__init__.py b/src/ctapipe_io_lst/__init__.py index ee03bc1f..08d084bb 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, EVBPreprocessing, 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 EVBPreprocessing.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) @@ -652,11 +651,18 @@ def _generator(self): self.fill_pointing_info(array_event) # apply low level corrections - if self.apply_drs4_corrections: - self.r0_r1_calibrator.apply_drs4_corrections(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 +670,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 +678,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..9b221d17 100644 --- a/src/ctapipe_io_lst/calibration.py +++ b/src/ctapipe_io_lst/calibration.py @@ -277,7 +277,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 +311,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 +345,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 c02169ad..22af518e 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,26 +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 index of this step in the ttype_pattern array. - Note - ---- + 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 - RESERVED1 = 4 # PPF4 # processing flags PEDESTAL_SUBTRACTION = 5 # PF0 PE_CALIBRATION = 6 # PF0 - RESERVED2 = 7 # PF0 - RESERVED3 = 8 # PF0 + + +class EVBPreprocessingFlag(IntFlag): + """ + IntFlag version of the EVBPreprocessing, as stored in Event.debug.tdp_action + """ + 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): 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