Merge pull request #216 from cta-observatory/ctar1_missing_module

Add support for data with missing pixels (offline DVR or missing modules) in EVBv6 / R1v1 format

Author: rlopezcoto

src/ctapipe_io_lst/__init__.py

@@ -479,27 +479,30 @@ def fill_from_cta_r1(self, array_event, zfits_event):
         offset = self.data_stream.waveform_offset
         pixel_id_map = self.camera_config.pixel_id_map
 
-        # FIXME: missing modules / pixels
-        # FIXME: DVR? should not happen in r1 but dl0, but our own converter uses the old R1
-
         # reorder to nominal pixel order
         pixel_status = _reorder_pixel_status(
-            zfits_event.pixel_status, pixel_id_map, set_dvr_bits=True
+            zfits_event.pixel_status, pixel_id_map, set_dvr_bits=not self.dvr_applied,
         )
 
         n_channels = zfits_event.num_channels
-        n_pixels = zfits_event.num_pixels
         n_samples = zfits_event.num_samples
 
+        if self.dvr_applied:
+            stored_pixels = (zfits_event.pixel_status & PixelStatus.DVR_STATUS) > 0
+            n_pixels = np.count_nonzero(stored_pixels)
+        else:
+            stored_pixels = slice(None)  # all pixels stored
+            n_pixels = zfits_event.num_pixels
+
         readout_shape = (n_channels, n_pixels, n_samples)
         raw_waveform = zfits_event.waveform.reshape(readout_shape)
         waveform = raw_waveform.astype(np.float32) / scale - offset
 
         reordered_waveform = np.full((n_channels, N_PIXELS, n_samples), 0.0, dtype=np.float32)
-        reordered_waveform[:, pixel_id_map] = waveform
+        reordered_waveform[:, pixel_id_map[stored_pixels]] = waveform
         waveform = reordered_waveform
 
-        # FIXME, check using evb_preprocessing and make ctapipe support 2 gains
+
         if zfits_event.num_channels == 2:
             selected_gain_channel = None
         else:
@@ -529,22 +532,23 @@ def fill_from_cta_r1(self, array_event, zfits_event):
 
         if DataLevel.R0 in self.datalevels:
             reordered_raw_waveform = np.full((n_channels, N_PIXELS, n_samples), 0, dtype=np.uint16)
-            reordered_raw_waveform[:, pixel_id_map] = raw_waveform
+            reordered_raw_waveform[:, pixel_id_map[stored_pixels]] = raw_waveform
             array_event.r0.tel[self.tel_id] = R0CameraContainer(
                 waveform=reordered_raw_waveform,
             )
 
     def fill_lst_from_ctar1(self, zfits_event):
+        pixel_status = _reorder_pixel_status(
+            zfits_event.pixel_status,
+            self.pixel_id_map,
+            set_dvr_bits=not self.dvr_applied,
+        )
         evt = LSTEventContainer(
-            pixel_status=zfits_event.pixel_status.copy(),
+            pixel_status=pixel_status,
             first_capacitor_id=zfits_event.first_cell_id,
             calibration_monitoring_id=zfits_event.calibration_monitoring_id,
             local_clock_counter=zfits_event.module_hires_local_clock_counter,
         )
-        # set bits for dvr if not already set
-        if not self.dvr_applied:
-            not_broken = get_channel_info(evt.pixel_status) != 0
-            evt.pixel_status[not_broken] |= PixelStatus.DVR_STATUS_0
 
         if zfits_event.debug is not None:
             debug = zfits_event.debug
@@ -971,7 +975,7 @@ def tag_flatfield_events(self, array_event):
         '''
         tel_id = self.tel_id
         waveform = array_event.r1.tel[tel_id].waveform
-        
+
         if waveform.ndim == 3:
             image = waveform[HIGH_GAIN].sum(axis=1)
         else:
@@ -981,7 +985,7 @@ def tag_flatfield_events(self, array_event):
         n_in_range = np.count_nonzero(in_range)
 
         looks_like_ff = n_in_range >= self.min_flatfield_pixel_fraction * image.size
-        
+
         if looks_like_ff:
             # Tag as FF only events with 2-gains waveforms: both gains are needed for calibration
             if waveform.ndim == 3:
@@ -1091,7 +1095,7 @@ def fill_r0r1_camera_container(self, zfits_event):
         if CTAPIPE_0_20:
             # reorder to nominal pixel order
             pixel_status = _reorder_pixel_status(
-                zfits_event.pixel_status, pixel_id_map, set_dvr_bits=True
+                zfits_event.pixel_status, pixel_id_map, set_dvr_bits=not self.dvr_applied
             )
             r1.pixel_status = pixel_status
             r1.event_time = cta_high_res_to_time(

src/ctapipe_io_lst/tests/test_lsteventsource.py

@@ -21,6 +21,7 @@
 test_r0_path_all_streams = test_r0_dir / 'LST-1.1.Run02008.0000_first50.fits.fz'
 
 test_missing_module_path = test_data / 'real/R0/20210215/LST-1.1.Run03669.0000_first50.fits.fz'
+test_missing_module_r1v1_path = test_data / 'real/R0/20240310/LST-1.1.Run17016.0000_first10.fits.fz'
 
 test_drive_report = test_data / 'real/monitoring/DrivePositioning/DrivePosition_log_20200218.txt'
 
@@ -142,6 +143,35 @@ def test_missing_modules():
         assert np.all(event.r0.tel[1].waveform[:, 514] == fill)
 
 
+def test_missing_modules_r1v1():
+    from ctapipe_io_lst import LSTEventSource
+    source = LSTEventSource(
+        test_missing_module_r1v1_path,
+        apply_drs4_corrections=False,
+        pointing_information=False,
+    )
+
+    assert source.lst_service.telescope_id == 1
+    assert source.lst_service.num_modules == 264
+
+    n_events = 0
+    for event in source:
+        n_events += 1
+        # one module missing, so 7 pixels
+        assert np.count_nonzero(event.mon.tel[1].pixel_status.hardware_failing_pixels) == N_PIXELS_MODULE * N_GAINS
+        assert np.count_nonzero(event.r0.tel[1].waveform == 0.0) == N_PIXELS_MODULE * N_SAMPLES * N_GAINS
+
+        missing_gain, missing_pixel = np.nonzero(event.mon.tel[1].pixel_status.hardware_failing_pixels)
+        # 514 is one of the missing pixels
+        for gain, pixel in zip(missing_gain, missing_pixel):
+            np.testing.assert_equal(event.r0.tel[1].waveform[gain, pixel], 0.0)
+
+        if CTAPIPE_0_20:
+            np.testing.assert_equal(event.lst.tel[1].evt.pixel_status, event.r1.tel[1].pixel_status)
+
+    assert n_events == 40
+
+
 def test_gain_selected():
     from ctapipe_io_lst import LSTEventSource