diff --git a/ctapipe/calib/camera/calibrator.py b/ctapipe/calib/camera/calibrator.py
index dfc63a67bbc..483e2c5052e 100644
--- a/ctapipe/calib/camera/calibrator.py
+++ b/ctapipe/calib/camera/calibrator.py
@@ -1,5 +1,6 @@
 import numpy as np
 from ctapipe.core import Component
+from ctapipe.calib.camera.gainselection import ManualGainSelector
 from ctapipe.image.reducer import NullDataVolumeReducer
 from ctapipe.image.extractor import NeighborPeakWindowSum
 import warnings
@@ -73,6 +74,7 @@ class CameraCalibrator(Component):
     the DL1 data level in the event container.
     """
     def __init__(self, config=None, parent=None,
+                 gain_selector=None,
                  data_volume_reducer=None,
                  image_extractor=None,
                  **kwargs):
@@ -87,6 +89,9 @@ def __init__(self, config=None, parent=None,
             Tool executable that is calling this component.
             Passes the correct logger to the component.
             Set to None if no Tool to pass.
+        gain_selector : ctapipe.calib.camera.gainselection.GainSelector
+            The GainSelector to use. If None, then ManualGainSelector will be
+            used, which by default selects the high/first gain channel.
         data_volume_reducer : ctapipe.image.reducer.DataVolumeReducer
             The DataVolumeReducer to use. If None, then
             NullDataVolumeReducer will be used by default, and waveforms
@@ -101,6 +106,10 @@ def __init__(self, config=None, parent=None,
         self._r1_empty_warn = False
         self._dl0_empty_warn = False
 
+        if gain_selector is None:
+            gain_selector = ManualGainSelector(parent=self)
+        self.gain_selector = gain_selector
+
         if data_volume_reducer is None:
             data_volume_reducer = NullDataVolumeReducer(parent=self)
         self.data_volume_reducer = data_volume_reducer
@@ -126,6 +135,7 @@ def _get_correction(self, event, telid):
         ndarray
         """
         try:
+            selected_gain_channel = event.r1.tel[telid].selected_gain_channel
             shift = self.image_extractor.window_shift
             width = self.image_extractor.window_width
             shape = event.mc.tel[telid].reference_pulse_shape
@@ -134,7 +144,8 @@ def _get_correction(self, event, telid):
             time_slice = event.mc.tel[telid].time_slice
             correction = integration_correction(n_chan, shape, step,
                                                 time_slice, width, shift)
-            return correction
+            pixel_correction = correction[selected_gain_channel]
+            return pixel_correction
         except (AttributeError, KeyError):
             # Don't apply correction when window_shift or window_width
             # does not exist in extractor, or when container does not have
@@ -165,15 +176,34 @@ def _calibrate_dl0(self, event, telid):
         waveforms = event.r1.tel[telid].waveform
         if self._check_r1_empty(waveforms):
             return
-        # TODO: Add gain selection
-        reduced_waveforms = self.data_volume_reducer(waveforms)
+
+        # Perform gain selection. This is typically not the responsibility of
+        # ctapipe; DL0 (and R1) waveforms are aleady gain selected and
+        # therefore single channel. However, the waveforms read from
+        # simtelarray do not have the gain selection applied, and so must be
+        # done as part of the calibration step to ensure the correct
+        # waveform dimensions.
+        waveforms_gs, selected_gain_channel = self.gain_selector(waveforms)
+        if selected_gain_channel is not None:
+            event.r1.tel[telid].selected_gain_channel = selected_gain_channel
+        else:
+            # If pixel_channel is None, then waveforms has already been
+            # pre-gainselected, and presumably the selected_gain_channel
+            # container is filled by the EventSource
+            if event.r1.tel[telid].selected_gain_channel is None:
+                raise ValueError(
+                    "EventSource is loading pre-gainselected waveforms "
+                    "without filling the selected_gain_channel container"
+                )
+
+        reduced_waveforms = self.data_volume_reducer(waveforms_gs)
         event.dl0.tel[telid].waveform = reduced_waveforms
 
     def _calibrate_dl1(self, event, telid):
         waveforms = event.dl0.tel[telid].waveform
         if self._check_dl0_empty(waveforms):
             return
-        n_samples = waveforms.shape[2]
+        n_pixels, n_samples = waveforms.shape
         if n_samples == 1:
             # To handle ASTRI and dst
             # TODO: Improved handling of ASTRI and dst
@@ -182,7 +212,7 @@ def _calibrate_dl1(self, event, telid):
             #   - Don't do anything if dl1 container already filled
             #   - Update on SST review decision
             corrected_charge = waveforms[..., 0]
-            pulse_time = np.zeros(waveforms.shape[0:2])
+            pulse_time = np.zeros(n_pixels)
         else:
             # TODO: pass camera to ImageExtractor.__init__
             if self.image_extractor.requires_neighbors():
@@ -192,7 +222,7 @@ def _calibrate_dl1(self, event, telid):
 
             # Apply integration correction
             # TODO: Remove integration correction
-            correction = self._get_correction(event, telid)[:, None]
+            correction = self._get_correction(event, telid)
             corrected_charge = charge * correction
 
         event.dl1.tel[telid].image = corrected_charge
diff --git a/ctapipe/calib/camera/gainselection.py b/ctapipe/calib/camera/gainselection.py
index 0a414cf6481..4dbaa723f54 100644
--- a/ctapipe/calib/camera/gainselection.py
+++ b/ctapipe/calib/camera/gainselection.py
@@ -45,17 +45,19 @@ def __call__(self, waveforms):
             Shape: (n_pix, n_samples)
         """
         if waveforms.ndim == 2:  # Return if already gain selected
-            pixel_channel = None  # Provided by EventSource
-            return waveforms, pixel_channel
+            selected_gain_channel = None  # Provided by EventSource
+            return waveforms, selected_gain_channel
         elif waveforms.ndim == 3:
             n_channels, n_pixels, _ = waveforms.shape
             if n_channels == 1:  # Reduce if already single channel
-                pixel_channel = np.zeros(n_pixels, dtype=int)
-                return waveforms[0], pixel_channel
+                selected_gain_channel = np.zeros(n_pixels, dtype=int)
+                return waveforms[0], selected_gain_channel
             else:
-                pixel_channel = self.select_channel(waveforms)
-                gain_selected = waveforms[pixel_channel, np.arange(n_pixels)]
-                return gain_selected, pixel_channel
+                selected_gain_channel = self.select_channel(waveforms)
+                selected_gain_waveforms = waveforms[
+                    selected_gain_channel, np.arange(n_pixels)
+                ]
+                return selected_gain_waveforms, selected_gain_channel
         else:
             raise ValueError(
                 f"Cannot handle waveform array of shape: {waveforms.ndim}"
@@ -77,7 +79,7 @@ def select_channel(self, waveforms):
 
         Returns
         -------
-        pixel_channel : ndarray
+        selected_gain_channel : ndarray
             Gain channel to use for each pixel
             Shape: n_pix
             Dtype: int
diff --git a/ctapipe/calib/camera/tests/test_calibrator.py b/ctapipe/calib/camera/tests/test_calibrator.py
index b1f624f09ac..6677d5505f2 100644
--- a/ctapipe/calib/camera/tests/test_calibrator.py
+++ b/ctapipe/calib/camera/tests/test_calibrator.py
@@ -2,9 +2,11 @@
     CameraCalibrator,
     integration_correction,
 )
+from ctapipe.io.containers import DataContainer
 from ctapipe.image.extractor import LocalPeakWindowSum
 from traitlets.config.configurable import Config
 import pytest
+import numpy as np
 
 
 def test_camera_calibrator(example_event):
@@ -12,7 +14,36 @@ def test_camera_calibrator(example_event):
     calibrator = CameraCalibrator()
     calibrator(example_event)
     image = example_event.dl1.tel[telid].image
+    pulse_time = example_event.dl1.tel[telid].pulse_time
     assert image is not None
+    assert pulse_time is not None
+    assert image.shape == (1764,)
+    assert pulse_time.shape == (1764,)
+
+
+def test_select_gain():
+    n_channels = 2
+    n_pixels = 2048
+    n_samples = 128
+    telid = 0
+
+    calibrator = CameraCalibrator()
+
+    event = DataContainer()
+    event.r1.tel[telid].waveform = np.ones((n_channels, n_pixels, n_samples))
+    calibrator._calibrate_dl0(event, telid)
+    assert event.dl0.tel[telid].waveform.shape == (n_pixels, n_samples)
+
+    event = DataContainer()
+    event.r1.tel[telid].waveform = np.ones((n_pixels, n_samples))
+    with pytest.raises(ValueError):
+        calibrator._calibrate_dl0(event, telid)
+
+    event = DataContainer()
+    event.r1.tel[telid].waveform = np.ones((n_pixels, n_samples))
+    event.r1.tel[telid].selected_gain_channel = np.zeros(n_pixels)
+    calibrator._calibrate_dl0(event, telid)
+    assert event.dl0.tel[telid].waveform.shape == (n_pixels, n_samples)
 
 
 def test_manual_extractor():
@@ -55,7 +86,7 @@ def test_integration_correction_no_ref_pulse(example_event):
     calibrator = CameraCalibrator()
     calibrator._calibrate_dl0(example_event, telid)
     correction = calibrator._get_correction(example_event, telid)
-    assert correction[0] == 1
+    assert (correction == 1).all()
 
 
 def test_check_r1_empty(example_event):
diff --git a/ctapipe/calib/camera/tests/test_gainselection.py b/ctapipe/calib/camera/tests/test_gainselection.py
index 9983a3bac00..bd18e99213d 100644
--- a/ctapipe/calib/camera/tests/test_gainselection.py
+++ b/ctapipe/calib/camera/tests/test_gainselection.py
@@ -14,9 +14,9 @@ def test_gain_selector():
     waveforms[1] *= 2
 
     gain_selector = TestGainSelector()
-    waveforms_gs, pixel_channel = gain_selector(waveforms)
+    waveforms_gs, selected_gain_channel = gain_selector(waveforms)
     np.testing.assert_equal(waveforms[GainChannel.HIGH], waveforms_gs)
-    np.testing.assert_equal(pixel_channel, 0)
+    np.testing.assert_equal(selected_gain_channel, 0)
 
 
 def test_pre_selected():
@@ -24,9 +24,9 @@ def test_pre_selected():
     waveforms = np.zeros(shape)
 
     gain_selector = TestGainSelector()
-    waveforms_gs, pixel_channel = gain_selector(waveforms)
+    waveforms_gs, selected_gain_channel = gain_selector(waveforms)
     assert waveforms.shape == waveforms_gs.shape
-    assert pixel_channel is None
+    assert selected_gain_channel is None
 
 
 def test_single_channel():
@@ -34,9 +34,9 @@ def test_single_channel():
     waveforms = np.zeros(shape)
 
     gain_selector = TestGainSelector()
-    waveforms_gs, pixel_channel = gain_selector(waveforms)
+    waveforms_gs, selected_gain_channel = gain_selector(waveforms)
     assert waveforms_gs.shape == (2048, 128)
-    assert (pixel_channel == 0).all()
+    assert (selected_gain_channel == 0).all()
 
 
 def test_manual_gain_selector():
@@ -45,14 +45,14 @@ def test_manual_gain_selector():
     waveforms[1] *= 2
 
     gs_high = ManualGainSelector(channel="HIGH")
-    waveforms_gs, pixel_channel = gs_high(waveforms)
+    waveforms_gs, selected_gain_channel = gs_high(waveforms)
     np.testing.assert_equal(waveforms[GainChannel.HIGH], waveforms_gs)
-    np.testing.assert_equal(pixel_channel, 0)
+    np.testing.assert_equal(selected_gain_channel, 0)
 
     gs_low = ManualGainSelector(channel="LOW")
-    waveforms_gs, pixel_channel = gs_low(waveforms)
+    waveforms_gs, selected_gain_channel = gs_low(waveforms)
     np.testing.assert_equal(waveforms[GainChannel.LOW], waveforms_gs)
-    np.testing.assert_equal(pixel_channel, 1)
+    np.testing.assert_equal(selected_gain_channel, 1)
 
 
 def test_threshold_gain_selector():
@@ -62,8 +62,8 @@ def test_threshold_gain_selector():
     waveforms[0, 0] = 100
 
     gain_selector = ThresholdGainSelector(threshold=50)
-    waveforms_gs, pixel_channel = gain_selector(waveforms)
+    waveforms_gs, selected_gain_channel = gain_selector(waveforms)
     assert (waveforms_gs[0] == 1).all()
     assert (waveforms_gs[np.arange(1, 2048)] == 0).all()
-    assert pixel_channel[0] == 1
-    assert (pixel_channel[np.arange(1, 2048)] == 0).all()
+    assert selected_gain_channel[0] == 1
+    assert (selected_gain_channel[np.arange(1, 2048)] == 0).all()
diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index e09fd11a367..68cb72f247d 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -22,6 +22,7 @@
 from traitlets import Int
 from ctapipe.core import Component
 from numba import njit, prange, guvectorize, float64, float32, int64, int32
+import warnings
 
 
 @guvectorize(
@@ -51,7 +52,7 @@ def sum_samples_around_peak(waveforms, peak_index, width, shift, ret):
     ----------
     waveforms : ndarray
         Waveforms stored in a numpy array.
-        Shape: (n_chan, n_pix, n_samples)
+        Shape: (n_pix, n_samples)
     peak_index : ndarray or int
         Peak index for each pixel.
     width : ndarray or int
@@ -65,7 +66,7 @@ def sum_samples_around_peak(waveforms, peak_index, width, shift, ret):
     -------
     charge : ndarray
         Extracted charge.
-        Shape: (n_chan, n_pix)
+        Shape: (n_pix)
 
     """
     n_samples = waveforms.size
@@ -78,12 +79,16 @@ def sum_samples_around_peak(waveforms, peak_index, width, shift, ret):
 
 
 @njit([
-    float64[:, :, :](float64[:, :, :], int64[:, :], int64),
-    float64[:, :, :](float64[:, :, :], int32[:, :], int64),
-    float64[:, :, :](float32[:, :, :], int32[:, :], int64),
-    float64[:, :, :](float32[:, :, :], int64[:, :], int64),
+    float64[:, :](float64[:, :], int64[:, :], int64),
+    float64[:, :](float64[:, :], int32[:, :], int64),
+    float64[:, :](float32[:, :], int32[:, :], int64),
+    float64[:, :](float32[:, :], int64[:, :], int64),
+    float64[:, :, :](float64[:, :, :], int64[:, :], int64),  # TODO: deprecate
+    float64[:, :, :](float64[:, :, :], int32[:, :], int64),  # TODO: deprecate
+    float64[:, :, :](float32[:, :, :], int32[:, :], int64),  # TODO: deprecate
+    float64[:, :, :](float32[:, :, :], int64[:, :], int64),  # TODO: deprecate
 ], parallel=True)
-def neighbor_average_waveform(waveforms, neighbors, lwt):
+def neighbor_average_waveform_jit(waveforms, neighbors, lwt):
     """
     Obtain the average waveform built from the neighbors of each pixel
 
@@ -91,7 +96,7 @@ def neighbor_average_waveform(waveforms, neighbors, lwt):
     ----------
     waveforms : ndarray
         Waveforms stored in a numpy array.
-        Shape: (n_chan, n_pix, n_samples)
+        Shape: (n_pix, n_samples)
     neighbors : ndarray
         2D array where each row is [pixel index, one neighbor of that pixel].
         Changes per telescope.
@@ -105,21 +110,45 @@ def neighbor_average_waveform(waveforms, neighbors, lwt):
     -------
     average_wf : ndarray
         Average of neighbor waveforms for each pixel.
-        Shape: (n_chan, n_pix, n_samples)
+        Shape: (n_pix, n_samples)
 
     """
+    # TODO: deprecate
+    if waveforms.ndim == 3:
+        n_neighbors = neighbors.shape[0]
+        sum_ = waveforms * lwt
+        n = np.zeros(waveforms.shape)
+        for i in prange(n_neighbors):
+            pixel = neighbors[i, 0]
+            neighbor = neighbors[i, 1]
+            for channel in range(waveforms.shape[0]):
+                sum_[channel, pixel] += waveforms[channel, neighbor]
+                n[channel, pixel] += 1
+        return sum_ / n
+
     n_neighbors = neighbors.shape[0]
     sum_ = waveforms * lwt
     n = np.zeros(waveforms.shape)
     for i in prange(n_neighbors):
         pixel = neighbors[i, 0]
         neighbor = neighbors[i, 1]
-        for channel in range(waveforms.shape[0]):
-            sum_[channel, pixel] += waveforms[channel, neighbor]
-            n[channel, pixel] += 1
+        sum_[pixel] += waveforms[neighbor]
+        n[pixel] += 1
     return sum_ / n
 
 
+def neighbor_average_waveform(waveforms, neighbors, lwt):
+    if waveforms.ndim == 3:
+        warnings.warn(
+            "A 3 dimensional waveforms array was passed "
+            "to neighbor_average_waveform. "
+            "Waveforms should be shape (n_pixels, n_samples). "
+            "This will raise an error in future versions.",
+            DeprecationWarning
+        )
+    return neighbor_average_waveform_jit(waveforms, neighbors, lwt)
+
+
 @guvectorize(
     [
         (float64[:], int64, int64, int64, float64[:]),
@@ -147,7 +176,7 @@ def extract_pulse_time_around_peak(waveforms, peak_index, width, shift, ret):
     ----------
     waveforms : ndarray
         Waveforms stored in a numpy array.
-        Shape: (n_chan, n_pix, n_samples)
+        Shape: (n_pix, n_samples)
     peak_index : ndarray or int
         Peak index in waveform for each pixel.
     width : ndarray or int
@@ -161,7 +190,7 @@ def extract_pulse_time_around_peak(waveforms, peak_index, width, shift, ret):
     -------
     pulse_time : ndarray
         Floating point pulse time in each pixel
-        Shape: (n_chan, n_pix)
+        Shape: (n_pix)
 
     """
     n_samples = waveforms.size
@@ -188,7 +217,7 @@ def subtract_baseline(waveforms, baseline_start, baseline_end):
     ----------
     waveforms : ndarray
         Waveforms stored in a numpy array.
-        Shape: (n_chan, n_pix, n_samples)
+        Shape: (n_pix, n_samples)
     baseline_start : int
         Sample where the baseline window starts
     baseline_end : int
@@ -200,7 +229,7 @@ def subtract_baseline(waveforms, baseline_start, baseline_end):
         Waveform with the baseline subtracted
     """
     baseline_corrected = waveforms - np.mean(
-        waveforms[..., baseline_start:baseline_end], axis=2
+        waveforms[..., baseline_start:baseline_end], axis=-1
     )[..., None]
 
     return baseline_corrected
@@ -274,16 +303,16 @@ def __call__(self, waveforms):
         ----------
         waveforms : ndarray
             Waveforms stored in a numpy array of shape
-            (n_chan, n_pix, n_samples).
+            (n_pix, n_samples).
 
         Returns
         -------
         charge : ndarray
             Extracted charge.
-            Shape: (n_chan, n_pix)
+            Shape: (n_pix)
         pulse_time : ndarray
             Floating point pulse time in each pixel.
-            Shape: (n_chan, n_pix)
+            Shape: (n_pix)
         """
 
 
@@ -293,7 +322,7 @@ class FullWaveformSum(ImageExtractor):
     """
 
     def __call__(self, waveforms):
-        charge = waveforms.sum(2)
+        charge = waveforms.sum(axis=-1)
         pulse_time = extract_pulse_time_around_peak(
             waveforms, 0, waveforms.shape[-1], 0
         )
@@ -314,7 +343,7 @@ class FixedWindowSum(ImageExtractor):
     def __call__(self, waveforms):
         start = self.window_start
         end = self.window_start + self.window_width
-        charge = waveforms[..., start:end].sum(2)
+        charge = waveforms[..., start:end].sum(axis=-1)
         pulse_time = extract_pulse_time_around_peak(
             waveforms, self.window_start, self.window_width, 0
         )
@@ -335,13 +364,22 @@ class GlobalPeakWindowSum(ImageExtractor):
     ).tag(config=True)
 
     def __call__(self, waveforms):
-        peak_index = waveforms.mean(1).argmax(1)
+        peak_index = waveforms.mean(axis=-2).argmax(axis=-1)
+        if waveforms.ndim == 3:
+            warnings.warn(
+                "A 3 dimensional waveforms array was passed "
+                "to GlobalPeakWindowSum. "
+                "Waveforms should be shape (n_pixels, n_samples). "
+                "This will raise an error in future versions.",
+                DeprecationWarning
+            )
+            peak_index = peak_index[:, np.newaxis]
         charge = sum_samples_around_peak(
-            waveforms, peak_index[:, np.newaxis],
+            waveforms, peak_index,
             self.window_width, self.window_shift
         )
         pulse_time = extract_pulse_time_around_peak(
-            waveforms, peak_index[:, np.newaxis],
+            waveforms, peak_index,
             self.window_width, self.window_shift
         )
         return charge, pulse_time
@@ -361,7 +399,7 @@ class LocalPeakWindowSum(ImageExtractor):
     ).tag(config=True)
 
     def __call__(self, waveforms):
-        peak_index = waveforms.argmax(2).astype(np.int)
+        peak_index = waveforms.argmax(axis=-1).astype(np.int)
         charge = sum_samples_around_peak(
             waveforms, peak_index, self.window_width, self.window_shift
         )
@@ -395,7 +433,7 @@ def __call__(self, waveforms):
         average_wfs = neighbor_average_waveform(
             waveforms, self.neighbors, self.lwt
         )
-        peak_index = average_wfs.argmax(2)
+        peak_index = average_wfs.argmax(axis=-1)
         charge = sum_samples_around_peak(
             waveforms, peak_index, self.window_width, self.window_shift
         )
diff --git a/ctapipe/image/muon/muon_diagnostic_plots.py b/ctapipe/image/muon/muon_diagnostic_plots.py
index 37ecc526c1f..5b995f76edd 100644
--- a/ctapipe/image/muon/muon_diagnostic_plots.py
+++ b/ctapipe/image/muon/muon_diagnostic_plots.py
@@ -113,7 +113,7 @@ def plot_muon_event(event, muonparams):
             # from the calibration
             ax1 = fig.add_subplot(1, npads, 1)
             plotter = CameraPlotter(event)
-            image = event.dl1.tel[tel_id].image[0]
+            image = event.dl1.tel[tel_id].image
             geom = event.inst.subarray.tel[tel_id].camera
 
             tailcuts = (5., 7.)
diff --git a/ctapipe/image/muon/muon_reco_functions.py b/ctapipe/image/muon/muon_reco_functions.py
index b85ebcecc42..55dc9d0d233 100644
--- a/ctapipe/image/muon/muon_reco_functions.py
+++ b/ctapipe/image/muon/muon_reco_functions.py
@@ -76,7 +76,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/reducer.py b/ctapipe/image/reducer.py
index 54b22bd4aaf..d494754cb6f 100644
--- a/ctapipe/image/reducer.py
+++ b/ctapipe/image/reducer.py
@@ -27,13 +27,13 @@ def __call__(self, waveforms):
         ----------
         waveforms : ndarray
             Waveforms stored in a numpy array of shape
-            (n_chan, n_pix, n_samples).
+            (n_pix, n_samples).
 
         Returns
         -------
         reduced_waveforms : ndarray
             Reduced waveforms stored in a numpy array of shape
-            (n_chan, n_pix, n_samples).
+            (n_pix, n_samples).
         """
 
 
diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 6045b6cdd7a..c27cdc317b5 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -22,6 +22,204 @@
 def camera_waveforms():
     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
+    t_pulse = random.uniform(mid - 10, mid + 10, n_pixels)[:, np.newaxis]
+
+    # Create pulses
+    y = norm.pdf(x, t_pulse, pulse_sigma)
+
+    # Randomize amplitudes
+    y *= random.uniform(100, 1000, n_pixels)[:, np.newaxis]
+
+    return y, camera
+
+
+def test_sum_samples_around_peak(camera_waveforms):
+    waveforms, _ = camera_waveforms
+    n_pixels, n_samples = waveforms.shape
+    rand = np.random.RandomState(1)
+    peak_index = rand.uniform(0, n_samples, n_pixels).astype(np.int)
+    charge = sum_samples_around_peak(waveforms, peak_index, 7, 3)
+    assert_allclose(charge[0], 146.022991, rtol=1e-3)
+
+
+def test_sum_samples_around_peak_expected(camera_waveforms):
+    waveforms, _ = camera_waveforms
+    waveforms = np.ones(waveforms.shape)
+    n_samples = waveforms.shape[-1]
+
+    peak_index = 0
+    width = 10
+    shift = 0
+    charge = sum_samples_around_peak(waveforms, peak_index, width, shift)
+    assert_equal(charge, 10)
+
+    peak_index = 0
+    width = 10
+    shift = 10
+    charge = sum_samples_around_peak(waveforms, peak_index, width, shift)
+    assert_equal(charge, 0)
+
+    peak_index = 0
+    width = 20
+    shift = 10
+    charge = sum_samples_around_peak(waveforms, peak_index, width, shift)
+    assert_equal(charge, 10)
+
+    peak_index = n_samples
+    width = 10
+    shift = 0
+    charge = sum_samples_around_peak(waveforms, peak_index, width, shift)
+    assert_equal(charge, 0)
+
+    peak_index = n_samples
+    width = 20
+    shift = 10
+    charge = sum_samples_around_peak(waveforms, peak_index, width, shift)
+    assert_equal(charge, 10)
+
+    peak_index = 0
+    width = n_samples*3
+    shift = n_samples
+    charge = sum_samples_around_peak(waveforms, peak_index, width, shift)
+    assert_equal(charge, n_samples)
+
+
+def test_neighbor_average_waveform(camera_waveforms):
+    waveforms, camera = camera_waveforms
+    nei = camera.neighbor_matrix_where
+    average_wf = neighbor_average_waveform(waveforms, nei, 0)
+    assert_allclose(average_wf[0, 48], 28.690154, rtol=1e-3)
+
+    average_wf = neighbor_average_waveform(waveforms, nei, 4)
+    assert_allclose(average_wf[0, 48], 98.565743, rtol=1e-3)
+
+
+def test_extract_pulse_time_around_peak(camera_waveforms):
+    x = np.arange(100)
+    y = norm.pdf(x, 41.2, 6)
+    pulse_time = extract_pulse_time_around_peak(
+        y[np.newaxis, :], 0, x.size, 0
+    )
+    assert_allclose(pulse_time[0], 41.2, rtol=1e-3)
+
+
+def test_baseline_subtractor(camera_waveforms):
+    waveforms, _ = camera_waveforms
+    n_pixels, _ = waveforms.shape
+    rand = np.random.RandomState(1)
+    offset = np.arange(n_pixels)[:, np.newaxis]
+    waveforms = rand.normal(0, 0.1, waveforms.shape) + offset
+    assert_allclose(waveforms[3].mean(), 3, rtol=1e-2)
+    baseline_subtracted = subtract_baseline(waveforms, 0, 10)
+    assert_allclose(baseline_subtracted.mean(), 0, atol=1e-3)
+
+
+def test_full_waveform_sum(camera_waveforms):
+    waveforms, _ = camera_waveforms
+    extractor = FullWaveformSum()
+    charge, pulse_time = extractor(waveforms)
+    assert_allclose(charge[0], 545.945, rtol=1e-3)
+    assert_allclose(pulse_time[0], 46.34044, rtol=1e-3)
+
+
+def test_fixed_window_sum(camera_waveforms):
+    waveforms, _ = camera_waveforms
+    extractor = FixedWindowSum(window_start=45)
+    charge, pulse_time = extractor(waveforms)
+    assert_allclose(charge[0], 232.559, rtol=1e-3)
+    assert_allclose(pulse_time[0], 47.823488, rtol=1e-3)
+
+
+def test_global_peak_window_sum(camera_waveforms):
+    waveforms, _ = camera_waveforms
+    extractor = GlobalPeakWindowSum()
+    charge, pulse_time = extractor(waveforms)
+    assert_allclose(charge[0], 232.559, rtol=1e-3)
+    assert_allclose(pulse_time[0], 47.823488, rtol=1e-3)
+
+
+def test_local_peak_window_sum(camera_waveforms):
+    waveforms, _ = camera_waveforms
+    extractor = LocalPeakWindowSum()
+    charge, pulse_time = extractor(waveforms)
+    assert_allclose(charge[0], 240.3, rtol=1e-3)
+    assert_allclose(pulse_time[0], 46.036266, rtol=1e-3)
+
+
+def test_neighbor_peak_window_sum(camera_waveforms):
+    waveforms, camera = camera_waveforms
+    nei = camera.neighbor_matrix_where
+    extractor = NeighborPeakWindowSum()
+    extractor.neighbors = nei
+    charge, pulse_time = extractor(waveforms)
+    assert_allclose(charge[0], 94.671, rtol=1e-3)
+    assert_allclose(pulse_time[0], 54.116092, rtol=1e-3)
+
+    extractor.lwt = 4
+    charge, pulse_time = extractor(waveforms)
+    assert_allclose(charge[0], 220.418657, rtol=1e-3)
+    assert_allclose(pulse_time[0], 48.717848, rtol=1e-3)
+
+
+def test_baseline_subtracted_neighbor_peak_window_sum(camera_waveforms):
+    waveforms, camera = camera_waveforms
+    nei = camera.neighbor_matrix_where
+    extractor = BaselineSubtractedNeighborPeakWindowSum()
+    extractor.neighbors = nei
+    charge, pulse_time = extractor(waveforms)
+    assert_allclose(charge[0], 94.671, rtol=1e-3)
+    assert_allclose(pulse_time[0], 54.116092, rtol=1e-3)
+
+
+def test_waveform_extractor_factory(camera_waveforms):
+    waveforms, _ = camera_waveforms
+    extractor = ImageExtractor.from_name('LocalPeakWindowSum')
+    extractor(waveforms)
+
+
+def test_waveform_extractor_factory_args():
+    """
+    Config is supposed to be created by a `Tool`
+    """
+    from traitlets.config.loader import Config
+    config = Config(
+        {
+            'ImageExtractor': {
+                'window_width': 20,
+                'window_shift': 3,
+            }
+        }
+    )
+
+    extractor = ImageExtractor.from_name(
+        'LocalPeakWindowSum',
+        config=config,
+    )
+    assert extractor.window_width == 20
+    assert extractor.window_shift == 3
+
+    with pytest.warns(UserWarning):
+        ImageExtractor.from_name(
+            'FullWaveformSum',
+            config=config,
+        )
+
+# TODO: Deperacate 2 channel waveforms
+
+
+@pytest.fixture(scope='module')
+def camera_waveforms_2chan():
+    camera = CameraGeometry.from_name("CHEC")
+
     n_pixels = camera.n_pixels
     n_samples = 96
     mid = n_samples // 2
@@ -48,8 +246,8 @@ def camera_waveforms():
     return y, camera
 
 
-def test_sum_samples_around_peak(camera_waveforms):
-    waveforms, _ = camera_waveforms
+def test_sum_samples_around_peak_2chan(camera_waveforms_2chan):
+    waveforms, _ = camera_waveforms_2chan
     _, n_pixels, n_samples = waveforms.shape
     rand = np.random.RandomState(1)
     peak_index = rand.uniform(0, n_samples, (2, n_pixels)).astype(np.int)
@@ -59,8 +257,8 @@ def test_sum_samples_around_peak(camera_waveforms):
     assert_allclose(charge[1][0], 22.393974, rtol=1e-3)
 
 
-def test_sum_samples_around_peak_expected(camera_waveforms):
-    waveforms, _ = camera_waveforms
+def test_sum_samples_around_peak_expected_2chan(camera_waveforms_2chan):
+    waveforms, _ = camera_waveforms_2chan
     waveforms = np.ones(waveforms.shape)
     n_samples = waveforms.shape[-1]
 
@@ -101,8 +299,8 @@ def test_sum_samples_around_peak_expected(camera_waveforms):
     assert_equal(charge, n_samples)
 
 
-def test_neighbor_average_waveform(camera_waveforms):
-    waveforms, camera = camera_waveforms
+def test_neighbor_average_waveform_2chan(camera_waveforms_2chan):
+    waveforms, camera = camera_waveforms_2chan
     nei = camera.neighbor_matrix_where
     average_wf = neighbor_average_waveform(waveforms, nei, 0)
 
@@ -115,7 +313,7 @@ def test_neighbor_average_waveform(camera_waveforms):
     assert_allclose(average_wf[1, 0, 48], 9.578896, rtol=1e-3)
 
 
-def test_extract_pulse_time_around_peak(camera_waveforms):
+def test_extract_pulse_time_around_peak_2chan(camera_waveforms_2chan):
     x = np.arange(100)
     y = norm.pdf(x, 41.2, 6)
     pulse_time = extract_pulse_time_around_peak(
@@ -125,8 +323,8 @@ def test_extract_pulse_time_around_peak(camera_waveforms):
     assert_allclose(pulse_time[0], 41.2, rtol=1e-3)
 
 
-def test_baseline_subtractor(camera_waveforms):
-    waveforms, _ = camera_waveforms
+def test_baseline_subtractor_2chan(camera_waveforms_2chan):
+    waveforms, _ = camera_waveforms_2chan
     n_chan, n_pixels, n_samples = waveforms.shape
     rand = np.random.RandomState(1)
     offset = np.arange(n_pixels)[np.newaxis, :, np.newaxis]
@@ -136,8 +334,8 @@ def test_baseline_subtractor(camera_waveforms):
     assert_allclose(baseline_subtracted.mean(), 0, atol=1e-3)
 
 
-def test_full_waveform_sum(camera_waveforms):
-    waveforms, _ = camera_waveforms
+def test_full_waveform_sum_2chan(camera_waveforms_2chan):
+    waveforms, _ = camera_waveforms_2chan
     extractor = FullWaveformSum()
     charge, pulse_time = extractor(waveforms)
 
@@ -147,8 +345,8 @@ def test_full_waveform_sum(camera_waveforms):
     assert_allclose(pulse_time[1][0], 62.359948, rtol=1e-3)
 
 
-def test_fixed_window_sum(camera_waveforms):
-    waveforms, _ = camera_waveforms
+def test_fixed_window_sum_2chan(camera_waveforms_2chan):
+    waveforms, _ = camera_waveforms_2chan
     extractor = FixedWindowSum(window_start=45)
     charge, pulse_time = extractor(waveforms)
 
@@ -158,8 +356,8 @@ def test_fixed_window_sum(camera_waveforms):
     assert_allclose(pulse_time[1][0], 49.370007, rtol=1e-3)
 
 
-def test_global_peak_window_sum(camera_waveforms):
-    waveforms, _ = camera_waveforms
+def test_global_peak_window_sum_2chan(camera_waveforms_2chan):
+    waveforms, _ = camera_waveforms_2chan
     extractor = GlobalPeakWindowSum()
     charge, pulse_time = extractor(waveforms)
 
@@ -169,8 +367,8 @@ def test_global_peak_window_sum(camera_waveforms):
     assert_allclose(pulse_time[1][0], 62.931829, rtol=1e-3)
 
 
-def test_local_peak_window_sum(camera_waveforms):
-    waveforms, _ = camera_waveforms
+def test_local_peak_window_sum_2chan(camera_waveforms_2chan):
+    waveforms, _ = camera_waveforms_2chan
     extractor = LocalPeakWindowSum()
     charge, pulse_time = extractor(waveforms)
 
@@ -180,8 +378,8 @@ def test_local_peak_window_sum(camera_waveforms):
     assert_allclose(pulse_time[1][0], 62.038344, rtol=1e-3)
 
 
-def test_neighbor_peak_window_sum(camera_waveforms):
-    waveforms, camera = camera_waveforms
+def test_neighbor_peak_window_sum_2chan(camera_waveforms_2chan):
+    waveforms, camera = camera_waveforms_2chan
     nei = camera.neighbor_matrix_where
     extractor = NeighborPeakWindowSum()
     extractor.neighbors = nei
@@ -201,8 +399,8 @@ def test_neighbor_peak_window_sum(camera_waveforms):
     assert_allclose(pulse_time[1][0], 62.038344, rtol=1e-3)
 
 
-def test_baseline_subtracted_neighbor_peak_window_sum(camera_waveforms):
-    waveforms, camera = camera_waveforms
+def test_baseline_subtracted_neighbor_pw_sum_2chan(camera_waveforms_2chan):
+    waveforms, camera = camera_waveforms_2chan
     nei = camera.neighbor_matrix_where
     extractor = BaselineSubtractedNeighborPeakWindowSum()
     extractor.neighbors = nei
@@ -212,37 +410,3 @@ def test_baseline_subtracted_neighbor_peak_window_sum(camera_waveforms):
     assert_allclose(charge[1][0], 426.887, rtol=1e-3)
     assert_allclose(pulse_time[0][0], 54.116092, rtol=1e-3)
     assert_allclose(pulse_time[1][0], 62.038344, rtol=1e-3)
-
-
-def test_waveform_extractor_factory(camera_waveforms):
-    waveforms, _ = camera_waveforms
-    extractor = ImageExtractor.from_name('LocalPeakWindowSum')
-    extractor(waveforms)
-
-
-def test_waveform_extractor_factory_args():
-    """
-    Config is supposed to be created by a `Tool`
-    """
-    from traitlets.config.loader import Config
-    config = Config(
-        {
-            'ImageExtractor': {
-                'window_width': 20,
-                'window_shift': 3,
-            }
-        }
-    )
-
-    extractor = ImageExtractor.from_name(
-        'LocalPeakWindowSum',
-        config=config,
-    )
-    assert extractor.window_width == 20
-    assert extractor.window_shift == 3
-
-    with pytest.warns(UserWarning):
-        ImageExtractor.from_name(
-            'FullWaveformSum',
-            config=config,
-        )
diff --git a/ctapipe/io/containers.py b/ctapipe/io/containers.py
index efd51ff7764..463cbe36383 100644
--- a/ctapipe/io/containers.py
+++ b/ctapipe/io/containers.py
@@ -94,19 +94,13 @@ class DL1CameraContainer(Container):
     image = Field(
         None,
         "Numpy array of camera image, after waveform extraction."
-        "Shape: (n_chan, n_pixel)"
+        "Shape: (n_pixel)"
     )
     pulse_time = Field(
         None,
         "Numpy array containing position of the pulse as determined by "
         "the extractor."
-        "Shape: (n_chan, n_pixel, n_samples)"
-    )
-    #TODO: Remove when gain selection added?
-    gain_channel = Field(
-        None,
-        "boolean numpy array of which gain channel was used for each pixel "
-        "in the image "
+        "Shape: (n_pixel, n_samples)"
     )
 
 class CameraCalibrationContainer(Container):
@@ -133,7 +127,7 @@ class R0CameraContainer(Container):
     trig_pix_id = Field(None, "pixels involved in the camera trigger")
     waveform = Field(None, (
         "numpy array containing ADC samples"
-        "(n_channels x n_pixels, n_samples)"
+        "(n_channels, n_pixels, n_samples)"
     ))
 
 
@@ -158,7 +152,11 @@ class R1CameraContainer(Container):
 
     waveform = Field(None, (
         "numpy array containing a set of images, one per ADC sample"
-        "(n_channels x n_pixels, n_samples)"
+        "Shape: (n_channels, n_pixels, n_samples)"
+    ))
+    selected_gain_channel = Field(None, (
+        "Numpy array containing the gain channel chosen for each pixel. "
+        "Shape: (n_pixels)"
     ))
 
 
diff --git a/ctapipe/plotting/bokeh_event_viewer.py b/ctapipe/plotting/bokeh_event_viewer.py
index 877bffaa4fa..7a6523fca0f 100644
--- a/ctapipe/plotting/bokeh_event_viewer.py
+++ b/ctapipe/plotting/bokeh_event_viewer.py
@@ -28,9 +28,9 @@ def __init__(self, event_viewer, fig=None):
         self._view_options = {
             'r0': lambda e, t, c, time: e.r0.tel[t].waveform[c, :, time],
             'r1': lambda e, t, c, time: e.r1.tel[t].waveform[c, :, time],
-            'dl0': lambda e, t, c, time: e.dl0.tel[t].waveform[c, :, time],
-            'dl1': lambda e, t, c, time: e.dl1.tel[t].image[c, :],
-            'pulse_time': lambda e, t, c, time: e.dl1.tel[t].pulse_time[c, :],
+            'dl0': lambda e, t, c, time: e.dl0.tel[t].waveform[:, time],
+            'dl1': lambda e, t, c, time: e.dl1.tel[t].image[:],
+            'pulse_time': lambda e, t, c, time: e.dl1.tel[t].pulse_time[:],
         }
 
         self.w_view = None
@@ -177,7 +177,7 @@ def __init__(self, event_viewer, fig=None):
         self._view_options = {
             'r0': lambda e, t, c, p: e.r0.tel[t].waveform[c, p],
             'r1': lambda e, t, c, p: e.r1.tel[t].waveform[c, p],
-            'dl0': lambda e, t, c, p: e.dl0.tel[t].waveform[c, p],
+            'dl0': lambda e, t, c, p: e.dl0.tel[t].waveform[p],
         }
 
         self.w_view = None
diff --git a/ctapipe/plotting/event_viewer.py b/ctapipe/plotting/event_viewer.py
index 7b32ff36484..048c18b2875 100644
--- a/ctapipe/plotting/event_viewer.py
+++ b/ctapipe/plotting/event_viewer.py
@@ -100,7 +100,7 @@ def draw_event(self, event, hillas_parameters=None):
             ax = plt.subplot(camera_grid[ii])
             geom = event.inst.subarray.tel[tel_id].camera
             self.get_camera_view(tel_id,
-                                 image=images.tel[tel_id].image[0],
+                                 image=images.tel[tel_id].image,
                                  axis=ax,
                                  geom=geom)
 
diff --git a/ctapipe/tools/display_dl1.py b/ctapipe/tools/display_dl1.py
index d9b1c9eb430..10c9525772d 100644
--- a/ctapipe/tools/display_dl1.py
+++ b/ctapipe/tools/display_dl1.py
@@ -68,9 +68,9 @@ def get_geometry(event, telid):
         return event.inst.subarray.tel[telid].camera
 
     def plot(self, event, telid):
-        chan = 0
-        image = event.dl1.tel[telid].image[chan]
-        pulse_time = event.dl1.tel[telid].pulse_time[chan]
+        image = event.dl1.tel[telid].image
+        pulse_time = event.dl1.tel[telid].pulse_time
+        print("plot", image.shape, pulse_time.shape)
 
         if self._current_tel != telid:
             self._current_tel = telid
@@ -83,7 +83,7 @@ def plot(self, event, telid):
             self.c_intensity = CameraDisplay(geom, ax=self.ax_intensity)
             self.c_pulse_time = CameraDisplay(geom, ax=self.ax_pulse_time)
 
-            tmaxmin = event.dl0.tel[telid].waveform.shape[2]
+            tmaxmin = event.dl0.tel[telid].waveform.shape[1]
             t_chargemax = pulse_time[image.argmax()]
             cmap_time = colors.LinearSegmentedColormap.from_list(
                 "cmap_t",
diff --git a/ctapipe/tools/display_events_single_tel.py b/ctapipe/tools/display_events_single_tel.py
index 7fde80cb196..417fd574e82 100755
--- a/ctapipe/tools/display_events_single_tel.py
+++ b/ctapipe/tools/display_events_single_tel.py
@@ -17,7 +17,7 @@
 from ctapipe.calib import CameraCalibrator
 from ctapipe.core import Tool
 from ctapipe.core.traits import Float, Dict, List
-from ctapipe.core.traits import Unicode, Int, Integer, Bool
+from ctapipe.core.traits import Unicode, Int, Bool
 from ctapipe.image import (
     tailcuts_clean, hillas_parameters, HillasParameterizationError
 )
@@ -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': 'SingleTelEventDisplay.infile',
         'tel': 'SingleTelEventDisplay.tel',
         'max-events': 'EventSource.max_events',
-        'channel': 'SingleTelEventDisplay.channel',
         'write': 'SingleTelEventDisplay.write',
         'clean': 'SingleTelEventDisplay.clean',
         'hillas': 'SingleTelEventDisplay.hillas',
@@ -118,7 +114,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)
@@ -132,7 +128,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/ctapipe/tools/display_integrator.py b/ctapipe/tools/display_integrator.py
index 74114d48cfd..d754c0b7d54 100644
--- a/ctapipe/tools/display_integrator.py
+++ b/ctapipe/tools/display_integrator.py
@@ -18,10 +18,10 @@
 
 def plot(event, telid, chan, extractor_name):
     # Extract required images
-    dl0 = event.dl0.tel[telid].waveform[chan]
+    dl0 = event.dl0.tel[telid].waveform
 
     t_pe = event.mc.tel[telid].photo_electron_image
-    dl1 = event.dl1.tel[telid].image[chan]
+    dl1 = event.dl1.tel[telid].image
     max_time = np.unravel_index(np.argmax(dl0), dl0.shape)[1]
     max_charges = np.max(dl0, axis=1)
     max_pix = int(np.argmax(max_charges))
diff --git a/ctapipe/tools/display_summed_images.py b/ctapipe/tools/display_summed_images.py
index b8016780083..926c69a1bb1 100644
--- a/ctapipe/tools/display_summed_images.py
+++ b/ctapipe/tools/display_summed_images.py
@@ -97,7 +97,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/ctapipe/tools/extract_charge_resolution.py b/ctapipe/tools/extract_charge_resolution.py
index 17b94026137..03685fc1c68 100644
--- a/ctapipe/tools/extract_charge_resolution.py
+++ b/ctapipe/tools/extract_charge_resolution.py
@@ -86,7 +86,7 @@ def start(self):
 
             for mc, dl1 in zip(event.mc.tel.values(), event.dl1.tel.values()):
                 true_charge = mc.photo_electron_image
-                measured_charge = dl1.image[0]
+                measured_charge = dl1.image
                 pixels = np.arange(measured_charge.size)
                 self.calculator.add(pixels, true_charge, measured_charge)
 
diff --git a/docs/tutorials/calibrated_data_exploration.ipynb b/docs/tutorials/calibrated_data_exploration.ipynb
index c4d7c72ef2a..43dca5e1439 100644
--- a/docs/tutorials/calibrated_data_exploration.ipynb
+++ b/docs/tutorials/calibrated_data_exploration.ipynb
@@ -160,7 +160,7 @@
     "tel_id = sorted(event.r1.tels_with_data)[1]\n",
     "sub = event.inst.subarray\n",
     "camera = sub.tel[tel_id].camera\n",
-    "image = event.dl1.tel[tel_id].image[0] "
+    "image = event.dl1.tel[tel_id].image"
    ]
   },
   {
@@ -293,7 +293,7 @@
     "image_sum = np.zeros_like(tel.camera.pix_x.value)  # just make an array of 0's in the same shape as the camera \n",
     "\n",
     "for tel_id in cams_in_event:\n",
-    "    image_sum += event.dl1.tel[tel_id].image[0]"
+    "    image_sum += event.dl1.tel[tel_id].image"
    ]
   },
   {
@@ -382,7 +382,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.2"
   }
  },
  "nbformat": 4,
diff --git a/docs/tutorials/ctapipe_handson.ipynb b/docs/tutorials/ctapipe_handson.ipynb
index becbd1efbbe..ae04b5bf960 100644
--- a/docs/tutorials/ctapipe_handson.ipynb
+++ b/docs/tutorials/ctapipe_handson.ipynb
@@ -398,7 +398,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "CameraDisplay(tel.camera, image=dl1tel.image[0])"
+    "CameraDisplay(tel.camera, image=dl1tel.image)"
    ]
   },
   {
@@ -407,7 +407,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "CameraDisplay(tel.camera, image=dl1tel.pulse_time[0])"
+    "CameraDisplay(tel.camera, image=dl1tel.pulse_time)"
    ]
   },
   {
@@ -432,7 +432,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "image = dl1tel.image[0]\n",
+    "image = dl1tel.image\n",
     "mask = tailcuts_clean(tel.camera, image, picture_thresh=10, boundary_thresh=5)\n",
     "mask"
    ]
@@ -557,8 +557,8 @@
     "    \n",
     "    for tel_id, tel_data in event.dl1.tel.items():\n",
     "        tel = event.inst.subarray.tel[tel_id]\n",
-    "        mask = tailcuts_clean(tel.camera, tel_data.image[0])\n",
-    "        params = hillas_parameters(tel.camera[mask], tel_data.image[0][mask])"
+    "        mask = tailcuts_clean(tel.camera, tel_data.image)\n",
+    "        params = hillas_parameters(tel.camera[mask], tel_data.image[mask])"
    ]
   },
   {
@@ -583,8 +583,8 @@
     "    \n",
     "        for tel_id, tel_data in event.dl1.tel.items():\n",
     "            tel = event.inst.subarray.tel[tel_id]\n",
-    "            mask = tailcuts_clean(tel.camera, tel_data.image[0])\n",
-    "            params = hillas_parameters(tel.camera[mask], tel_data.image[0][mask])\n",
+    "            mask = tailcuts_clean(tel.camera, tel_data.image)\n",
+    "            params = hillas_parameters(tel.camera[mask], tel_data.image[mask])\n",
     "            writer.write(\"hillas\", params)"
    ]
   },
@@ -649,7 +649,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.2"
   }
  },
  "nbformat": 4,
diff --git a/docs/tutorials/lst_analysis_bootcamp_2018.ipynb b/docs/tutorials/lst_analysis_bootcamp_2018.ipynb
index 4d5ba92feb4..a971ccc2cf3 100644
--- a/docs/tutorials/lst_analysis_bootcamp_2018.ipynb
+++ b/docs/tutorials/lst_analysis_bootcamp_2018.ipynb
@@ -356,7 +356,7 @@
     "\n",
     "# right now, there might be one image per gain channel.\n",
     "# This will change as soon as \n",
-    "display.image = dl1.image[0]\n",
+    "display.image = dl1.image\n",
     "display.add_colorbar()"
    ]
   },
@@ -401,7 +401,7 @@
     "\n",
     "clean = tailcuts_clean(\n",
     "    camera, \n",
-    "    dl1.image[0],\n",
+    "    dl1.image,\n",
     "    boundary_thresh=boundary,\n",
     "    picture_thresh=picture,\n",
     "    min_number_picture_neighbors=min_neighbors\n",
@@ -420,11 +420,11 @@
     "d2 = CameraDisplay(camera, ax=ax2)\n",
     "\n",
     "ax1.set_title('Image')\n",
-    "d1.image = dl1.image[0]\n",
+    "d1.image = dl1.image\n",
     "d1.add_colorbar(ax=ax1)\n",
     "\n",
     "ax2.set_title('Pulse Time')\n",
-    "d2.image = dl1.pulse_time[0] - np.average(dl1.pulse_time[0], weights=dl1.image[0])\n",
+    "d2.image = dl1.pulse_time - np.average(dl1.pulse_time, weights=dl1.image)\n",
     "d2.cmap = 'RdBu_r'\n",
     "d2.add_colorbar(ax=ax2)\n",
     "d2.set_limits_minmax(-20,20)\n",
@@ -457,7 +457,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hillas = hillas_parameters(camera[clean], dl1.image[0][clean])\n",
+    "hillas = hillas_parameters(camera[clean], dl1.image[clean])\n",
     "\n",
     "print(hillas)"
    ]
@@ -471,7 +471,7 @@
     "display = CameraDisplay(camera)\n",
     "\n",
     "# set \"unclean\" pixels to 0\n",
-    "cleaned = dl1.image[0].copy()\n",
+    "cleaned = dl1.image.copy()\n",
     "cleaned[~clean] = 0.0\n",
     "\n",
     "display.image = cleaned\n",
@@ -488,8 +488,8 @@
    "source": [
     "timing = timing_parameters(\n",
     "    camera[clean],\n",
-    "    dl1.image[0][clean],\n",
-    "    dl1.pulse_time[0][clean],\n",
+    "    dl1.image[clean],\n",
+    "    dl1.pulse_time[clean],\n",
     "    hillas,\n",
     ")\n",
     "\n",
@@ -506,7 +506,7 @@
     "    camera.pix_x, camera.pix_y,hillas.x, hillas.y, hillas.psi\n",
     ")\n",
     "\n",
-    "plt.plot(long[clean], dl1.pulse_time[0][clean], 'o')\n",
+    "plt.plot(long[clean], dl1.pulse_time[clean], 'o')\n",
     "plt.plot(long[clean], timing.slope * long[clean] + timing.intercept)"
    ]
   },
@@ -516,7 +516,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "l = leakage(camera, dl1.image[0], clean)\n",
+    "l = leakage(camera, dl1.image, clean)\n",
     "print(l)"
    ]
   },
@@ -527,7 +527,7 @@
    "outputs": [],
    "source": [
     "disp = CameraDisplay(camera)\n",
-    "disp.image = dl1.image[0]\n",
+    "disp.image = dl1.image\n",
     "disp.highlight_pixels(camera.get_border_pixel_mask(1), linewidth=2, color='xkcd:red')"
    ]
   },
@@ -548,7 +548,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "conc = concentration(camera, dl1.image[0], hillas)\n",
+    "conc = concentration(camera, dl1.image, hillas)\n",
     "print(conc)"
    ]
   },
@@ -616,8 +616,8 @@
     "\n",
     "        for telescope_id, dl1 in event.dl1.tel.items():\n",
     "            camera = event.inst.subarray.tels[telescope_id].camera\n",
-    "            image = dl1.image[0]\n",
-    "            pulse_time = dl1.pulse_time[0]\n",
+    "            image = dl1.image\n",
+    "            pulse_time = dl1.pulse_time\n",
     "\n",
     "            boundary, picture, min_neighbors = cleaning_level[camera.cam_id]\n",
     "\n",
@@ -790,8 +790,8 @@
     "        for telescope_id, dl1 in event.dl1.tel.items():      \n",
     "\n",
     "            camera = event.inst.subarray.tels[telescope_id].camera\n",
-    "            image = dl1.image[0]\n",
-    "            pulse_time = dl1.pulse_time[0]\n",
+    "            image = dl1.image\n",
+    "            pulse_time = dl1.pulse_time\n",
     "\n",
     "            boundary, picture, min_neighbors = cleaning_level[camera.cam_id]\n",
     "\n",
@@ -1095,7 +1095,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.2"
   }
  },
  "nbformat": 4,
diff --git a/docs/tutorials/theta_square.ipynb b/docs/tutorials/theta_square.ipynb
index c42dd31cb52..a5f19744307 100644
--- a/docs/tutorials/theta_square.ipynb
+++ b/docs/tutorials/theta_square.ipynb
@@ -116,7 +116,7 @@
     "        camgeom = subarray.tel[tel_id].camera\n",
     "\n",
     "        # note the [0] is for channel 0 which is high-gain channel\n",
-    "        image = event.dl1.tel[tel_id].image[0]\n",
+    "        image = event.dl1.tel[tel_id].image\n",
     "\n",
     "        # Cleaning  of the image\n",
     "        cleaned_image = image\n",
@@ -235,7 +235,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.2"
   },
   "toc": {
    "nav_menu": {
diff --git a/examples/plot_array_hillas.py b/examples/plot_array_hillas.py
index f90afd665b7..f1df9fec12e 100644
--- a/examples/plot_array_hillas.py
+++ b/examples/plot_array_hillas.py
@@ -94,7 +94,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.copy()
diff --git a/examples/plot_showers_in_nominal.py b/examples/plot_showers_in_nominal.py
index 2be3f112e57..32e2cf46a88 100644
--- a/examples/plot_showers_in_nominal.py
+++ b/examples/plot_showers_in_nominal.py
@@ -38,7 +38,7 @@
         for tel_id, dl1 in event.dl1.tel.items():
             camera = event.inst.subarray.tels[tel_id].camera
             focal_length = event.inst.subarray.tels[tel_id].optics.equivalent_focal_length
-            image = dl1.image[0]
+            image = dl1.image
 
             # telescope mc info
             mc_tel = event.mc.tel[tel_id]
diff --git a/examples/plot_theta_square.py b/examples/plot_theta_square.py
index 051814d64b8..b0790ad7c63 100644
--- a/examples/plot_theta_square.py
+++ b/examples/plot_theta_square.py
@@ -53,7 +53,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_event_writer.py b/examples/simple_event_writer.py
index 3d0f939d7c1..7fbff5470ab 100755
--- a/examples/simple_event_writer.py
+++ b/examples/simple_event_writer.py
@@ -91,7 +91,7 @@ def start(self):
                 dl1_tel = event.dl1.tel[tel_id]
 
                 # Image cleaning
-                image = dl1_tel.image[0]  # Waiting for automatic gain selection
+                image = dl1_tel.image  # Waiting for automatic gain selection
                 mask = tailcuts_clean(camera, image, picture_thresh=10, boundary_thresh=5)
                 cleaned = image.copy()
                 cleaned[~mask] = 0
diff --git a/examples/stereo_reconstruction.py b/examples/stereo_reconstruction.py
index 7aa78cbeb3e..bcb953c03ea 100644
--- a/examples/stereo_reconstruction.py
+++ b/examples/stereo_reconstruction.py
@@ -41,8 +41,8 @@
     for telescope_id, dl1 in event.dl1.tel.items():
         camera = event.inst.subarray.tels[telescope_id].camera
 
-        image = dl1.image[0]
-        peakpos = dl1.peakpos[0]
+        image = dl1.image
+        peakpos = dl1.pulse_time
 
         # cleaning
         boundary, picture, min_neighbors = cleaning_level[camera.cam_id]