diff --git a/ctapipe/calib/camera/calibrator.py b/ctapipe/calib/camera/calibrator.py
index 008ed8efc8c..8f79df960ae 100644
--- a/ctapipe/calib/camera/calibrator.py
+++ b/ctapipe/calib/camera/calibrator.py
@@ -32,8 +32,7 @@ class CameraCalibrator(Component):
         extractor_t0='ChargeExtractorFactory.t0',
         window_width='ChargeExtractorFactory.window_width',
         window_shift='ChargeExtractorFactory.window_shift',
-        sig_amp_cut_HG='ChargeExtractorFactory.sig_amp_cut_HG',
-        sig_amp_cut_LG='ChargeExtractorFactory.sig_amp_cut_LG',
+        peak_detection_threshold='ChargeExtractorFactory.peak_detection_threshold',
         lwt='ChargeExtractorFactory.lwt',
         clip_amplitude='CameraDL1Calibrator.clip_amplitude',
         radius='CameraDL1Calibrator.radius',
diff --git a/ctapipe/calib/camera/dl0.py b/ctapipe/calib/camera/dl0.py
index 91c16311531..7d8fe2fcf61 100644
--- a/ctapipe/calib/camera/dl0.py
+++ b/ctapipe/calib/camera/dl0.py
@@ -81,11 +81,16 @@ def reduce(self, event):
             A `ctapipe` event container
         """
         tels = event.r1.tels_with_data
-        for telid in tels:
-            r1 = event.r1.tel[telid].waveform
+        for telid, r1_tel in tels.items():
             if self.check_r1_exists(event, telid):
+                dl0_tel = event.dl0.tel[telid]
+                r1_waveform = r1_tel.waveform
                 if self._reducer is None:
-                    event.dl0.tel[telid].waveform = r1
+                    dl0_tel.waveform = r1_waveform
                 else:
-                    reduction = self._reducer.reduce_waveforms(r1)
-                    event.dl0.tel[telid].waveform = reduction
+                    reduction = self._reducer.reduce_waveforms(r1_waveform)
+                    dl0_tel.waveform = reduction
+
+                dl0_tel.gain_channel = r1_tel.gain_channel
+                dl0_tel.trigger_time = r1_tel.trigger_time
+                dl0_tel.trigger_type = r1_tel.trigger_type
diff --git a/ctapipe/calib/camera/dl1.py b/ctapipe/calib/camera/dl1.py
index a083b33a648..1b069dc2060 100644
--- a/ctapipe/calib/camera/dl1.py
+++ b/ctapipe/calib/camera/dl1.py
@@ -194,13 +194,14 @@ def calibrate(self, event):
         for telid in event.dl0.tels_with_data:
 
             if self.check_dl0_exists(event, telid):
+                dl0_tel = event.dl0.tel[telid]
                 waveforms = event.dl0.tel[telid].waveform
-                n_samples = waveforms.shape[2]
+                n_pixels, n_samples = waveforms.shape
                 if n_samples == 1:
                     # To handle ASTRI and dst
-                    corrected = waveforms[..., 0]
+                    corrected = waveforms[:, 0]
                     window = np.ones(waveforms.shape)
-                    peakpos = np.zeros(waveforms.shape[0:2])
+                    peakpos = np.zeros(n_pixels)
                     cleaned = waveforms
                 else:
                     # Clean waveforms
@@ -215,7 +216,7 @@ def calibrate(self, event):
                     charge, peakpos, window = extract(cleaned)
 
                     # Apply integration correction
-                    correction = self.get_correction(event, telid)[:, None]
+                    correction = self.get_correction(event, telid)[0]
                     corrected = charge * correction
 
                 # Clip amplitude
@@ -227,4 +228,5 @@ def calibrate(self, event):
                 event.dl1.tel[telid].image = corrected
                 event.dl1.tel[telid].extracted_samples = window
                 event.dl1.tel[telid].peakpos = peakpos
-                event.dl1.tel[telid].cleaned = cleaned
+                event.dl1.tel[telid].waveform = cleaned
+                event.dl1.tel[telid].gain_channel = dl0_tel.gain_channel
diff --git a/ctapipe/calib/camera/gainselection.py b/ctapipe/calib/camera/gainselection.py
index 13a8a211765..dd021923b07 100644
--- a/ctapipe/calib/camera/gainselection.py
+++ b/ctapipe/calib/camera/gainselection.py
@@ -1,7 +1,8 @@
 """
 Algorithms to select correct gain channel
 """
-from abc import ABCMeta, abstractclassmethod
+from abc import ABCMeta, abstractmethod
+from enum import IntEnum
 
 import numpy as np
 
@@ -11,7 +12,13 @@
 __all__ = ['GainSelectorFactory',
            'ThresholdGainSelector',
            'SimpleGainSelector',
-           'pick_gain_channel']
+           'pick_gain_channel',
+           'Channel']
+
+class Channel(IntEnum):
+    """ book-keeping for gain channel index """
+    HIGH = 0
+    LOW = 1
 
 
 def pick_gain_channel(waveforms, threshold, select_by_sample=False):
@@ -37,26 +44,25 @@ def pick_gain_channel(waveforms, threshold, select_by_sample=False):
 
     # if we have 2 channels:
     if waveforms.shape[0] == 2:
-        waveforms = np.squeeze(waveforms)
-        new_waveforms = waveforms[0].copy()
+        new_waveforms = waveforms[Channel.HIGH].copy()
 
         if select_by_sample:
             # replace any samples that are above threshold with low-gain ones:
-            gain_mask = waveforms[0] > threshold
-            new_waveforms[gain_mask] = waveforms[1][gain_mask]
+            gain_mask = waveforms[Channel.HIGH] > threshold
+            new_waveforms[gain_mask] = waveforms[Channel.LOW][gain_mask]
         else:
             # use entire low-gain waveform if any sample of high-gain
             # waveform is above threshold
-            gain_mask = (waveforms[0] > threshold).any(axis=1)
-            new_waveforms[gain_mask] = waveforms[1][gain_mask]
+            gain_mask = (waveforms[Channel.HIGH] > threshold).any(axis=1)
+            new_waveforms[gain_mask] = waveforms[Channel.LOW][gain_mask]
 
     elif waveforms.shape[0] == 1:
-        new_waveforms = np.squeeze(waveforms)
+        new_waveforms = waveforms[Channel.HIGH]
         gain_mask = np.zeros_like(new_waveforms).astype(bool)
 
     else:
-        raise ValueError("input waveforms has shape %s. not sure what to do "
-                         "with that.", waveforms.shape)
+        raise ValueError("input waveforms has shape {}. not sure what to do "
+                         "with that.".format(waveforms.shape))
 
     return new_waveforms, gain_mask
 
@@ -67,7 +73,10 @@ class GainSelector(Component, metaclass=ABCMeta):
     single waveform.
     """
 
-    @abstractclassmethod
+    def __init__(self, parent=None, config=None, **kwargs):
+        super().__init__(parent=parent, config=config, **kwargs)
+
+    @abstractmethod
     def select_gains(self, cam_id, multi_gain_waveform):
         """
         Takes an input waveform and cam_id  and performs gain selection
@@ -96,14 +105,17 @@ class SimpleGainSelector(GainSelector):
     Simply choose a single gain channel always.
     """
 
-    channel = traits.Int(default_value=0, help="which gain channel to "
-                                               "retain").tag(config=True)
+    channel = traits.CaselessStrEnum(
+        [x.name for x in Channel],
+        default_value=Channel.HIGH.name,
+        help="which gain channel to retain"
+    ).tag(config=True)
 
     def select_gains(self, cam_id, multi_gain_waveform):
+        chan = Channel[self.channel].value
         return (
-            multi_gain_waveform[self.channel],
-            (np.ones(multi_gain_waveform.shape[1]) * self.channel).astype(
-                np.bool)
+            multi_gain_waveform[chan],
+            np.full(multi_gain_waveform.shape[1], self.channel, dtype=np.bool)
         )
 
 
@@ -149,7 +161,13 @@ def __init__(self, config=None, parent=None, **kwargs):
             self.threshold_table_name,
             role='dl0.tel.svc.gain_thresholds'
         )
-        self.thresholds = dict(zip(tab['cam_id'], tab['gain_threshold_pe']))
+
+        if 'gain_threshold_dc' in tab:
+            thresh = tab['gain_threshold_dc']
+        else:
+            thresh = tab['gain_threshold_pe']
+
+        self.thresholds = dict(zip(tab['cam_id'], thresh))
         self.log.debug("Loaded threshold table: \n %s", tab)
 
     def __str__(self):
diff --git a/ctapipe/calib/camera/r1.py b/ctapipe/calib/camera/r1.py
index 6603ea0bffb..4e08afe3698 100644
--- a/ctapipe/calib/camera/r1.py
+++ b/ctapipe/calib/camera/r1.py
@@ -14,9 +14,10 @@
 of the data.
 """
 from abc import abstractmethod
+from .gainselection import ThresholdGainSelector, SimpleGainSelector
 import numpy as np
 from ...core import Component, Factory
-from ...core.traits import Unicode
+from ...core.traits import Unicode, Float
 from ...io import EventSource
 
 __all__ = [
@@ -48,7 +49,7 @@ class CameraR1Calibrator(Component):
     kwargs
     """
 
-    def __init__(self, config=None, tool=None, **kwargs):
+    def __init__(self, config=None, tool=None, gain_selector=None, **kwargs):
         """
         Parent class for the r1 calibrators. Fills the r1 container.
 
@@ -67,6 +68,8 @@ def __init__(self, config=None, tool=None, **kwargs):
         super().__init__(config=config, parent=tool, **kwargs)
         self._r0_empty_warn = False
 
+        self.gain_selector = gain_selector or SimpleGainSelector(tool, config)
+
     @abstractmethod
     def calibrate(self, event):
         """
@@ -114,7 +117,8 @@ def check_r0_exists(self, event, telid):
 class NullR1Calibrator(CameraR1Calibrator):
     """
     A dummy R1 calibrator that simply fills the r1 container with the samples
-    from the r0 container.
+    from the r0 container.  To do the R1 gain selection, it uses a
+    SimpleGainSelector, which selects by default the high-gain channel.
 
     Parameters
     ----------
@@ -129,16 +133,21 @@ class NullR1Calibrator(CameraR1Calibrator):
     kwargs
     """
 
-    def __init__(self, config=None, tool=None, **kwargs):
-        super().__init__(config, tool, **kwargs)
+    def __init__(self, config=None, tool=None, gain_selector=None, **kwargs):
+        super().__init__(config, tool, gain_selector, **kwargs)
         self.log.info("Using NullR1Calibrator, if event source is at "
                       "the R0 level, then r1 samples will equal r0 samples")
-
+        self.gain_selector = gain_selector or SimpleGainSelector(parent=tool,
+                                                                 config=config)
+        
     def calibrate(self, event):
         for telid in event.r0.tels_with_data:
             if self.check_r0_exists(event, telid):
                 samples = event.r0.tel[telid].waveform
-                event.r1.tel[telid].waveform = samples.astype('float32')
+                cam_id = event.inst.subarray.tel[telid].camera.cam_id
+                wf, mask = self.gain_selector.select_gains(cam_id, samples)
+                event.r1.tel[telid].waveform = wf.astype('float32')
+                event.r1.tel[telid].gain_channel = mask
 
 
 class HESSIOR1Calibrator(CameraR1Calibrator):
@@ -158,23 +167,29 @@ class HESSIOR1Calibrator(CameraR1Calibrator):
         Tool executable that is calling this component.
         Passes the correct logger to the component.
         Set to None if no Tool to pass.
-    kwargs
-    """
 
-    calib_scale = 1.05
+    Attributes
+    ----------
+    calib_scale: float
+        the factor needed to transform from mean p.e. units to
+        units of the single-p.e. peak: Depends on the collection efficiency,
+        the asymmetry of the single p.e. amplitude  distribution and the
+        electronic noise added to the signals. Default value of 1.05 is for
+        GCT, but should be 0.92 for e.g. prod3 SimTelarray simulations.
     """
-    CALIB_SCALE is only relevant for MC calibration.
 
-    CALIB_SCALE is the factor needed to transform from mean p.e. units to 
-    units of the single-p.e. peak: Depends on the collection efficiency, 
-    the asymmetry of the single p.e. amplitude  distribution and the 
-    electronic noise added to the signals. Default value is for GCT.
+    calib_scale = Float(
+        1.05,
+        help="factor to transform from mean PE units to units single PE "
+             "peak. (0.92 for standard MCs)",
+    ).tag(config=True)
 
-    To correctly calibrate to number of photoelectron, a fresh SPE calibration
-    should be applied using a SPE sim_telarray run with an 
-    artificial light source.
-    """
-    # TODO: Handle calib_scale differently per simlated telescope
+    def __init__(self, config=None, tool=None, gain_selector=None, **kwargs):
+        if gain_selector is None:
+            gain_selector = ThresholdGainSelector(config, tool)
+
+        super().__init__(config=config, tool=tool,
+                         gain_selector=gain_selector, **kwargs)
 
     def calibrate(self, event):
         if event.meta['origin'] != 'hessio':
@@ -183,16 +198,22 @@ def calibrate(self, event):
 
         for telid in event.r0.tels_with_data:
             if self.check_r0_exists(event, telid):
-                samples = event.r0.tel[telid].waveform
-                n_samples = samples.shape[2]
+                # todo: KPK  do gain selection first, then need to multiply
+                # correct pedestals with a mask
+                waveform = event.r0.tel[telid].waveform
+                n_samples = waveform.shape[2]
                 ped = event.mc.tel[telid].pedestal / n_samples
                 gain = event.mc.tel[telid].dc_to_pe * self.calib_scale
-                calibrated = (samples - ped[..., None]) * gain[..., None]
-                event.r1.tel[telid].waveform = calibrated
+                calibrated = (waveform - ped[..., None]) * gain[..., None]
 
+                cam_id = event.inst.subarray.tel[telid].camera.cam_id
+                wf, mask = self.gain_selector.select_gains(cam_id, calibrated)
+                event.r1.tel[telid].waveform_full = calibrated
+                event.r1.tel[telid].waveform = wf
+                event.r1.tel[telid].gain_channel = mask
 
-class TargetIOR1Calibrator(CameraR1Calibrator):
 
+class TargetIOR1Calibrator(CameraR1Calibrator):
     pedestal_path = Unicode(
         '',
         allow_none=True,
@@ -288,7 +309,13 @@ def fake_calibrate(self, event):
 
         if self.check_r0_exists(event, self.telid):
             samples = event.r0.tel[self.telid].waveform
-            event.r1.tel[self.telid].waveform = samples.astype('float32')
+
+            cam_id = event.inst.subarray.tel[self.telid].camera.cam_id
+            waveform, mask = self.gain_selector.select_gains(cam_id,
+                                                             samples)
+            event.r1.tel[self.telid].waveform_full = samples.astype('float32')
+            event.r1.tel[self.telid].waveform = waveform.astype('float32')
+            event.r1.tel[self.telid].gain_channel = mask
 
     def real_calibrate(self, event):
         """
@@ -305,12 +332,23 @@ def real_calibrate(self, event):
 
         if self.check_r0_exists(event, self.telid):
             samples = event.r0.tel[self.telid].waveform
-            if self._r1_wf is None:
-                self._r1_wf = np.zeros(samples.shape, dtype=np.float32)
+            cam_id = event.inst.subarray.tel[self.telid].camera.cam_id
+            waveform, mask = self.gain_selector.select_gains(cam_id,
+                                                             samples)
+            event.r1.tel[self.telid].waveform_full = samples.astype('float32')
+            event.r1.tel[self.telid].gain_channel = mask
+
             fci = event.targetio.tel[self.telid].first_cell_ids
-            r1 = event.r1.tel[self.telid].waveform[0]
-            self.calibrator.ApplyEvent(samples[0], fci, self._r1_wf[0])
-            event.r1.tel[self.telid].waveform = self._r1_wf
+            r1 = event.r1.tel[self.telid].waveform
+            self.calibrator.ApplyEvent(waveform, fci, r1)
+
+            # KK: replace prev code with this
+            # if self._r1_wf is None:
+            #     self._r1_wf = np.zeros(samples.shape, dtype=np.float32)
+            # fci = event.targetio.tel[self.telid].first_cell_ids
+            # r1 = event.r1.tel[self.telid].waveform[0]
+            # self.calibrator.ApplyEvent(samples[0], fci, self._r1_wf[0])
+            # event.r1.tel[self.telid].waveform = self._r1_wf
 
 
 class CameraR1CalibratorFactory(Factory):
diff --git a/ctapipe/calib/camera/tests/test_calibrator.py b/ctapipe/calib/camera/tests/test_calibrator.py
index 28dca237c5d..e7c3cf69f8e 100644
--- a/ctapipe/calib/camera/tests/test_calibrator.py
+++ b/ctapipe/calib/camera/tests/test_calibrator.py
@@ -1,5 +1,7 @@
 from copy import deepcopy
+
 from numpy.testing import assert_allclose
+
 from ctapipe.calib.camera import (
     CameraCalibrator,
     HESSIOR1Calibrator,
@@ -11,14 +13,14 @@
 
 
 def test_camera_calibrator(test_event):
-    event = deepcopy(test_event) # so we don't modify the test event
+    event = deepcopy(test_event)  # so we don't modify the test event
     telid = 11
 
     calibrator = CameraCalibrator(r1_product="HESSIOR1Calibrator")
 
     calibrator.calibrate(event)
     image = event.dl1.tel[telid].image
-    assert_allclose(image[0, 0], -2.216, 1e-3)
+    assert_allclose(image[0], -2.216, 1e-3)
 
 
 def test_manual_r1():
diff --git a/ctapipe/calib/camera/tests/test_dl0.py b/ctapipe/calib/camera/tests/test_dl0.py
index 40138fad9d3..49c05a9fdf5 100644
--- a/ctapipe/calib/camera/tests/test_dl0.py
+++ b/ctapipe/calib/camera/tests/test_dl0.py
@@ -19,7 +19,7 @@ def test_camera_dl0_reducer(test_event):
     reducer = CameraDL0Reducer()
     reducer.reduce(event)
     waveforms = event.dl0.tel[telid].waveform
-    assert_almost_equal(waveforms[0, 0, 0], -0.091, 3)
+    assert_almost_equal(waveforms[0, 0], -0.091, 3)
 
 
 def test_check_r1_exists(test_event):
diff --git a/ctapipe/calib/camera/tests/test_dl1.py b/ctapipe/calib/camera/tests/test_dl1.py
index 99c1f8d4231..758b1ac4bc5 100644
--- a/ctapipe/calib/camera/tests/test_dl1.py
+++ b/ctapipe/calib/camera/tests/test_dl1.py
@@ -51,7 +51,13 @@ def test_camera_dl1_calibrator(test_event):
 
     calibrator.calibrate(event)
     image = event.dl1.tel[telid].image
-    assert_allclose(image[0, 0], -2.216, 1e-3)
+    waveform = event.dl1.tel[telid].waveform
+    assert_allclose(image[0], -2.216, 1e-3)
+    assert image.ndim == 1
+    assert len(image) > 1  # make sure image has more than one pixel
+    assert waveform.ndim == 2
+    assert event.dl1.tel[telid].gain_channel.shape == waveform.shape
+    assert event.dl1.tel[telid]
 
 
 def test_check_dl0_exists(test_event):
@@ -59,6 +65,6 @@ def test_check_dl0_exists(test_event):
     telid = 11
     previous_calibration(event)
     calibrator = CameraDL1Calibrator()
-    assert(calibrator.check_dl0_exists(event, telid) is True)
+    assert calibrator.check_dl0_exists(event, telid) is True
     event.dl0.tel[telid].waveform = None
-    assert(calibrator.check_dl0_exists(event, telid) is False)
+    assert calibrator.check_dl0_exists(event, telid) is False
diff --git a/ctapipe/calib/camera/tests/test_gainselection.py b/ctapipe/calib/camera/tests/test_gainselection.py
index 8c90d344d6e..af88da24e1d 100644
--- a/ctapipe/calib/camera/tests/test_gainselection.py
+++ b/ctapipe/calib/camera/tests/test_gainselection.py
@@ -1,8 +1,8 @@
 import numpy as np
 import pytest
 
+from ctapipe.calib.camera.gainselection import SimpleGainSelector, Channel
 from ctapipe.calib.camera.gainselection import ThresholdGainSelector
-from ctapipe.calib.camera.gainselection import SimpleGainSelector
 from ctapipe.calib.camera.gainselection import pick_gain_channel
 
 
@@ -97,22 +97,34 @@ def test_threshold_gain_selector():
     with pytest.raises(ValueError):
         selector.select_gains("NectarCam", np.ones((3, 1000, 30)))
 
-    # 1-gain channel input:
-    wf0 = np.ones((1, 1000, 1))
+    # 1-gain channel input with many samples:
+    wf0 = np.ones((1, 1000, 30))
     wf1, gm = selector.select_gains("ASTRICam", wf0)
-    assert wf1.shape == (1000,)
-    assert gm.shape == (1000,)
+    assert wf1.shape == (1000, 30)
+    assert gm.shape == (1000, 30)
+
+    # 2-gain channel input with no samples:
+    wf0 = np.random.uniform(10, size=(2, 2368, 1))
+    wf1, gm = selector.select_gains("ASTRICam", wf0)
+    assert wf1.shape == (2368, 1)
+    assert gm.shape == (2368, 1)
+
+    # 1-gain channel input with no samples:
+    wf0 = np.random.uniform(10, size=(1, 2368, 1))
+    wf1, gm = selector.select_gains("ASTRICam", wf0)
+    assert wf1.shape == (2368, 1)
+    assert gm.shape == (2368, 1)
 
 
 def test_simple_gain_selector():
     gs = SimpleGainSelector()
 
-    for chan in [0, 1]:
-        gs.channel = chan
+    for chan in [Channel.HIGH, Channel.LOW]:
+        gs.channel = chan.name
 
         waveforms_2g = np.random.normal(size=(2, 1000, 30))
         waveforms_1g, mask = gs.select_gains("NectarCam", waveforms_2g)
 
         assert waveforms_1g.shape == (1000, 30)
-        assert (waveforms_1g == waveforms_2g[chan]).all()
+        assert (waveforms_1g == waveforms_2g[chan.value]).all()
         assert mask.shape == (1000,)
diff --git a/ctapipe/calib/camera/tests/test_r1.py b/ctapipe/calib/camera/tests/test_r1.py
index 8dabfd92d8a..d9eff018a94 100644
--- a/ctapipe/calib/camera/tests/test_r1.py
+++ b/ctapipe/calib/camera/tests/test_r1.py
@@ -20,7 +20,7 @@ def test_hessio_r1_calibrator(test_event):
     calibrator = HESSIOR1Calibrator()
     calibrator.calibrate(event)
     r1 = event.r1.tel[telid].waveform
-    assert_almost_equal(r1[0, 0, 0], -0.091, 3)
+    assert_almost_equal(r1[0, 0], -0.091, 3)
 
 
 def test_null_r1_calibrator(test_event):
@@ -28,7 +28,7 @@ def test_null_r1_calibrator(test_event):
     event = deepcopy(test_event)
     calibrator = NullR1Calibrator()
     calibrator.calibrate(event)
-    r0 = event.r0.tel[telid].waveform
+    r0 = event.r0.tel[telid].waveform[0]
     r1 = event.r1.tel[telid].waveform
     assert_array_equal(r0, r1)
 
@@ -48,7 +48,7 @@ def test_targetio_calibrator():
     event_r1 = source_r1._get_event_by_index(0)
 
     r1c.calibrate(event_r0)
-    assert_array_equal(event_r0.r0.tel[0].waveform,
+    assert_array_equal(event_r0.r0.tel[0].waveform[0],
                        event_r0.r1.tel[0].waveform)
 
     r1c = CameraR1CalibratorFactory.produce(
diff --git a/ctapipe/image/charge_extractors.py b/ctapipe/image/charge_extractors.py
index d23d6221a52..ed7f6910616 100644
--- a/ctapipe/image/charge_extractors.py
+++ b/ctapipe/image/charge_extractors.py
@@ -6,14 +6,26 @@
            'GlobalPeakIntegrator', 'LocalPeakIntegrator',
            'NeighbourPeakIntegrator', 'AverageWfPeakIntegrator']
 
-
 from abc import abstractmethod
+
 import numpy as np
-from traitlets import Int, CaselessStrEnum, Float
+from traitlets import Int, Float
+
 from ctapipe.core import Component, Factory
 from ctapipe.utils.neighbour_sum import get_sum_array
 
 
+class FullWaveAxes:
+    channel = 0
+    pixel = 1
+    sample = 2
+
+
+class WaveAxes:
+    pixel = 0
+    sample = 1
+
+
 class ChargeExtractor(Component):
 
     def __init__(self, config=None, tool=None, **kwargs):
@@ -230,7 +242,7 @@ def get_start_and_width(self, waveforms, peakpos):
         """
         w_start = self._get_window_start(waveforms, peakpos)
         w_width = self._get_window_width(waveforms)
-        n_samples = waveforms.shape[2]
+        n_samples = waveforms.shape[WaveAxes.sample]
         self.check_window_width_and_start(n_samples, w_start, w_width)
         return w_start, w_width
 
@@ -262,7 +274,7 @@ def get_window(waveforms, start, width):
         end = start + width
 
         # Obtain integration window using the indices of the waveforms array
-        ind = np.indices(waveforms.shape)[2]
+        ind = np.indices(waveforms.shape)[WaveAxes.sample]
         integration_window = (ind >= start[..., None]) & (ind < end[..., None])
         return integration_window
 
@@ -288,7 +300,7 @@ def extract_from_window(waveforms, window):
             Extracted charge stored in a numpy array of shape (n_chan, n_pix).
         """
         windowed = np.ma.array(waveforms, mask=~window)
-        charge = windowed.sum(2).data
+        charge = windowed.sum(axis=WaveAxes.sample).data
         return charge
 
     def get_window_from_waveforms(self, waveforms):
@@ -350,16 +362,16 @@ def __init__(self, config=None, tool=None, **kwargs):
         super().__init__(config=config, tool=tool, **kwargs)
 
     def _get_window_start(self, waveforms, peakpos):
-        nchan, npix, nsamples = waveforms.shape
-        return np.zeros((nchan, npix), dtype=np.intp)
+        npix, nsamples = waveforms.shape
+        return np.zeros(npix, dtype=np.intp)
 
     def _get_window_width(self, waveforms):
-        nchan, npix, nsamples = waveforms.shape
-        return np.full((nchan, npix), nsamples, dtype=np.intp)
+        npix, nsamples = waveforms.shape
+        return np.full(npix, nsamples, dtype=np.intp)
 
     def get_peakpos(self, waveforms):
-        nchan, npix, nsamples = waveforms.shape
-        return np.zeros((nchan, npix), dtype=np.intp)
+        npix, nsamples = waveforms.shape
+        return np.zeros(npix, dtype=np.intp)
 
 
 class WindowIntegrator(Integrator):
@@ -400,8 +412,8 @@ def _get_window_start(self, waveforms, peakpos):
         return peakpos - self.window_shift
 
     def _get_window_width(self, waveforms):
-        nchan, npix, nsamples = waveforms.shape
-        return np.full((nchan, npix), self.window_width, dtype=np.intp)
+        npix, nsamples = waveforms.shape
+        return np.full(npix, self.window_width, dtype=np.intp)
 
     @abstractmethod
     def _obtain_peak_position(self, waveforms):
@@ -451,19 +463,16 @@ def __init__(self, config=None, tool=None, **kwargs):
         super().__init__(config=config, tool=tool, **kwargs)
 
     def _obtain_peak_position(self, waveforms):
-        nchan, npix, nsamples = waveforms.shape
-        return np.full((nchan, npix), self.t0, dtype=np.intp)
+        npix, nsamples = waveforms.shape
+        return np.full(npix, self.t0, dtype=np.intp)
 
 
 class PeakFindingIntegrator(WindowIntegrator):
-    sig_amp_cut_HG = Float(None, allow_none=True,
-                           help='Define the cut above which a sample is '
-                                'considered as significant for PeakFinding '
-                                'in the HG channel').tag(config=True)
-    sig_amp_cut_LG = Float(None, allow_none=True,
-                           help='Define the cut above which a sample is '
-                                'considered as significant for PeakFinding '
-                                'in the LG channel').tag(config=True)
+    peak_detection_threshold = Float(
+        None, allow_none=True,
+        help='Define the cut above which a sample is considered as '
+             'significant for PeakFinding in the HG channel'
+    ).tag(config=True)
 
     def __init__(self, config=None, tool=None, **kwargs):
         """
@@ -488,7 +497,6 @@ def __init__(self, config=None, tool=None, **kwargs):
         kwargs
         """
         super().__init__(config=config, tool=tool, **kwargs)
-        self._sig_channel = None
         self._sig_pixels = None
 
     # Extract significant entries
@@ -509,21 +517,15 @@ def _extract_significant_entries(self, waveforms):
             masked.
 
         """
-        nchan, npix, nsamples = waveforms.shape
-        if self.sig_amp_cut_LG or self.sig_amp_cut_HG:
+        npix, nsamples = waveforms.shape
+        if self.peak_detection_threshold:
             sig_entries = np.ones(waveforms.shape, dtype=bool)
-            if self.sig_amp_cut_HG:
-                sig_entries[0] = waveforms[0] > self.sig_amp_cut_HG
-            if nchan > 1 and self.sig_amp_cut_LG:
-                sig_entries[1] = waveforms[1] > self.sig_amp_cut_LG
-            self._sig_pixels = np.any(sig_entries, axis=2)
-            self._sig_channel = np.any(self._sig_pixels, axis=1)
-            if not self._sig_channel[0]:
-                self.log.error("sigamp excludes all values in HG channel")
+            if self.peak_detection_threshold:
+                sig_entries = waveforms > self.peak_detection_threshold
+            self._sig_pixels = np.any(sig_entries, axis=WaveAxes.sample)
             return np.ma.array(waveforms, mask=~sig_entries)
         else:
-            self._sig_channel = np.ones(nchan, dtype=bool)
-            self._sig_pixels = np.ones((nchan, npix), dtype=bool)
+            self._sig_pixels = np.ones(npix, dtype=bool)
             return waveforms
 
     @abstractmethod
@@ -555,25 +557,17 @@ class GlobalPeakIntegrator(PeakFindingIntegrator):
     """
 
     def __init__(self, config=None, tool=None, **kwargs):
-
         super().__init__(config=config, tool=tool, **kwargs)
 
     def _obtain_peak_position(self, waveforms):
-        nchan, npix, nsamples = waveforms.shape
+        npix, nsamples = waveforms.shape
         significant_samples = self._extract_significant_entries(waveforms)
-        max_t = significant_samples.argmax(2)
-        max_s = significant_samples.max(2)
-
-        peakpos = np.zeros((nchan, npix), dtype=np.int)
-        peakpos[0, :] = np.round(np.average(max_t[0], weights=max_s[0]))
-        if nchan > 1:
-            if self._sig_channel[1]:
-                peakpos[1, :] = np.round(
-                    np.average(max_t[1], weights=max_s[1]))
-            else:
-                self.log.info("LG not significant, using HG for peak finding "
-                              "instead")
-                peakpos[1, :] = peakpos[0]
+        max_t = significant_samples.argmax(axis=WaveAxes.sample)
+        max_s = significant_samples.max(axis=WaveAxes.sample)
+
+        peakpos = np.zeros(npix, dtype=np.int)
+        peakpos[:] = np.round(np.average(max_t, weights=max_s))
+
         return peakpos
 
 
@@ -604,14 +598,13 @@ def __init__(self, config=None, tool=None, **kwargs):
         super().__init__(config=config, tool=tool, **kwargs)
 
     def _obtain_peak_position(self, waveforms):
-        nchan, npix, nsamples = waveforms.shape
+        npix, nsamples = waveforms.shape
         significant_samples = self._extract_significant_entries(waveforms)
-        peakpos = np.full((nchan, npix), significant_samples.argmax(2),
-                          dtype=np.int)
-        sig_pix = self._sig_pixels
-        if nchan > 1:  # If the LG is not significant, use the HG peakpos
-            peakpos[1] = np.where(sig_pix[1] < sig_pix[0],
-                                  peakpos[0], peakpos[1])
+        peakpos = np.full(
+            npix,
+            significant_samples.argmax(WaveAxes.sample),
+            dtype=np.int
+        )
         return peakpos
 
 
@@ -649,13 +642,16 @@ def requires_neighbours():
         return True
 
     def _obtain_peak_position(self, waveforms):
-        shape = waveforms.shape
+        npix, nsamp = waveforms.shape
         significant_samples = self._extract_significant_entries(waveforms)
         sig_sam = significant_samples.astype(np.float32)
         sum_data = np.zeros_like(sig_sam)
-        n = self.neighbours.astype(np.uint16)
-        get_sum_array(sig_sam, sum_data, *shape, n, n.shape[0], self.lwt)
-        return sum_data.argmax(2).astype(np.int)
+        neighbors = self.neighbours.astype(np.uint16)
+        neighbors_length = neighbors.shape[0]
+        get_sum_array(sig_sam, sum_data,
+                      npix, nsamp,
+                      neighbors, neighbors_length, self.lwt)
+        return sum_data.argmax(WaveAxes.sample).astype(np.int)
 
 
 class AverageWfPeakIntegrator(PeakFindingIntegrator):
@@ -685,11 +681,11 @@ def __init__(self, config=None, tool=None, **kwargs):
         super().__init__(config=config, tool=tool, **kwargs)
 
     def _obtain_peak_position(self, waveforms):
-        nchan, npix, nsamples = waveforms.shape
+        npix, nsamples = waveforms.shape
         significant_samples = self._extract_significant_entries(waveforms)
-        peakpos = np.zeros((nchan, npix), dtype=np.int)
-        avg_wf = np.mean(significant_samples, axis=1)
-        peakpos += np.argmax(avg_wf, axis=1)[:, None]
+        peakpos = np.zeros(npix, dtype=np.int)
+        avg_wf = np.mean(significant_samples, axis=WaveAxes.pixel)
+        peakpos += np.argmax(avg_wf, axis=WaveAxes.pixel)
         return peakpos
 
 
diff --git a/ctapipe/image/muon/muon_reco_functions.py b/ctapipe/image/muon/muon_reco_functions.py
index 6eb0df48ed9..e5429664334 100644
--- a/ctapipe/image/muon/muon_reco_functions.py
+++ b/ctapipe/image/muon/muon_reco_functions.py
@@ -75,7 +75,7 @@ def analyze_muon_event(event):
     for telid in event.dl0.tels_with_data:
 
         logger.debug("Analysing muon event for tel %d", telid)
-        image = event.dl1.tel[telid].image[0]
+        image = event.dl1.tel[telid].image
 
         # Get geometry
         teldes = event.inst.subarray.tel[telid]
diff --git a/ctapipe/image/tests/test_charge_extraction.py b/ctapipe/image/tests/test_charge_extraction.py
index 69faa084b9e..245eec8ee89 100644
--- a/ctapipe/image/tests/test_charge_extraction.py
+++ b/ctapipe/image/tests/test_charge_extraction.py
@@ -12,115 +12,78 @@
                                              AverageWfPeakIntegrator)
 
 
-def test_full_integration(test_event):
-    telid = 11
-    event = deepcopy(test_event)
+def dummy_calib(event, telid=11):
     data = event.r0.tel[telid].waveform
     nsamples = data.shape[2]
     ped = event.mc.tel[telid].pedestal
     data_ped = data - np.atleast_3d(ped / nsamples)
     data_ped = np.array([data_ped[0], data_ped[0]])  # Test LG functionality
+    return data_ped[0]
+
+
+def test_full_integration(test_event):
+    data_ped = dummy_calib(deepcopy(test_event))
 
     integrator = FullIntegrator()
     integration, peakpos, window = integrator.extract_charge(data_ped)
 
-    assert_almost_equal(integration[0][0], 149, 0)
-    assert_almost_equal(integration[1][0], 149, 0)
-    assert peakpos[0][0] == 0
-    assert peakpos[1][0] == 0
+    assert_almost_equal(integration[0], 149, 0)
+
+    assert peakpos[0] == 0
 
 
 def test_simple_integration(test_event):
-    telid = 11
-    event = deepcopy(test_event)
-    data = event.r0.tel[telid].waveform
-    nsamples = data.shape[2]
-    ped = event.mc.tel[telid].pedestal
-    data_ped = data - np.atleast_3d(ped / nsamples)
-    data_ped = np.array([data_ped[0], data_ped[0]])  # Test LG functionality
+    data_ped = dummy_calib(deepcopy(test_event), telid=11)
 
     integrator = SimpleIntegrator()
     integration, peakpos, window = integrator.extract_charge(data_ped)
 
-    assert_almost_equal(integration[0][0], 74, 0)
-    assert_almost_equal(integration[1][0], 74, 0)
-    assert peakpos[0][0] == 0
-    assert peakpos[1][0] == 0
+    assert_almost_equal(integration[0], 74, 0)
+    assert peakpos[0] == 0
 
 
 def test_global_peak_integration(test_event):
-    telid = 11
-    event = deepcopy(test_event)
-    data = event.r0.tel[telid].waveform
-    nsamples = data.shape[2]
-    ped = event.mc.tel[telid].pedestal
-    data_ped = data - np.atleast_3d(ped / nsamples)
-    data_ped = np.array([data_ped[0], data_ped[0]])  # Test LG functionality
+    data_ped = dummy_calib(deepcopy(test_event), telid=11)
 
     integrator = GlobalPeakIntegrator()
     integration, peakpos, window = integrator.extract_charge(data_ped)
 
-    assert_almost_equal(integration[0][0], 58, 0)
-    assert_almost_equal(integration[1][0], 58, 0)
-    assert peakpos[0][0] == 14
-    assert peakpos[1][0] == 14
+    assert_almost_equal(integration[0], 58, 0)
+    assert peakpos[0] == 14
 
 
 def test_local_peak_integration(test_event):
-    telid = 11
-    event = deepcopy(test_event)
-    data = event.r0.tel[telid].waveform
-    nsamples = data.shape[2]
-    ped = event.mc.tel[telid].pedestal
-    data_ped = data - np.atleast_3d(ped / nsamples)
-    data_ped = np.array([data_ped[0], data_ped[0]])  # Test LG functionality
+    data_ped = dummy_calib(deepcopy(test_event), telid=11)
 
     integrator = LocalPeakIntegrator()
     integration, peakpos, window = integrator.extract_charge(data_ped)
 
-    assert_almost_equal(integration[0][0], 76, 0)
-    assert_almost_equal(integration[1][0], 76, 0)
-    assert peakpos[0][0] == 13
-    assert peakpos[1][0] == 13
+    assert_almost_equal(integration[0], 76, 0)
+    assert peakpos[0] == 13
 
 
 def test_nb_peak_integration(test_event):
-    telid = 11
-    event = deepcopy(test_event)
-    data = event.r0.tel[telid].waveform
-    nsamples = data.shape[2]
-    ped = event.mc.tel[telid].pedestal
-    data_ped = data - np.atleast_3d(ped / nsamples)
-    data_ped = np.array([data_ped[0], data_ped[0]])  # Test LG functionality
-    geom = event.inst.subarray.tel[telid].camera
+    data_ped = dummy_calib(deepcopy(test_event), telid=11)
+
+    geom = test_event.inst.subarray.tel[11].camera
     nei = geom.neighbor_matrix_where
 
     integrator = NeighbourPeakIntegrator()
     integrator.neighbours = nei
     integration, peakpos, window = integrator.extract_charge(data_ped)
 
-    assert_almost_equal(integration[0][0], -64, 0)
-    assert_almost_equal(integration[1][0], -64, 0)
-    assert peakpos[0][0] == 20
-    assert peakpos[1][0] == 20
+    assert_almost_equal(integration[0], -64, 0)
+    assert peakpos[0] == 20
 
 
 def test_averagewf_peak_integration(test_event):
-    telid = 11
-    event = deepcopy(test_event)
-    data = event.r0.tel[telid].waveform
-    nsamples = data.shape[2]
-    ped = event.mc.tel[telid].pedestal
-    data_ped = data - np.atleast_3d(ped / nsamples)
-    data_ped = np.array([data_ped[0], data_ped[0]])  # Test LG functionality
+    data_ped = dummy_calib(deepcopy(test_event), telid=11)
 
     integrator = AverageWfPeakIntegrator()
     integration, peakpos, window = integrator.extract_charge(data_ped)
 
-    assert_almost_equal(integration[0][0], 73, 0)
-    assert_almost_equal(integration[1][0], 73, 0)
-    assert peakpos[0][0] == 10
-    assert peakpos[1][0] == 10
+    assert_almost_equal(integration[0], 73, 0)
+    assert peakpos[0] == 10
 
 
 def test_charge_extractor_factory(test_event):
@@ -129,13 +92,8 @@ def test_charge_extractor_factory(test_event):
         product='LocalPeakIntegrator'
     )
 
-    telid = 11
-    event = deepcopy(test_event)
-    data = event.r0.tel[telid].waveform
-    nsamples = data.shape[2]
-    ped = event.mc.tel[telid].pedestal
-    data_ped = data - np.atleast_3d(ped / nsamples)
+    data_ped = dummy_calib(deepcopy(test_event), telid=11)
 
     integration, peakpos, window = extractor.extract_charge(data_ped)
 
-    assert_almost_equal(integration[0][0], 76, 0)
+    assert_almost_equal(integration[0], 76, 0)
diff --git a/ctapipe/image/waveform_cleaning.py b/ctapipe/image/waveform_cleaning.py
index 41491ecf196..6e56b5109cd 100644
--- a/ctapipe/image/waveform_cleaning.py
+++ b/ctapipe/image/waveform_cleaning.py
@@ -122,20 +122,18 @@ def get_extractor(self):
         """
 
     def apply(self, waveforms):
-        samples = waveforms[0]
-
         # Subtract initial baseline
-        baseline_sub = samples - np.mean(samples[:, :32], axis=1)[:, None]
+        baseline_sub = waveforms - np.mean(waveforms[:, :32], axis=1)[:, None]
 
         # Obtain waveform with pulse masked
         baseline_sub_b = baseline_sub[None, ...]
         window, _ = self.extractor.get_window_from_waveforms(waveforms)
-        windowed = np.ma.array(baseline_sub_b, mask=window[0])
+        windowed = np.ma.array(baseline_sub_b, mask=window)
         no_pulse = np.ma.filled(windowed, 0)[0]
 
         # Get smooth baseline (no pulse)
         smooth_flat = np.convolve(no_pulse.ravel(), self.kernel, "same")
-        smooth_baseline = np.reshape(smooth_flat, samples.shape)
+        smooth_baseline = np.reshape(smooth_flat, waveforms.shape)
         no_pulse_std = np.std(no_pulse, axis=1)
         smooth_baseline_std = np.std(smooth_baseline, axis=1)
         with np.errstate(divide='ignore', invalid='ignore'):
@@ -148,7 +146,7 @@ def apply(self, waveforms):
         # Subtract smooth baseline
         cleaned = smooth_wf - smooth_baseline
 
-        self.stages['0: raw'] = samples
+        self.stages['0: raw'] = waveforms
         self.stages['1: baseline_sub'] = baseline_sub
         self.stages['window'] = window
         self.stages['2: no_pulse'] = no_pulse
diff --git a/ctapipe/io/containers.py b/ctapipe/io/containers.py
index 71361348844..14964db70ec 100644
--- a/ctapipe/io/containers.py
+++ b/ctapipe/io/containers.py
@@ -103,8 +103,9 @@ class DL1CameraContainer(Container):
         "numpy array containing position of the peak as determined by "
         "the peak-finding algorithm for each pixel"
     )
-    cleaned = Field(
-        None, "numpy array containing the waveform after cleaning"
+
+    waveform = Field(
+        None, "numpy array containing the waveform after waveform-cleaning"
     )
 
 
@@ -160,10 +161,13 @@ class R1CameraContainer(Container):
                                "readout")
     trigger_type = Field(0o0, "camera trigger type")
 
-    waveform = Field(None, (
-        "numpy array containing a set of images, one per ADC sample"
-        "(n_channels x n_pixels, n_samples)"
-    ))
+    waveform_full = Field(None, "Calibrated, but not gain-selected waveforms "
+                                "per pixel (n_chan, n_pixels, n_samples")
+
+    waveform = Field(None, "calibrated images images, one per ADC sample ("
+                           "n_pixels, n_samples)")
+
+    gain_channel = Field(None, "Mask of which gain channel was used")
 
 
 class R1Container(Container):
@@ -192,6 +196,10 @@ class DL0CameraContainer(Container):
         "if pixels or time slices are zero-suppressed"
     ))
 
+    gain_channel = Field(None, "boolean numpy array of which gain channel was "
+                               "used for each pixel in the image ")
+
+
 
 class DL0Container(Container):
     """
diff --git a/ctapipe/io/targetioeventsource.py b/ctapipe/io/targetioeventsource.py
index 55b8d067b57..6c619ec00c6 100644
--- a/ctapipe/io/targetioeventsource.py
+++ b/ctapipe/io/targetioeventsource.py
@@ -103,14 +103,14 @@ def __init__(self, config=None, tool=None, **kwargs):
                 dtype=np.float32
             )
             self._get_tio_event = self._reader.GetR1Event
-            self._samples = self._r1_samples[0]
+            self._samples = self._r1_samples
         else:
             self._r0_samples = np.zeros(
                 (1, n_pix, n_samples),
                 dtype=np.uint16
             )
             self._get_tio_event = self._reader.GetR0Event
-            self._samples = self._r0_samples[0]
+            self._samples = self._r0_samples
 
         self._init_container()
 
diff --git a/ctapipe/tools/extract_charge_resolution.py b/ctapipe/tools/extract_charge_resolution.py
index de5d2b10393..927c346d29b 100644
--- a/ctapipe/tools/extract_charge_resolution.py
+++ b/ctapipe/tools/extract_charge_resolution.py
@@ -35,8 +35,8 @@ class ChargeResolutionGenerator(Tool):
                         window_width='ChargeExtractorFactory.window_width',
                         t0='ChargeExtractorFactory.t0',
                         window_shift='ChargeExtractorFactory.window_shift',
-                        sig_amp_cut_HG='ChargeExtractorFactory.sig_amp_cut_HG',
-                        sig_amp_cut_LG='ChargeExtractorFactory.sig_amp_cut_LG',
+                        peak_detection_threshold='ChargeExtractorFactory'
+                                     '.peak_detection_threshold',
                         lwt='ChargeExtractorFactory.lwt',
                         clip_amplitude='CameraDL1Calibrator.clip_amplitude',
                         radius='CameraDL1Calibrator.radius',
diff --git a/ctapipe/utils/neighbour_sum.py b/ctapipe/utils/neighbour_sum.py
index b1ae90f26f1..1529e255817 100644
--- a/ctapipe/utils/neighbour_sum.py
+++ b/ctapipe/utils/neighbour_sum.py
@@ -3,7 +3,7 @@
 neighbour_sum_c.cc through ctypes.
 """
 
-__all__ = []
+__all__ = ['get_sum_array']
 
 import numpy as np
 import ctypes
@@ -16,7 +16,7 @@
     get_sum_array.restype = None
     get_sum_array.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"),
                               ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"),
-                              ctypes.c_size_t, ctypes.c_size_t, ctypes.c_size_t,
+                              ctypes.c_size_t, ctypes.c_size_t,
                               ndpointer(ctypes.c_uint16, flags="C_CONTIGUOUS"),
                               ctypes.c_size_t,
                               ctypes.c_int]
diff --git a/ctapipe/utils/neighbour_sum_c.cc b/ctapipe/utils/neighbour_sum_c.cc
index 3e5f93a8285..e4f6aafb46f 100644
--- a/ctapipe/utils/neighbour_sum_c.cc
+++ b/ctapipe/utils/neighbour_sum_c.cc
@@ -7,30 +7,30 @@ ctapipe.image.charge_extractors.NeighbourPeakIntegrator.
 #include <stdint.h>
 #include <stdio.h>
 
-extern "C" void get_sum_array(const float* waveforms, float* sum_array, size_t n_chan, size_t n_pix, size_t n_samples, const uint16_t* nei, size_t nei_length, int lwt)
+extern "C" void get_sum_array(const float* waveforms, float* sum_array,
+                              size_t n_pix, size_t n_samples,
+                              const uint16_t* nei, size_t nei_length, int lwt)
 {
-    for (size_t c = 0; c < n_chan; ++c) {
-        if (lwt > 0){
-            for (size_t p = 0; p < n_pix; ++p) {
-                int index = c * n_pix * n_samples + p * n_samples;
-                const float* wf = waveforms + index;
-                float* sum = sum_array + index;
-                for (size_t s = 0; s < n_samples; ++s) {
-                    sum[s] += wf[s] * lwt;
-                }
-            }
-        }
-        for (size_t ni = 0; ni < nei_length; ++ni) {
-            const uint16_t* nei_cur = nei + ni * 2;
-            int p = nei_cur[0];
-            int n = nei_cur[1];
-            int index = c * n_pix * n_samples + p * n_samples;
-            int nei_index = c * n_pix * n_samples + n * n_samples;
-            const float* wfn = waveforms + nei_index;
+    if (lwt > 0){
+        for (size_t p = 0; p < n_pix; ++p) {
+            int index = p * n_samples;
+            const float* wf = waveforms + index;
             float* sum = sum_array + index;
             for (size_t s = 0; s < n_samples; ++s) {
-                sum[s] += wfn[s];
+                sum[s] += wf[s] * lwt;
             }
         }
     }
+    for (size_t ni = 0; ni < nei_length; ++ni) {
+        const uint16_t* nei_cur = nei + ni * 2;
+        int p = nei_cur[0];
+        int n = nei_cur[1];
+        int index = p * n_samples;
+        int nei_index = n * n_samples;
+        const float* wfn = waveforms + nei_index;
+        float* sum = sum_array + index;
+        for (size_t s = 0; s < n_samples; ++s) {
+            sum[s] += wfn[s];
+        }
+    }
 }
diff --git a/examples/calibration_pipeline.py b/examples/calibration_pipeline.py
index 7c07bdfe518..4a01c453add 100644
--- a/examples/calibration_pipeline.py
+++ b/examples/calibration_pipeline.py
@@ -13,15 +13,15 @@
 class ImagePlotter(Component):
     display = Bool(
         False,
-        help='Display the photoelectron images on-screen as they '
-        'are produced.'
+        help='Display the photoelectron images on-screen as they are produced.'
     ).tag(config=True)
     output_path = Unicode(
         None,
         allow_none=True,
-        help='Output path for the pdf containing all the '
-        'images. Set to None for no saved '
-        'output.'
+        help=(
+            'Output path for the pdf containing all the images. Set to None '
+            'for no saved output.'
+        )
     ).tag(config=True)
 
     def __init__(self, config=None, tool=None, **kwargs):
@@ -78,7 +78,7 @@ def plot(self, event, telid):
             self.c_intensity = CameraDisplay(geom, ax=self.ax_intensity)
             self.c_peakpos = CameraDisplay(geom, ax=self.ax_peakpos)
 
-            tmaxmin = event.dl0.tel[telid].waveform.shape[2]
+            tmaxmin = event.dl0.tel[telid].waveform.shape[1]
             t_chargemax = peakpos[image.argmax()]
             cmap_time = colors.LinearSegmentedColormap.from_list(
                 'cmap_t',
@@ -128,8 +128,7 @@ def finish(self):
 
 class DisplayDL1Calib(Tool):
     name = "DisplayDL1Calib"
-    description = "Calibrate dl0 data to dl1, and plot the photoelectron " \
-                  "images."
+    description = "Calibrate dl0 data to dl1, and plot the photoelectron images."
 
     telescope = Int(
         None,
@@ -145,8 +144,7 @@ class DisplayDL1Calib(Tool):
             window_width='ChargeExtractorFactory.window_width',
             t0='ChargeExtractorFactory.t0',
             window_shift='ChargeExtractorFactory.window_shift',
-            sig_amp_cut_HG='ChargeExtractorFactory.sig_amp_cut_HG',
-            sig_amp_cut_LG='ChargeExtractorFactory.sig_amp_cut_LG',
+            peak_detection_threshold='ChargeExtractorFactory.peak_detection_threshold',
             lwt='ChargeExtractorFactory.lwt',
             clip_amplitude='CameraDL1Calibrator.clip_amplitude',
             T='DisplayDL1Calib.telescope',
diff --git a/examples/display_events_single_tel.py b/examples/display_events_single_tel.py
index 5fa5fe2366b..8dded78a9ed 100755
--- a/examples/display_events_single_tel.py
+++ b/examples/display_events_single_tel.py
@@ -31,9 +31,6 @@ class SingleTelEventDisplay(Tool):
 
     infile = Unicode(help="input file to read", default='').tag(config=True)
     tel = Int(help='Telescope ID to display', default=0).tag(config=True)
-    channel = Integer(
-        help="channel number to display", min=0, max=1
-    ).tag(config=True)
     write = Bool(
         help="Write out images to PNG files", default=False
     ).tag(config=True)
@@ -56,7 +53,6 @@ class SingleTelEventDisplay(Tool):
         'infile': 'EventSourceFactory.input_url',
         'tel': 'SingleTelEventDisplay.tel',
         'max-events': 'EventSourceFactory.max_events',
-        'channel': 'SingleTelEventDisplay.channel',
         'write': 'SingleTelEventDisplay.write',
         'clean': 'SingleTelEventDisplay.clean',
         'hillas': 'SingleTelEventDisplay.hillas',
@@ -117,7 +113,7 @@ def start(self):
 
             if self.samples:
                 # display time-varying event
-                data = event.dl0.tel[self.tel].waveform[self.channel]
+                data = event.dl0.tel[self.tel].waveform
                 for ii in range(data.shape[1]):
                     disp.image = data[:, ii]
                     disp.set_limits_percent(70)
@@ -131,7 +127,7 @@ def start(self):
                         )
             else:
                 # display integrated event:
-                im = event.dl1.tel[self.tel].image[self.channel]
+                im = event.dl1.tel[self.tel].image
 
                 if self.clean:
                     mask = tailcuts_clean(
diff --git a/examples/display_integrator.py b/examples/display_integrator.py
index 81e3ca5f344..b71284d0d2b 100644
--- a/examples/display_integrator.py
+++ b/examples/display_integrator.py
@@ -245,24 +245,24 @@ def plot(event, telid, chan, extractor_name):
 
 class DisplayIntegrator(Tool):
     name = "DisplayIntegrator"
-    description = "Calibrate dl0 data to dl1, and plot the various camera " \
-                  "images that characterise the event and calibration. Also " \
-                  "plot some examples of waveforms with the " \
-                  "integration window."
+    description = ("Calibrate dl0 data to dl1, and plot the various camera " 
+                   "images that characterise the event and calibration. Also " 
+                   "plot some examples of waveforms with the " 
+                   "integration window.")
 
     event_index = Int(0, help='Event index to view.').tag(config=True)
     use_event_id = Bool(
         False,
-        help='event_index will obtain an event using'
-        'event_id instead of '
-        'index.'
+        help='event_index will obtain an event using event_id instead of index.'
     ).tag(config=True)
+
     telescope = Int(
         None,
         allow_none=True,
         help='Telescope to view. Set to None to display the first'
         'telescope with data.'
     ).tag(config=True)
+
     channel = Enum([0, 1], 0, help='Channel to view').tag(config=True)
 
     aliases = Dict(
@@ -273,8 +273,8 @@ class DisplayIntegrator(Tool):
             extractor='ChargeExtractorFactory.product',
             window_width='ChargeExtractorFactory.window_width',
             window_shift='ChargeExtractorFactory.window_shift',
-            sig_amp_cut_HG='ChargeExtractorFactory.sig_amp_cut_HG',
-            sig_amp_cut_LG='ChargeExtractorFactory.sig_amp_cut_LG',
+            peak_detection_threshold='ChargeExtractorFactory'
+            '.peak_detection_threshold',
             lwt='ChargeExtractorFactory.lwt',
             clip_amplitude='CameraDL1Calibrator.clip_amplitude',
             radius='CameraDL1Calibrator.radius',
@@ -341,10 +341,8 @@ def start(self):
         if telid is None:
             telid = tels[0]
         if telid not in tels:
-            self.log.error(
-                "[event] please specify one of the following "
-                "telescopes for this event: {}".format(tels)
-            )
+            self.log.error("please specify one of the following"
+                           " telescopes for this event: %s", tels)
             exit()
 
         extractor_name = self.extractor.__class__.__name__
diff --git a/examples/display_summed_images.py b/examples/display_summed_images.py
index a0b3e0078bd..cd8dc24b29f 100644
--- a/examples/display_summed_images.py
+++ b/examples/display_summed_images.py
@@ -98,7 +98,7 @@ def start(self):
 
             imsum[:] = 0
             for telid in event.dl0.tels_with_data:
-                imsum += event.dl1.tel[telid].image[0]
+                imsum += event.dl1.tel[telid].image
 
             self.log.info(
                 "event={} ntels={} energy={}".format(
diff --git a/examples/plot_theta_square.py b/examples/plot_theta_square.py
index a3216a651f0..93d8eb33dfd 100644
--- a/examples/plot_theta_square.py
+++ b/examples/plot_theta_square.py
@@ -7,7 +7,8 @@
 import matplotlib.pyplot as plt
 import numpy as np
 from astropy import units as u
-from astropy.coordinates.angle_utilities import angular_separation
+import sys
+from tqdm import tqdm
 
 from ctapipe.calib import CameraCalibrator
 from ctapipe.image import hillas_parameters
@@ -16,14 +17,12 @@
 from ctapipe.reco import HillasReconstructor
 from ctapipe.utils import datasets
 
-
 if len(sys.argv) >= 2:
     filename = sys.argv[1]
 else:
     # importing data from avaiable datasets in ctapipe
     filename = datasets.get_dataset_path("gamma_test_large.simtel.gz")
 
-
 # reading the Monte Carlo file for LST
 source = event_source(filename, allowed_tels={1, 2, 3, 4})
 
@@ -31,7 +30,7 @@
 calib = CameraCalibrator(r1_product="HESSIOR1Calibrator")
 off_angles = []
 
-for event in source:
+for event in tqdm(source):
 
     # calibrating the event
     calib.calibrate(event)
@@ -53,7 +52,7 @@
         camgeom = subarray.tel[tel_id].camera
 
         # note the [0] is for channel 0 which is high-gain channel
-        image = event.dl1.tel[tel_id].image[0]
+        image = event.dl1.tel[tel_id].image
 
         # Cleaning  of the image
         cleaned_image = image
diff --git a/examples/simple_pipeline.py b/examples/simple_pipeline.py
index 3a15bab7bb1..570a49eb144 100644
--- a/examples/simple_pipeline.py
+++ b/examples/simple_pipeline.py
@@ -11,8 +11,12 @@
 if __name__ == '__main__':
 
     filename = sys.argv[1]
+    try:
+        max_events = int(sys.argv[2])
+    except IndexError:
+        max_events = None
 
-    source = event_source(filename, max_events=None)
+    source = event_source(filename, max_events=max_events)
 
     cal = CameraCalibrator(r1_product="HESSIOR1Calibrator")