From 76b7f620065b08de995c392193acdc17f019ed93 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20N=C3=B6the?= <maximilian.noethe@tu-dortmund.de>
Date: Mon, 23 May 2022 10:39:42 +0200
Subject: [PATCH 01/35] Add basic structure for flashcam extractor

---
 ctapipe/image/extractor.py            | 16 ++++++++++++++++
 ctapipe/image/tests/test_extractor.py | 14 +++++++++++++-
 ctapipe/resources/base_config.yaml    |  7 ++++++-
 3 files changed, 35 insertions(+), 2 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 0e078582ef8..32ed15d67c8 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1285,3 +1285,19 @@ def __call__(
             peak_time=pulse_time2.astype("float32"),
             is_valid=is_valid,
         )
+
+
+@lru_cache
+def exponential_scale_from_pulse_shape(camera: CameraDescription):
+    pulse_shape = camera.readout.reference_pulse_shape
+    sample_width = camera.readout.reference_pulse_sample_width
+
+
+class FlashCamExtractor(ImageExtractor):
+    def __call__(self, waveforms, telid, selected_gain_channel) -> DL1CameraContainer:
+
+        return DL1CameraContainer(
+            image=np.zeros(...),
+            peak_time=...,
+            is_valid=True,
+        )
diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 9246dddfda1..2a6cd7e8240 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -23,7 +23,12 @@
     subtract_baseline,
 )
 from ctapipe.image.toymodel import SkewedGaussian, WaveformModel, obtain_time_image
-from ctapipe.instrument import SubarrayDescription
+from ctapipe.instrument import SubarrayDescription, TelescopeDescription
+from numpy.testing import assert_allclose, assert_equal
+from scipy.stats import norm
+from traitlets.config.loader import Config
+from traitlets.traitlets import TraitError
+from ctapipe.io import EventSource
 
 extractors = non_abstract_children(ImageExtractor)
 # FixedWindowSum has no peak finding and need to be set manually
@@ -617,3 +622,10 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray):
 
     expected = np.average([29, 30, 31], weights=[5, 10, 3])
     assert np.allclose(dl1.peak_time[bright_pixels], expected / sample_rate)
+
+
+def test_flashcam_extractor(tmp_path):
+    path = "gamma_20deg_0deg_run2___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10-100evts.simtel.zst"
+    source = EventSource(f"dataset://{path}")
+    subarray = source.subarray
+    print(subarray.get_tel_ids_for_type("MST_MST_FlashCam"))
diff --git a/ctapipe/resources/base_config.yaml b/ctapipe/resources/base_config.yaml
index 539bfb0fa84..48ce249d7b1 100644
--- a/ctapipe/resources/base_config.yaml
+++ b/ctapipe/resources/base_config.yaml
@@ -39,7 +39,12 @@ CameraCalibrator:
   #
   # Note this is a telescope-wise parameter, so can be specified per telescope
   # if necessary (see below for an example)
-  image_extractor_type: NeighborPeakWindowSum
+  image_extractor_type:
+    - ['type', '*', 'NeighborPeakWindowSum']
+    - ['type', '*FlashCam', 'FlashCamExtractor']
+
+  FlashCamExtractor:
+    exponential_scale: 7.5
 
 # The ImageProcessor performs the DL1a-> DL1b (image parameters) transition. It
 # is run only if the parameters `DataWriter.write_image_parameters=True` and the

From b685a002dceb0da5c97ff5288b04791c9ff6354d Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Thu, 22 Sep 2022 10:23:58 +0200
Subject: [PATCH 02/35] Add draft implementation of FlashCamExtractor and
 utility functions. Charge test still fails because of neighbour jitter.

---
 ctapipe/image/extractor.py            | 131 ++++++++++++++++++++++++--
 ctapipe/image/tests/test_extractor.py |  40 +++++++-
 2 files changed, 160 insertions(+), 11 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 32ed15d67c8..b627aff6d4d 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -27,6 +27,7 @@
 import numpy as np
 from numba import float32, float64, guvectorize, int64, njit, prange
 from scipy.ndimage import convolve1d
+from scipy.signal import filtfilt
 from traitlets import Bool, Int
 
 from ctapipe.containers import DL1CameraContainer
@@ -198,7 +199,8 @@ def extract_sliding_window(waveforms, width, sampling_rate_ghz, sum_, peak_time)
     peak_time[0] /= sampling_rate_ghz
 
 
-@njit(cache=True)
+# FIXME JIT has been disabled because it asserts when passing floating-point waveforms (e.g., FlashCamExtractor, at least on arm64)
+# @njit(cache=True)
 def neighbor_average_maximum(
     waveforms, neighbors_indices, neighbors_indptr, local_weight, broken_pixels
 ):
@@ -1286,18 +1288,127 @@ def __call__(
             is_valid=is_valid,
         )
 
-
 @lru_cache
-def exponential_scale_from_pulse_shape(camera: CameraDescription):
-    pulse_shape = camera.readout.reference_pulse_shape
-    sample_width = camera.readout.reference_pulse_sample_width
+def deconvolution_parameters(camera: CameraDescription, upsampling: int, window_width: int, window_shift: int):
+    assert(upsampling > 0)
+    assert(window_width > 0)
+
+    ref_pulse_shapes = camera.readout.reference_pulse_shape
+    ref_sample_width = camera.readout.reference_pulse_sample_width.to_value("ns")
+    camera_sample_width = 1.0 / camera.readout.sampling_rate.to_value("GHz")
+
+    assert(camera_sample_width >= ref_sample_width)
+    avg_step = int(camera_sample_width / ref_sample_width + 0.5)
+
+    pzs = []
+    for ref_pulse_shape in ref_pulse_shapes:
+        phase_pzs = []
+        for phase in range(avg_step):
+            x = ref_pulse_shape[phase::avg_step]
+            i_min = np.argmin(np.diff(x)) + 1
+            if i_min < x.size and x[i_min] != 0:
+                phase_pzs.append(x[i_min + 1] / x[i_min])
+
+        assert(len(phase_pzs))
+        pzs.append(np.mean(phase_pzs))
+
+    gains, shifts = [], []
+    for pz, ref_pulse_shape in zip(pzs, ref_pulse_shapes):
+        integral = ref_pulse_shape.sum() * ref_sample_width / camera_sample_width
+        phase_gains, phase_shifts = [], []
+        for phase in range(avg_step):
+            x = ref_pulse_shape[phase::avg_step]
+            y = deconvolve(np.atleast_2d(x), 0.0, upsampling, pz)[0]
+            i_max = y.argmax()
+            start = i_max - window_shift  # TODO should use extract_around_peak here
+            stop = start + window_width
+            if start >= 0 and stop <= y.size:
+                phase_shifts.append((i_max / upsampling * avg_step - ref_pulse_shape.argmax()) * ref_sample_width)
+                phase_gains.append(y[start:stop].sum() / integral)
+
+        assert(len(phase_gains))
+        gains.append(np.mean(phase_gains))
+        shifts.append(np.mean(phase_shifts))
+
+    return pzs, gains, shifts
+
+
+def deconvolve(waveforms, bls, up : int, pz : float):
+    assert(len(waveforms.shape) == 2)
+    bls = np.atleast_2d(bls).T
+    y = waveforms - bls
+    y[:,1:] -= pz * y[:,:-1]
+    y[:,0] = 0
+    if up > 1:
+        return filtfilt(np.ones(up), up, np.repeat(y, up, axis=-1))
+    return y
 
 
 class FlashCamExtractor(ImageExtractor):
-    def __call__(self, waveforms, telid, selected_gain_channel) -> DL1CameraContainer:
+    upsampling = IntTelescopeParameter(
+        default_value=4, help="Define the upsampling factor for waveforms"
+    ).tag(config=True)
 
-        return DL1CameraContainer(
-            image=np.zeros(...),
-            peak_time=...,
-            is_valid=True,
+    window_width = IntTelescopeParameter(
+        default_value=4, help="Define the width of the integration window"
+    ).tag(config=True)
+
+    window_shift = IntTelescopeParameter(
+        default_value=2,
+        help="Define the shift of the integration window from the peak_index "
+        "(peak_index - shift)",
+    ).tag(config=True)
+
+    apply_integration_correction = BoolTelescopeParameter(
+        default_value=True, help="Apply the integration window correction"
+    ).tag(config=True)
+
+    local_weight = IntTelescopeParameter(
+        default_value=0,
+        help="Weight of the local pixel (0: peak from neighbors only, "
+        "1: local pixel counts as much as any neighbor)",
+    ).tag(config=True)
+
+    def __init__(self, subarray, **kwargs):
+        super().__init__(subarray=subarray, **kwargs)
+
+        self.sampling_rate_ghz = {
+            tel_id: telescope.camera.readout.sampling_rate.to_value("GHz")
+            for tel_id, telescope in subarray.tel.items()
+        }
+
+    def __call__(self, waveforms, tel_id, _, broken_pixels) -> DL1CameraContainer:
+        upsampling = max(1, self.upsampling.tel[tel_id])
+        window_width = max(1, self.window_width.tel[tel_id])
+        window_shift = self.window_shift.tel[tel_id]
+        pzs, gains, shifts = deconvolution_parameters(self.subarray.tel[tel_id].camera, upsampling, window_width, window_shift)
+        pz, gain, shift = pzs[0], gains[0], shifts[0]
+
+        bls = waveforms[:,0]  # FIXME fetch from NSB estimates
+        waveforms = deconvolve(waveforms, bls, upsampling, pz)
+
+        # FIXME near-duplicate of neighbour peak sum for now
+        neighbors = self.subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse
+        peak_index = neighbor_average_maximum(
+            waveforms,
+            neighbors_indices=neighbors.indices,
+            neighbors_indptr=neighbors.indptr,
+            local_weight=self.local_weight.tel[tel_id],
+            broken_pixels=broken_pixels,
+        )
+
+        charge, peak_time = extract_around_peak(
+            waveforms,
+            peak_index,
+            window_width,
+            window_shift,
+            self.sampling_rate_ghz[tel_id] * upsampling,
         )
+
+        if gain != 0:
+            charge /= gain
+
+        if shift != 0:
+            peak_time -= shift
+
+        return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True)
diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 2a6cd7e8240..d55240d6067 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -11,6 +11,7 @@
 from ctapipe.core import non_abstract_children
 from ctapipe.image.extractor import (
     FixedWindowSum,
+    FlashCamExtractor,
     FullWaveformSum,
     ImageExtractor,
     NeighborPeakWindowSum,
@@ -33,6 +34,7 @@
 extractors = non_abstract_children(ImageExtractor)
 # FixedWindowSum has no peak finding and need to be set manually
 extractors.remove(FixedWindowSum)
+extractors.remove(FlashCamExtractor)
 
 
 @pytest.fixture(scope="module")
@@ -73,6 +75,16 @@ def subarray_1_LST(prod3_lst):
     return subarray
 
 
+@pytest.fixture(scope="module")
+def subarray_1_MST_FC(prod5_mst_flashcam):
+    subarray = SubarrayDescription(
+        "One MST with FlashCam",
+        tel_positions={1: np.zeros(3) * u.m},
+        tel_descriptions={1: prod5_mst_flashcam},
+    )
+    return subarray
+
+
 def get_test_toymodel(subarray, minCharge=100, maxCharge=1000):
     tel_id = list(subarray.tel.keys())[0]
     n_pixels = subarray.tel[tel_id].camera.geometry.n_pixels
@@ -97,6 +109,11 @@ def toymodel(subarray):
     return get_test_toymodel(subarray)
 
 
+@pytest.fixture(scope="module")
+def toymodel_1_MST_FC(subarray_1_MST_FC):
+    return get_test_toymodel(subarray_1_MST_FC)
+
+
 def test_extract_around_peak(toymodel):
     waveforms, _, _, _, _, _ = toymodel
     n_pixels, n_samples = waveforms.shape
@@ -624,8 +641,29 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray):
     assert np.allclose(dl1.peak_time[bright_pixels], expected / sample_rate)
 
 
-def test_flashcam_extractor(tmp_path):
+def test_flashcam_extractor(toymodel_1_MST_FC):
+    (
+        waveforms,
+        subarray,
+        tel_id,
+        selected_gain_channel,
+        true_charge,
+        true_time,
+    ) = toymodel_1_MST_FC
+
+    assert(subarray.tel[1].camera.name == "FlashCam")
+    extractor = FlashCamExtractor(subarray=subarray)
+
+    # Test toymodel
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+    dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
+    assert_allclose(dl1.image, true_charge, rtol=0.1)
+    assert_allclose(dl1.peak_time, true_time, atol=2.5)  # FIXME could be much closer for large pulses...
+    assert dl1.is_valid == True
+
+    # Test prod5 simulations
     path = "gamma_20deg_0deg_run2___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10-100evts.simtel.zst"
     source = EventSource(f"dataset://{path}")
     subarray = source.subarray
     print(subarray.get_tel_ids_for_type("MST_MST_FlashCam"))
+    # TODO Test!

From ac70558e25838a386c971af8677b5885f721befd Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Thu, 22 Sep 2022 11:31:06 +0200
Subject: [PATCH 03/35] Add FlashCamExtractor to module export. Fix naming of
 selected_gain_channel arg.

---
 ctapipe/image/extractor.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index b627aff6d4d..530f66dcfee 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -4,6 +4,7 @@
 
 __all__ = [
     "ImageExtractor",
+    "FlashCamExtractor",
     "FullWaveformSum",
     "FixedWindowSum",
     "GlobalPeakWindowSum",
@@ -1377,7 +1378,7 @@ def __init__(self, subarray, **kwargs):
             for tel_id, telescope in subarray.tel.items()
         }
 
-    def __call__(self, waveforms, tel_id, _, broken_pixels) -> DL1CameraContainer:
+    def __call__(self, waveforms, tel_id, selected_gain_channel, broken_pixels) -> DL1CameraContainer:
         upsampling = max(1, self.upsampling.tel[tel_id])
         window_width = max(1, self.window_width.tel[tel_id])
         window_shift = self.window_shift.tel[tel_id]

From 8febc8ed4eb31594b667b093a91d03dfdc6f86b3 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Thu, 22 Sep 2022 11:37:27 +0200
Subject: [PATCH 04/35] Format according to pre-commit hooks.

---
 ctapipe/image/extractor.py | 41 ++++++++++++++++++++++++--------------
 1 file changed, 26 insertions(+), 15 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 530f66dcfee..eb0abefd0f7 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -26,7 +26,7 @@
 from typing import Tuple
 
 import numpy as np
-from numba import float32, float64, guvectorize, int64, njit, prange
+from numba import float32, float64, guvectorize, int64, prange
 from scipy.ndimage import convolve1d
 from scipy.signal import filtfilt
 from traitlets import Bool, Int
@@ -39,6 +39,7 @@
     FloatTelescopeParameter,
     IntTelescopeParameter,
 )
+from ctapipe.instrument import CameraDescription
 
 from .cleaning import tailcuts_clean
 from .hillas import camera_to_shower_coordinates, hillas_parameters
@@ -1289,16 +1290,19 @@ def __call__(
             is_valid=is_valid,
         )
 
+
 @lru_cache
-def deconvolution_parameters(camera: CameraDescription, upsampling: int, window_width: int, window_shift: int):
-    assert(upsampling > 0)
-    assert(window_width > 0)
+def deconvolution_parameters(
+    camera: CameraDescription, upsampling: int, window_width: int, window_shift: int
+):
+    assert upsampling > 0
+    assert window_width > 0
 
     ref_pulse_shapes = camera.readout.reference_pulse_shape
     ref_sample_width = camera.readout.reference_pulse_sample_width.to_value("ns")
     camera_sample_width = 1.0 / camera.readout.sampling_rate.to_value("GHz")
 
-    assert(camera_sample_width >= ref_sample_width)
+    assert camera_sample_width >= ref_sample_width
     avg_step = int(camera_sample_width / ref_sample_width + 0.5)
 
     pzs = []
@@ -1310,7 +1314,7 @@ def deconvolution_parameters(camera: CameraDescription, upsampling: int, window_
             if i_min < x.size and x[i_min] != 0:
                 phase_pzs.append(x[i_min + 1] / x[i_min])
 
-        assert(len(phase_pzs))
+        assert len(phase_pzs)
         pzs.append(np.mean(phase_pzs))
 
     gains, shifts = [], []
@@ -1324,22 +1328,25 @@ def deconvolution_parameters(camera: CameraDescription, upsampling: int, window_
             start = i_max - window_shift  # TODO should use extract_around_peak here
             stop = start + window_width
             if start >= 0 and stop <= y.size:
-                phase_shifts.append((i_max / upsampling * avg_step - ref_pulse_shape.argmax()) * ref_sample_width)
+                phase_shifts.append(
+                    (i_max / upsampling * avg_step - ref_pulse_shape.argmax())
+                    * ref_sample_width
+                )
                 phase_gains.append(y[start:stop].sum() / integral)
 
-        assert(len(phase_gains))
+        assert len(phase_gains)
         gains.append(np.mean(phase_gains))
         shifts.append(np.mean(phase_shifts))
 
     return pzs, gains, shifts
 
 
-def deconvolve(waveforms, bls, up : int, pz : float):
-    assert(len(waveforms.shape) == 2)
+def deconvolve(waveforms, bls, up: int, pz: float):
+    assert len(waveforms.shape) == 2
     bls = np.atleast_2d(bls).T
     y = waveforms - bls
-    y[:,1:] -= pz * y[:,:-1]
-    y[:,0] = 0
+    y[:, 1:] -= pz * y[:, :-1]
+    y[:, 0] = 0
     if up > 1:
         return filtfilt(np.ones(up), up, np.repeat(y, up, axis=-1))
     return y
@@ -1378,14 +1385,18 @@ def __init__(self, subarray, **kwargs):
             for tel_id, telescope in subarray.tel.items()
         }
 
-    def __call__(self, waveforms, tel_id, selected_gain_channel, broken_pixels) -> DL1CameraContainer:
+    def __call__(
+        self, waveforms, tel_id, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
         upsampling = max(1, self.upsampling.tel[tel_id])
         window_width = max(1, self.window_width.tel[tel_id])
         window_shift = self.window_shift.tel[tel_id]
-        pzs, gains, shifts = deconvolution_parameters(self.subarray.tel[tel_id].camera, upsampling, window_width, window_shift)
+        pzs, gains, shifts = deconvolution_parameters(
+            self.subarray.tel[tel_id].camera, upsampling, window_width, window_shift
+        )
         pz, gain, shift = pzs[0], gains[0], shifts[0]
 
-        bls = waveforms[:,0]  # FIXME fetch from NSB estimates
+        bls = waveforms[:, 0]  # FIXME fetch from NSB estimates
         waveforms = deconvolve(waveforms, bls, upsampling, pz)
 
         # FIXME near-duplicate of neighbour peak sum for now

From 92a12b59ead30764ed54e5c7d3ecc9c724fa80d8 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Thu, 22 Sep 2022 11:49:48 +0200
Subject: [PATCH 05/35] Reactivate njit.

---
 ctapipe/image/extractor.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index eb0abefd0f7..a15de0f8c86 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -26,7 +26,7 @@
 from typing import Tuple
 
 import numpy as np
-from numba import float32, float64, guvectorize, int64, prange
+from numba import float32, float64, guvectorize, int64, njit, prange
 from scipy.ndimage import convolve1d
 from scipy.signal import filtfilt
 from traitlets import Bool, Int
@@ -201,8 +201,7 @@ def extract_sliding_window(waveforms, width, sampling_rate_ghz, sum_, peak_time)
     peak_time[0] /= sampling_rate_ghz
 
 
-# FIXME JIT has been disabled because it asserts when passing floating-point waveforms (e.g., FlashCamExtractor, at least on arm64)
-# @njit(cache=True)
+@njit(cache=True)
 def neighbor_average_maximum(
     waveforms, neighbors_indices, neighbors_indptr, local_weight, broken_pixels
 ):

From 1193167bcbadee1d54d9efc4e729271cabddaddb Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Thu, 22 Sep 2022 12:06:08 +0200
Subject: [PATCH 06/35] Throw exceptions instead of asserting.

---
 ctapipe/image/extractor.py | 23 +++++++++++++++++------
 1 file changed, 17 insertions(+), 6 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index a15de0f8c86..423cf451774 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1294,14 +1294,19 @@ def __call__(
 def deconvolution_parameters(
     camera: CameraDescription, upsampling: int, window_width: int, window_shift: int
 ):
-    assert upsampling > 0
-    assert window_width > 0
+    if upsampling < 1:
+        raise ValueError(f"upsampling must be > 0, got {upsampling}")
+    if window_width < 1:
+        raise ValueError(f"window_width must be > 0, got {window_width}")
 
     ref_pulse_shapes = camera.readout.reference_pulse_shape
     ref_sample_width = camera.readout.reference_pulse_sample_width.to_value("ns")
     camera_sample_width = 1.0 / camera.readout.sampling_rate.to_value("GHz")
 
-    assert camera_sample_width >= ref_sample_width
+    if camera_sample_width < ref_sample_width:
+        raise ValueError(
+            f"ref_sample_width must be equal to or shorter than camera_sample_width, got ref_sample_width={ref_sample_width} and camera_sample_width={camera_sample_width}"
+        )
     avg_step = int(camera_sample_width / ref_sample_width + 0.5)
 
     pzs = []
@@ -1313,7 +1318,10 @@ def deconvolution_parameters(
             if i_min < x.size and x[i_min] != 0:
                 phase_pzs.append(x[i_min + 1] / x[i_min])
 
-        assert len(phase_pzs)
+        if len(phase_pzs) == 0:
+            raise ValueError(
+                "ref_pulse_shape is malformed - cannot find deconvolution scale"
+            )
         pzs.append(np.mean(phase_pzs))
 
     gains, shifts = [], []
@@ -1333,7 +1341,10 @@ def deconvolution_parameters(
                 )
                 phase_gains.append(y[start:stop].sum() / integral)
 
-        assert len(phase_gains)
+        if len(phase_gains) == 0:
+            raise ValueError(
+                "ref_pulse_shape is malformed - peak is not well contained"
+            )
         gains.append(np.mean(phase_gains))
         shifts.append(np.mean(phase_shifts))
 
@@ -1341,7 +1352,7 @@ def deconvolution_parameters(
 
 
 def deconvolve(waveforms, bls, up: int, pz: float):
-    assert len(waveforms.shape) == 2
+    waveforms = np.atleast_2d(waveforms)
     bls = np.atleast_2d(bls).T
     y = waveforms - bls
     y[:, 1:] -= pz * y[:, :-1]

From 85f5a26c8812b6227e5f5e7205084e28f8b22ab5 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Thu, 22 Sep 2022 12:21:39 +0200
Subject: [PATCH 07/35] Move parameter estimation into constructor.

---
 ctapipe/image/extractor.py | 19 ++++++++++++++-----
 1 file changed, 14 insertions(+), 5 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 423cf451774..2d9b3e889a7 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1395,15 +1395,24 @@ def __init__(self, subarray, **kwargs):
             for tel_id, telescope in subarray.tel.items()
         }
 
+        self.deconvolution_pars = {
+            tel_id: deconvolution_parameters(
+                tel.camera,
+                self.upsampling.tel[tel_id],
+                self.window_width.tel[tel_id],
+                self.window_shift.tel[tel_id],
+            )
+            for tel_id, tel in subarray.tel.items()
+        }
+
     def __call__(
         self, waveforms, tel_id, selected_gain_channel, broken_pixels
     ) -> DL1CameraContainer:
-        upsampling = max(1, self.upsampling.tel[tel_id])
-        window_width = max(1, self.window_width.tel[tel_id])
+        upsampling = self.upsampling.tel[tel_id]
+        window_width = self.window_width.tel[tel_id]
         window_shift = self.window_shift.tel[tel_id]
-        pzs, gains, shifts = deconvolution_parameters(
-            self.subarray.tel[tel_id].camera, upsampling, window_width, window_shift
-        )
+
+        pzs, gains, shifts = self.deconvolution_pars[tel_id]
         pz, gain, shift = pzs[0], gains[0], shifts[0]
 
         bls = waveforms[:, 0]  # FIXME fetch from NSB estimates

From dbd006ba86d4e8b515b2cc66d678ea2cc2a84cd8 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Thu, 22 Sep 2022 12:30:07 +0200
Subject: [PATCH 08/35] Enforce positivity of upsampling and window_width.

---
 ctapipe/image/extractor.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 2d9b3e889a7..1c0ccb0c660 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1364,11 +1364,11 @@ def deconvolve(waveforms, bls, up: int, pz: float):
 
 class FlashCamExtractor(ImageExtractor):
     upsampling = IntTelescopeParameter(
-        default_value=4, help="Define the upsampling factor for waveforms"
+        default_value=4, min=1, help="Define the upsampling factor for waveforms"
     ).tag(config=True)
 
     window_width = IntTelescopeParameter(
-        default_value=4, help="Define the width of the integration window"
+        default_value=4, min=1, help="Define the width of the integration window"
     ).tag(config=True)
 
     window_shift = IntTelescopeParameter(

From 2b52429b790f493aedd92fdba022a9f845344f1a Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Thu, 22 Sep 2022 13:59:23 +0200
Subject: [PATCH 09/35] Refine error message.

---
 ctapipe/image/extractor.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 1c0ccb0c660..1cd9419c47f 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1305,7 +1305,7 @@ def deconvolution_parameters(
 
     if camera_sample_width < ref_sample_width:
         raise ValueError(
-            f"ref_sample_width must be equal to or shorter than camera_sample_width, got ref_sample_width={ref_sample_width} and camera_sample_width={camera_sample_width}"
+            f"ref_sample_width (got {ref_sample_width}) must be equal to or shorter than camera_sample_width (got {camera_sample_width}); need a reference single p.e. pulse shape with finer sampling!"
         )
     avg_step = int(camera_sample_width / ref_sample_width + 0.5)
 

From 789ccc79fadce50e60a4c90663e3f9b1601ee442 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Fri, 23 Sep 2022 14:23:30 +0200
Subject: [PATCH 10/35] Improve compatibility with traitlets >4.0.

---
 ctapipe/image/extractor.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 1cd9419c47f..cf03be2debb 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1364,12 +1364,12 @@ def deconvolve(waveforms, bls, up: int, pz: float):
 
 class FlashCamExtractor(ImageExtractor):
     upsampling = IntTelescopeParameter(
-        default_value=4, min=1, help="Define the upsampling factor for waveforms"
-    ).tag(config=True)
+        default_value=4, help="Define the upsampling factor for waveforms"
+    ).tag(config=True, min=1)
 
     window_width = IntTelescopeParameter(
-        default_value=4, min=1, help="Define the width of the integration window"
-    ).tag(config=True)
+        default_value=4, help="Define the width of the integration window"
+    ).tag(config=True, min=1)
 
     window_shift = IntTelescopeParameter(
         default_value=2,

From b2a072341f2814d5c4e4ce82e51565fe5bff1339 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Fri, 23 Sep 2022 14:25:31 +0200
Subject: [PATCH 11/35] Test FlashCamExtractor on prod5 simulation set. Extend
 error bounds to make tests pass.

---
 ctapipe/image/tests/test_extractor.py | 40 ++++++++++++++++++---------
 1 file changed, 27 insertions(+), 13 deletions(-)

diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index d55240d6067..714d53c27b6 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -641,7 +641,8 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray):
     assert np.allclose(dl1.peak_time[bright_pixels], expected / sample_rate)
 
 
-def test_flashcam_extractor(toymodel_1_MST_FC):
+def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
+    # Test on toy model
     (
         waveforms,
         subarray,
@@ -650,20 +651,33 @@ def test_flashcam_extractor(toymodel_1_MST_FC):
         true_charge,
         true_time,
     ) = toymodel_1_MST_FC
-
-    assert(subarray.tel[1].camera.name == "FlashCam")
     extractor = FlashCamExtractor(subarray=subarray)
-
-    # Test toymodel
     broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
     dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
-    assert_allclose(dl1.image, true_charge, rtol=0.1)
-    assert_allclose(dl1.peak_time, true_time, atol=2.5)  # FIXME could be much closer for large pulses...
     assert dl1.is_valid == True
 
-    # Test prod5 simulations
-    path = "gamma_20deg_0deg_run2___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10-100evts.simtel.zst"
-    source = EventSource(f"dataset://{path}")
-    subarray = source.subarray
-    print(subarray.get_tel_ids_for_type("MST_MST_FlashCam"))
-    # TODO Test!
+    # TODO Shrink tolerances after optimising the extractor parameters
+    assert_allclose(dl1.image, true_charge, rtol=0.2)
+    assert_allclose(dl1.peak_time, true_time, atol=2.5)
+
+    # Test on prod5 simulations
+    with EventSource(prod5_gamma_simtel_path) as source:
+        subarray = source.subarray
+        extractor = extractor = FlashCamExtractor(subarray)
+        for event in source:
+            for tel_id in subarray.get_tel_ids_for_type("MST_MST_FlashCam"):
+                true_charge = event.simulation.tel[tel_id].true_image
+                if true_charge is None:
+                    continue  # telescope did not trigger
+
+                waveforms = event.r1.tel[tel_id].waveform
+                selected_gain_channel = np.zeros(waveforms.shape[0], dtype=np.int64)
+                broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+
+                dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
+                assert dl1.is_valid == True
+
+                bright_pixels = true_charge > 30
+                assert_allclose(
+                    dl1.image[bright_pixels], true_charge[bright_pixels], rtol=0.33
+                )

From 89415b83bd6c47a934cfa04dfa1aa20abc3f69dd Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Fri, 23 Sep 2022 15:24:22 +0200
Subject: [PATCH 12/35] Add docstrings and type annotations.

---
 ctapipe/image/extractor.py | 83 +++++++++++++++++++++++++++++++-------
 1 file changed, 69 insertions(+), 14 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index cf03be2debb..23e9810d824 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -23,9 +23,10 @@
 
 from abc import abstractmethod
 from functools import lru_cache
-from typing import Tuple
+from typing import List, Tuple
 
 import numpy as np
+import numpy.typing as npt
 from numba import float32, float64, guvectorize, int64, njit, prange
 from scipy.ndimage import convolve1d
 from scipy.signal import filtfilt
@@ -1293,7 +1294,32 @@ def __call__(
 @lru_cache
 def deconvolution_parameters(
     camera: CameraDescription, upsampling: int, window_width: int, window_shift: int
-):
+) -> Tuple[List[float], List[float], List[float]]:
+    """
+    Estimates deconvolution and recalibration parameters from the camera's reference
+    single-p.e. pulse shape for the given configuration of FlashCamExtractor.
+
+    Parameters
+    ----------
+    camera : CameraDescription
+        Description of the target camera.
+    upsampling : int
+        Upsampling factor (>= 1); see also `deconvolve(...)`.
+    window_width : int
+        Integration window width (>= 1); see also `extract_around_peak(...)`.
+    window_shift : int
+        Shift of the integration window relative to the peak; see also
+        `extract_around_peak(...)`.
+
+    Returns
+    -------
+    pole_zeros : list of floats
+        Pole-zero parameter for each channel to be passed to `deconvolve(...)`.
+    gain_losses : list of floats
+        Gain loss of each channel that needs to be corrected after deconvolution.
+    time_shift_nsec : list of floats
+        Timing shift of each channel that needs to be corrected after deconvolution.
+    """
     if upsampling < 1:
         raise ValueError(f"upsampling must be > 0, got {upsampling}")
     if window_width < 1:
@@ -1309,7 +1335,7 @@ def deconvolution_parameters(
         )
     avg_step = int(camera_sample_width / ref_sample_width + 0.5)
 
-    pzs = []
+    pzs = []  # avg. pole-zero deconvolution parameters
     for ref_pulse_shape in ref_pulse_shapes:
         phase_pzs = []
         for phase in range(avg_step):
@@ -1320,11 +1346,11 @@ def deconvolution_parameters(
 
         if len(phase_pzs) == 0:
             raise ValueError(
-                "ref_pulse_shape is malformed - cannot find deconvolution scale"
+                "ref_pulse_shape is malformed - cannot find deconvolution parameter"
             )
         pzs.append(np.mean(phase_pzs))
 
-    gains, shifts = [], []
+    gains, shifts = [], []  # avg. gains and timing shifts of the deconvolved pulses
     for pz, ref_pulse_shape in zip(pzs, ref_pulse_shapes):
         integral = ref_pulse_shape.sum() * ref_sample_width / camera_sample_width
         phase_gains, phase_shifts = [], []
@@ -1351,15 +1377,44 @@ def deconvolution_parameters(
     return pzs, gains, shifts
 
 
-def deconvolve(waveforms, bls, up: int, pz: float):
-    waveforms = np.atleast_2d(waveforms)
-    bls = np.atleast_2d(bls).T
-    y = waveforms - bls
-    y[:, 1:] -= pz * y[:, :-1]
-    y[:, 0] = 0
-    if up > 1:
-        return filtfilt(np.ones(up), up, np.repeat(y, up, axis=-1))
-    return y
+def deconvolve(
+    waveforms: npt.ArrayLike,
+    baselines: npt.ArrayLike,
+    upsampling: int,
+    pole_zero: float,
+) -> np.ndarray:
+    """
+    Applies pole-zero deconvolution and upsampling to pixel waveforms.
+
+    Parameters
+    ----------
+    waveforms : ndarray
+        Waveforms stored in a numpy array.
+        Shape: (n_pix, n_samples)
+    baselines : ndarray or float
+        Baseline estimates for each pixel.
+        Shape: (n_pix, ) or scalar
+    upsampling : int
+        Upsampling factor to use (>= 1); the input waveforms are resampled at upsampling times their original sampling rate.
+    pole_zero : float
+        Deconvolution factor obtained from `deconvolution_parameters(...)` applied to the reference single p.e. pulse shape.
+
+    Returns
+    -------
+    deconvolved_waveforms : ndarray
+        Deconvolved and upsampled waveforms stored in a numpy array.
+        Shape: (n_pix, upsampling * n_samples)
+    """
+    deconvolved_waveforms = np.atleast_2d(waveforms) - np.atleast_2d(baselines).T
+    deconvolved_waveforms[:, 1:] -= pole_zero * deconvolved_waveforms[:, :-1]
+    deconvolved_waveforms[:, 0] = 0
+    if upsampling > 1:
+        return filtfilt(
+            np.ones(upsampling),
+            upsampling,
+            np.repeat(deconvolved_waveforms, upsampling, axis=-1),
+        )
+    return deconvolved_waveforms
 
 
 class FlashCamExtractor(ImageExtractor):

From 56100632c0fca4a28a0e52f1eb088c5e5a5c397f Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Fri, 23 Sep 2022 15:38:26 +0200
Subject: [PATCH 13/35] Polish docstrings.

---
 ctapipe/image/extractor.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 23e9810d824..1feb781f06f 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1384,7 +1384,9 @@ def deconvolve(
     pole_zero: float,
 ) -> np.ndarray:
     """
-    Applies pole-zero deconvolution and upsampling to pixel waveforms.
+    Applies pole-zero deconvolution and upsampling to pixel waveforms. Use
+    `deconvolution_parameters(...)` to estimate the required `pole_zero` parameter
+    for the specific camera model.
 
     Parameters
     ----------
@@ -1392,12 +1394,14 @@ def deconvolve(
         Waveforms stored in a numpy array.
         Shape: (n_pix, n_samples)
     baselines : ndarray or float
-        Baseline estimates for each pixel.
+        Baseline estimates for each pixel that are subtracted from the waveforms
+        before deconvolution.
         Shape: (n_pix, ) or scalar
     upsampling : int
-        Upsampling factor to use (>= 1); the input waveforms are resampled at upsampling times their original sampling rate.
+        Upsampling factor to use (>= 1); if > 1, the input waveforms are resampled
+        at upsampling times their original sampling rate.
     pole_zero : float
-        Deconvolution factor obtained from `deconvolution_parameters(...)` applied to the reference single p.e. pulse shape.
+        Deconvolution parameter obtained from `deconvolution_parameters(...)`.
 
     Returns
     -------

From 84f18359f1f6a0f91f883d0267c8f9717ea0a4f1 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Fri, 23 Sep 2022 16:37:03 +0200
Subject: [PATCH 14/35] Improve consistency in variable naming and error
 messages.

---
 ctapipe/image/extractor.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 1feb781f06f..b58ae9a6289 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1317,7 +1317,7 @@ def deconvolution_parameters(
         Pole-zero parameter for each channel to be passed to `deconvolve(...)`.
     gain_losses : list of floats
         Gain loss of each channel that needs to be corrected after deconvolution.
-    time_shift_nsec : list of floats
+    time_shifts_nsec : list of floats
         Timing shift of each channel that needs to be corrected after deconvolution.
     """
     if upsampling < 1:
@@ -1326,14 +1326,14 @@ def deconvolution_parameters(
         raise ValueError(f"window_width must be > 0, got {window_width}")
 
     ref_pulse_shapes = camera.readout.reference_pulse_shape
-    ref_sample_width = camera.readout.reference_pulse_sample_width.to_value("ns")
-    camera_sample_width = 1.0 / camera.readout.sampling_rate.to_value("GHz")
+    ref_sample_width_nsec = camera.readout.reference_pulse_sample_width.to_value("ns")
+    camera_sample_width_nsec = 1.0 / camera.readout.sampling_rate.to_value("GHz")
 
-    if camera_sample_width < ref_sample_width:
+    if camera_sample_width_nsec < ref_sample_width_nsec:
         raise ValueError(
-            f"ref_sample_width (got {ref_sample_width}) must be equal to or shorter than camera_sample_width (got {camera_sample_width}); need a reference single p.e. pulse shape with finer sampling!"
+            f"Reference pulse sampling time (got {ref_sample_width_nsec} ns) must be equal to or shorter than the camera sampling time (got {camera_sample_width_nsec} ns); need a reference single p.e. pulse shape with finer sampling!"
         )
-    avg_step = int(camera_sample_width / ref_sample_width + 0.5)
+    avg_step = int(camera_sample_width_nsec / ref_sample_width_nsec + 0.5)
 
     pzs = []  # avg. pole-zero deconvolution parameters
     for ref_pulse_shape in ref_pulse_shapes:
@@ -1352,7 +1352,9 @@ def deconvolution_parameters(
 
     gains, shifts = [], []  # avg. gains and timing shifts of the deconvolved pulses
     for pz, ref_pulse_shape in zip(pzs, ref_pulse_shapes):
-        integral = ref_pulse_shape.sum() * ref_sample_width / camera_sample_width
+        integral = (
+            ref_pulse_shape.sum() * ref_sample_width_nsec / camera_sample_width_nsec
+        )
         phase_gains, phase_shifts = [], []
         for phase in range(avg_step):
             x = ref_pulse_shape[phase::avg_step]
@@ -1363,7 +1365,7 @@ def deconvolution_parameters(
             if start >= 0 and stop <= y.size:
                 phase_shifts.append(
                     (i_max / upsampling * avg_step - ref_pulse_shape.argmax())
-                    * ref_sample_width
+                    * ref_sample_width_nsec
                 )
                 phase_gains.append(y[start:stop].sum() / integral)
 

From f0b08ebd812cbb391fbea8791b1775cc06b6f77f Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Fri, 23 Sep 2022 16:49:46 +0200
Subject: [PATCH 15/35] Remove obsolete parameter.

---
 ctapipe/resources/base_config.yaml | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/ctapipe/resources/base_config.yaml b/ctapipe/resources/base_config.yaml
index 48ce249d7b1..eca7da71e87 100644
--- a/ctapipe/resources/base_config.yaml
+++ b/ctapipe/resources/base_config.yaml
@@ -43,9 +43,6 @@ CameraCalibrator:
     - ['type', '*', 'NeighborPeakWindowSum']
     - ['type', '*FlashCam', 'FlashCamExtractor']
 
-  FlashCamExtractor:
-    exponential_scale: 7.5
-
 # The ImageProcessor performs the DL1a-> DL1b (image parameters) transition. It
 # is run only if the parameters `DataWriter.write_image_parameters=True` and the
 # parameters don't already exist in the input file (or if the user forces them

From f4a64d3e1a2198e054b573ee8d1016560a3bbad7 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Mon, 26 Sep 2022 15:22:08 +0200
Subject: [PATCH 16/35] Add configurable compensation for effective pulse time
 profile.

---
 ctapipe/image/extractor.py | 33 +++++++++++++++++++++++++++++++--
 1 file changed, 31 insertions(+), 2 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index b58ae9a6289..b22a15b77f1 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -23,10 +23,11 @@
 
 from abc import abstractmethod
 from functools import lru_cache
-from typing import List, Tuple
+from typing import Callable, List, Optional, Tuple
 
 import numpy as np
 import numpy.typing as npt
+import scipy.stats
 from numba import float32, float64, guvectorize, int64, njit, prange
 from scipy.ndimage import convolve1d
 from scipy.signal import filtfilt
@@ -1293,7 +1294,11 @@ def __call__(
 
 @lru_cache
 def deconvolution_parameters(
-    camera: CameraDescription, upsampling: int, window_width: int, window_shift: int
+    camera: CameraDescription,
+    upsampling: int,
+    window_width: int,
+    window_shift: int,
+    time_profile_pdf: Optional[Callable[[npt.ArrayLike], npt.ArrayLike]] = None,
 ) -> Tuple[List[float], List[float], List[float]]:
     """
     Estimates deconvolution and recalibration parameters from the camera's reference
@@ -1310,6 +1315,11 @@ def deconvolution_parameters(
     window_shift : int
         Shift of the integration window relative to the peak; see also
         `extract_around_peak(...)`.
+    time_profile_pdf : callable or None
+        PDF of the assumed effective Cherenkov time profile to assume when
+        calculating the gain loss; takes nanoseconds as arguments and returns
+        probability density (with mode at ~0 ns); default: None (assume
+        instantaneous pulse).
 
     Returns
     -------
@@ -1352,6 +1362,13 @@ def deconvolution_parameters(
 
     gains, shifts = [], []  # avg. gains and timing shifts of the deconvolved pulses
     for pz, ref_pulse_shape in zip(pzs, ref_pulse_shapes):
+        if time_profile_pdf:  # convolve ref_pulse_shape with time profile PDF
+            t = (
+                np.arange(ref_pulse_shape.size) - ref_pulse_shape.size / 2
+            ) * ref_sample_width_nsec
+            time_profile = time_profile_pdf(t)
+            ref_pulse_shape = np.convolve(ref_pulse_shape, time_profile, "same")
+
         integral = (
             ref_pulse_shape.sum() * ref_sample_width_nsec / camera_sample_width_nsec
         )
@@ -1448,6 +1465,12 @@ class FlashCamExtractor(ImageExtractor):
         "1: local pixel counts as much as any neighbor)",
     ).tag(config=True)
 
+    effective_time_profile_std = FloatTelescopeParameter(
+        default_value=0.0,
+        help="Effective Cherenkov time profile std. dev. (in nanoseconds) to "
+        "assume for calculating the gain correction",
+    ).tag(config=True, min=0)
+
     def __init__(self, subarray, **kwargs):
         super().__init__(subarray=subarray, **kwargs)
 
@@ -1456,12 +1479,18 @@ def __init__(self, subarray, **kwargs):
             for tel_id, telescope in subarray.tel.items()
         }
 
+        def time_profile_pdf_gen(std_dev: float):
+            if std_dev == 0:
+                return None
+            return scipy.stats.norm(0.0, std_dev).pdf
+
         self.deconvolution_pars = {
             tel_id: deconvolution_parameters(
                 tel.camera,
                 self.upsampling.tel[tel_id],
                 self.window_width.tel[tel_id],
                 self.window_shift.tel[tel_id],
+                time_profile_pdf_gen(self.effective_time_profile_std.tel[tel_id]),
             )
             for tel_id, tel in subarray.tel.items()
         }

From 67261469682b60692bc3364c76026e53fc7955fc Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Tue, 4 Oct 2022 15:57:12 +0200
Subject: [PATCH 17/35] Make baseline estimator configurable and disable in
 default (example) config (may be removed in future).

---
 ctapipe/image/extractor.py         | 30 +++++++++++++++++++++++++-----
 ctapipe/resources/base_config.yaml |  3 +++
 2 files changed, 28 insertions(+), 5 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index b22a15b77f1..851de1aea68 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1471,6 +1471,16 @@ class FlashCamExtractor(ImageExtractor):
         "assume for calculating the gain correction",
     ).tag(config=True, min=0)
 
+    baseline_window_start = IntTelescopeParameter(
+        default_value=0,
+        help="Define the first sample to use for baseline estimation",
+    ).tag(config=True, min=0)
+
+    baseline_window_width = IntTelescopeParameter(
+        default_value=1,
+        help="Define the number of samples to use for baseline estimation",
+    ).tag(config=True, min=0)
+
     def __init__(self, subarray, **kwargs):
         super().__init__(subarray=subarray, **kwargs)
 
@@ -1499,13 +1509,23 @@ def __call__(
         self, waveforms, tel_id, selected_gain_channel, broken_pixels
     ) -> DL1CameraContainer:
         upsampling = self.upsampling.tel[tel_id]
-        window_width = self.window_width.tel[tel_id]
-        window_shift = self.window_shift.tel[tel_id]
+        baseline_window_start = self.baseline_window_start.tel[tel_id]
+        baseline_window_width = self.baseline_window_width.tel[tel_id]
+        integration_window_width = self.window_width.tel[tel_id]
+        integration_window_shift = self.window_shift.tel[tel_id]
 
         pzs, gains, shifts = self.deconvolution_pars[tel_id]
         pz, gain, shift = pzs[0], gains[0], shifts[0]
 
-        bls = waveforms[:, 0]  # FIXME fetch from NSB estimates
+        if (
+            baseline_window_width == 0
+        ):  # assume baseline has been subtracted properly in R1 pre-calibration step
+            bls = 0.0
+        else:
+            bls = waveforms[
+                :, baseline_window_start : baseline_window_start + baseline_window_width
+            ].mean(1)
+
         waveforms = deconvolve(waveforms, bls, upsampling, pz)
 
         # FIXME near-duplicate of neighbour peak sum for now
@@ -1521,8 +1541,8 @@ def __call__(
         charge, peak_time = extract_around_peak(
             waveforms,
             peak_index,
-            window_width,
-            window_shift,
+            integration_window_width,
+            integration_window_shift,
             self.sampling_rate_ghz[tel_id] * upsampling,
         )
 
diff --git a/ctapipe/resources/base_config.yaml b/ctapipe/resources/base_config.yaml
index eca7da71e87..990f60a97f0 100644
--- a/ctapipe/resources/base_config.yaml
+++ b/ctapipe/resources/base_config.yaml
@@ -43,6 +43,9 @@ CameraCalibrator:
     - ['type', '*', 'NeighborPeakWindowSum']
     - ['type', '*FlashCam', 'FlashCamExtractor']
 
+FlashCamExtractor:
+  baseline_window_width: 0
+
 # The ImageProcessor performs the DL1a-> DL1b (image parameters) transition. It
 # is run only if the parameters `DataWriter.write_image_parameters=True` and the
 # parameters don't already exist in the input file (or if the user forces them

From 937ecc921bcee7928bd110bafa48a43326644b96 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Tue, 4 Oct 2022 16:23:54 +0200
Subject: [PATCH 18/35] Assume 2 ns as default Cherenkov time profile width.

---
 ctapipe/resources/base_config.yaml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/ctapipe/resources/base_config.yaml b/ctapipe/resources/base_config.yaml
index 990f60a97f0..a4482d31143 100644
--- a/ctapipe/resources/base_config.yaml
+++ b/ctapipe/resources/base_config.yaml
@@ -45,6 +45,7 @@ CameraCalibrator:
 
 FlashCamExtractor:
   baseline_window_width: 0
+  effective_time_profile_std: 2.0
 
 # The ImageProcessor performs the DL1a-> DL1b (image parameters) transition. It
 # is run only if the parameters `DataWriter.write_image_parameters=True` and the

From cd0c2e7994a3eee6edbd2f529a3248b6abb601bb Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Fri, 7 Oct 2022 17:28:30 +0200
Subject: [PATCH 19/35] Add preliminary option to define neighbour-based
 windows based on the integrated signals.

---
 ctapipe/image/extractor.py | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 851de1aea68..ada508bc871 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1481,6 +1481,12 @@ class FlashCamExtractor(ImageExtractor):
         help="Define the number of samples to use for baseline estimation",
     ).tag(config=True, min=0)
 
+    integral_based_windows = BoolTelescopeParameter(
+        default_value=False,
+        help="Whether to define the integration windows based on the neighbour averages"
+        "of the integrated waveforms",
+    ).tag(config=True)
+
     def __init__(self, subarray, **kwargs):
         super().__init__(subarray=subarray, **kwargs)
 
@@ -1528,10 +1534,17 @@ def __call__(
 
         waveforms = deconvolve(waveforms, bls, upsampling, pz)
 
+        if self.integral_based_windows.tel[tel_id]:
+            nn_waveforms = scipy.signal.convolve(
+                waveforms, np.ones((1, integration_window_width)), "same"
+            )
+        else:
+            nn_waveforms = waveforms
+
         # FIXME near-duplicate of neighbour peak sum for now
         neighbors = self.subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse
         peak_index = neighbor_average_maximum(
-            waveforms,
+            nn_waveforms,
             neighbors_indices=neighbors.indices,
             neighbors_indptr=neighbors.indptr,
             local_weight=self.local_weight.tel[tel_id],

From 306c269adc0f29658a4e32b355e799dc5c39ada5 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Mon, 17 Oct 2022 07:32:29 +0200
Subject: [PATCH 20/35] Cleanup, update defaults and add option to clip inputs
 to neighbour summation.

---
 ctapipe/image/extractor.py            | 51 +++++++++------------------
 ctapipe/image/tests/test_extractor.py | 13 +++----
 ctapipe/resources/base_config.yaml    |  4 ---
 3 files changed, 24 insertions(+), 44 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index ada508bc871..322c82691be 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1446,11 +1446,11 @@ class FlashCamExtractor(ImageExtractor):
     ).tag(config=True, min=1)
 
     window_width = IntTelescopeParameter(
-        default_value=4, help="Define the width of the integration window"
+        default_value=7, help="Define the width of the integration window"
     ).tag(config=True, min=1)
 
     window_shift = IntTelescopeParameter(
-        default_value=2,
+        default_value=3,
         help="Define the shift of the integration window from the peak_index "
         "(peak_index - shift)",
     ).tag(config=True)
@@ -1466,25 +1466,15 @@ class FlashCamExtractor(ImageExtractor):
     ).tag(config=True)
 
     effective_time_profile_std = FloatTelescopeParameter(
-        default_value=0.0,
+        default_value=2.0,
         help="Effective Cherenkov time profile std. dev. (in nanoseconds) to "
         "assume for calculating the gain correction",
     ).tag(config=True, min=0)
 
-    baseline_window_start = IntTelescopeParameter(
-        default_value=0,
-        help="Define the first sample to use for baseline estimation",
-    ).tag(config=True, min=0)
-
-    baseline_window_width = IntTelescopeParameter(
-        default_value=1,
-        help="Define the number of samples to use for baseline estimation",
-    ).tag(config=True, min=0)
-
-    integral_based_windows = BoolTelescopeParameter(
-        default_value=False,
-        help="Whether to define the integration windows based on the neighbour averages"
-        "of the integrated waveforms",
+    neighbour_sum_clipping = FloatTelescopeParameter(
+        default_value=5.0,
+        help="(Soft) clipping level of a pixel's contribution to a neighbour sum "
+        "(set to 0 or inf to disable clipping)",
     ).tag(config=True)
 
     def __init__(self, subarray, **kwargs):
@@ -1511,35 +1501,28 @@ def time_profile_pdf_gen(std_dev: float):
             for tel_id, tel in subarray.tel.items()
         }
 
+    @staticmethod
+    def clip(x, lo=0.0, hi=np.inf):
+        """Applies soft clipping to ±1 and then hard clipping to (lo, hi)."""
+        return np.clip(x / (1.0 + np.abs(x)), lo, hi)
+
     def __call__(
         self, waveforms, tel_id, selected_gain_channel, broken_pixels
     ) -> DL1CameraContainer:
         upsampling = self.upsampling.tel[tel_id]
-        baseline_window_start = self.baseline_window_start.tel[tel_id]
-        baseline_window_width = self.baseline_window_width.tel[tel_id]
         integration_window_width = self.window_width.tel[tel_id]
         integration_window_shift = self.window_shift.tel[tel_id]
+        neighbour_sum_clipping = self.neighbour_sum_clipping.tel[tel_id]
 
         pzs, gains, shifts = self.deconvolution_pars[tel_id]
         pz, gain, shift = pzs[0], gains[0], shifts[0]
 
-        if (
-            baseline_window_width == 0
-        ):  # assume baseline has been subtracted properly in R1 pre-calibration step
-            bls = 0.0
-        else:
-            bls = waveforms[
-                :, baseline_window_start : baseline_window_start + baseline_window_width
-            ].mean(1)
-
-        waveforms = deconvolve(waveforms, bls, upsampling, pz)
+        waveforms = deconvolve(waveforms, 0.0, upsampling, pz)
 
-        if self.integral_based_windows.tel[tel_id]:
-            nn_waveforms = scipy.signal.convolve(
-                waveforms, np.ones((1, integration_window_width)), "same"
-            )
-        else:
+        if neighbour_sum_clipping == 0.0 or np.isinf(neighbour_sum_clipping):
             nn_waveforms = waveforms
+        else:
+            nn_waveforms = self.clip(waveforms / neighbour_sum_clipping)
 
         # FIXME near-duplicate of neighbour peak sum for now
         neighbors = self.subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse
diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 714d53c27b6..2053623bc7a 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -656,14 +656,13 @@ def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
     dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
     assert dl1.is_valid == True
 
-    # TODO Shrink tolerances after optimising the extractor parameters
-    assert_allclose(dl1.image, true_charge, rtol=0.2)
+    assert_allclose(dl1.image, true_charge, rtol=0.15)
     assert_allclose(dl1.peak_time, true_time, atol=2.5)
 
     # Test on prod5 simulations
     with EventSource(prod5_gamma_simtel_path) as source:
         subarray = source.subarray
-        extractor = extractor = FlashCamExtractor(subarray)
+        extractor = FlashCamExtractor(subarray)
         for event in source:
             for tel_id in subarray.get_tel_ids_for_type("MST_MST_FlashCam"):
                 true_charge = event.simulation.tel[tel_id].true_image
@@ -672,12 +671,14 @@ def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
 
                 waveforms = event.r1.tel[tel_id].waveform
                 selected_gain_channel = np.zeros(waveforms.shape[0], dtype=np.int64)
-                broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+                broken_pixels = event.mon.tel[
+                    tel_id
+                ].pixel_status.hardware_failing_pixels[0]
 
                 dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
                 assert dl1.is_valid == True
 
-                bright_pixels = true_charge > 30
+                bright_pixels = (true_charge > 30) * ~broken_pixels
                 assert_allclose(
-                    dl1.image[bright_pixels], true_charge[bright_pixels], rtol=0.33
+                    dl1.image[bright_pixels], true_charge[bright_pixels], rtol=0.35
                 )
diff --git a/ctapipe/resources/base_config.yaml b/ctapipe/resources/base_config.yaml
index a4482d31143..eca7da71e87 100644
--- a/ctapipe/resources/base_config.yaml
+++ b/ctapipe/resources/base_config.yaml
@@ -43,10 +43,6 @@ CameraCalibrator:
     - ['type', '*', 'NeighborPeakWindowSum']
     - ['type', '*FlashCam', 'FlashCamExtractor']
 
-FlashCamExtractor:
-  baseline_window_width: 0
-  effective_time_profile_std: 2.0
-
 # The ImageProcessor performs the DL1a-> DL1b (image parameters) transition. It
 # is run only if the parameters `DataWriter.write_image_parameters=True` and the
 # parameters don't already exist in the input file (or if the user forces them

From 3471f9ae6c8b3d6bb009ec6a0d6660ea73d17878 Mon Sep 17 00:00:00 2001
From: Felix Werner <felix.werner@mpi-hd.mpg.de>
Date: Fri, 21 Oct 2022 17:45:10 +0200
Subject: [PATCH 21/35] Implement leading edge timing.

---
 ctapipe/image/extractor.py            | 129 ++++++++++++++++++++++++--
 ctapipe/image/tests/test_extractor.py |   4 +-
 2 files changed, 122 insertions(+), 11 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 322c82691be..b4f5101f3d0 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1298,6 +1298,8 @@ def deconvolution_parameters(
     upsampling: int,
     window_width: int,
     window_shift: int,
+    leading_edge_timing: bool,
+    leading_edge_rel_descend_limit: float,
     time_profile_pdf: Optional[Callable[[npt.ArrayLike], npt.ArrayLike]] = None,
 ) -> Tuple[List[float], List[float], List[float]]:
     """
@@ -1315,6 +1317,10 @@ def deconvolution_parameters(
     window_shift : int
         Shift of the integration window relative to the peak; see also
         `extract_around_peak(...)`.
+    leading_edge_timing : bool
+        Whether time calculation will be done on the leading edge.
+    leading_edge_rel_descend_limit : float
+        If leading edge timing is used, the fraction of the peak value down to which samples are accumulated in the centroid calculation.
     time_profile_pdf : callable or None
         PDF of the assumed effective Cherenkov time profile to assume when
         calculating the gain loss; takes nanoseconds as arguments and returns
@@ -1329,6 +1335,8 @@ def deconvolution_parameters(
         Gain loss of each channel that needs to be corrected after deconvolution.
     time_shifts_nsec : list of floats
         Timing shift of each channel that needs to be corrected after deconvolution.
+    leading_edge_shifts : list of floats
+        Offset of the leading edge peak w.r.t. the deconvolved peak.
     """
     if upsampling < 1:
         raise ValueError(f"upsampling must be > 0, got {upsampling}")
@@ -1360,7 +1368,7 @@ def deconvolution_parameters(
             )
         pzs.append(np.mean(phase_pzs))
 
-    gains, shifts = [], []  # avg. gains and timing shifts of the deconvolved pulses
+    gains, shifts, pz2d_shifts = [], [], []  # avg. gains and timing shifts
     for pz, ref_pulse_shape in zip(pzs, ref_pulse_shapes):
         if time_profile_pdf:  # convolve ref_pulse_shape with time profile PDF
             t = (
@@ -1372,16 +1380,26 @@ def deconvolution_parameters(
         integral = (
             ref_pulse_shape.sum() * ref_sample_width_nsec / camera_sample_width_nsec
         )
-        phase_gains, phase_shifts = [], []
+        phase_gains, phase_shifts, phase_pz2d_shifts = [], [], []
         for phase in range(avg_step):
-            x = ref_pulse_shape[phase::avg_step]
-            y = deconvolve(np.atleast_2d(x), 0.0, upsampling, pz)[0]
+            x = np.atleast_2d(ref_pulse_shape[phase::avg_step])
+            y = deconvolve(x, 0.0, upsampling, pz)[0]
             i_max = y.argmax()
             start = i_max - window_shift  # TODO should use extract_around_peak here
             stop = start + window_width
             if start >= 0 and stop <= y.size:
+                if leading_edge_timing:
+                    d = deconvolve(x, 0.0, upsampling, 1)[0]
+                    d_pk_idx = d.argmax()
+                    phase_pz2d_shifts.append(i_max - d_pk_idx)
+                    time = adaptive_centroid(
+                        d, d_pk_idx, leading_edge_rel_descend_limit
+                    )
+                else:
+                    time = i_max
+
                 phase_shifts.append(
-                    (i_max / upsampling * avg_step - ref_pulse_shape.argmax())
+                    (time / upsampling * avg_step - ref_pulse_shape.argmax())
                     * ref_sample_width_nsec
                 )
                 phase_gains.append(y[start:stop].sum() / integral)
@@ -1392,8 +1410,9 @@ def deconvolution_parameters(
             )
         gains.append(np.mean(phase_gains))
         shifts.append(np.mean(phase_shifts))
+        pz2d_shifts.append(np.mean(phase_pz2d_shifts))
 
-    return pzs, gains, shifts
+    return pzs, gains, shifts, pz2d_shifts
 
 
 def deconvolve(
@@ -1440,6 +1459,75 @@ def deconvolve(
     return deconvolved_waveforms
 
 
+@guvectorize(
+    [
+        (float64[:], int64, float64, float32[:]),
+        (float32[:], int64, float64, float32[:]),
+    ],
+    "(s),(),()->()",
+    nopython=True,
+    cache=True,
+)
+def adaptive_centroid(waveforms, peak_index, rel_descend_limit, cog):
+    """
+    Calculates the pulse centroid for all samples down to rel_descend_limit * peak_amplitude.
+
+    The ret argument is required by numpy to create the numpy array which is
+    returned. It can be ignored when calling this function.
+
+    Parameters
+    ----------
+    waveforms : ndarray
+        Waveforms stored in a numpy array.
+        Shape: (n_pix, n_samples)
+    peak_index : ndarray or int
+        Peak index for each pixel.
+    rel_descend_limit : ndarray or float
+        Fraction of the peak value down to which samples are accumulated in the centroid calculation.
+    cog : ndarray
+        Return argument for ufunc (ignore)
+        Returns the peak centroid in units "samples"
+
+    Returns
+    -------
+    cog : ndarray
+        Peak centroid in units "samples"
+    """
+    cog[0] = peak_index  # preload in case of errors
+
+    nsamples = waveforms.size
+    if nsamples == 0:
+        return
+
+    peak_index = max(0, min(peak_index, nsamples - 1))
+    peak_amplitude = waveforms[peak_index]
+    if peak_amplitude < 0.0:
+        return
+
+    descend_limit = rel_descend_limit * peak_amplitude
+
+    sum = float64(0.0)
+    jsum = float64(0.0)
+    j = peak_index
+    while j >= 0 and waveforms[j] > descend_limit:
+        sum += waveforms[j]
+        jsum += j * waveforms[j]
+        j -= 1
+        if waveforms[j] > peak_amplitude:
+            descend_limit = rel_descend_limit * peak_amplitude
+
+    j = peak_index + 1
+    while j < nsamples and waveforms[j] > descend_limit:
+        sum += waveforms[j]
+        jsum += j * waveforms[j]
+        j += 1
+        if waveforms[j] > peak_amplitude:
+            descend_limit = rel_descend_limit * peak_amplitude
+
+    if sum != 0.0:
+        cog[0] = jsum / sum
+
+
 class FlashCamExtractor(ImageExtractor):
     upsampling = IntTelescopeParameter(
         default_value=4, help="Define the upsampling factor for waveforms"
@@ -1477,6 +1565,16 @@ class FlashCamExtractor(ImageExtractor):
         "(set to 0 or inf to disable clipping)",
     ).tag(config=True)
 
+    leading_edge_timing = BoolTelescopeParameter(
+        default_value=True, help="Calculate leading edge time instead of peak time"
+    ).tag(config=False)
+
+    leading_edge_rel_descend_limit = FloatTelescopeParameter(
+        default_value=0.05,
+        help="Fraction of the peak value down to which samples are accumulated "
+        "in the leading edge centroid calculation",
+    ).tag(config=True, min=0.0, max=1.0)
+
     def __init__(self, subarray, **kwargs):
         super().__init__(subarray=subarray, **kwargs)
 
@@ -1496,6 +1594,8 @@ def time_profile_pdf_gen(std_dev: float):
                 self.upsampling.tel[tel_id],
                 self.window_width.tel[tel_id],
                 self.window_shift.tel[tel_id],
+                self.leading_edge_timing.tel[tel_id],
+                self.leading_edge_rel_descend_limit.tel[tel_id],
                 time_profile_pdf_gen(self.effective_time_profile_std.tel[tel_id]),
             )
             for tel_id, tel in subarray.tel.items()
@@ -1513,9 +1613,14 @@ def __call__(
         integration_window_width = self.window_width.tel[tel_id]
         integration_window_shift = self.window_shift.tel[tel_id]
         neighbour_sum_clipping = self.neighbour_sum_clipping.tel[tel_id]
+        leading_edge_timing = self.leading_edge_timing.tel[tel_id]
+        leading_edge_rel_descend_limit = self.leading_edge_rel_descend_limit.tel[tel_id]
+
+        pzs, gains, shifts, pz2ds = self.deconvolution_pars[tel_id]
+        pz, gain, shift, pz2d = pzs[0], gains[0], shifts[0], pz2ds[0]
 
-        pzs, gains, shifts = self.deconvolution_pars[tel_id]
-        pz, gain, shift = pzs[0], gains[0], shifts[0]
+        if leading_edge_timing:
+            d_waveforms = deconvolve(waveforms, 0.0, upsampling, 1)
 
         waveforms = deconvolve(waveforms, 0.0, upsampling, pz)
 
@@ -1524,7 +1629,6 @@ def __call__(
         else:
             nn_waveforms = self.clip(waveforms / neighbour_sum_clipping)
 
-        # FIXME near-duplicate of neighbour peak sum for now
         neighbors = self.subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse
         peak_index = neighbor_average_maximum(
             nn_waveforms,
@@ -1542,6 +1646,13 @@ def __call__(
             self.sampling_rate_ghz[tel_id] * upsampling,
         )
 
+        if leading_edge_timing:
+            peak_time = adaptive_centroid(
+                d_waveforms,
+                np.round(peak_index - pz2d).astype(int),
+                leading_edge_rel_descend_limit,
+            )
+
         if gain != 0:
             charge /= gain
 
diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 2053623bc7a..9687af67211 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -657,7 +657,7 @@ def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
     assert dl1.is_valid == True
 
     assert_allclose(dl1.image, true_charge, rtol=0.15)
-    assert_allclose(dl1.peak_time, true_time, atol=2.5)
+    assert_allclose(dl1.peak_time, true_time, atol=1.0)
 
     # Test on prod5 simulations
     with EventSource(prod5_gamma_simtel_path) as source:
@@ -678,7 +678,7 @@ def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
                 dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
                 assert dl1.is_valid == True
 
-                bright_pixels = (true_charge > 30) * ~broken_pixels
+                bright_pixels = (true_charge > 30) & ~broken_pixels
                 assert_allclose(
                     dl1.image[bright_pixels], true_charge[bright_pixels], rtol=0.35
                 )

From 28bc35edc3230a0c44a0b6b7deebab5176abeba0 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Fri, 6 Jan 2023 14:53:45 +0100
Subject: [PATCH 22/35] select only trigger telescopes

---
 ctapipe/image/tests/test_extractor.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 9687af67211..9cb6d1feec3 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -640,7 +640,6 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray):
     expected = np.average([29, 30, 31], weights=[5, 10, 3])
     assert np.allclose(dl1.peak_time[bright_pixels], expected / sample_rate)
 
-
 def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
     # Test on toy model
     (
@@ -663,11 +662,13 @@ def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
     with EventSource(prod5_gamma_simtel_path) as source:
         subarray = source.subarray
         extractor = FlashCamExtractor(subarray)
-        for event in source:
-            for tel_id in subarray.get_tel_ids_for_type("MST_MST_FlashCam"):
+        
+        def is_flashcam(tel_id):
+            return subarray.tel[tel_id].camera.name == "FlashCam"
+
+        for event in source:    
+            for tel_id in filter(is_flashcam, event.trigger.tels_with_trigger):
                 true_charge = event.simulation.tel[tel_id].true_image
-                if true_charge is None:
-                    continue  # telescope did not trigger
 
                 waveforms = event.r1.tel[tel_id].waveform
                 selected_gain_channel = np.zeros(waveforms.shape[0], dtype=np.int64)

From 2beadf6f576e702a246722f8ede5ea7fee7b8913 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Sat, 14 Jan 2023 13:42:26 +0100
Subject: [PATCH 23/35] implemented some of the suggested changes

---
 ctapipe/image/extractor.py            | 77 ++++++++++++++++-----------
 ctapipe/image/tests/test_extractor.py |  5 +-
 2 files changed, 48 insertions(+), 34 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index b4f5101f3d0..b45388d5c4e 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -48,7 +48,7 @@
 from .invalid_pixels import InvalidPixelHandler
 from .morphology import brightest_island, number_of_islands
 from .timing import timing_parameters
-
+import astropy.units as u
 
 @guvectorize(
     [
@@ -1292,7 +1292,6 @@ def __call__(
         )
 
 
-@lru_cache
 def deconvolution_parameters(
     camera: CameraDescription,
     upsampling: int,
@@ -1344,32 +1343,31 @@ def deconvolution_parameters(
         raise ValueError(f"window_width must be > 0, got {window_width}")
 
     ref_pulse_shapes = camera.readout.reference_pulse_shape
-    ref_sample_width_nsec = camera.readout.reference_pulse_sample_width.to_value("ns")
-    camera_sample_width_nsec = 1.0 / camera.readout.sampling_rate.to_value("GHz")
+    ref_sample_width_nsec = camera.readout.reference_pulse_sample_width.to_value(u.ns)
+    camera_sample_width_nsec = 1.0 / camera.readout.sampling_rate.to_value(u.GHz)
 
     if camera_sample_width_nsec < ref_sample_width_nsec:
         raise ValueError(
             f"Reference pulse sampling time (got {ref_sample_width_nsec} ns) must be equal to or shorter than the camera sampling time (got {camera_sample_width_nsec} ns); need a reference single p.e. pulse shape with finer sampling!"
         )
-    avg_step = int(camera_sample_width_nsec / ref_sample_width_nsec + 0.5)
+    avg_step = int(round(camera_sample_width_nsec / ref_sample_width_nsec))
 
-    pzs = []  # avg. pole-zero deconvolution parameters
+    pole_zeros = []  # avg. pole-zero deconvolution parameters
     for ref_pulse_shape in ref_pulse_shapes:
         phase_pzs = []
         for phase in range(avg_step):
             x = ref_pulse_shape[phase::avg_step]
             i_min = np.argmin(np.diff(x)) + 1
-            if i_min < x.size and x[i_min] != 0:
-                phase_pzs.append(x[i_min + 1] / x[i_min])
+            phase_pzs.append(x[i_min + 1] / x[i_min])
 
         if len(phase_pzs) == 0:
             raise ValueError(
                 "ref_pulse_shape is malformed - cannot find deconvolution parameter"
             )
-        pzs.append(np.mean(phase_pzs))
+        pole_zeros.append(np.mean(phase_pzs))
 
     gains, shifts, pz2d_shifts = [], [], []  # avg. gains and timing shifts
-    for pz, ref_pulse_shape in zip(pzs, ref_pulse_shapes):
+    for pz, ref_pulse_shape in zip(pole_zeros, ref_pulse_shapes):
         if time_profile_pdf:  # convolve ref_pulse_shape with time profile PDF
             t = (
                 np.arange(ref_pulse_shape.size) - ref_pulse_shape.size / 2
@@ -1385,7 +1383,7 @@ def deconvolution_parameters(
             x = np.atleast_2d(ref_pulse_shape[phase::avg_step])
             y = deconvolve(x, 0.0, upsampling, pz)[0]
             i_max = y.argmax()
-            start = i_max - window_shift  # TODO should use extract_around_peak here
+            start = i_max - window_shift
             stop = start + window_width
             if start >= 0 and stop <= y.size:
                 if leading_edge_timing:
@@ -1412,7 +1410,7 @@ def deconvolution_parameters(
         shifts.append(np.mean(phase_shifts))
         pz2d_shifts.append(np.mean(phase_pz2d_shifts))
 
-    return pzs, gains, shifts, pz2d_shifts
+    return pole_zeros, gains, shifts, pz2d_shifts
 
 
 def deconvolve(
@@ -1468,7 +1466,7 @@ def deconvolve(
     nopython=True,
     cache=True,
 )
-def adaptive_centroid(waveforms, peak_index, rel_descend_limit, cog):
+def adaptive_centroid(waveforms, peak_index, rel_descend_limit, centroids):
     """
     Calculates the pulse centroid for all samples down to rel_descend_limit * peak_amplitude.
 
@@ -1484,22 +1482,24 @@ def adaptive_centroid(waveforms, peak_index, rel_descend_limit, cog):
         Peak index for each pixel.
     rel_descend_limit : ndarray or float
         Fraction of the peak value down to which samples are accumulated in the centroid calculation.
-    cog : ndarray
+    centroids : ndarray
         Return argument for ufunc (ignore)
         Returns the peak centroid in units "samples"
 
     Returns
     -------
-    cog : ndarray
+    centroids : ndarray
         Peak centroid in units "samples"
     """
-    cog[0] = peak_index  # preload in case of errors
+    centroids[0] = peak_index  # preload in case of errors
 
-    nsamples = waveforms.size
-    if nsamples == 0:
+    n_samples = waveforms.size
+    if n_samples == 0:
         return
 
-    peak_index = max(0, min(peak_index, nsamples - 1))
+    if (peak_index > (n_samples - 1)) or (peak_index < 0):
+        raise ValueError("peak_index must be within the waveform limits")
+
     peak_amplitude = waveforms[peak_index]
     if peak_amplitude < 0.0:
         return
@@ -1508,6 +1508,7 @@ def adaptive_centroid(waveforms, peak_index, rel_descend_limit, cog):
 
     sum = float64(0.0)
     jsum = float64(0.0)
+
     j = peak_index
     while j >= 0 and waveforms[j] > descend_limit:
         sum += waveforms[j]
@@ -1517,7 +1518,7 @@ def adaptive_centroid(waveforms, peak_index, rel_descend_limit, cog):
             descend_limit = rel_descend_limit * peak_amplitude
 
     j = peak_index + 1
-    while j < nsamples and waveforms[j] > descend_limit:
+    while j < n_samples and waveforms[j] > descend_limit:
         sum += waveforms[j]
         jsum += j * waveforms[j]
         j += 1
@@ -1525,10 +1526,26 @@ def adaptive_centroid(waveforms, peak_index, rel_descend_limit, cog):
             descend_limit = rel_descend_limit * peak_amplitude
 
     if sum != 0.0:
-        cog[0] = jsum / sum
+        centroids[0] = jsum / sum
 
 
 class FlashCamExtractor(ImageExtractor):
+    """
+
+    The waveforms are first upsampled to achieve one nanosecond sampling (as a default, for the FlashCam). 
+    A pole-zero deconvolution [1] is then performed to the waveforms to recover the original impulse or narrow 
+    the resulting pulse due to convolution with detector response. The modified waveform is integrated in a 
+    window around a peak, which is defined by the neighbors of the pixel. The waveforms are clipped in
+    order to reduce the contribution from the afterpulses in the neighbor sum. If leading_edge_timing is
+    set to True, the so-called leading edge time is found (with the adaptive_centroid function) instead of the peak time.
+    
+    This extractor has been optimized for the FlashCam [2].
+    
+    [1] Smith, S. W. 1997, The Scientist and Engineer’s Guide to Digital Signal Processing (California Technical Publishing)
+    [2] FlashCam: a novel Cherenkov telescope camera with continuous signal digitization. CTA Consortium. 
+    A. Gadola (Zurich U.) et al. DOI: 10.1088/1748-0221/10/01/C01014. Published in: JINST 10 (2015) 01, C01014
+
+    """
     upsampling = IntTelescopeParameter(
         default_value=4, help="Define the upsampling factor for waveforms"
     ).tag(config=True, min=1)
@@ -1616,18 +1633,15 @@ def __call__(
         leading_edge_timing = self.leading_edge_timing.tel[tel_id]
         leading_edge_rel_descend_limit = self.leading_edge_rel_descend_limit.tel[tel_id]
 
-        pzs, gains, shifts, pz2ds = self.deconvolution_pars[tel_id]
-        pz, gain, shift, pz2d = pzs[0], gains[0], shifts[0], pz2ds[0]
-
-        if leading_edge_timing:
-            d_waveforms = deconvolve(waveforms, 0.0, upsampling, 1)
+        pole_zeros, gains, shifts, pz2ds = self.deconvolution_pars[tel_id]
+        pz, gain, shift, pz2d = pole_zeros[0], gains[0], shifts[0], pz2ds[0]
 
-        waveforms = deconvolve(waveforms, 0.0, upsampling, pz)
+        t_waveforms = deconvolve(waveforms, 0.0, upsampling, pz)
 
         if neighbour_sum_clipping == 0.0 or np.isinf(neighbour_sum_clipping):
-            nn_waveforms = waveforms
+            nn_waveforms = t_waveforms
         else:
-            nn_waveforms = self.clip(waveforms / neighbour_sum_clipping)
+            nn_waveforms = self.clip(t_waveforms / neighbour_sum_clipping)
 
         neighbors = self.subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse
         peak_index = neighbor_average_maximum(
@@ -1639,7 +1653,7 @@ def __call__(
         )
 
         charge, peak_time = extract_around_peak(
-            waveforms,
+            t_waveforms,
             peak_index,
             integration_window_width,
             integration_window_shift,
@@ -1647,9 +1661,10 @@ def __call__(
         )
 
         if leading_edge_timing:
+            d_waveforms = deconvolve(waveforms, 0.0, upsampling, 1)
             peak_time = adaptive_centroid(
                 d_waveforms,
-                np.round(peak_index - pz2d).astype(int),
+                max(0, min(np.round(peak_index - pz2d).astype(int))),
                 leading_edge_rel_descend_limit,
             )
 
diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 9cb6d1feec3..3d28562abf6 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -653,10 +653,9 @@ def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
     extractor = FlashCamExtractor(subarray=subarray)
     broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
     dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
-    assert dl1.is_valid == True
-
     assert_allclose(dl1.image, true_charge, rtol=0.15)
     assert_allclose(dl1.peak_time, true_time, atol=1.0)
+    assert dl1.is_valid == True
 
     # Test on prod5 simulations
     with EventSource(prod5_gamma_simtel_path) as source:
@@ -679,7 +678,7 @@ def is_flashcam(tel_id):
                 dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
                 assert dl1.is_valid == True
 
-                bright_pixels = (true_charge > 30) & ~broken_pixels
+                bright_pixels = (true_charge > 30) & (true_charge < 3000) & ~broken_pixels
                 assert_allclose(
                     dl1.image[bright_pixels], true_charge[bright_pixels], rtol=0.35
                 )

From 3f5e287d7dd24773bfa1b812b9456287a70fed81 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Wed, 15 Mar 2023 18:23:39 +0100
Subject: [PATCH 24/35] code style

---
 ctapipe/image/extractor.py | 25 +++++++++++++------------
 1 file changed, 13 insertions(+), 12 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index b45388d5c4e..624d9069116 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -25,14 +25,10 @@
 from functools import lru_cache
 from typing import Callable, List, Optional, Tuple
 
+import astropy.units as u
 import numpy as np
 import numpy.typing as npt
 import scipy.stats
-from numba import float32, float64, guvectorize, int64, njit, prange
-from scipy.ndimage import convolve1d
-from scipy.signal import filtfilt
-from traitlets import Bool, Int
-
 from ctapipe.containers import DL1CameraContainer
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import (
@@ -42,13 +38,17 @@
     IntTelescopeParameter,
 )
 from ctapipe.instrument import CameraDescription
+from numba import float32, float64, guvectorize, int64, njit, prange
+from scipy.ndimage import convolve1d
+from scipy.signal import filtfilt
+from traitlets import Bool, Int
 
 from .cleaning import tailcuts_clean
 from .hillas import camera_to_shower_coordinates, hillas_parameters
 from .invalid_pixels import InvalidPixelHandler
 from .morphology import brightest_island, number_of_islands
 from .timing import timing_parameters
-import astropy.units as u
+
 
 @guvectorize(
     [
@@ -1532,20 +1532,21 @@ def adaptive_centroid(waveforms, peak_index, rel_descend_limit, centroids):
 class FlashCamExtractor(ImageExtractor):
     """
 
-    The waveforms are first upsampled to achieve one nanosecond sampling (as a default, for the FlashCam). 
-    A pole-zero deconvolution [1] is then performed to the waveforms to recover the original impulse or narrow 
-    the resulting pulse due to convolution with detector response. The modified waveform is integrated in a 
+    The waveforms are first upsampled to achieve one nanosecond sampling (as a default, for the FlashCam).
+    A pole-zero deconvolution [1] is then performed to the waveforms to recover the original impulse or narrow
+    the resulting pulse due to convolution with detector response. The modified waveform is integrated in a
     window around a peak, which is defined by the neighbors of the pixel. The waveforms are clipped in
     order to reduce the contribution from the afterpulses in the neighbor sum. If leading_edge_timing is
     set to True, the so-called leading edge time is found (with the adaptive_centroid function) instead of the peak time.
-    
+
     This extractor has been optimized for the FlashCam [2].
-    
+
     [1] Smith, S. W. 1997, The Scientist and Engineer’s Guide to Digital Signal Processing (California Technical Publishing)
-    [2] FlashCam: a novel Cherenkov telescope camera with continuous signal digitization. CTA Consortium. 
+    [2] FlashCam: a novel Cherenkov telescope camera with continuous signal digitization. CTA Consortium.
     A. Gadola (Zurich U.) et al. DOI: 10.1088/1748-0221/10/01/C01014. Published in: JINST 10 (2015) 01, C01014
 
     """
+
     upsampling = IntTelescopeParameter(
         default_value=4, help="Define the upsampling factor for waveforms"
     ).tag(config=True, min=1)

From 137c7abfbd57c4b21b6b832a56be64df475d5d5c Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Wed, 15 Mar 2023 18:37:39 +0100
Subject: [PATCH 25/35] code style

---
 ctapipe/image/extractor.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 624d9069116..faf0a4f85b5 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -29,6 +29,11 @@
 import numpy as np
 import numpy.typing as npt
 import scipy.stats
+from numba import float32, float64, guvectorize, int64, njit, prange
+from scipy.ndimage import convolve1d
+from scipy.signal import filtfilt
+from traitlets import Bool, Int
+
 from ctapipe.containers import DL1CameraContainer
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import (
@@ -38,10 +43,6 @@
     IntTelescopeParameter,
 )
 from ctapipe.instrument import CameraDescription
-from numba import float32, float64, guvectorize, int64, njit, prange
-from scipy.ndimage import convolve1d
-from scipy.signal import filtfilt
-from traitlets import Bool, Int
 
 from .cleaning import tailcuts_clean
 from .hillas import camera_to_shower_coordinates, hillas_parameters

From b94836b0494ccd3c75095b438c9628ffd2bd1b78 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Wed, 15 Mar 2023 18:48:52 +0100
Subject: [PATCH 26/35] test code style

---
 ctapipe/image/tests/test_extractor.py | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 3d28562abf6..5690ba9e0e2 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -24,11 +24,7 @@
     subtract_baseline,
 )
 from ctapipe.image.toymodel import SkewedGaussian, WaveformModel, obtain_time_image
-from ctapipe.instrument import SubarrayDescription, TelescopeDescription
-from numpy.testing import assert_allclose, assert_equal
-from scipy.stats import norm
-from traitlets.config.loader import Config
-from traitlets.traitlets import TraitError
+from ctapipe.instrument import SubarrayDescription
 from ctapipe.io import EventSource
 
 extractors = non_abstract_children(ImageExtractor)
@@ -640,6 +636,7 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray):
     expected = np.average([29, 30, 31], weights=[5, 10, 3])
     assert np.allclose(dl1.peak_time[bright_pixels], expected / sample_rate)
 
+
 def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
     # Test on toy model
     (
@@ -661,11 +658,11 @@ def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
     with EventSource(prod5_gamma_simtel_path) as source:
         subarray = source.subarray
         extractor = FlashCamExtractor(subarray)
-        
+
         def is_flashcam(tel_id):
             return subarray.tel[tel_id].camera.name == "FlashCam"
 
-        for event in source:    
+        for event in source:
             for tel_id in filter(is_flashcam, event.trigger.tels_with_trigger):
                 true_charge = event.simulation.tel[tel_id].true_image
 
@@ -678,7 +675,9 @@ def is_flashcam(tel_id):
                 dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
                 assert dl1.is_valid == True
 
-                bright_pixels = (true_charge > 30) & (true_charge < 3000) & ~broken_pixels
+                bright_pixels = (
+                    (true_charge > 30) & (true_charge < 3000) & ~broken_pixels
+                )
                 assert_allclose(
                     dl1.image[bright_pixels], true_charge[bright_pixels], rtol=0.35
                 )

From e5d8c8699071f7044aa7e9260e9189c2ff950435 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Wed, 15 Mar 2023 19:43:11 +0100
Subject: [PATCH 27/35] add missing config

---
 ctapipe/image/extractor.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index faf0a4f85b5..64988013837 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1586,7 +1586,7 @@ class FlashCamExtractor(ImageExtractor):
 
     leading_edge_timing = BoolTelescopeParameter(
         default_value=True, help="Calculate leading edge time instead of peak time"
-    ).tag(config=False)
+    ).tag(config=True)
 
     leading_edge_rel_descend_limit = FloatTelescopeParameter(
         default_value=0.05,

From f6aa071dac9bcae37ff1529a8bfb12e06fb9cc7d Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Wed, 15 Mar 2023 21:40:28 +0100
Subject: [PATCH 28/35] fix codacy warnings

---
 ctapipe/image/extractor.py            | 23 ++++++++++++++---------
 ctapipe/image/tests/test_extractor.py | 12 ++++++------
 2 files changed, 20 insertions(+), 15 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 64988013837..9ecf8f67aca 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1320,7 +1320,8 @@ def deconvolution_parameters(
     leading_edge_timing : bool
         Whether time calculation will be done on the leading edge.
     leading_edge_rel_descend_limit : float
-        If leading edge timing is used, the fraction of the peak value down to which samples are accumulated in the centroid calculation.
+        If leading edge timing is used, the fraction of the peak value down to which samples are accumulated in the
+        centroid calculation.
     time_profile_pdf : callable or None
         PDF of the assumed effective Cherenkov time profile to assume when
         calculating the gain loss; takes nanoseconds as arguments and returns
@@ -1349,7 +1350,9 @@ def deconvolution_parameters(
 
     if camera_sample_width_nsec < ref_sample_width_nsec:
         raise ValueError(
-            f"Reference pulse sampling time (got {ref_sample_width_nsec} ns) must be equal to or shorter than the camera sampling time (got {camera_sample_width_nsec} ns); need a reference single p.e. pulse shape with finer sampling!"
+            f"Reference pulse sampling time (got {ref_sample_width_nsec} ns) must be equal to or shorter than the "
+            f"camera sampling time (got {camera_sample_width_nsec} ns); need a reference single p.e. pulse shape with "
+            "finer sampling!"
         )
     avg_step = int(round(camera_sample_width_nsec / ref_sample_width_nsec))
 
@@ -1507,12 +1510,12 @@ def adaptive_centroid(waveforms, peak_index, rel_descend_limit, centroids):
 
     descend_limit = rel_descend_limit * peak_amplitude
 
-    sum = float64(0.0)
+    sum_ = float64(0.0)
     jsum = float64(0.0)
 
     j = peak_index
     while j >= 0 and waveforms[j] > descend_limit:
-        sum += waveforms[j]
+        sum_ += waveforms[j]
         jsum += j * waveforms[j]
         j -= 1
         if waveforms[j] > peak_amplitude:
@@ -1520,14 +1523,14 @@ def adaptive_centroid(waveforms, peak_index, rel_descend_limit, centroids):
 
     j = peak_index + 1
     while j < n_samples and waveforms[j] > descend_limit:
-        sum += waveforms[j]
+        sum_ += waveforms[j]
         jsum += j * waveforms[j]
         j += 1
         if waveforms[j] > peak_amplitude:
             descend_limit = rel_descend_limit * peak_amplitude
 
-    if sum != 0.0:
-        centroids[0] = jsum / sum
+    if sum_ != 0.0:
+        centroids[0] = jsum / sum_
 
 
 class FlashCamExtractor(ImageExtractor):
@@ -1538,11 +1541,13 @@ class FlashCamExtractor(ImageExtractor):
     the resulting pulse due to convolution with detector response. The modified waveform is integrated in a
     window around a peak, which is defined by the neighbors of the pixel. The waveforms are clipped in
     order to reduce the contribution from the afterpulses in the neighbor sum. If leading_edge_timing is
-    set to True, the so-called leading edge time is found (with the adaptive_centroid function) instead of the peak time.
+    set to True, the so-called leading edge time is found (with the adaptive_centroid function) instead of the peak
+    time.
 
     This extractor has been optimized for the FlashCam [2].
 
-    [1] Smith, S. W. 1997, The Scientist and Engineer’s Guide to Digital Signal Processing (California Technical Publishing)
+    [1] Smith, S. W. 1997, The Scientist and Engineer’s Guide to Digital Signal Processing (California Technical
+    Publishing)
     [2] FlashCam: a novel Cherenkov telescope camera with continuous signal digitization. CTA Consortium.
     A. Gadola (Zurich U.) et al. DOI: 10.1088/1748-0221/10/01/C01014. Published in: JINST 10 (2015) 01, C01014
 
diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 5690ba9e0e2..2f0ac2791e9 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -72,7 +72,7 @@ def subarray_1_LST(prod3_lst):
 
 
 @pytest.fixture(scope="module")
-def subarray_1_MST_FC(prod5_mst_flashcam):
+def subarray_mst_fc(prod5_mst_flashcam):
     subarray = SubarrayDescription(
         "One MST with FlashCam",
         tel_positions={1: np.zeros(3) * u.m},
@@ -106,8 +106,8 @@ def toymodel(subarray):
 
 
 @pytest.fixture(scope="module")
-def toymodel_1_MST_FC(subarray_1_MST_FC):
-    return get_test_toymodel(subarray_1_MST_FC)
+def toymodel_mst_fc(subarray_mst_fc: object) -> object:
+    return get_test_toymodel(subarray_mst_fc)
 
 
 def test_extract_around_peak(toymodel):
@@ -637,7 +637,7 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray):
     assert np.allclose(dl1.peak_time[bright_pixels], expected / sample_rate)
 
 
-def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
+def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path):
     # Test on toy model
     (
         waveforms,
@@ -646,7 +646,7 @@ def test_flashcam_extractor(toymodel_1_MST_FC, prod5_gamma_simtel_path):
         selected_gain_channel,
         true_charge,
         true_time,
-    ) = toymodel_1_MST_FC
+    ) = toymodel_mst_fc
     extractor = FlashCamExtractor(subarray=subarray)
     broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
     dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
@@ -676,7 +676,7 @@ def is_flashcam(tel_id):
                 assert dl1.is_valid == True
 
                 bright_pixels = (
-                    (true_charge > 30) & (true_charge < 3000) & ~broken_pixels
+                    (true_charge > 30) & (true_charge < 3000) & (~broken_pixels)
                 )
                 assert_allclose(
                     dl1.image[bright_pixels], true_charge[bright_pixels], rtol=0.35

From 657db5166de145e842a1f6cb49d6fef1c3c0d3f9 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Wed, 15 Mar 2023 21:58:41 +0100
Subject: [PATCH 29/35] fix codacy warnings

---
 docs/changes/2188.feature.rst | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 docs/changes/2188.feature.rst

diff --git a/docs/changes/2188.feature.rst b/docs/changes/2188.feature.rst
new file mode 100644
index 00000000000..c5e1cf5f977
--- /dev/null
+++ b/docs/changes/2188.feature.rst
@@ -0,0 +1 @@
+Add signal extraction algorithm for the FlashCam
\ No newline at end of file

From 8e48d5f5c8a4437fb9940e848585e731fb446019 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Tue, 21 Mar 2023 16:58:20 +0100
Subject: [PATCH 30/35] substituted filtfilt upsampling

---
 ctapipe/image/extractor.py | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 9ecf8f67aca..91ea9f86bfc 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -31,7 +31,6 @@
 import scipy.stats
 from numba import float32, float64, guvectorize, int64, njit, prange
 from scipy.ndimage import convolve1d
-from scipy.signal import filtfilt
 from traitlets import Bool, Int
 
 from ctapipe.containers import DL1CameraContainer
@@ -703,7 +702,6 @@ def _calculate_correction(self, tel_id):
         n_channels = len(readout.reference_pulse_shape)
         correction = np.ones(n_channels, dtype=np.float64)
         for ichannel, pulse_shape in enumerate(readout.reference_pulse_shape):
-
             # apply the same method as sliding window to find the highest sum
             cwf = np.cumsum(pulse_shape)
             # add zero at the begining so it is easier to substract the two arrays later
@@ -1452,12 +1450,18 @@ def deconvolve(
     deconvolved_waveforms = np.atleast_2d(waveforms) - np.atleast_2d(baselines).T
     deconvolved_waveforms[:, 1:] -= pole_zero * deconvolved_waveforms[:, :-1]
     deconvolved_waveforms[:, 0] = 0
+
+    def filtfilt_custom(signal, filt):
+        forward = convolve1d(signal, filt, axis=-1, mode="nearest")
+        backward = convolve1d(forward[::-1], filt, axis=-1, mode="nearest")
+        return backward[::-1]
+
     if upsampling > 1:
-        return filtfilt(
-            np.ones(upsampling),
-            upsampling,
-            np.repeat(deconvolved_waveforms, upsampling, axis=-1),
-        )
+        filt = np.ones(upsampling)
+        filt_weighted = filt / upsampling
+        signal = np.repeat(deconvolved_waveforms, upsampling, axis=-1)
+        return filtfilt_custom(signal, filt_weighted)
+
     return deconvolved_waveforms
 
 

From e4408da19fa21df3b8f2aadbc85b731829f599b3 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Wed, 22 Mar 2023 12:34:42 +0100
Subject: [PATCH 31/35] added private fast filtfilt function and changed doc

---
 ctapipe/image/extractor.py | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 91ea9f86bfc..eb4a5f515e3 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1415,6 +1415,16 @@ def deconvolution_parameters(
     return pole_zeros, gains, shifts, pz2d_shifts
 
 
+def __filtfilt_fast(signal, filt):
+    """
+    Apply a linear filter forward and backward to a signal, based on scipy.signal.filtfilt.
+    filtfilt has some speed issues (https://github.com/scipy/scipy/issues/17080)
+    """
+    forward = convolve1d(signal, filt, axis=-1, mode="nearest")
+    backward = convolve1d(forward[::-1], filt, axis=-1, mode="nearest")
+    return backward[::-1]
+
+
 def deconvolve(
     waveforms: npt.ArrayLike,
     baselines: npt.ArrayLike,
@@ -1451,16 +1461,11 @@ def deconvolve(
     deconvolved_waveforms[:, 1:] -= pole_zero * deconvolved_waveforms[:, :-1]
     deconvolved_waveforms[:, 0] = 0
 
-    def filtfilt_custom(signal, filt):
-        forward = convolve1d(signal, filt, axis=-1, mode="nearest")
-        backward = convolve1d(forward[::-1], filt, axis=-1, mode="nearest")
-        return backward[::-1]
-
     if upsampling > 1:
         filt = np.ones(upsampling)
         filt_weighted = filt / upsampling
         signal = np.repeat(deconvolved_waveforms, upsampling, axis=-1)
-        return filtfilt_custom(signal, filt_weighted)
+        return __filtfilt_fast(signal, filt_weighted)
 
     return deconvolved_waveforms
 
@@ -1539,6 +1544,7 @@ def adaptive_centroid(waveforms, peak_index, rel_descend_limit, centroids):
 
 class FlashCamExtractor(ImageExtractor):
     """
+    Image extractor applying signal preprocessing for FlashCam
 
     The waveforms are first upsampled to achieve one nanosecond sampling (as a default, for the FlashCam).
     A pole-zero deconvolution [1] is then performed to the waveforms to recover the original impulse or narrow

From 75f7b9ec69adf745fd4cc94cb2dc07c0b3e56b14 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Mon, 27 Mar 2023 12:19:40 +0200
Subject: [PATCH 32/35] added test for adaptive_centroid

---
 ctapipe/image/tests/test_extractor.py | 42 +++++++++++++++++++++++++--
 1 file changed, 40 insertions(+), 2 deletions(-)

diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 2f0ac2791e9..0f8c1160d6c 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -17,6 +17,7 @@
     NeighborPeakWindowSum,
     SlidingWindowMaxSum,
     TwoPassWindowSum,
+    adaptive_centroid,
     extract_around_peak,
     extract_sliding_window,
     integration_correction,
@@ -408,7 +409,6 @@ def test_neighbor_peak_window_sum_local_weight(toymodel):
 
 
 def test_Two_pass_window_sum_no_noise(subarray_1_LST):
-
     rng = np.random.default_rng(0)
 
     subarray = subarray_1_LST
@@ -585,7 +585,6 @@ def test_extractor_tel_param(toymodel):
 
 @pytest.mark.parametrize("Extractor", non_abstract_children(ImageExtractor))
 def test_dtype(Extractor, subarray):
-
     tel_id = 1
     n_pixels = subarray.tel[tel_id].camera.geometry.n_pixels
     selected_gain_channel = np.zeros(n_pixels, dtype=int)
@@ -637,6 +636,45 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray):
     assert np.allclose(dl1.peak_time[bright_pixels], expected / sample_rate)
 
 
+def test_adaptive_centroid(toymodel_mst_fc):
+    (
+        waveforms,
+        subarray,
+        tel_id,
+        selected_gain_channel,
+        true_charge,
+        true_time,
+    ) = toymodel_mst_fc
+
+    neighbors = subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse
+    broken_pixels = np.zeros(waveforms.shape[0], dtype=bool)
+
+    trig_time = np.argmax(waveforms, axis=-1)
+    peak_time = adaptive_centroid(
+        waveforms,
+        trig_time,
+        1,
+    )
+
+    assert (peak_time == trig_time).all()
+
+    waveforms = waveforms[np.min(waveforms, axis=-1) > 0.0]
+    peak_pos = neighbor_average_maximum(
+        waveforms,
+        neighbors_indices=neighbors.indices,
+        neighbors_indptr=neighbors.indptr,
+        local_weight=0,
+        broken_pixels=broken_pixels,
+    )
+
+    peak_time = adaptive_centroid(
+        waveforms,
+        peak_pos,
+        0.0,
+    )
+    assert (peak_pos == peak_time).all()
+
+
 def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path):
     # Test on toy model
     (

From 5c5eca8ba660b75ce13cf3f8a2eb991e2b97600a Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Mon, 27 Mar 2023 17:44:53 +0200
Subject: [PATCH 33/35] added tests for deconvolve and upsampling

---
 ctapipe/image/tests/test_extractor.py | 48 +++++++++++++++++++++++++++
 1 file changed, 48 insertions(+)

diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index 0f8c1160d6c..d0ba1f1acb9 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -4,6 +4,7 @@
 import numpy as np
 import pytest
 from numpy.testing import assert_allclose, assert_equal
+from scipy.signal import filtfilt
 from scipy.stats import norm
 from traitlets.config.loader import Config
 from traitlets.traitlets import TraitError
@@ -17,7 +18,9 @@
     NeighborPeakWindowSum,
     SlidingWindowMaxSum,
     TwoPassWindowSum,
+    __filtfilt_fast,
     adaptive_centroid,
+    deconvolve,
     extract_around_peak,
     extract_sliding_window,
     integration_correction,
@@ -675,6 +678,51 @@ def test_adaptive_centroid(toymodel_mst_fc):
     assert (peak_pos == peak_time).all()
 
 
+def test_deconvolve(toymodel_mst_fc):
+    (
+        waveforms,
+        subarray,
+        tel_id,
+        selected_gain_channel,
+        true_charge,
+        true_time,
+    ) = toymodel_mst_fc
+
+    deconvolved_waveforms_0 = deconvolve(waveforms, 0, 0, 0.0)
+
+    assert (deconvolved_waveforms_0[:, 1:] == waveforms[:, 1:]).all()
+
+    deconvolved_waveforms_1 = deconvolve(waveforms, 0, 0, 1.0)
+
+    assert (deconvolved_waveforms_1[:, 1:] == np.diff(waveforms, axis=-1)).all()
+
+
+def test_upsampling(toymodel_mst_fc):
+    (
+        waveforms,
+        subarray,
+        tel_id,
+        selected_gain_channel,
+        true_charge,
+        true_time,
+    ) = toymodel_mst_fc
+
+    upsampling = 4
+    filt = np.ones(upsampling)
+    filt_weighted = filt / upsampling
+    signal = np.repeat(waveforms, upsampling, axis=-1)
+    up_waveforms = __filtfilt_fast(signal, filt_weighted)
+
+    assert (
+        up_waveforms
+        - filtfilt(
+            np.ones(upsampling),
+            upsampling,
+            np.repeat(waveforms, upsampling, axis=-1),
+        )
+    ).all() < 1e-10
+
+
 def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path):
     # Test on toy model
     (

From 2e0e69a23bc958e6bcf8fb8be6ef3fa8a5077707 Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Mon, 27 Mar 2023 19:43:01 +0200
Subject: [PATCH 34/35] modified test for upsampling and solved bug in
 upsampling function

---
 ctapipe/image/extractor.py            |  4 +--
 ctapipe/image/tests/test_extractor.py | 46 +++++++++++++++++++--------
 2 files changed, 34 insertions(+), 16 deletions(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index eb4a5f515e3..9813ba01a4b 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1421,8 +1421,8 @@ def __filtfilt_fast(signal, filt):
     filtfilt has some speed issues (https://github.com/scipy/scipy/issues/17080)
     """
     forward = convolve1d(signal, filt, axis=-1, mode="nearest")
-    backward = convolve1d(forward[::-1], filt, axis=-1, mode="nearest")
-    return backward[::-1]
+    backward = convolve1d(forward[:, ::-1], filt, axis=-1, mode="nearest")
+    return backward[:, ::-1]
 
 
 def deconvolve(
diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py
index d0ba1f1acb9..768329bd07c 100644
--- a/ctapipe/image/tests/test_extractor.py
+++ b/ctapipe/image/tests/test_extractor.py
@@ -706,21 +706,39 @@ def test_upsampling(toymodel_mst_fc):
         true_charge,
         true_time,
     ) = toymodel_mst_fc
+    upsampling_even = 4
+    upsampling_odd = 3
+    filt_even = np.ones(upsampling_even)
+    filt_weighted_even = filt_even / upsampling_even
+    signal_even = np.repeat(waveforms, upsampling_even, axis=-1)
+    up_waveforms_even = __filtfilt_fast(signal_even, filt_weighted_even)
+
+    np.testing.assert_allclose(
+        up_waveforms_even,
+        filtfilt(
+            np.ones(upsampling_even),
+            upsampling_even,
+            np.repeat(waveforms, upsampling_even, axis=-1),
+        ),
+        rtol=1e-4,
+        atol=1e-4,
+    )
 
-    upsampling = 4
-    filt = np.ones(upsampling)
-    filt_weighted = filt / upsampling
-    signal = np.repeat(waveforms, upsampling, axis=-1)
-    up_waveforms = __filtfilt_fast(signal, filt_weighted)
-
-    assert (
-        up_waveforms
-        - filtfilt(
-            np.ones(upsampling),
-            upsampling,
-            np.repeat(waveforms, upsampling, axis=-1),
-        )
-    ).all() < 1e-10
+    filt_odd = np.ones(upsampling_odd)
+    filt_weighted_odd = filt_odd / upsampling_odd
+    signal_odd = np.repeat(waveforms, upsampling_odd, axis=-1)
+    up_waveforms_odd = __filtfilt_fast(signal_odd, filt_weighted_odd)
+
+    np.testing.assert_allclose(
+        up_waveforms_odd,
+        filtfilt(
+            np.ones(upsampling_odd),
+            upsampling_odd,
+            np.repeat(waveforms, upsampling_odd, axis=-1),
+        ),
+        rtol=1e-4,
+        atol=1e-4,
+    )
 
 
 def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path):

From 3c06ad7b2046bb4892bba1d4a7928d0fb32772da Mon Sep 17 00:00:00 2001
From: nieves <clara.escanuela@mpi-hd.mpg.de>
Date: Tue, 28 Mar 2023 17:02:10 +0200
Subject: [PATCH 35/35] changed order of guvectorize input

---
 ctapipe/image/extractor.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py
index 9813ba01a4b..5893ebf1dd0 100644
--- a/ctapipe/image/extractor.py
+++ b/ctapipe/image/extractor.py
@@ -1472,8 +1472,8 @@ def deconvolve(
 
 @guvectorize(
     [
-        (float64[:], int64, float64, float32[:]),
         (float32[:], int64, float64, float32[:]),
+        (float64[:], int64, float64, float32[:]),
     ],
     "(s),(),()->()",
     nopython=True,