diff --git a/docs/changes/2636.feature.rst b/docs/changes/2636.feature.rst new file mode 100644 index 00000000000..d8bdc488322 --- /dev/null +++ b/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 diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index 853ba3f7da8..3c7683da24d 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/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: diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index b5a3875f62e..b320ae75191 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/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)