Merge pull request #179 from cta-observatory/offline_dvr_r1

Add support for reading DVRed data using the offline dvr tool

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

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

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)

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

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

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

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

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: