R1- and DL0-waveform data model

src/ctapipe_io_nectarcam/__init__.py

@@ -45,7 +45,7 @@
 
 from .anyarray_dtypes import CDTS_AFTER_37201_DTYPE, CDTS_BEFORE_37201_DTYPE, TIB_DTYPE
 from .calibration import NectarCAMR0Corrections
-from .constants import N_GAINS, N_PIXELS, N_SAMPLES, nectarcam_location
+from .constants import N_GAINS, N_MODULES, N_PIXELS, N_SAMPLES, nectarcam_location
 from .containers import (
     NectarCAMDataContainer,
     NectarCAMDataStreamContainer,
@@ -653,13 +653,13 @@ def is_simulation(self):
 
     @property
     def datalevels(self):
-        '''
+        """
         TODO: TO BE COMPLETED
         Needs to be updated to take into account different cases
         if select_gain: (R0,R1)
         if input data file is already gain selected: (R1)
         otherwise: (R0)
-        '''
+        """
         if self.r0_r1_calibrator.select_gain:
             return (DataLevel.R0, DataLevel.R1)
         if self.r0_r1_calibrator.calibration_path is not None:
@@ -962,9 +962,21 @@ def fill_nectarcam_event_container_from_zfile(self, array_event, event):
             # Above info not in Event for v6. Put None instead
 
         event_container.event_id = event.event_id
-        event_container.pixel_status = event.pixel_status
-
-        event_container.module_status = nectarcam_data.module_status
+        zfits_pixel_status = event.pixel_status
+        event_container.pixel_status = np.ones(N_PIXELS, dtype=zfits_pixel_status.dtype)
+        # R1 data should have the bit 0 at 1 according to the R1 data
+        event_container.pixel_status[
+            self.nectarcam_service.pixel_ids
+        ] = zfits_pixel_status
+
+        zfits_module_status = nectarcam_data.module_status
+        event_container.module_status = np.zeros(
+            N_MODULES, dtype=zfits_module_status.dtype
+        )
+        # module_status: 1 = present, 0 = absent
+        event_container.module_status[
+            self.nectarcam_service.module_ids
+        ] = zfits_module_status
         event_container.extdevices_presence = nectarcam_data.extdevices_presence
         event_container.swat_data = nectarcam_data.swat_data
         event_container.counters = nectarcam_data.counters
@@ -1017,8 +1029,9 @@ def fill_nectarcam_event_container_from_zfile(self, array_event, event):
             self.unpack_feb_data(event_container, event, nectarcam_data)
 
         # Fill information of the trigger mask in the pixel_status if needed
-        # Note that FEB data must have been loaded
-        if self._missing_trig_pat:
+        # Note that FEB data must have been loaded.
+        # Should it be forced if _missing_trig_pat is at False ?
+        if self._missing_trig_pat and event_container.trigger_pattern is not None:
             tpat = np.any(event_container.trigger_pattern, axis=0)
 
             event_container.pixel_status = (
@@ -1250,6 +1263,9 @@ def fill_r0r1_camera_container(self, zfits_event):
                 (N_PIXELS, N_SAMPLES), fill, dtype=dtype
             )  # VIM : Replace full by empty ?
             reordered_waveform[expected_pixels] = waveform
+            # R1- and DL0-waveform data shape should be (n_gains, n_pixels,
+            # n_samples) from ctapipe v0.21 onwards:
+            reordered_waveform = reordered_waveform[np.newaxis, ...]
 
             reordered_selected_gain = np.full(N_PIXELS, -1, dtype=np.int8)
             reordered_selected_gain[expected_pixels] = selected_gain

src/ctapipe_io_nectarcam/calibration.py

@@ -1,22 +1,22 @@
 import numpy as np
 import tables
 from ctapipe.calib.camera.gainselection import ThresholdGainSelector
-from ctapipe.containers import MonitoringContainer, MonitoringCameraContainer, FlatFieldContainer, WaveformCalibrationContainer
-from ctapipe.core import TelescopeComponent
-from ctapipe.core.traits import (
-    Path, FloatTelescopeParameter, Bool, Float
+from ctapipe.containers import (
+    FlatFieldContainer,
+    MonitoringCameraContainer,
+    MonitoringContainer,
+    WaveformCalibrationContainer,
 )
+from ctapipe.core import TelescopeComponent
+from ctapipe.core.traits import Bool, Float, FloatTelescopeParameter, Path
 from ctapipe.io import HDF5TableReader
 from pkg_resources import resource_filename
 
-from .constants import (
-    N_GAINS, N_PIXELS, HIGH_GAIN, LOW_GAIN,
-    PIXEL_INDEX,
-)
+from .constants import HIGH_GAIN, LOW_GAIN, N_GAINS, N_PIXELS, PIXEL_INDEX
 from .containers import NectarCAMDataContainer
 
 __all__ = [
-    'NectarCAMR0Corrections',
+    "NectarCAMR0Corrections",
 ]
 
 
@@ -30,36 +30,33 @@ class NectarCAMR0Corrections(TelescopeComponent):
 
     calibration_path = Path(
         default_value=resource_filename(
-            'ctapipe_io_nectarcam',
-            'resources/calibrationfile_run3255_pedrun3255_gainrun3155_ctapipe018.hdf5'
+            "ctapipe_io_nectarcam",
+            "resources/calibrationfile_run3255_pedrun3255_gainrun3155_ctapipe018.hdf5",
         ),
-        exists=True, directory_ok=False, allow_none=True,
-        help='Path to calibration file',
+        exists=True,
+        directory_ok=False,
+        allow_none=True,
+        help="Path to calibration file",
     ).tag(config=True)
 
     calib_scale_high_gain = FloatTelescopeParameter(
-        default_value=1.0,
-        help='High gain waveform is multiplied by this number'
+        default_value=1.0, help="High gain waveform is multiplied by this number"
     ).tag(config=True)
 
     calib_scale_low_gain = FloatTelescopeParameter(
-        default_value=1.0,
-        help='Low gain waveform is multiplied by this number'
+        default_value=1.0, help="Low gain waveform is multiplied by this number"
     ).tag(config=True)
 
     select_gain = Bool(
-        default_value=False,
-        help='Set to False to keep both gains.'
+        default_value=False, help="Set to False to keep both gains."
     ).tag(config=True)
 
     gain_selection_threshold = Float(
-        default_value=3500.,
-        help='Threshold for the ThresholdGainSelector.'
+        default_value=3500.0, help="Threshold for the ThresholdGainSelector."
     ).tag(config=True)
 
     apply_flatfield = Bool(
-        default_value=True,
-        help='Apply flatfielding coefficients.'
+        default_value=True, help="Apply flatfielding coefficients."
     ).tag(config=True)
 
     def __init__(self, subarray, config=None, parent=None, **kwargs):
@@ -69,16 +66,13 @@ def __init__(self, subarray, config=None, parent=None, **kwargs):
         Parameters
         ----------
         """
-        super().__init__(
-            subarray=subarray, config=config, parent=parent, **kwargs
-        )
+        super().__init__(subarray=subarray, config=config, parent=parent, **kwargs)
 
         self.mon_data = None
 
         if self.select_gain:
             self.gain_selector = ThresholdGainSelector(
-                threshold=self.gain_selection_threshold,
-                parent=self
+                threshold=self.gain_selection_threshold, parent=self
             )
         else:
             self.gain_selector = None
@@ -92,7 +86,7 @@ def calibrate(self, event: NectarCAMDataContainer):
             # check if waveform is already filled
             if r1.waveform is None:
                 r1.waveform = event.r0.tel[tel_id].waveform.copy()
-                # Force a copy as v6 have waveform as float 
+                # Force a copy as v6 have waveform as float
                 # so in this case r0 == r1 (in memory)
 
             r1.waveform = r1.waveform.astype(np.float32, copy=False)
@@ -111,35 +105,42 @@ def calibrate(self, event: NectarCAMDataContainer):
                 convert_to_pe(
                     waveform=r1.waveform,
                     calibration=calibration,
-                    selected_gain_channel=r1.selected_gain_channel
+                    selected_gain_channel=r1.selected_gain_channel,
                 )
                 # flatfielding
                 if self.apply_flatfield:
                     flatfield = self.mon_data.tel[tel_id].flatfield
                     apply_flatfield(
                         waveform=r1.waveform,
                         flatfield=flatfield,
-                        selected_gain_channel=r1.selected_gain_channel
+                        selected_gain_channel=r1.selected_gain_channel,
                     )
 
             # do not use pixels that are broken in this run or have unusable calibration
             unusable_pixels = np.logical_and(
                 event.mon.tel[tel_id].pixel_status.hardware_failing_pixels,
-                self.mon_data.tel[tel_id].calibration.unusable_pixels
+                self.mon_data.tel[tel_id].calibration.unusable_pixels,
             )
             # set r1 waveforms to 0 for broken pixels
             if r1.selected_gain_channel is None:
                 r1.waveform[unusable_pixels] = 0.0
+                r1.pixel_status = np.uint8(unusable_pixels)
             else:
-                r1.waveform[unusable_pixels[r1.selected_gain_channel, PIXEL_INDEX]] = 0.0
+                r1.waveform[
+                    unusable_pixels[r1.selected_gain_channel, PIXEL_INDEX]
+                ] = 0.0
+                r1.pixel_status = np.uint8(unusable_pixels[r1.selected_gain_channel])
+                r1.waveform = r1.waveform[np.newaxis, ...]
 
-            # needed for charge scaling in ctpaipe dl1 calib
+            # 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]
+                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]
@@ -158,21 +159,27 @@ def _read_calibration_file(path):
 
         with tables.open_file(path) as f:
             tel_ids = [
-                int(key[4:]) for key in f.root._v_children.keys()
-                if key.startswith('tel_')
+                int(key[4:])
+                for key in f.root._v_children.keys()
+                if key.startswith("tel_")
             ]
         for tel_id in tel_ids:
             with HDF5TableReader(path) as h5_table:
                 mon.tel[tel_id] = MonitoringCameraContainer(
-                    calibration=next(h5_table.read(f'/tel_{tel_id}/calibration', WaveformCalibrationContainer)),
-                    flatfield=next(h5_table.read(f'/tel_{tel_id}/flatfield', FlatFieldContainer)),
+                    calibration=next(
+                        h5_table.read(
+                            f"/tel_{tel_id}/calibration", WaveformCalibrationContainer
+                        )
+                    ),
+                    flatfield=next(
+                        h5_table.read(f"/tel_{tel_id}/flatfield", FlatFieldContainer)
+                    ),
                 )
 
         return mon
 
 
 def convert_to_pe(waveform, calibration, selected_gain_channel):
-
     if selected_gain_channel is None:
         waveform -= calibration.pedestal_per_sample
         waveform *= calibration.dc_to_pe[:, :, np.newaxis]
@@ -186,4 +193,5 @@ def apply_flatfield(waveform, flatfield, selected_gain_channel):
         waveform *= flatfield.relative_gain_mean[:, :, np.newaxis]
     else:
         waveform *= flatfield.relative_gain_mean[
-            selected_gain_channel, PIXEL_INDEX, np.newaxis]
+            selected_gain_channel, PIXEL_INDEX, np.newaxis
+        ]

src/ctapipe_io_nectarcam/tests/test_nectarcameventsource.py

@@ -249,7 +249,7 @@ def test_r1_waveforms():
                 input_url=example_file_path, max_events=n_events, config=config
             )
 
-            waveform_shape = (N_PIXELS, N_SAMPLES)
+            waveform_shape = (1, N_PIXELS, N_SAMPLES)
             for event in inputfile_reader:
                 for telid in event.trigger.tels_with_trigger:
                     assert event.r1.tel[telid].waveform.shape == waveform_shape