Add application of DL1 charge calibration to CameraCalibrator

ctapipe/calib/camera/calibrator.py

@@ -211,24 +211,29 @@ def _calibrate_dl1(self, event, telid):
             #   - Read into dl1 container directly?
             #   - Don't do anything if dl1 container already filled
             #   - Update on SST review decision
-            corrected_charge = waveforms[..., 0]
+            charge = waveforms[..., 0]
             pulse_time = np.zeros(n_pixels)
         else:
             # TODO: pass camera to ImageExtractor.__init__
             if self.image_extractor.requires_neighbors():
                 camera = event.inst.subarray.tel[telid].camera
                 self.image_extractor.neighbors = camera.neighbor_matrix_where
+            # TODO: apply timing correction to waveforms before charge extraction
             charge, pulse_time = self.image_extractor(waveforms)
 
             # Apply integration correction
             # TODO: Remove integration correction
             correction = self._get_correction(event, telid)
-            corrected_charge = charge * correction
+            charge = charge * correction
 
-        event.dl1.tel[telid].image = corrected_charge
-        event.dl1.tel[telid].pulse_time = pulse_time
+        # Calibrate extracted charge
+        pedestal = event.calibration.tel[telid].dl1.pedestal_offset
+        absolute = event.calibration.tel[telid].dl1.absolute_factor
+        relative = event.calibration.tel[telid].dl1.relative_factor
+        charge = (charge - pedestal) * relative / absolute
 
-        # TODO: Add charge calibration
+        event.dl1.tel[telid].image = charge
+        event.dl1.tel[telid].pulse_time = pulse_time
 
     def __call__(self, event):
         """

ctapipe/calib/camera/tests/test_calibrator.py

@@ -2,11 +2,13 @@
     CameraCalibrator,
     integration_correction,
 )
-from ctapipe.io.containers import DataContainer
-from ctapipe.image.extractor import LocalPeakWindowSum
+from ctapipe.instrument import CameraGeometry
+from ctapipe.io.containers import DataContainer, EventAndMonDataContainer
+from ctapipe.image.extractor import LocalPeakWindowSum, FullWaveformSum
 from traitlets.config.configurable import Config
 import pytest
 import numpy as np
+from scipy.stats import norm
 
 
 def test_camera_calibrator(example_event):
@@ -114,3 +116,50 @@ def test_check_dl0_empty(example_event):
 
     assert calibrator._check_dl0_empty(None) is True
     assert calibrator._check_dl0_empty(waveform) is False
+
+
+def test_dl1_charge_calib():
+    camera = CameraGeometry.from_name("CHEC")
+    n_pixels = camera.n_pixels
+    n_samples = 96
+    mid = n_samples // 2
+    pulse_sigma = 6
+    random = np.random.RandomState(1)
+    x = np.arange(n_samples)
+
+    # Randomize times and create pulses
+    time_offset = random.uniform(mid - 10, mid + 10, n_pixels)[:, np.newaxis]
+    y = norm.pdf(x, time_offset, pulse_sigma)
+
+    # Define absolute calibration coefficients
+    absolute = random.uniform(100, 1000, n_pixels)
+    y *= absolute[:, np.newaxis]
+
+    # Define relative coefficients
+    relative = random.normal(1, 0.01, n_pixels)
+    y /= relative[:, np.newaxis]
+
+    # Define pedestal
+    pedestal = random.uniform(-4, 4, n_pixels)
+    y += pedestal[:, np.newaxis]
+
+    event = DataContainer()
+    telid = 0
+    event.r1.tel[telid].waveform = y[np.newaxis, ...]
+
+    # Test default
+    calibrator = CameraCalibrator(image_extractor=FullWaveformSum())
+    calibrator(event)
+    np.testing.assert_allclose(event.dl1.tel[telid].image, y.sum(1))
+
+    event.calibration.tel[telid].dl1.time_shift = time_offset
+    event.calibration.tel[telid].dl1.pedestal_offset = pedestal * n_samples
+    event.calibration.tel[telid].dl1.absolute_factor = absolute
+    event.calibration.tel[telid].dl1.relative_factor = relative
+
+    # Test without need for timing corrections
+    calibrator = CameraCalibrator(image_extractor=FullWaveformSum())
+    calibrator(event)
+    np.testing.assert_allclose(event.dl1.tel[telid].image, 1)
+
+    # TODO: Test with timing corrections

ctapipe/io/containers.py

@@ -20,12 +20,14 @@
     "DL0CameraContainer",
     "DL1Container",
     "DL1CameraContainer",
+    "EventCameraCalibrationContainer",
+    "EventCalibrationContainer",
     "SST1MContainer",
     "SST1MCameraContainer",
     "MCEventContainer",
     "MCHeaderContainer",
     "MCCameraEventContainer",
-    "CameraCalibrationContainer",
+    "DL1CameraCalibrationContainer",
     "CentralTriggerContainer",
     "ReconstructedContainer",
     "ReconstructedShowerContainer",
@@ -126,21 +128,40 @@ class DL1CameraContainer(Container):
     )
 
 
-class CameraCalibrationContainer(Container):
-    """
-    Storage of externally calculated calibration parameters (not per-event)
-    """
-
-    dc_to_pe = Field(None, "DC/PE calibration arrays from MC file")
-    pedestal = Field(None, "pedestal calibration arrays from MC file")
-
-
 class DL1Container(Container):
     """ DL1 Calibrated Camera Images and associated data"""
 
     tel = Field(Map(DL1CameraContainer), "map of tel_id to DL1CameraContainer")
 
 
+class DL1CameraCalibrationContainer(Container):
+    """
+    Storage of DL1 calibration parameters for the current event
+    """
+
+    pedestal_offset = Field(
+        0,
+        "Additive coefficients for the pedestal calibration of extracted charge "
+        "for each pixel"
+    )
+    absolute_factor = Field(
+        1,
+        "Multiplicative coefficients for the absolute calibration of extracted charge into "
+        "physical units (e.g. photoelectrons or photons) for each pixel"
+    )
+    relative_factor = Field(
+        1,
+        "Multiplicative Coefficients for the relative correction between pixels to achieve a "
+        "uniform charge response (post absolute calibration) from a "
+        "uniform illumination."
+    )
+    time_shift = Field(
+        0,
+        "Additive coefficients for the timing correction before charge extraction "
+        "for each pixel"
+    )
+
+
 class R0CameraContainer(Container):
     """
     Storage of raw data from a single telescope
@@ -472,6 +493,29 @@ class TelescopePointingContainer(Container):
     altitude = Field(nan * u.rad, "Altitude", unit=u.rad)
 
 
+class EventCameraCalibrationContainer(Container):
+    """
+    Container for the calibration coefficients for the current event and camera
+    """
+    dl1 = Field(
+        DL1CameraCalibrationContainer(), "Container for DL1 calibration coefficients"
+    )
+
+
+class EventCalibrationContainer(Container):
+    """
+    Container for calibration coefficients for the current event
+    """
+
+    tels_with_data = Field([], "list of telescopes with data")
+
+    # create the camera container
+    tel = Field(
+        Map(EventCameraCalibrationContainer),
+        "map of tel_id to EventCameraCalibrationContainer"
+    )
+
+
 class DataContainer(Container):
     """ Top-level container for all event information """
 
@@ -492,6 +536,10 @@ class DataContainer(Container):
         reason="will be separated from event structure in future version",
     )
     pointing = Field(Map(TelescopePointingContainer), "Telescope pointing positions")
+    calibration = Field(
+        EventCalibrationContainer(),
+        "Container for calibration coefficients for the current event"
+    )
 
 
 class SST1MDataContainer(DataContainer):
@@ -818,7 +866,6 @@ class WaveformCalibrationContainer(Container):
     """
     Container for the pixel calibration coefficients
     """
-
     time = Field(0, "Time associated to the calibration event", unit=u.s)
     time_range = Field(
         [],