Add support for calibration scale and shift in R1

ctapipe/io/simteleventsource.py

@@ -22,7 +22,7 @@
     TelescopeTriggerContainer,
 )
 from ..coordinates import CameraFrame
-from ..core.traits import Bool, CaselessStrEnum, create_class_enum_trait
+from ..core.traits import Bool, Float, CaselessStrEnum, create_class_enum_trait
 from ..instrument import (
     CameraDescription,
     CameraGeometry,
@@ -100,7 +100,7 @@ def build_camera(cam_settings, pixel_settings, telescope, frame):
     )
 
 
-def apply_simtel_r1_calibration(r0_waveforms, pedestal, dc_to_pe, gain_selector):
+def apply_simtel_r1_calibration(r0_waveforms, pedestal, dc_to_pe, gain_selector, calib_scale=1.0, calib_shift=0.0):
     """
     Perform the R1 calibration for R0 simtel waveforms. This includes:
         - Gain selection
@@ -123,6 +123,14 @@ def apply_simtel_r1_calibration(r0_waveforms, pedestal, dc_to_pe, gain_selector)
         simtel file for each gain channel
         Shape: (n_channels, n_pixels)
     gain_selector : ctapipe.calib.camera.gainselection.GainSelector
+    calib_scale : float
+        Extra global scale factor for calibration.
+        Conversion factor to transform the integrated charges
+        (in ADC counts) into number of photoelectrons on top of dc_to_pe.
+        Defaults to no scaling.
+    calib_shift: float
+        Shift the resulting R1 output in p.e. for simulating miscalibration.
+        Defaults to no shift.
 
     Returns
     -------
@@ -135,8 +143,9 @@ def apply_simtel_r1_calibration(r0_waveforms, pedestal, dc_to_pe, gain_selector)
     """
     n_channels, n_pixels, n_samples = r0_waveforms.shape
     ped = pedestal[..., np.newaxis]
-    gain = dc_to_pe[..., np.newaxis]
-    r1_waveforms = (r0_waveforms - ped) * gain
+    DC_to_PHE = dc_to_pe[..., np.newaxis]
+    gain = DC_to_PHE * calib_scale
+    r1_waveforms = (r0_waveforms - ped) * gain + calib_shift
     if n_channels == 1:
         selected_gain_channel = np.zeros(n_pixels, dtype=np.int8)
         r1_waveforms = r1_waveforms[0]
@@ -178,6 +187,22 @@ class SimTelEventSource(EventSource):
         base_class=GainSelector, default_value="ThresholdGainSelector"
     ).tag(config=True)
 
+    calib_scale = Float(
+        default_value=1.0,
+        help=(
+            "Factor to transform ADC counts into number of photoelectrons."
+            " Corrects the DC_to_PHE factor."
+        )
+    ).tag(config=True)
+
+    calib_shift = Float(
+        default_value=0.0,
+        help=(
+            "Factor to shift the R1 photoelectron samples. "
+            "Can be used to simulate mis-calibration."
+        )
+    ).tag(config=True)
+
     def __init__(self, input_url=None, config=None, parent=None, **kwargs):
         """
         EventSource for simtelarray files using the pyeventio library.
@@ -422,7 +447,12 @@ def _generate_events(self):
                 mon.calibration.pedestal_per_sample = pedestal
 
                 r1.waveform, r1.selected_gain_channel = apply_simtel_r1_calibration(
-                    adc_samples, pedestal, dc_to_pe, self.gain_selector
+                    adc_samples,
+                    pedestal,
+                    dc_to_pe,
+                    self.gain_selector,
+                    self.calib_scale,
+                    self.calib_shift
                 )
 
                 # get time_shift from laser calibration

ctapipe/io/tests/test_simteleventsource.py

@@ -254,7 +254,10 @@ def test_apply_simtel_r1_calibration_1_channel():
 
     gain_selector = ThresholdGainSelector(threshold=90)
     r1_waveforms, selected_gain_channel = apply_simtel_r1_calibration(
-        r0_waveforms, pedestal, dc_to_pe, gain_selector
+        r0_waveforms,
+        pedestal,
+        dc_to_pe,
+        gain_selector
     )
 
     assert (selected_gain_channel == 0).all()
@@ -285,7 +288,10 @@ def test_apply_simtel_r1_calibration_2_channel():
 
     gain_selector = ThresholdGainSelector(threshold=90)
     r1_waveforms, selected_gain_channel = apply_simtel_r1_calibration(
-        r0_waveforms, pedestal, dc_to_pe, gain_selector
+        r0_waveforms,
+        pedestal,
+        dc_to_pe,
+        gain_selector
     )
 
     assert selected_gain_channel[0] == 1
@@ -331,3 +337,41 @@ def test_only_config():
 
     s = SimTelEventSource(config=config)
     assert s.input_url == Path(gamma_test_large_path).absolute()
+
+
+def test_calibscale_and_calibshift(prod5_gamma_simtel_path):
+
+    telid = 25
+
+    with SimTelEventSource(
+        input_url=prod5_gamma_simtel_path, max_events=1
+    ) as source:
+
+        for event in source:
+            pass
+
+    calib_scale = 2.0
+
+    with SimTelEventSource(
+        input_url=prod5_gamma_simtel_path, max_events=1, calib_scale=calib_scale
+    ) as source:
+
+        for event_scaled in source:
+            pass
+
+    np.testing.assert_allclose(event.r1.tel[telid].waveform[0],
+                               event_scaled.r1.tel[telid].waveform[0] / calib_scale,
+                               rtol=0.1)
+
+    calib_shift = 2.0  # p.e.
+
+    with SimTelEventSource(
+        input_url=prod5_gamma_simtel_path, max_events=1, calib_shift=calib_shift
+    ) as source:
+
+        for event_shifted in source:
+            pass
+
+    np.testing.assert_allclose(event.r1.tel[telid].waveform[0],
+                               event_shifted.r1.tel[telid].waveform[0] - calib_shift,
+                               rtol=0.1)