Apply CALIB_SCALE directly from ctapipe.io.SimtelEventSource

protopipe/pipeline/event_preparer.py

@@ -145,13 +145,6 @@ def __init__(
     ):
         """Initiliaze an EventPreparer object."""
 
-        # Calibscale
-        try:
-            self.calibscale = config["Calibration"]["calibscale"]
-        except KeyError:
-            # defaults for no calibscale applied
-            self.calibscale = 1.0
-
         # Cleaning for reconstruction
         self.cleaner_reco = ImageCleaner(  # for reconstruction
             config=config["ImageCleaning"]["biggest"],
@@ -467,7 +460,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
 
                 # use ctapipe's functionality to get the calibrated image
                 # and scale the reconstructed values if required
-                pmt_signal = event.dl1.tel[tel_id].image / self.calibscale
+                pmt_signal = event.dl1.tel[tel_id].image
 
                 # If required...
                 if save_images is True:

protopipe/pipeline/temp.py

@@ -17,7 +17,13 @@
 import astropy.units as u
 from astropy.coordinates import SkyCoord
 
+from ctapipe.io import SimTelEventSource
+from ctapipe.containers import (ArrayEventContainer,
+                                SimulatedCameraContainer,
+                                SimulatedEventContainer,
+                                EventType)
 from ctapipe.core import Container, Field
+from ctapipe.core.traits import Float
 from ctapipe.coordinates import CameraFrame, TelescopeFrame
 from ctapipe.instrument import CameraGeometry
 from ctapipe.reco.hillas_reconstructor import HillasReconstructor
@@ -26,6 +32,171 @@
 logger = logging.getLogger(__name__)
 
 
+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
+        - Pedestal subtraction
+        - Conversion of samples into units proportional to photoelectrons
+          (If the full signal in the waveform was integrated, then the resulting
+          value would be in photoelectrons.)
+          (Also applies flat-fielding)
+    Parameters
+    ----------
+    r0_waveforms : ndarray
+        Raw ADC waveforms from a simtel file. All gain channels available.
+        Shape: (n_channels, n_pixels, n_samples)
+    pedestal : ndarray
+        Pedestal stored in the simtel file for each gain channel
+        Shape: (n_channels, n_pixels)
+    dc_to_pe : ndarray
+        Conversion factor between R0 waveform samples and ~p.e., stored in the
+        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
+    -------
+    r1_waveforms : ndarray
+        Calibrated waveforms
+        Shape: (n_pixels, n_samples)
+    selected_gain_channel : ndarray
+        The gain channel selected for each pixel
+        Shape: (n_pixels)
+    """
+    n_channels, n_pixels, n_samples = r0_waveforms.shape
+    ped = pedestal[..., np.newaxis]
+    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]
+    else:
+        selected_gain_channel = gain_selector(r0_waveforms)
+        r1_waveforms = r1_waveforms[selected_gain_channel, np.arange(n_pixels)]
+    return r1_waveforms, selected_gain_channel
+
+
+class MySimTelEventSource(SimTelEventSource):
+    """ Read events from a SimTelArray data file (in EventIO format)."""
+
+    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 _generate_events(self):
+        data = ArrayEventContainer()
+        data.simulation = SimulatedEventContainer()
+        data.meta["origin"] = "hessio"
+        data.meta["input_url"] = self.input_url
+        data.meta["max_events"] = self.max_events
+
+        self._fill_array_pointing(data)
+
+        for counter, array_event in enumerate(self.file_):
+
+            event_id = array_event.get("event_id", -1)
+            obs_id = self.file_.header["run"]
+            data.count = counter
+            data.index.obs_id = obs_id
+            data.index.event_id = event_id
+
+            self._fill_trigger_info(data, array_event)
+
+            if data.trigger.event_type == EventType.SUBARRAY:
+                self._fill_simulated_event_information(data, array_event)
+
+            # this should be done in a nicer way to not re-allocate the
+            # data each time (right now it's just deleted and garbage
+            # collected)
+            data.r0.tel.clear()
+            data.r1.tel.clear()
+            data.dl0.tel.clear()
+            data.dl1.tel.clear()
+            data.pointing.tel.clear()
+            data.simulation.tel.clear()
+
+            telescope_events = array_event["telescope_events"]
+            tracking_positions = array_event["tracking_positions"]
+
+            for tel_id, telescope_event in telescope_events.items():
+                adc_samples = telescope_event.get("adc_samples")
+                if adc_samples is None:
+                    adc_samples = telescope_event["adc_sums"][:, :, np.newaxis]
+
+                n_gains, n_pixels, n_samples = adc_samples.shape
+                true_image = (
+                    array_event.get("photoelectrons", {})
+                    .get(tel_id - 1, {})
+                    .get("photoelectrons", None)
+                )
+
+                data.simulation.tel[tel_id] = SimulatedCameraContainer(
+                    true_image=true_image
+                )
+
+                data.pointing.tel[tel_id] = self._fill_event_pointing(
+                    tracking_positions[tel_id]
+                )
+
+                r0 = data.r0.tel[tel_id]
+                r1 = data.r1.tel[tel_id]
+                r0.waveform = adc_samples
+
+                cam_mon = array_event["camera_monitorings"][tel_id]
+                pedestal = cam_mon["pedestal"] / cam_mon["n_ped_slices"]
+                dc_to_pe = array_event["laser_calibrations"][tel_id]["calib"]
+
+                # fill dc_to_pe and pedestal_per_sample info into monitoring
+                # container
+                mon = data.mon.tel[tel_id]
+                mon.calibration.dc_to_pe = dc_to_pe
+                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,
+                    self.calib_scale,
+                    self.calib_shift
+                )
+
+                # get time_shift from laser calibration
+                time_calib = array_event["laser_calibrations"][tel_id]["tm_calib"]
+                pix_index = np.arange(n_pixels)
+
+                dl1_calib = data.calibration.tel[tel_id].dl1
+                dl1_calib.time_shift = time_calib[r1.selected_gain_channel, pix_index]
+
+            yield data
+
+
 class HillasParametersTelescopeFrameContainer(Container):
     container_prefix = "hillas"
 

protopipe/scripts/data_training.py

@@ -9,9 +9,10 @@
 import signal
 import tables as tb
 import pandas as pd
+from traitlets.config import Config
 
 from ctapipe.utils import CutFlow
-from ctapipe.io import EventSource
+from protopipe.pipeline.temp import MySimTelEventSource
 
 from protopipe.pipeline import EventPreparer
 from protopipe.pipeline.utils import (
@@ -222,13 +223,27 @@ def main():
     outfile = tb.open_file(args.outfile, mode="w")
     outTable = {}
     outData = {}
+
+    # Configuration options for SimTelEventSource
+    # Readout window integration correction
+    try:
+        calib_scale = cfg["Calibration"]["calib_scale"]
+    except KeyError:
+        # defaults for no calibscale applied
+        calib_scale = 1.0
+
+    cfg_SimTelEventSource = Config()
+    cfg_SimTelEventSource.MySimTelEventSource.calib_scale = calib_scale
 
     for i, filename in enumerate(filenamelist):
 
         print("file: {} filename = {}".format(i, filename))
 
-        source = EventSource(
-            input_url=filename, allowed_tels=allowed_tels, max_events=args.max_events
+        source = MySimTelEventSource(
+            input_url=filename,
+            config=cfg_SimTelEventSource,
+            allowed_tels=allowed_tels,
+            max_events=args.max_events
         )
 
         # loop that cleans and parametrises the images and performs the

protopipe/scripts/write_dl2.py

@@ -8,12 +8,13 @@
 from astropy.coordinates.angle_utilities import angular_separation
 import tables as tb
 import astropy.units as u
+from traitlets.config import Config
 
 # ctapipe
-from ctapipe.io import EventSource
 from ctapipe.utils import CutFlow
 
 # Utilities
+from protopipe.pipeline.temp import MySimTelEventSource
 from protopipe.pipeline import EventPreparer
 from protopipe.pipeline.utils import (
     bcolors,
@@ -272,11 +273,24 @@ class RecoEvent(tb.IsDescription):
         images_outfile = tb.open_file("images.h5", mode="w")
         images_table = {}
         images_phe = {}
+
+    # Configuration options for MySimTelEventSource
+    try:
+        calib_scale = cfg["Calibration"]["calib_scale"]
+    except KeyError:
+        # defaults for no calibscale applied
+        calib_scale = 1.0
+
+    cfg_SimTelEventSource = Config()
+    cfg_SimTelEventSource.MySimTelEventSource.calib_scale = calib_scale
 
     for i, filename in enumerate(filenamelist):
 
-        source = EventSource(
-            input_url=filename, allowed_tels=allowed_tels, max_events=args.max_events
+        source = MySimTelEventSource(
+            input_url=filename,
+            config=cfg_SimTelEventSource,
+            allowed_tels=allowed_tels,
+            max_events=args.max_events
         )
         # loop that cleans and parametrises the images and performs the reconstruction
         for (