Variance calibration

docs/changes/2636.feature.rst

@@ -0,0 +1,3 @@
+Update ``CameraCalibrator`` in ``ctapipe.calib.camera.calibrator`` allowing it to correctly calibrate variance images generated with the ``VarianceExtractor``.
+  - If the VarianceExtractor is used for the ``CameraCalibrator`` the element-wise square of the relative and absolute gain calibration factors are applied to the image;
+  - For other image extractors the plain factors are still applied.

src/ctapipe/calib/camera/calibrator.py

@@ -283,7 +283,11 @@ def _calibrate_dl1(self, event, tel_id):
             )
 
             # correct non-integer remainder of the shift if given
-            if self.apply_peak_time_shift.tel[tel_id] and remaining_shift is not None:
+            if (
+                self.apply_peak_time_shift.tel[tel_id]
+                and remaining_shift is not None
+                and dl1.peak_time is not None
+            ):
                 dl1.peak_time -= remaining_shift
 
         # Calibrate extracted charge
@@ -292,13 +296,21 @@ def _calibrate_dl1(self, event, tel_id):
             and dl1_calib.absolute_factor is not None
         ):
             if selected_gain_channel is None:
-                dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor
+                if extractor.__class__.__name__ == "VarianceExtractor":
+                    dl1.image *= np.square(
+                        dl1_calib.relative_factor / dl1_calib.absolute_factor
+                    )
+                else:
+                    dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor
             else:
                 corr = (
                     dl1_calib.relative_factor[selected_gain_channel, pixel_index]
                     / dl1_calib.absolute_factor[selected_gain_channel, pixel_index]
                 )
-                dl1.image *= corr
+                if extractor.__class__.__name__ == "VarianceExtractor":
+                    dl1.image *= np.square(corr)
+                else:
+                    dl1.image *= corr
 
         # handle invalid pixels
         if self.invalid_pixel_handler is not None:

src/ctapipe/calib/camera/tests/test_calibrator.py

@@ -16,6 +16,7 @@
     GlobalPeakWindowSum,
     LocalPeakWindowSum,
     NeighborPeakWindowSum,
+    VarianceExtractor,
 )
 from ctapipe.image.reducer import NullDataVolumeReducer, TailCutsDataVolumeReducer
 
@@ -130,6 +131,66 @@ def test_check_dl0_empty(example_event, example_subarray):
     assert (event.dl1.tel[tel_id].image == 2).all()
 
 
+def test_dl1_variance_calib(example_event, example_subarray):
+    # test the calibration of variance images
+    tel_id = list(example_event.r0.tel)[0]
+    calibrator = CameraCalibrator(
+        subarray=example_subarray,
+        image_extractor=VarianceExtractor(subarray=example_subarray),
+        apply_waveform_time_shift=False,
+    )
+    calibrator(example_event)
+    image = example_event.dl1.tel[tel_id].image
+    assert image is not None
+    assert image.shape == (
+        1,
+        1764,
+    )
+
+
+def test_calib_LST_camera(example_subarray):
+    n_channels = 2
+    n_pixels = 1855  # number of pixels in LSTcam
+    n_samples = 100
+
+    random = np.random.default_rng(1)
+
+    tel_id = 1
+    y = random.normal(0, 6, (n_channels, n_pixels, n_samples))
+
+    absolute = random.uniform(100, 1000, (n_channels, n_pixels)).astype("float32")
+    y *= absolute[..., np.newaxis]
+
+    relative = random.normal(1, 0.01, (n_channels, n_pixels))
+    y /= relative[..., np.newaxis]
+
+    pedestal = random.uniform(-4, 4, (n_channels, n_pixels))
+    y += pedestal[..., np.newaxis]
+
+    event = ArrayEventContainer()
+    event.dl0.tel[tel_id].waveform = y
+    event.calibration.tel[tel_id].dl1.pedestal_offset = pedestal
+    event.calibration.tel[tel_id].dl1.absolute_factor = absolute
+    event.calibration.tel[tel_id].dl1.relative_factor = relative
+    event.dl0.tel[tel_id].selected_gain_channel = None
+    event.r1.tel[tel_id].selected_gain_channel = None
+
+    calibrator = CameraCalibrator(
+        subarray=example_subarray,
+        image_extractor=VarianceExtractor(subarray=example_subarray),
+        apply_waveform_time_shift=False,
+    )
+    calibrator(event)
+
+    image = event.dl1.tel[tel_id].image
+
+    assert image is not None
+    assert image.shape == (
+        2,
+        1855,
+    )
+
+
 def test_dl1_charge_calib(example_subarray):
     # copy because we mutate the camera, should not affect other tests
     subarray = deepcopy(example_subarray)