Merge pull request #2338 from cta-observatory/ctao_dl0

Add fields defined in CTAO DL0 and R1 data models

ctapipe/calib/camera/calibrator.py

@@ -9,7 +9,7 @@
 import numpy as np
 from numba import float32, float64, guvectorize, int64
 
-from ctapipe.containers import DL1CameraContainer
+from ctapipe.containers import DL0CameraContainer, DL1CameraContainer, PixelStatus
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import (
     BoolTelescopeParameter,
@@ -178,19 +178,35 @@ def _check_dl0_empty(self, waveforms):
             return False
 
     def _calibrate_dl0(self, event, tel_id):
-        waveforms = event.r1.tel[tel_id].waveform
-        selected_gain_channel = event.r1.tel[tel_id].selected_gain_channel
-        if self._check_r1_empty(waveforms):
+        r1 = event.r1.tel[tel_id]
+
+        if self._check_r1_empty(r1.waveform):
             return
 
-        reduced_waveforms_mask = self.data_volume_reducer(
-            waveforms, tel_id=tel_id, selected_gain_channel=selected_gain_channel
+        signal_pixels = self.data_volume_reducer(
+            r1.waveform,
+            tel_id=tel_id,
+            selected_gain_channel=r1.selected_gain_channel,
         )
 
-        waveforms_copy = waveforms.copy()
-        waveforms_copy[~reduced_waveforms_mask] = 0
-        event.dl0.tel[tel_id].waveform = waveforms_copy
-        event.dl0.tel[tel_id].selected_gain_channel = selected_gain_channel
+        dl0_waveform = r1.waveform.copy()
+        dl0_waveform[~signal_pixels] = 0
+
+        dl0_pixel_status = r1.pixel_status.copy()
+        # set dvr pixel bit in pixel_status for pixels kept by DVR
+        dl0_pixel_status[signal_pixels] |= PixelStatus.DVR_STORED_AS_SIGNAL
+        # unset dvr bits for removed pixels
+        dl0_pixel_status[~signal_pixels] &= ~np.uint8(PixelStatus.DVR_STATUS)
+
+        event.dl0.tel[tel_id] = DL0CameraContainer(
+            event_type=r1.event_type,
+            event_time=r1.event_time,
+            waveform=dl0_waveform,
+            selected_gain_channel=r1.selected_gain_channel,
+            pixel_status=dl0_pixel_status,
+            first_cell_id=r1.first_cell_id,
+            calibration_monitoring_id=r1.calibration_monitoring_id,
+        )
 
     def _calibrate_dl1(self, event, tel_id):
         waveforms = event.dl0.tel[tel_id].waveform

ctapipe/containers.py

@@ -165,6 +165,70 @@ class EventType(enum.Enum):
     UNKNOWN = 255
 
 
+class PixelStatus(enum.IntFlag):
+    """
+    Pixel status information
+
+    See DL0 Data Model specification:
+    https://redmine.cta-observatory.org/dmsf/files/17552/view
+    """
+
+    DVR_STORED_AS_SIGNAL = enum.auto()
+    DVR_STORED_NO_SIGNAL = enum.auto()
+    HIGH_GAIN_STORED = enum.auto()
+    LOW_GAIN_STORED = enum.auto()
+    SATURATED = enum.auto()
+    PIXEL_TRIGGER_0 = enum.auto()
+    PIXEL_TRIGGER_1 = enum.auto()
+    PIXEL_TRIGGER_2 = enum.auto()
+
+    #: DVR status uses two bits
+    #: 0 = not stored, 1 = identified as signal, 2 = stored, not identified as signal
+    DVR_STATUS = DVR_STORED_AS_SIGNAL | DVR_STORED_NO_SIGNAL
+
+    #: Pixel trigger information, TBD
+    PIXEL_TRIGGER = PIXEL_TRIGGER_0 | PIXEL_TRIGGER_1 | PIXEL_TRIGGER_2
+
+    @staticmethod
+    def get_dvr_status(pixel_status):
+        """
+        Return only the bits corresponding to the DVR_STATUS
+
+        Returns
+        -------
+        dvr_status: int or array[uint8]
+            0 = pixel not stored
+            1 = pixel was identified as signal pixel and stored
+            2 = pixel was stored, but not identified as signal
+        """
+        return pixel_status & PixelStatus.DVR_STATUS
+
+    @staticmethod
+    def get_channel_info(pixel_status):
+        """
+        Return only the bits corresponding to the channel info (high/low gain stored)
+
+        Returns
+        -------
+        channel_info: int or array[uint8]
+            0 = pixel broken / disabled
+            1 = only high gain read out
+            2 = only low gain read out
+            3 = both gains read out
+        """
+        gain_bits = PixelStatus.HIGH_GAIN_STORED | PixelStatus.LOW_GAIN_STORED
+        return (pixel_status & gain_bits) >> 2
+
+    @staticmethod
+    def is_invalid(pixel_status):
+        """Return if pixel values are marked as invalid
+
+        This is encoded in the data model as neither high gain nor low gain marked as stored
+        """
+        gain_bits = PixelStatus.HIGH_GAIN_STORED | PixelStatus.LOW_GAIN_STORED
+        return (pixel_status & gain_bits) == 0
+
+
 class EventIndexContainer(Container):
     """index columns to include in event lists, common to all data levels"""
 
@@ -503,17 +567,55 @@ class R1CameraContainer(Container):
     Storage of r1 calibrated data from a single telescope
     """
 
+    event_type = Field(EventType.UNKNOWN, "type of event", type=EventType)
+    event_time = Field(NAN_TIME, "event timestamp")
+
     waveform = Field(
         None,
         (
             "numpy array containing a set of images, one per ADC sample"
             "Shape: (n_pixels, n_samples)"
         ),
     )
+
+    pixel_status = Field(
+        None,
+        "Array of pixel status values, see PixelStatus for definition of the values",
+        ndim=1,
+        dtype=np.uint8,
+    )
+
+    first_cell_id = Field(
+        None,
+        "Array of first cell ids of the readout chips. Only used by LST and SST.",
+        dtype=np.uint16,
+        ndim=1,
+    )
+
+    module_hires_local_clock_counter = Field(
+        None,
+        "Clock counter values of the camera modules. See R1 data model for details.",
+        dtype=np.uint64,
+        ndim=1,
+    )
+
+    pedestal_intensity = Field(
+        None,
+        "Pedestal intensity in each pixel in DC",
+        dtype=np.float32,
+        ndim=1,
+    )
+
+    calibration_monitoring_id = Field(
+        None,
+        "ID of the CalibrationMonitoringSet containing the applied pre-calibration parameters",
+    )
+
     selected_gain_channel = Field(
         None,
         (
             "Numpy array containing the gain channel chosen for each pixel. "
+            "Note: should be replaced by using ``pixel_status`` "
             "Shape: (n_pixels)"
         ),
     )
@@ -532,9 +634,15 @@ class R1Container(Container):
 
 class DL0CameraContainer(Container):
     """
-    Storage of data volume reduced dl0 data from a single telescope
+    Storage of data volume reduced dl0 data from a single telescope.
+
+    See DL0 Data Model specification:
+    https://redmine.cta-observatory.org/dmsf/files/17552/view
     """
 
+    event_type = Field(EventType.UNKNOWN, "type of event", type=EventType)
+    event_time = Field(NAN_TIME, "event timestamp")
+
     waveform = Field(
         None,
         (
@@ -545,10 +653,30 @@ class DL0CameraContainer(Container):
         ),
     )
 
+    pixel_status = Field(
+        None,
+        "Array of pixel status values, see PixelStatus for definition of the values",
+        dtype=np.uint8,
+        ndim=1,
+    )
+
+    first_cell_id = Field(
+        None,
+        "Array of first cell ids of the readout chips. Only used by LST and SST.",
+        dtype=np.uint16,
+        ndim=1,
+    )
+
+    calibration_monitoring_id = Field(
+        None,
+        "ID of the CalibrationMonitoringSet containing the applied pre-calibration parameters",
+    )
+
     selected_gain_channel = Field(
         None,
         (
             "Numpy array containing the gain channel chosen for each pixel. "
+            "Note: this should be replaced by only using ``pixel_status`` "
             "Shape: (n_pixels)"
         ),
     )

ctapipe/io/simteleventsource.py

@@ -20,7 +20,7 @@
     FiveLayerAtmosphereDensityProfile,
     TableAtmosphereDensityProfile,
 )
-from ..calib.camera.gainselection import GainSelector
+from ..calib.camera.gainselection import GainChannel, GainSelector
 from ..containers import (
     ArrayEventContainer,
     CoordinateFrameType,
@@ -29,6 +29,7 @@
     ObservationBlockContainer,
     ObservationBlockState,
     ObservingMode,
+    PixelStatus,
     PixelStatusContainer,
     PointingContainer,
     PointingMode,
@@ -274,17 +275,20 @@ def apply_simtel_r1_calibration(
         The gain channel selected for each pixel
         Shape: (n_pixels)
     """
-    n_channels, n_pixels, n_samples = r0_waveforms.shape
+    n_channels, n_pixels, _ = r0_waveforms.shape
     ped = pedestal[..., np.newaxis]
     DC_to_PHE = dc_to_pe[..., np.newaxis]
     gain = DC_to_PHE * calib_scale
+
     r1_waveforms = (r0_waveforms - ped) * gain + calib_shift
+
     if n_channels == 1:
         selected_gain_channel = np.zeros(n_pixels, dtype=np.int8)
         r1_waveforms = r1_waveforms[0]
     else:
         selected_gain_channel = gain_selector(r0_waveforms)
         r1_waveforms = r1_waveforms[selected_gain_channel, np.arange(n_pixels)]
+
     return r1_waveforms, selected_gain_channel
 
 
@@ -809,7 +813,7 @@ def _generate_events(self):
                 mon = data.mon.tel[tel_id]
                 mon.calibration.dc_to_pe = dc_to_pe
                 mon.calibration.pedestal_per_sample = pedestal
-                mon.pixel_status = self._get_pixels_status(tel_id)
+                mon.pixel_status = self._fill_mon_pixels_status(tel_id)
 
                 r1_waveform, selected_gain_channel = apply_simtel_r1_calibration(
                     adc_samples,
@@ -819,9 +823,16 @@ def _generate_events(self):
                     self.calib_scale,
                     self.calib_shift,
                 )
+
+                pixel_status = self._get_r1_pixel_status(
+                    tel_id=tel_id,
+                    selected_gain_channel=selected_gain_channel,
+                )
                 data.r1.tel[tel_id] = R1CameraContainer(
+                    event_type=trigger.event_type,
                     waveform=r1_waveform,
                     selected_gain_channel=selected_gain_channel,
+                    pixel_status=pixel_status,
                 )
 
                 # get time_shift from laser calibration
@@ -833,7 +844,26 @@ def _generate_events(self):
 
             yield data
 
-    def _get_pixels_status(self, tel_id):
+    def _get_r1_pixel_status(self, tel_id, selected_gain_channel):
+        tel_desc = self.file_.telescope_descriptions[tel_id]
+        n_pixels = tel_desc["camera_organization"]["n_pixels"]
+        pixel_status = np.zeros(n_pixels, dtype=np.uint8)
+
+        high_gain_stored = selected_gain_channel == GainChannel.HIGH
+        low_gain_stored = selected_gain_channel == GainChannel.LOW
+
+        # set gain bits
+        pixel_status[high_gain_stored] |= PixelStatus.HIGH_GAIN_STORED
+        pixel_status[low_gain_stored] |= PixelStatus.LOW_GAIN_STORED
+
+        # reset gain bits for completely disabled pixels
+        disabled = tel_desc["disabled_pixels"]["HV_disabled"]
+        channel_bits = PixelStatus.HIGH_GAIN_STORED | PixelStatus.LOW_GAIN_STORED
+        pixel_status[disabled] &= ~np.uint8(channel_bits)
+
+        return pixel_status
+
+    def _fill_mon_pixels_status(self, tel_id):
         tel = self.file_.telescope_descriptions[tel_id]
         n_pixels = tel["camera_organization"]["n_pixels"]
         n_gains = tel["camera_organization"]["n_gains"]

ctapipe/tests/test_datamodel.py

@@ -0,0 +1,18 @@
+"""Tests for the data structures defined in ctapipe.containers"""
+import numpy as np
+from numpy.testing import assert_equal
+
+
+def test_pixel_status():
+    """Test methods of the PixelStatus enum on numpy arrays"""
+    from ctapipe.containers import PixelStatus
+
+    pixel_status = np.array([0b1101, 0b1000, 0b1110, 0b1101, 0b0000], dtype=np.uint8)
+
+    assert_equal(
+        PixelStatus.is_invalid(pixel_status), [False, False, False, False, True]
+    )
+    assert_equal(PixelStatus.get_dvr_status(pixel_status), [1, 0, 2, 1, 0])
+    assert_equal(
+        PixelStatus.get_channel_info(pixel_status), [0b11, 0b10, 0b11, 0b11, 0b00]
+    )

docs/changes/2338.feature.rst

@@ -0,0 +1 @@
+Added missing fields defined in the CTAO R1 and DL0 data models to the corresponding containers.