Merge pull request #1182 from cta-observatory/drs4_evb6

Change dtype of drs4 baseline values to int16 to make it work with pre-corrected EVBv6 data.

lstchain/scripts/lstchain_convert_drs4_pedestal_to_evb.py

@@ -59,8 +59,8 @@ def main():
 
     pedestal = read_pedestal(args.drs4_baseline_file, args.tel_id)
 
-    if pedestal.dtype != np.uint16:
-        print(f'Input pedestal must have dtype uint16, got {pedestal.dtype}.', file=sys.stderr)
+    if not (pedestal.dtype == np.uint16 or pedestal.dtype == np.int16):
+        print(f'Input pedestal must have dtype (u)int16, got {pedestal.dtype}.', file=sys.stderr)
         print('Make sure the pedestal file was generated by lstchain > 0.8.2', file=sys.stderr)
         sys.exit(1)
 

lstchain/tools/lstchain_create_drs4_pedestal_file.py

@@ -14,7 +14,7 @@
     flag,
 )
 from ctapipe.io.hdf5tableio import HDF5TableWriter
-from ctapipe_io_lst import LSTEventSource
+from ctapipe_io_lst import LSTEventSource, EVBPreprocessing, TriggerBits
 from ctapipe_io_lst.calibration import get_spike_A_positions
 from ctapipe_io_lst.constants import (
     N_CAPACITORS_PIXEL,
@@ -31,13 +31,13 @@ class DRS4CalibrationContainer(Container):
     baseline_mean = Field(
         None,
         "Mean baseline of each capacitor, shape (N_GAINS, N_PIXELS, N_CAPACITORS_PIXEL)",
-        dtype=np.uint16,
+        dtype=np.int16,
         ndim=3,
     )
     baseline_std = Field(
         None,
         "Std Dev. of the baseline calculation, shape (N_GAINS, N_PIXELS, N_CAPACITORS_PIXEL)",
-        dtype=np.uint16,
+        dtype=np.int16,
         ndim=3,
     )
     baseline_counts = Field(
@@ -51,14 +51,16 @@ class DRS4CalibrationContainer(Container):
         None,
         "Mean spike height for each pixel, shape (N_GAINS, N_PIXELS, 3)",
         ndim=3,
-        dtype=np.uint16,
+        dtype=np.int16,
     )
 
 
-def convert_to_uint16(array):
-    '''Convert an array to uint16, rounding and clipping before to avoid under/overflow'''
-    array = np.clip(np.round(array), 0, np.iinfo(np.uint16).max)
-    return array.astype(np.uint16)
+def convert_to_int16(array):
+    '''Convert an array to int16, rounding and clipping before to avoid under/overflow'''
+    dtype = np.int16
+    info = np.iinfo(dtype)
+    array = np.clip(np.round(array), info.min, info.max)
+    return array.astype(dtype)
 
 
 @numba.njit(cache=True, inline='always')
@@ -257,6 +259,65 @@ def mean_spike_height(self):
 
         return spike_heights
 
+    def check_baseline_values(self, baseline_mean):
+        """
+        Check the values of the drs4 baseline for issues.
+
+        In case of no EVBv5 data or no pre-calibration by EVBv6,
+        we expect no negative and no small values.
+
+        In case of EVBv6 applied baseline corrections, values should
+        be close to 0, but negative values are ok.
+        """
+
+        check_large_values = False
+        check_negative = True
+        check_small = True
+
+        if self.source.cta_r1:
+            preprocessing = self.source.evb_preprocessing[TriggerBits.MONO]
+            if EVBPreprocessing.BASELINE_SUBTRACTION in preprocessing:
+                check_large_values = True
+                check_negative = False
+                check_small = False
+
+
+        if check_negative:
+            negative = baseline_mean < 0
+            n_negative = np.count_nonzero(negative)
+            if n_negative > 0:
+                gain, pixel, capacitor = np.nonzero(negative)
+                self.log.critical(f'{n_negative} baseline values are smaller than 0')
+                self.log.info("Gain | Pixel | Capacitor | Baseline ")
+                for g, p, c in zip(gain, pixel, capacitor):
+                    self.log.info(
+                        f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}"
+                    )
+
+        if check_small:
+            small = baseline_mean < 25
+            n_small = np.count_nonzero(small)
+            if n_small > 0:
+                gain, pixel, capacitor = np.nonzero(small)
+                self.log.warning(f'{n_small} baseline values are smaller than 25')
+                self.log.info("Gain | Pixel | Capacitor | Baseline ")
+                for g, p, c in zip(gain, pixel, capacitor):
+                    self.log.info(
+                        f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}"
+                    )
+
+        if check_large_values:
+            large = np.abs(baseline_mean) > 25
+            n_large = np.count_nonzero(large)
+            if n_large > 0:
+                gain, pixel, capacitor = np.nonzero(large)
+                self.log.warning(f'{n_large} baseline values have an abs value >= 25')
+                self.log.info("Gain | Pixel | Capacitor | Baseline ")
+                for g, p, c in zip(gain, pixel, capacitor):
+                    self.log.info(
+                        f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}"
+                    )
+
     def finish(self):
         tel_id = self.source.tel_id
         self.log.info('Writing output to %s', self.output_path)
@@ -267,28 +328,13 @@ def finish(self):
         baseline_std = self.baseline_stats.std.reshape(shape)
         baseline_counts = self.baseline_stats.counts.reshape(shape).astype(np.uint16)
 
-        n_negative = np.count_nonzero(baseline_mean < 0)
-        if n_negative > 0:
-            gain, pixel, capacitor = np.nonzero(baseline_mean < 0)
-            self.log.critical(f'{n_negative} baseline values are smaller than 0')
-            self.log.info("Gain | Pixel | Capacitor | Baseline ")
-            for g, p, c in zip(gain, pixel, capacitor):
-                self.log.info(f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}")
-
-        n_small = np.count_nonzero(baseline_mean < 25)
-        if n_small > 0:
-            gain, pixel, capacitor = np.nonzero(baseline_mean < 25)
-            self.log.warning(f'{n_small} baseline values are smaller than 25')
-            self.log.info("Gain | Pixel | Capacitor | Baseline ")
-            for g, p, c in zip(gain, pixel, capacitor):
-                self.log.info(f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}")
-
+        self.check_baseline_values(baseline_mean)
 
         # Convert baseline mean and spike heights to uint16, handle missing
         # values and values smaller 0, larger maxint
-        baseline_mean = convert_to_uint16(baseline_mean)
-        baseline_std = convert_to_uint16(baseline_std)
-        spike_height = convert_to_uint16(self.mean_spike_height())
+        baseline_mean = convert_to_int16(baseline_mean)
+        baseline_std = convert_to_int16(baseline_std)
+        spike_height = convert_to_int16(self.mean_spike_height())
 
         with HDF5TableWriter(self.output_path) as writer:
             Provenance().add_output_file(str(self.output_path))

lstchain/tools/tests/test_create_drs4_pedestal_file.py

@@ -35,12 +35,12 @@ def test_create_drs4_pedestal_file(tmp_path):
     baseline_mean = drs4_data['baseline_mean']
     baseline_counts = drs4_data['baseline_std']
 
-    assert baseline_mean.dtype == np.uint16
+    assert baseline_mean.dtype == np.int16
     assert baseline_mean.shape == (N_GAINS, N_PIXELS, N_CAPACITORS_PIXEL)
     assert np.isclose(np.average(baseline_mean[baseline_counts > 0], weights=baseline_counts[baseline_counts > 0]), 400, rtol=0.05)
 
     spike_height = drs4_data['spike_height']
-    assert spike_height.dtype == np.uint16
+    assert spike_height.dtype == np.int16
     mean_spike_height = np.nanmean(spike_height, axis=(0, 1))
 
     # these are the expected spike heights, but due to the low statistics,