Merge pull request #2636 from cta-observatory/Variance_calibration

Variance calibration

docs/changes/2636.feature.rst

@@ -0,0 +1,4 @@
+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.
+  - The ``VarianceExtractor`` provides no peak time and the calibrator will skip shifting the peak time for extractors like the ``VarianceExtractor`` that similarly do not provide a peak time

src/ctapipe/calib/camera/calibrator.py

@@ -16,7 +16,7 @@
     ComponentName,
     TelescopeParameter,
 )
-from ctapipe.image.extractor import ImageExtractor
+from ctapipe.image.extractor import ImageExtractor, VarianceExtractor
 from ctapipe.image.invalid_pixels import InvalidPixelHandler
 from ctapipe.image.reducer import DataVolumeReducer
 
@@ -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 (
+                dl1.peak_time is not None
+                and self.apply_peak_time_shift.tel[tel_id]
+                and remaining_shift is not None
+            ):
                 dl1.peak_time -= remaining_shift
 
         # Calibrate extracted charge
@@ -292,13 +296,17 @@ 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
+                calibration = dl1_calib.relative_factor / dl1_calib.absolute_factor
             else:
-                corr = (
+                calibration = (
                     dl1_calib.relative_factor[selected_gain_channel, pixel_index]
                     / dl1_calib.absolute_factor[selected_gain_channel, pixel_index]
                 )
-                dl1.image *= corr
+
+            if isinstance(extractor, VarianceExtractor):
+                calibration = calibration**2
+
+            dl1.image *= calibration
 
         # 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,53 @@ def test_check_dl0_empty(example_event, example_subarray):
     assert (event.dl1.tel[tel_id].image == 2).all()
 
 
+def test_dl1_variance_calib(example_subarray):
+    calibrator = CameraCalibrator(
+        subarray=example_subarray,
+        image_extractor=VarianceExtractor(subarray=example_subarray),
+        apply_waveform_time_shift=False,
+    )
+    n_samples = 100
+
+    event = ArrayEventContainer()
+
+    for tel_type in example_subarray.telescope_types:
+        tel_id = example_subarray.get_tel_ids_for_type(tel_type)[0]
+        n_pixels = example_subarray.tel[tel_id].camera.geometry.n_pixels
+        n_channels = example_subarray.tel[tel_id].camera.readout.n_channels
+
+        random = np.random.default_rng(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.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(event)
+
+    for tel_type in example_subarray.telescope_types:
+        tel_id = example_subarray.get_tel_ids_for_type(tel_type)[0]
+        image = event.dl1.tel[tel_id].image
+        camera = example_subarray.tel[tel_id].camera
+        assert image is not None
+        assert image.shape == (
+            camera.readout.n_channels,
+            example_subarray.tel[tel_id].camera.geometry.n_pixels,
+        )
+
+
 def test_dl1_charge_calib(example_subarray):
     # copy because we mutate the camera, should not affect other tests
     subarray = deepcopy(example_subarray)