Merge pull request #209 from cta-observatory/automatic_evb_corrections

Fix definition of evb preprocessing flags, apply calibration only when needed

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,26 +653,37 @@ 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)
 
             if self.pedestal_ids is not None:
                 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
                     or array_event.trigger.event_type not in {EventType.FLATFIELD, EventType.SKY_PEDESTAL}
                 ):
                     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

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):

src/ctapipe_io_lst/evb_preprocessing.py

@@ -1,47 +1,69 @@
-from enum import IntEnum
+from enum import IntEnum, IntFlag
 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.
+    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):
     """
     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

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