From 040d4ad666ff34a8d06562d5f880a806e5d366af Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 25 Feb 2026 16:27:33 +0000
Subject: [PATCH 01/49] add sipm selection function using median and std

---
 invisible_cities/reco/wfm_functions.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 4ecced9e2..ff463d70f 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -3,6 +3,8 @@
 authors: J.J. Gomez-Cadenas, G. Martinez
 """
 import numpy as np
+from typing import Optional
+from typing  import Callable
 
 from .. core.core_functions import define_window
 from .. calib               import calib_sensors_functions as csf
@@ -129,3 +131,23 @@ def compare_cwf_blr(cwf, pmtblr, event_list, window_size=500):
             DIFF.append(diff)
 
     return np.array(DIFF)
+
+
+def median_std_method(wfs    : np.ndarray,  
+                      nsigma : Optional[float] = 3.) -> np.ndarray:
+    """
+    Computes the median and standard deviation of the time summed SiPM 
+    waveforms and selects the SiPMs that are nsigma over the median.
+
+    Parameters
+    ----------
+    wfs    : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
+    nsigma : Number of standard deviations above the median, default 3.
+    
+    Returns
+    -------
+    Boolean numpy array of shape (n_sipms,) where True indicates that the SiPM is selected.
+    """
+    charges = np.sum(wfs, axis=1)
+    threshold = np.median(charges) + nsigma * np.std(charges)
+    return charges >= threshold
\ No newline at end of file

From eecf819fa275c163050cb37bfdf1f1b4df395313 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 25 Feb 2026 16:31:38 +0000
Subject: [PATCH 02/49] add sipm selection method using a charge threshold

---
 invisible_cities/reco/wfm_functions.py | 20 +++++++++++++++++++-
 1 file changed, 19 insertions(+), 1 deletion(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index ff463d70f..9e0456bc5 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -150,4 +150,22 @@ def median_std_method(wfs    : np.ndarray,
     """
     charges = np.sum(wfs, axis=1)
     threshold = np.median(charges) + nsigma * np.std(charges)
-    return charges >= threshold
\ No newline at end of file
+    return charges >= threshold
+
+
+def charge_threshold_method(wfs       : np.ndarray,
+                            threshold : Optional[float] = 5.) -> np.ndarray:
+    """
+    Selects the SiPMs whose time summed waveforms are above a threshold.
+    
+    Parameters
+    ----------
+    wfs       : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
+    threshold : Charge threshold in PE, default 5.
+    
+    Returns
+    -------
+    Boolean numpy array of shape (n_sipms,) where True indicates that the SiPM is selected.
+    """
+    charges = np.sum(wfs, axis=1)
+    return charges >= threshold

From 38364ffd5eb2a3bd4d7b584731beaddae839b5a2 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 25 Feb 2026 16:32:32 +0000
Subject: [PATCH 03/49] add sipm selection function using top n most energetic
 sipms

---
 invisible_cities/reco/wfm_functions.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 9e0456bc5..62468b6ae 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -169,3 +169,25 @@ def charge_threshold_method(wfs       : np.ndarray,
     """
     charges = np.sum(wfs, axis=1)
     return charges >= threshold
+
+
+def top_n_method(wfs : np.ndarray, 
+                 n   : Optional[int] = 10) -> np.ndarray:
+    """
+    Selects the SiPMs with the top n highest time summed waveforms.
+        
+    Parameters
+    ----------
+    wfs       : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
+    n         : Number of most energeticSiPMs to select, default 10.
+    
+    Returns
+    -------
+    Boolean numpy array of shape (n_sipms,) where True indicates that the SiPM is selected.
+    """
+    charges = np.sum(wfs, axis=1)
+    idx = np.argsort(charges)[-n:]
+
+    selected_ids = np.zeros_like(charges, dtype=bool)
+    selected_ids[idx] = True
+    return selected_ids

From 57955f4dfec57f0311921deb469e76fe32940274 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 25 Feb 2026 16:46:58 +0000
Subject: [PATCH 04/49] add function that kills isolated sipms

---
 invisible_cities/reco/wfm_functions.py | 35 ++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 62468b6ae..625a337b3 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -191,3 +191,38 @@ def top_n_method(wfs : np.ndarray,
     selected_ids = np.zeros_like(charges, dtype=bool)
     selected_ids[idx] = True
     return selected_ids
+
+
+def kill_isolated_sipms(selected_ids        : np.ndarray,
+                        sipm_x              : np.ndarray, 
+                        sipm_y              : np.ndarray, 
+                        proximity_threshold : float) -> np.ndarray:
+    """
+    For the SiPMs that have passed the previous selection, scans through the SiPMs to check if they 
+    have neighbouring SiPMs - i.e., within the proximity_threshold - that have also passed the selection. 
+    If no neighbours are found, the SiPMs are classed as isolated, and are removed.
+
+    Parameters
+    ----------
+    selected_ids        : Boolean array of shape (n_sipms,) indicating which SiPMs passed the previous selection.
+    sipm_x              : 1D array of shape (n_sipms,) containing the x positions of the SiPMs.
+    sipm_y              : 1D array of shape (n_sipms,) containing the y positions of the SiPMs.
+    proximity_threshold : Distance threshold in mm used to identify isolated SiPMs.
+
+    Returns
+    -------
+    selected_ids_no_isolated : Boolean array of shape (n_sipms,) where True indicates that the SiPM is selected.
+    """
+    selected_ids_no_isolated = selected_ids.copy()
+
+    for i in np.where(selected_ids)[0]:
+        x, y = sipm_x[i], sipm_y[i]
+
+        distances = np.sqrt((sipm_x - x)**2 + (sipm_y - y)**2)
+
+        n_neighbors = np.sum((distances < proximity_threshold) & selected_ids)
+
+        if n_neighbors <= 1:
+            selected_ids_no_isolated[i] = False
+            
+    return selected_ids_no_isolated
\ No newline at end of file

From 34ec2ef7eb2d12e26d2a803d7bdce673ba1e9b1e Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 25 Feb 2026 16:50:48 +0000
Subject: [PATCH 05/49] add function that creates circular padding around
 selected sipms

---
 invisible_cities/reco/wfm_functions.py | 31 +++++++++++++++++++++++++-
 1 file changed, 30 insertions(+), 1 deletion(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 625a337b3..399bf1b0f 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -225,4 +225,33 @@ def kill_isolated_sipms(selected_ids        : np.ndarray,
         if n_neighbors <= 1:
             selected_ids_no_isolated[i] = False
             
-    return selected_ids_no_isolated
\ No newline at end of file
+    return selected_ids_no_isolated
+
+
+def apply_circular_padding(selected_ids_no_isolated : np.ndarray,
+                           sipm_x                   : np.ndarray, 
+                           sipm_y                   : np.ndarray, 
+                           padding_radius           : float) -> np.ndarray:
+    """
+    For the SiPMs that pass the previous selection, creates circular padding of radius padding_radius, 
+    selecting all SiPMs within that radius. Stores the union of all selected SiPMs.
+
+    Parameters
+    ----------
+    selected_ids_no_isolated : Boolean array of shape (n_sipms,) indicating which SiPMs passed the previous selection.
+    sipm_x                   : 1D array of shape (n_sipms,) containing the x positions of the SiPMs.
+    sipm_y                   : 1D array of shape (n_sipms,) containing the y positions of the SiPMs.
+    padding_radius           : Distance threshold in mm used to create circular padding around selected SiPMs.
+
+    Returns
+    -------
+    sipm_ids_with_signal : Boolean array of shape (n_sipms,) where True indicates that the SiPM is selected.
+    """
+    sipm_ids_with_signal = np.zeros_like(selected_ids_no_isolated, dtype=bool)
+
+    for i in np.where(selected_ids_no_isolated)[0]:
+        x, y = sipm_x[i], sipm_y[i]
+        distances = np.sqrt((sipm_x - x)**2 + (sipm_y - y)**2)
+        sipm_ids_with_signal |= distances < padding_radius
+
+    return sipm_ids_with_signal
\ No newline at end of file

From a55132e670c91f560fb4195b633befa15aa9f4f2 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 25 Feb 2026 16:51:50 +0000
Subject: [PATCH 06/49] add function that creates the overall sipm masks

---
 invisible_cities/reco/wfm_functions.py | 46 +++++++++++++++++++++++++-
 1 file changed, 45 insertions(+), 1 deletion(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 399bf1b0f..a39ed2520 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -254,4 +254,48 @@ def apply_circular_padding(selected_ids_no_isolated : np.ndarray,
         distances = np.sqrt((sipm_x - x)**2 + (sipm_y - y)**2)
         sipm_ids_with_signal |= distances < padding_radius
 
-    return sipm_ids_with_signal
\ No newline at end of file
+    return sipm_ids_with_signal
+
+
+def make_sipm_selection(wfs                 : np.ndarray, 
+                        selection_func      : Callable, 
+                        selection_kwargs    : dict, 
+                        proximity_threshold : float, 
+                        padding_radius      : float) -> np.ndarray: 
+    """
+    SiPM selection pipeline, applies SiPM cuts based on user input.
+    A first selection of SiPMs is made, isolated SiPMs are removed
+    and padding is added around the SiPMs that are left.
+
+    Parameters
+    ----------
+    wfs                 : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
+    selection_func      : Function returning a boolean mask of energetic SiPMs.
+    selection_kwargs    : Dictionary of arguments passed to selection_func.
+    proximity_threshold : Threshold used to identify isolated SiPMs.
+    padding_radius      : Radial padding added to each SiPM that passes the selections.
+
+    Returns
+    -------
+    sipm_ids_with_signal : Array of shape (n_sipms,) with boolean values indicating which SiPMs are selected.
+    """
+
+    sipm_x = np.array(detector_info.X)
+    sipm_y = np.array(detector_info.Y)
+
+    selected_ids = selection_func(wfs, **selection_kwargs)
+
+    selected_ids_no_isolated = kill_isolated_sipms(
+        selected_ids,
+        sipm_x,
+        sipm_y,
+        proximity_threshold
+    )
+
+    sipm_ids_with_signal = apply_circular_padding(
+        selected_ids_no_isolated,
+        sipm_x,
+        sipm_y,
+        padding_radius
+    )
+    return sipm_ids_with_signal

From f137e230a8000a6c0cfb5fbc5db147ba2f446b50 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 25 Feb 2026 16:54:08 +0000
Subject: [PATCH 07/49] load detector geometry in

---
 invisible_cities/reco/wfm_functions.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index a39ed2520..fb4388966 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -3,12 +3,13 @@
 authors: J.J. Gomez-Cadenas, G. Martinez
 """
 import numpy as np
-from typing import Optional
+from typing  import Optional
 from typing  import Callable
 
 from .. core.core_functions import define_window
 from .. calib               import calib_sensors_functions as csf
 from .. sierpe              import blr
+from .. database            import load_db
 
 def to_adc(wfs, adc_to_pes):
     """
@@ -279,7 +280,7 @@ def make_sipm_selection(wfs                 : np.ndarray,
     -------
     sipm_ids_with_signal : Array of shape (n_sipms,) with boolean values indicating which SiPMs are selected.
     """
-
+    detector_info = load_db.DataSiPM('next100', 0)
     sipm_x = np.array(detector_info.X)
     sipm_y = np.array(detector_info.Y)
 

From 62b56aa4f283bcbb0bc889b6c684f40cdf428979 Mon Sep 17 00:00:00 2001
From: casper <johnwaiton@hotmail.co.uk>
Date: Wed, 25 Feb 2026 18:33:44 +0000
Subject: [PATCH 08/49] remove thresholds from `calibrate_sipms()`

---
 invisible_cities/calib/calib_sensors_functions.py | 10 +++++-----
 invisible_cities/cities/components.py             |  3 +--
 invisible_cities/cities/irene.py                  |  2 +-
 3 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/invisible_cities/calib/calib_sensors_functions.py b/invisible_cities/calib/calib_sensors_functions.py
index f8db736c3..d5387bdb0 100644
--- a/invisible_cities/calib/calib_sensors_functions.py
+++ b/invisible_cities/calib/calib_sensors_functions.py
@@ -136,15 +136,15 @@ def pmt_subtract_maw(cwfs, n_maw=100):
     return cwfs - maw
 
 
-def calibrate_sipms(sipm_wfs, adc_to_pes, thr, *, bls_mode=BlsMode.mode):
+def calibrate_sipms(sipm_wfs, adc_to_pes, *, bls_mode=BlsMode.mode):
     """
-    Subtracts the baseline, calibrates waveforms to pes
-    and suppresses values below `thr` (in pes).
+    Subtract baseline, and calibrates waveforms to pes.
     """
-    thr  = to_col_vector(np.full(sipm_wfs.shape[0], thr))
+    #thr  = to_col_vector(np.full(sipm_wfs.shape[0], thr))
     bls  = subtract_baseline(sipm_wfs, bls_mode=bls_mode)
     cwfs = calibrate_wfs(bls, adc_to_pes)
-    return np.where(cwfs > thr, cwfs, 0)
+    return cwfs
+    #return np.where(cwfs > thr, cwfs, 0)
 
 
 def subtract_mean  (wfs): return subtract_baseline(wfs, bls_mode=BlsMode.mean  )
diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index e5945bd7d..e78a8b1e8 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -785,14 +785,13 @@ def calibrate_pmts(cwf):# -> CCwfs:
     return calibrate_pmts
 
 
-def calibrate_sipms(dbfile, run_number, thr_sipm):
+def calibrate_sipms(dbfile, run_number):
     DataSiPM   = load_db.DataSiPM(dbfile, run_number)
     adc_to_pes = np.abs(DataSiPM.adc_to_pes.values)
 
     def calibrate_sipms(rwf):
         return csf.calibrate_sipms(rwf,
                                    adc_to_pes = adc_to_pes,
-                                   thr        = thr_sipm,
                                    bls_mode   = BlsMode.mode)
 
     return calibrate_sipms
diff --git a/invisible_cities/cities/irene.py b/invisible_cities/cities/irene.py
index 1c3e91ec8..7f1462045 100644
--- a/invisible_cities/cities/irene.py
+++ b/invisible_cities/cities/irene.py
@@ -90,7 +90,7 @@ def irene( files_in        : OneOrManyFiles
                               out  = ("s1_indices", "s2_indices", "s2_energies"))
 
     # Remove baseline and calibrate SiPMs
-    sipm_rwf_to_cal  = fl.map(calibrate_sipms(detector_db, run_number, sipm_thr),
+    sipm_rwf_to_cal  = fl.map(calibrate_sipms(detector_db, run_number),
                               item = "sipm")
 
     event_count_in  = fl.spy_count()

From 3e2f3339027643aec479f1173ab955144bae748a Mon Sep 17 00:00:00 2001
From: casper <johnwaiton@hotmail.co.uk>
Date: Wed, 25 Feb 2026 18:34:09 +0000
Subject: [PATCH 09/49] create `CutAlgo()` class

---
 invisible_cities/types/symbols.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/invisible_cities/types/symbols.py b/invisible_cities/types/symbols.py
index e5916f326..123679ff1 100644
--- a/invisible_cities/types/symbols.py
+++ b/invisible_cities/types/symbols.py
@@ -130,6 +130,8 @@ class XYReco(AutoNameEnumBase):
     barycenter = auto()
     corona     = auto()
 
+class CutAlgo(AutoNameEnumBase):
+    threshold  = auto()
 
 class WfType(AutoNameEnumBase):
     rwf  = auto()

From ff9346b6bdf85d2aa3db16adafaf3773969435c1 Mon Sep 17 00:00:00 2001
From: casper <johnwaiton@hotmail.co.uk>
Date: Wed, 25 Feb 2026 18:39:13 +0000
Subject: [PATCH 10/49] modify `charge_threshold_method()`

Done for backwards compatibility, this function now mirrors that
of `calibrate_sipms()` before thresholding was removed.
---
 invisible_cities/reco/wfm_functions.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index fb4388966..16b2dbc06 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -6,11 +6,12 @@
 from typing  import Optional
 from typing  import Callable
 
-from .. core.core_functions import define_window
+from .. core.core_functions import define_window, to_col_vector
 from .. calib               import calib_sensors_functions as csf
 from .. sierpe              import blr
 from .. database            import load_db
 
+
 def to_adc(wfs, adc_to_pes):
     """
     Convert waveform in pes to adc.
@@ -166,10 +167,12 @@ def charge_threshold_method(wfs       : np.ndarray,
     
     Returns
     -------
-    Boolean numpy array of shape (n_sipms,) where True indicates that the SiPM is selected.
+    2D array of shape (n_sipms, n_time_bins) containing the waveforms of
+    each SiPM with all values below threshold set to zero.
     """
-    charges = np.sum(wfs, axis=1)
-    return charges >= threshold
+
+    thr = to_col_vector(np.full(wfs.shape[0], threshold))
+    return np.where(wfs > thr, wfs, 0)
 
 
 def top_n_method(wfs : np.ndarray, 

From ea10adeb19855ebcc96621e28c22e17a04244da0 Mon Sep 17 00:00:00 2001
From: casper <johnwaiton@hotmail.co.uk>
Date: Wed, 25 Feb 2026 18:41:59 +0000
Subject: [PATCH 11/49] Implement `threshold_sipm_selection()`

---
 invisible_cities/cities/components.py | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index e78a8b1e8..7259a576e 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -797,6 +797,23 @@ def calibrate_sipms(rwf):
     return calibrate_sipms
 
 
+
+
+def threshold_sipm_selection(thr_sipm_type, thr_sipm, detector_db, run_number):
+    '''
+    Function that applies thresholding to the sipms in standard irene manner
+    '''
+
+    # extract sipm threshold
+    sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number)
+
+    def threshold_sipm_selection(wfs):
+        return wfm.charge_threshold_method(wfs, threshold = sipm_thr)
+
+    return threshold_sipm_selection
+
+
+
 def calibrate_with_mean(dbfile, run_number):
     DataSiPM   = load_db.DataSiPM(dbfile, run_number)
     adc_to_pes = np.abs(DataSiPM.adc_to_pes.values)

From 78dda04f99543f2992c93d9d200f0f5001ea56c3 Mon Sep 17 00:00:00 2001
From: casper <johnwaiton@hotmail.co.uk>
Date: Wed, 25 Feb 2026 18:42:30 +0000
Subject: [PATCH 12/49] Include `apply_cutting_function()`

This is the main controlling function, that decides which type of cut
will be applied.
---
 invisible_cities/cities/components.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index 7259a576e..a456978eb 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -83,6 +83,7 @@
 from .. types  .ic_types          import        types_dict_summary
 from .. types  .ic_types          import         types_dict_tracks
 from .. types  .symbols           import                    WfType
+from .. types  .symbols           import                   CutAlgo
 from .. types  .symbols           import               RebinMethod
 from .. types  .symbols           import                SiPMCharge
 from .. types  .symbols           import                   BlsMode
@@ -797,6 +798,18 @@ def calibrate_sipms(rwf):
     return calibrate_sipms
 
 
+def apply_cutting_function(algo, **cutting_params):
+
+    if algo is CutAlgo.threshold:
+        func = threshold_sipm_selection(**cutting_params)
+    else:
+        # temporary solution, think of a nicer method
+        func = threshold_sipm_selection(**cutting_params)
+
+    def apply_cutting_function(wfs):
+        return func(wfs)
+
+    return apply_cutting_function
 
 
 def threshold_sipm_selection(thr_sipm_type, thr_sipm, detector_db, run_number):

From 9c9172216e62aea47cf3e47617de9f726ec0037f Mon Sep 17 00:00:00 2001
From: casper <johnwaiton@hotmail.co.uk>
Date: Wed, 25 Feb 2026 18:44:47 +0000
Subject: [PATCH 13/49] Implement `apply_cut()` into irene

---
 invisible_cities/cities/irene.py | 68 +++++++++++++++++++++-----------
 1 file changed, 44 insertions(+), 24 deletions(-)

diff --git a/invisible_cities/cities/irene.py b/invisible_cities/cities/irene.py
index 7f1462045..6fa3d34d0 100644
--- a/invisible_cities/cities/irene.py
+++ b/invisible_cities/cities/irene.py
@@ -25,6 +25,7 @@
 from .. io   .run_and_event_io import run_and_event_writer
 from .. io   .trigger_io       import       trigger_writer
 from .. types.symbols          import WfType
+from .. types.symbols          import CutAlgo
 from .. types.symbols          import SiPMThreshold
 
 from .. dataflow            import dataflow as fl
@@ -43,34 +44,49 @@
 from .  components import wf_from_files
 from .  components import get_number_of_active_pmts
 from .  components import compute_and_write_pmaps
-from .  components import get_actual_sipm_thr
+from .  components import apply_cutting_function
 
+from typing import Dict
+from typing import Any
 
 @city
-def irene( files_in        : OneOrManyFiles
-         , file_out        : str
-         , compression     : str
-         , event_range     : EventRangeType
-         , print_mod       : int
-         , detector_db     : str
-         , run_number      : int
-         , n_baseline      : int
-         , n_maw           : int
-         , thr_maw         : float
-         , thr_sipm        : float
-         , thr_sipm_type   : SiPMThreshold
-         , s1_lmin         : int  , s1_lmax      : int
-         , s1_tmin         : float, s1_tmax      : float
-         , s1_rebin_stride : int  , s1_stride    : int
-         , thr_csum_s1     : float
-         , s2_lmin         : int  , s2_lmax      : int
-         , s2_tmin         : float, s2_tmax      : float
-         , s2_rebin_stride : int  , s2_stride    : int
-         , thr_csum_s2     : float, thr_sipm_s2  : float
-         , pmt_samp_wid    : float, sipm_samp_wid: float
+def irene( files_in         : OneOrManyFiles
+         , file_out         : str
+         , compression      : str
+         , event_range      : EventRangeType
+         , print_mod        : int
+         , detector_db      : str
+         , run_number       : int
+         , n_baseline       : int
+         , n_maw            : int
+         , thr_maw          : float
+         , thr_sipm         : float
+         , thr_sipm_type    : SiPMThreshold
+         , s1_lmin          : int  , s1_lmax      : int
+         , s1_tmin          : float, s1_tmax      : float
+         , s1_rebin_stride  : int  , s1_stride    : int
+         , thr_csum_s1      : float
+         , s2_lmin          : int  , s2_lmax      : int
+         , s2_tmin          : float, s2_tmax      : float
+         , s2_rebin_stride  : int  , s2_stride    : int
+         , thr_csum_s2      : float, thr_sipm_s2  : float
+         , pmt_samp_wid     : float, sipm_samp_wid: float
+         , cutting_function : CutAlgo
+         , cutting_params   : Dict[str, Any]
          ):
+    '''
+    `cutting_function` is defined within components.py, and can vary, resulting
+    in the need for `cutting_params`, which are defined as such to allow for
+    the prior function to run. Currently implemented are `threshold_sipm_selection`,
+    with differing selections methods added soon.
 
-    sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number)
+
+    params for `threshold_sipm_selection`:
+        thr_sipm_type : SiPMThreshold
+        thr_sipm      : float
+        detector_db   : str
+        run_number    : int
+    '''
 
     #### Define data transformations
 
@@ -89,9 +105,13 @@ def irene( files_in        : OneOrManyFiles
                               args = ("cwf_sum", "cwf_sum_maw"),
                               out  = ("s1_indices", "s2_indices", "s2_energies"))
 
+
     # Remove baseline and calibrate SiPMs
     sipm_rwf_to_cal  = fl.map(calibrate_sipms(detector_db, run_number),
                               item = "sipm")
+    # apply function depending on user input, from provided list of functions
+    apply_cut        = fl.map(apply_cutting_function(cutting_function, **cutting_params),
+                              item = "sipm")
 
     event_count_in  = fl.spy_count()
     event_count_out = fl.spy_count()
@@ -115,7 +135,7 @@ def irene( files_in        : OneOrManyFiles
                                          s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin,
                                          s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin,
                                          thr_sipm_s2,
-                                         h5out, sipm_rwf_to_cal)
+                                         h5out, apply_cut, sipm_rwf_to_cal)
 
         result = push(source = wf_from_files(files_in, WfType.rwf),
                       pipe   = pipe(fl.slice(*event_range, close_all=True),

From b268b328b15ec69715987f1b3c5fe45a9d10f6b2 Mon Sep 17 00:00:00 2001
From: casper <johnwaiton@hotmail.co.uk>
Date: Wed, 25 Feb 2026 18:45:20 +0000
Subject: [PATCH 14/49] Implement `apply_cut()` into
 `compute_and_write_pmaps()`

This included the adjustment of the cuts, applying
`select_wfs_above_time_integrated_thr()` and the zeroing integration cut
before the waveform rebinning.
---
 invisible_cities/cities/components.py   | 10 ++++----
 invisible_cities/cities/irene.py        |  4 ++--
 invisible_cities/reco/peak_functions.py | 32 +++++++++++++++----------
 3 files changed, 27 insertions(+), 19 deletions(-)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index a456978eb..bef70b6b8 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -747,7 +747,7 @@ def sensor_data(path, wf_type):
 
 def build_pmap(detector_db, run_number, pmt_samp_wid, sipm_samp_wid,
                s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin,
-               s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2):
+               s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, apply_cut):
     s1_params = dict(time        = minmax(min = s1_tmin,
                                           max = s1_tmax),
                     length       = minmax(min = s1_lmin,
@@ -767,8 +767,8 @@ def build_pmap(detector_db, run_number, pmt_samp_wid, sipm_samp_wid,
 
     def build_pmap(ccwf, s1_indx, s2_indx, sipmzs): # -> PMap
         return pkf.get_pmap(ccwf, s1_indx, s2_indx, sipmzs,
-                            s1_params, s2_params, thr_sipm_s2, pmt_ids,
-                            pmt_samp_wid, sipm_samp_wid)
+                            s1_params, s2_params, pmt_ids,
+                            pmt_samp_wid, sipm_samp_wid, apply_cut)
 
     return build_pmap
 
@@ -1243,7 +1243,7 @@ def integrate_wfs(wfs):
 def compute_and_write_pmaps(detector_db, run_number, pmt_samp_wid, sipm_samp_wid,
                   s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin,
                   s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2,
-                  h5out, sipm_rwf_to_cal=None):
+                  h5out, apply_cut, sipm_rwf_to_cal=None):
 
     # Filter events without signal over threshold
     indices_pass    = fl.map(check_nonempty_indices,
@@ -1254,7 +1254,7 @@ def compute_and_write_pmaps(detector_db, run_number, pmt_samp_wid, sipm_samp_wid
     # Build the PMap
     compute_pmap     = fl.map(build_pmap(detector_db, run_number, pmt_samp_wid, sipm_samp_wid,
                                          s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin,
-                                         s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2),
+                                         s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, apply_cut),
                               args = ("ccwfs", "s1_indices", "s2_indices", "sipm"),
                               out  = "pmap")
 
diff --git a/invisible_cities/cities/irene.py b/invisible_cities/cities/irene.py
index 6fa3d34d0..31a6b8041 100644
--- a/invisible_cities/cities/irene.py
+++ b/invisible_cities/cities/irene.py
@@ -110,8 +110,8 @@ def irene( files_in         : OneOrManyFiles
     sipm_rwf_to_cal  = fl.map(calibrate_sipms(detector_db, run_number),
                               item = "sipm")
     # apply function depending on user input, from provided list of functions
-    apply_cut        = fl.map(apply_cutting_function(cutting_function, **cutting_params),
-                              item = "sipm")
+    apply_cut        = apply_cutting_function(cutting_function, **cutting_params)
+
 
     event_count_in  = fl.spy_count()
     event_count_out = fl.spy_count()
diff --git a/invisible_cities/reco/peak_functions.py b/invisible_cities/reco/peak_functions.py
index 686b16c87..2aa9ef186 100644
--- a/invisible_cities/reco/peak_functions.py
+++ b/invisible_cities/reco/peak_functions.py
@@ -76,13 +76,20 @@ def build_pmt_responses(indices, times, widths, ccwf,
 
 
 def build_sipm_responses(indices, times, widths,
-                         sipm_wfs, rebin_stride, thr_sipm_s2):
-    _, _, sipm_wfs_ = pick_slice_and_rebin(indices , times, widths,
+                         sipm_wfs, rebin_stride, apply_cut):
+
+    if apply_cut is not None:
+        # apply cut before slicing and rebinning
+        (sipm_ids,
+         sipm_wfs)   = apply_cut(sipm_wfs)
+    else:
+        # give all sipm ids as index if no cut is applied
+        sipm_ids     = np.arange(sipm_wfs.shape[0])
+    # rebin
+    _, _, sipm_wfs = pick_slice_and_rebin(indices , times, widths,
                                            sipm_wfs, rebin_stride,
                                            pad_zeros = False)
-    (sipm_ids,
-     sipm_wfs)   = select_wfs_above_time_integrated_thr(sipm_wfs_,
-                                                        thr_sipm_s2)
+
     return SiPMResponses(sipm_ids, sipm_wfs)
 
 
@@ -93,7 +100,7 @@ def build_peak(indices, times,
                pmt_samp_wid  = 25 * units.ns,
                sipm_samp_wid =  1 * units.mus,
                sipm_wfs      = None,
-               thr_sipm_s2   = 0):
+               apply_cut     = None):
     sipm_pmt_bin_ratio = int(sipm_samp_wid/pmt_samp_wid)
     (pk_times ,
      pk_widths,
@@ -107,7 +114,7 @@ def build_peak(indices, times,
                                       widths * sipm_pmt_bin_ratio,
                                       sipm_wfs,
                                       rebin_stride // sipm_pmt_bin_ratio,
-                                      thr_sipm_s2)
+                                      apply_cut)
     else:
         sipm_r = SiPMResponses.build_empty_instance()
 
@@ -120,7 +127,8 @@ def find_peaks(ccwfs, index,
                Pk, pmt_ids,
                pmt_samp_wid = 25*units.ns,
                sipm_samp_wid = 1*units.mus,
-               sipm_wfs=None, thr_sipm_s2=0):
+               sipm_wfs=None, apply_cut = None):
+
     ccwfs = np.array(ccwfs, ndmin=2)
 
     peaks           = []
@@ -136,20 +144,20 @@ def find_peaks(ccwfs, index,
                         rebin_stride,
                         with_sipms, Pk,
                         pmt_samp_wid, sipm_samp_wid,
-                        sipm_wfs, thr_sipm_s2)
+                        sipm_wfs, apply_cut)
         peaks.append(pk)
     return peaks
 
 
 def get_pmap(ccwf, s1_indx, s2_indx, sipm_zs_wf,
-             s1_params, s2_params, thr_sipm_s2, pmt_ids,
-             pmt_samp_wid, sipm_samp_wid):
+             s1_params, s2_params, pmt_ids,
+             pmt_samp_wid, sipm_samp_wid, apply_cut = None):
     return PMap(find_peaks(ccwf, s1_indx, Pk=S1, pmt_ids=pmt_ids,
                            pmt_samp_wid=pmt_samp_wid,
                            **s1_params),
                 find_peaks(ccwf, s2_indx, Pk=S2, pmt_ids=pmt_ids,
                            sipm_wfs      = sipm_zs_wf,
-                           thr_sipm_s2   = thr_sipm_s2,
+                           apply_cut     = apply_cut,
                            pmt_samp_wid  = pmt_samp_wid,
                            sipm_samp_wid = sipm_samp_wid,
                            **s2_params))

From 1ead58b52b6b448f5d5f71403a96ec6cd9093a79 Mon Sep 17 00:00:00 2001
From: casper <johnwaiton@hotmail.co.uk>
Date: Wed, 25 Feb 2026 18:45:52 +0000
Subject: [PATCH 15/49] Update default irene config

---
 invisible_cities/config/irene.conf | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/invisible_cities/config/irene.conf b/invisible_cities/config/irene.conf
index 9b536c7c1..f62f962c4 100644
--- a/invisible_cities/config/irene.conf
+++ b/invisible_cities/config/irene.conf
@@ -52,3 +52,9 @@ thr_sipm_s2 = 10 * pes  # Threshold for the full sipm waveform
 
 pmt_samp_wid  = 25 * ns
 sipm_samp_wid = 1 * mus
+
+cutting_function = threshold 
+cutting_params   = dict(  thr_sipm_type = thr_sipm_type 
+                        , thr_sipm      = thr_sipm 
+                        , detector_db   = detector_db 
+                        , run_number    = run_number)

From d0c1c2433d362dfa5ece6db7747ca26c85b6f027 Mon Sep 17 00:00:00 2001
From: casper <johnwaiton@hotmail.co.uk>
Date: Wed, 25 Feb 2026 18:46:10 +0000
Subject: [PATCH 16/49] Remove whitespaces, add newline

---
 invisible_cities/reco/wfm_functions.py | 43 +++++++++++++-------------
 1 file changed, 22 insertions(+), 21 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 16b2dbc06..a66297d9c 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -135,17 +135,17 @@ def compare_cwf_blr(cwf, pmtblr, event_list, window_size=500):
     return np.array(DIFF)
 
 
-def median_std_method(wfs    : np.ndarray,  
+def median_std_method(wfs    : np.ndarray,
                       nsigma : Optional[float] = 3.) -> np.ndarray:
     """
-    Computes the median and standard deviation of the time summed SiPM 
+    Computes the median and standard deviation of the time summed SiPM
     waveforms and selects the SiPMs that are nsigma over the median.
 
     Parameters
     ----------
     wfs    : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
     nsigma : Number of standard deviations above the median, default 3.
-    
+
     Returns
     -------
     Boolean numpy array of shape (n_sipms,) where True indicates that the SiPM is selected.
@@ -159,12 +159,12 @@ def charge_threshold_method(wfs       : np.ndarray,
                             threshold : Optional[float] = 5.) -> np.ndarray:
     """
     Selects the SiPMs whose time summed waveforms are above a threshold.
-    
+
     Parameters
     ----------
     wfs       : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
     threshold : Charge threshold in PE, default 5.
-    
+
     Returns
     -------
     2D array of shape (n_sipms, n_time_bins) containing the waveforms of
@@ -175,16 +175,16 @@ def charge_threshold_method(wfs       : np.ndarray,
     return np.where(wfs > thr, wfs, 0)
 
 
-def top_n_method(wfs : np.ndarray, 
+def top_n_method(wfs : np.ndarray,
                  n   : Optional[int] = 10) -> np.ndarray:
     """
     Selects the SiPMs with the top n highest time summed waveforms.
-        
+
     Parameters
     ----------
     wfs       : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
     n         : Number of most energeticSiPMs to select, default 10.
-    
+
     Returns
     -------
     Boolean numpy array of shape (n_sipms,) where True indicates that the SiPM is selected.
@@ -198,12 +198,12 @@ def top_n_method(wfs : np.ndarray,
 
 
 def kill_isolated_sipms(selected_ids        : np.ndarray,
-                        sipm_x              : np.ndarray, 
-                        sipm_y              : np.ndarray, 
+                        sipm_x              : np.ndarray,
+                        sipm_y              : np.ndarray,
                         proximity_threshold : float) -> np.ndarray:
     """
-    For the SiPMs that have passed the previous selection, scans through the SiPMs to check if they 
-    have neighbouring SiPMs - i.e., within the proximity_threshold - that have also passed the selection. 
+    For the SiPMs that have passed the previous selection, scans through the SiPMs to check if they
+    have neighbouring SiPMs - i.e., within the proximity_threshold - that have also passed the selection.
     If no neighbours are found, the SiPMs are classed as isolated, and are removed.
 
     Parameters
@@ -228,16 +228,16 @@ def kill_isolated_sipms(selected_ids        : np.ndarray,
 
         if n_neighbors <= 1:
             selected_ids_no_isolated[i] = False
-            
+
     return selected_ids_no_isolated
 
 
 def apply_circular_padding(selected_ids_no_isolated : np.ndarray,
-                           sipm_x                   : np.ndarray, 
-                           sipm_y                   : np.ndarray, 
+                           sipm_x                   : np.ndarray,
+                           sipm_y                   : np.ndarray,
                            padding_radius           : float) -> np.ndarray:
     """
-    For the SiPMs that pass the previous selection, creates circular padding of radius padding_radius, 
+    For the SiPMs that pass the previous selection, creates circular padding of radius padding_radius,
     selecting all SiPMs within that radius. Stores the union of all selected SiPMs.
 
     Parameters
@@ -261,11 +261,11 @@ def apply_circular_padding(selected_ids_no_isolated : np.ndarray,
     return sipm_ids_with_signal
 
 
-def make_sipm_selection(wfs                 : np.ndarray, 
-                        selection_func      : Callable, 
-                        selection_kwargs    : dict, 
-                        proximity_threshold : float, 
-                        padding_radius      : float) -> np.ndarray: 
+def make_sipm_selection(wfs                 : np.ndarray,
+                        selection_func      : Callable,
+                        selection_kwargs    : dict,
+                        proximity_threshold : float,
+                        padding_radius      : float) -> np.ndarray:
     """
     SiPM selection pipeline, applies SiPM cuts based on user input.
     A first selection of SiPMs is made, isolated SiPMs are removed
@@ -303,3 +303,4 @@ def make_sipm_selection(wfs                 : np.ndarray,
         padding_radius
     )
     return sipm_ids_with_signal
+

From a675fcb13581b1d41d282e8382bba3e96a19a587 Mon Sep 17 00:00:00 2001
From: jwaiton <john.waiton@postgrad.manchester.ac.uk>
Date: Fri, 27 Feb 2026 11:06:35 +0000
Subject: [PATCH 17/49] Allow `calibrate_sipms()` to apply threshold.

This is to avoid possible situations where `calibrate_sipms()` is used
outwith irene. `thr` is now a keyword argument to ensure that if a
threshold is applied it is done so with intent.
---
 invisible_cities/calib/calib_sensors_functions.py     | 11 +++++++----
 .../calib/calib_sensors_functions_test.py             |  6 +++---
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/invisible_cities/calib/calib_sensors_functions.py b/invisible_cities/calib/calib_sensors_functions.py
index d5387bdb0..ba2c8854e 100644
--- a/invisible_cities/calib/calib_sensors_functions.py
+++ b/invisible_cities/calib/calib_sensors_functions.py
@@ -136,15 +136,18 @@ def pmt_subtract_maw(cwfs, n_maw=100):
     return cwfs - maw
 
 
-def calibrate_sipms(sipm_wfs, adc_to_pes, *, bls_mode=BlsMode.mode):
+def calibrate_sipms(sipm_wfs, adc_to_pes, *, thr = None, bls_mode=BlsMode.mode):
     """
     Subtract baseline, and calibrates waveforms to pes.
     """
-    #thr  = to_col_vector(np.full(sipm_wfs.shape[0], thr))
     bls  = subtract_baseline(sipm_wfs, bls_mode=bls_mode)
     cwfs = calibrate_wfs(bls, adc_to_pes)
-    return cwfs
-    #return np.where(cwfs > thr, cwfs, 0)
+    if thr is None:
+        return cwfs
+    else:
+        # if you apply a threshold here, apply as usual
+        thr  = to_col_vector(np.full(sipm_wfs.shape[0], thr))
+        return np.where(cwfs > thr, cwfs, 0)
 
 
 def subtract_mean  (wfs): return subtract_baseline(wfs, bls_mode=BlsMode.mean  )
diff --git a/invisible_cities/calib/calib_sensors_functions_test.py b/invisible_cities/calib/calib_sensors_functions_test.py
index 452d22d90..bdbfc673e 100644
--- a/invisible_cities/calib/calib_sensors_functions_test.py
+++ b/invisible_cities/calib/calib_sensors_functions_test.py
@@ -207,7 +207,7 @@ def test_calibrate_sipms_stat(oscillating_waveform_with_baseline,
      baseline )              = oscillating_waveform_with_baseline
     #n_maw                    = n_samples // 500
 
-    ccwfs = csf.calibrate_sipms(wfs, adc_to_pes, nsigma * noise_sigma, bls_mode=BlsMode.mode)
+    ccwfs = csf.calibrate_sipms(wfs, adc_to_pes, thr=nsigma * noise_sigma, bls_mode=BlsMode.mode)
 
     number_of_zeros = np.count_nonzero(ccwfs == 0)
     assert number_of_zeros > fraction * ccwfs.size
@@ -219,7 +219,7 @@ def test_calibrate_sipms_common_threshold(toy_sipm_signal):
      common_threshold, _) = toy_sipm_signal
 
     zs_wf = csf.calibrate_sipms(signal_adc, adc_to_pes,
-                                common_threshold, bls_mode=BlsMode.mode)
+                                thr=common_threshold, bls_mode=BlsMode.mode)
 
     for actual, expected in zip(zs_wf, signal_zs_common_threshold):
         assert actual == approx(expected)
@@ -232,7 +232,7 @@ def test_calibrate_sipms_individual_thresholds(toy_sipm_signal):
 
 
     zs_wf = csf.calibrate_sipms(signal_adc, adc_to_pes,
-                                individual_thresholds,
+                                thr=individual_thresholds,
                                 bls_mode=BlsMode.mode)
     for actual, expected in zip(zs_wf, signal_zs_individual_thresholds):
         assert actual == approx(expected)

From bdb026e9caa9f76f738eb83fc67c2b6a8a01ac0d Mon Sep 17 00:00:00 2001
From: jwaiton <john.waiton@postgrad.manchester.ac.uk>
Date: Fri, 27 Feb 2026 11:20:10 +0000
Subject: [PATCH 18/49] Include `thr_sipm_s2` into `cutting_params`

---
 invisible_cities/config/irene.conf | 1 +
 1 file changed, 1 insertion(+)

diff --git a/invisible_cities/config/irene.conf b/invisible_cities/config/irene.conf
index f62f962c4..1f2312a78 100644
--- a/invisible_cities/config/irene.conf
+++ b/invisible_cities/config/irene.conf
@@ -56,5 +56,6 @@ sipm_samp_wid = 1 * mus
 cutting_function = threshold 
 cutting_params   = dict(  thr_sipm_type = thr_sipm_type 
                         , thr_sipm      = thr_sipm 
+                        , thr_sipm_s2   = thr_sipm_s2
                         , detector_db   = detector_db 
                         , run_number    = run_number)

From 1e5e88faecdad21dab5549985b138fed8148d178 Mon Sep 17 00:00:00 2001
From: jwaiton <john.waiton@postgrad.manchester.ac.uk>
Date: Fri, 27 Feb 2026 11:33:08 +0000
Subject: [PATCH 19/49] Implement `thr_sipm_s2` in `charge_threshold_method()`

---
 invisible_cities/cities/components.py  | 22 +++++++++++++++-------
 invisible_cities/reco/wfm_functions.py | 15 +++++++++++----
 2 files changed, 26 insertions(+), 11 deletions(-)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index bef70b6b8..0aadbc4b3 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -812,16 +812,24 @@ def apply_cutting_function(wfs):
     return apply_cutting_function
 
 
-def threshold_sipm_selection(thr_sipm_type, thr_sipm, detector_db, run_number):
+def threshold_sipm_selection(thr_sipm_type
+                            , thr_sipm
+                            , thr_sipm_s2
+                            , run_number
+                            , detector_db = None):
     '''
-    Function that applies thresholding to the sipms in standard irene manner
+    Function that applies thresholding to the sipms in standard irene manner,
+    by zeroing all waveform values below a threshold.
     '''
+    # assume that if the detector_db is None, you return sipm threshold as the number provided
+    if detector_db is None:
+        sipm_thr = thr_sipm
+    else:
+        # extract sipm threshold
+        sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number)
 
-    # extract sipm threshold
-    sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number)
-
-    def threshold_sipm_selection(wfs):
-        return wfm.charge_threshold_method(wfs, threshold = sipm_thr)
+    def threshold_sipm_selection(wfs, indices):
+        return wfm.charge_threshold_method(wfs, indices, zeroing_thr = sipm_thr, integration_thr=thr_sipm_s2)
 
     return threshold_sipm_selection
 
diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index a66297d9c..9426651ab 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -155,10 +155,12 @@ def median_std_method(wfs    : np.ndarray,
     return charges >= threshold
 
 
-def charge_threshold_method(wfs       : np.ndarray,
-                            threshold : Optional[float] = 5.) -> np.ndarray:
+def charge_threshold_method(wfs             : np.ndarray,
+                            indices         : np.ndarray,
+                            zeroing_thr     : Optional[float] = 2.,
+                            integration_thr : Optional[float] = 5.) -> np.ndarray:
     """
-    Selects the SiPMs whose time summed waveforms are above a threshold.
+    Selects the SiPMs whose time summed waveforms within the s2 windows are above a threshold.
 
     Parameters
     ----------
@@ -168,8 +170,13 @@ def charge_threshold_method(wfs       : np.ndarray,
     Returns
     -------
     2D array of shape (n_sipms, n_time_bins) containing the waveforms of
-    each SiPM with all values below threshold set to zero.
+    each SiPM with all values below threshold set to zero, with a mask applied depending
+    on summed waveforms over certain indices.
     """
+    thr = to_col_vector(np.full(wfs.shape[0], zeroing_thr))
+
+    # zero entries below threshold
+    zwfs = np.where(wfs > thr, wfs, 0)
 
     thr = to_col_vector(np.full(wfs.shape[0], threshold))
     return np.where(wfs > thr, wfs, 0)

From 78cb788c7dc59558fe46343a2508de92a786aa0a Mon Sep 17 00:00:00 2001
From: jwaiton <john.waiton@postgrad.manchester.ac.uk>
Date: Fri, 27 Feb 2026 11:35:39 +0000
Subject: [PATCH 20/49] Implement integration cut based on slices

Initially this was done over the entire waveforms due to the reordering
of the cuts when applied in `build_sipm_responses()`. As suggested by
@Ian0sborne, it would be wise to move these cuts out of
`compute_and_write_pmaps()`
---
 invisible_cities/cities/components.py   |  4 ++--
 invisible_cities/reco/peak_functions.py | 17 ++++++++++++++++-
 invisible_cities/reco/wfm_functions.py  |  6 +++---
 3 files changed, 21 insertions(+), 6 deletions(-)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index 0aadbc4b3..3b9447b14 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -806,8 +806,8 @@ def apply_cutting_function(algo, **cutting_params):
         # temporary solution, think of a nicer method
         func = threshold_sipm_selection(**cutting_params)
 
-    def apply_cutting_function(wfs):
-        return func(wfs)
+    def apply_cutting_function(wfs, indices):
+        return func(wfs, indices)
 
     return apply_cutting_function
 
diff --git a/invisible_cities/reco/peak_functions.py b/invisible_cities/reco/peak_functions.py
index 2aa9ef186..1a0d829ef 100644
--- a/invisible_cities/reco/peak_functions.py
+++ b/invisible_cities/reco/peak_functions.py
@@ -25,9 +25,24 @@ def indices_and_wf_above_threshold(wf, thr):
     return ZsWf(indices_above_thr, wf_above_thr)
 
 
+def select_wf_slices_above_time_integrated_thr(wfs, indices, thr):
+    '''
+    function that integrates over certain time slices, and passes the waveform
+    based on the aforementioned slices passing the threshold
+    '''
+    slice_ = slice(indices[0], indices[-1] + 1)
+    wfs_   = wfs[:, slice_]
+
+    selected_ids = np.where(np.sum(wfs_, axis = 1) >= thr)[0]
+    selected_wfs = wfs[selected_ids]
+
+    return selected_ids, selected_wfs
+
+
 def select_wfs_above_time_integrated_thr(wfs, thr):
     selected_ids = np.where(np.sum(wfs, axis=1) >= thr)[0]
     selected_wfs = wfs[selected_ids]
+
     return selected_ids, selected_wfs
 
 
@@ -81,7 +96,7 @@ def build_sipm_responses(indices, times, widths,
     if apply_cut is not None:
         # apply cut before slicing and rebinning
         (sipm_ids,
-         sipm_wfs)   = apply_cut(sipm_wfs)
+         sipm_wfs)   = apply_cut(sipm_wfs, indices)
     else:
         # give all sipm ids as index if no cut is applied
         sipm_ids     = np.arange(sipm_wfs.shape[0])
diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 9426651ab..82015d0b7 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -10,7 +10,7 @@
 from .. calib               import calib_sensors_functions as csf
 from .. sierpe              import blr
 from .. database            import load_db
-
+from .. reco.peak_functions import select_wf_slices_above_time_integrated_thr
 
 def to_adc(wfs, adc_to_pes):
     """
@@ -178,8 +178,8 @@ def charge_threshold_method(wfs             : np.ndarray,
     # zero entries below threshold
     zwfs = np.where(wfs > thr, wfs, 0)
 
-    thr = to_col_vector(np.full(wfs.shape[0], threshold))
-    return np.where(wfs > thr, wfs, 0)
+    # returns selected ids and waveforms above integral
+    return select_wf_slices_above_time_integrated_thr(zwfs, indices, integration_thr)
 
 
 def top_n_method(wfs : np.ndarray,

From 621bc7977c793c69dcd8ac4da38c7553dd1db541 Mon Sep 17 00:00:00 2001
From: jwaiton <john.waiton@postgrad.manchester.ac.uk>
Date: Fri, 27 Feb 2026 11:38:15 +0000
Subject: [PATCH 21/49] Adjust tests to account for `cut_params`

---
 invisible_cities/reco/peak_functions_test.py | 21 ++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/invisible_cities/reco/peak_functions_test.py b/invisible_cities/reco/peak_functions_test.py
index 16efa9bcf..bdde0a2b3 100644
--- a/invisible_cities/reco/peak_functions_test.py
+++ b/invisible_cities/reco/peak_functions_test.py
@@ -13,6 +13,7 @@
 from hypothesis.strategies    import integers
 from hypothesis.extra.numpy   import arrays
 
+from ..cities.components      import apply_cutting_function
 from ..core.testing_utils   import exactly
 from ..core.testing_utils   import previous_float
 from ..core.testing_utils   import assert_Peak_equality
@@ -26,6 +27,8 @@
 from ..evm .pmaps           import PMap
 from ..io  .pmaps_io        import load_pmaps
 from ..types.ic_types       import minmax
+from ..types.symbols        import SiPMThreshold
+from ..types.symbols        import CutAlgo
 from .                      import peak_functions as pf
 
 
@@ -350,10 +353,18 @@ def test_build_sipm_responses(wf_with_indices):
     wfs_slice       = wfs[:, indices]
     peak_integrals  = wfs_slice.sum(axis=1)
     below_thr_index = np.argmin (peak_integrals)
-    # next_float doesn't work here
     thr             = peak_integrals[below_thr_index] * 1.000001
+
+    # next_float doesn't work here
+    cut_params = dict(detector_db   = 'None',
+                  thr_sipm_s2   = thr,
+                  thr_sipm      = 0,
+                  thr_sipm_type = SiPMThreshold.common,
+                  run_number    = 0)
+    apply_cut  = apply_cutting_function(CutAlgo.threshold, **cut_params)
+
     sipm_r          = pf.build_sipm_responses(indices, times, widths,
-                                              wfs, 1, thr)
+                                              wfs, 1, apply_cut)
 
     expected_ids = np.delete(      ids, below_thr_index)
     expected_wfs = np.delete(wfs_slice, below_thr_index, axis=0)
@@ -396,7 +407,7 @@ def test_build_peak_development(pmt_and_sipm_wfs_with_indices,
                          with_sipms   = with_sipms,
                          Pk           = Pk,
                          sipm_wfs     = sipm_wfs,
-                         thr_sipm_s2  = -1)
+                         )
 
     assert_Peak_equality(peak, expected_peak)
 
@@ -473,8 +484,7 @@ def test_find_peaks_s2_style(pmt_and_sipm_wfs_with_indices):
                           time_range, length_range,
                           stride, rebin_stride,
                           S2, pmt_ids,
-                          sipm_wfs    = sipm_wfs,
-                          thr_sipm_s2 = -1)
+                          sipm_wfs    = sipm_wfs)
 
     (rebinned_times,
      rebinned_widths,
@@ -502,7 +512,6 @@ def test_get_pmap(s1_and_s2_with_indices):
     sipm_samp_wid = 1  * units.mus
     pmap = pf.get_pmap(pmt_wfs, s1_indx, s2_indx, sipm_wfs,
                        s1_params, s2_params,
-                       thr_sipm_s2 = -1,
                        pmt_ids     = pmt_ids,
                        pmt_samp_wid  = pmt_samp_wid ,
                        sipm_samp_wid = sipm_samp_wid)

From dac33dd581dad547a9068183727bfd6e909854b1 Mon Sep 17 00:00:00 2001
From: jwaiton <john.waiton@postgrad.manchester.ac.uk>
Date: Fri, 27 Feb 2026 11:38:35 +0000
Subject: [PATCH 22/49] Implement modifiable cuts in hypathia

---
 invisible_cities/cities/hypathia.py   | 61 ++++++++++++++++-----------
 invisible_cities/config/hypathia.conf |  7 +++
 2 files changed, 43 insertions(+), 25 deletions(-)

diff --git a/invisible_cities/cities/hypathia.py b/invisible_cities/cities/hypathia.py
index 6a125b542..08bbd3be8 100644
--- a/invisible_cities/cities/hypathia.py
+++ b/invisible_cities/cities/hypathia.py
@@ -30,6 +30,7 @@
 from .. io   .run_and_event_io import run_and_event_writer
 from .. io   .      trigger_io import       trigger_writer
 from .. types.symbols          import WfType
+from .. types.symbols          import CutAlgo
 from .. types.symbols          import SiPMThreshold
 
 from .. dataflow            import dataflow as fl
@@ -49,32 +50,38 @@
 from .  components import simulate_sipm_response
 from .  components import calibrate_sipms
 from .  components import get_actual_sipm_thr
+from .  components import apply_cutting_function
+
+from typing import Dict
+from typing import Any
 
 
 @city
-def hypathia( files_in        : OneOrManyFiles
-            , file_out        : str
-            , compression     : str
-            , event_range     : EventRangeType
-            , print_mod       : int
-            , detector_db     : str
-            , run_number      : int
-            , sipm_noise_cut  : float
-            , filter_padding  : int
-            , thr_sipm        : float
-            , thr_sipm_type   : SiPMThreshold
-            , pmt_wfs_rebin   : int
-            , pmt_pe_rms      : float
-            , s1_lmin         : int  , s1_lmax     : int
-            , s1_tmin         : float, s1_tmax     : float
-            , s1_rebin_stride : int  , s1_stride   : int
-            , thr_csum_s1     : float
-            , s2_lmin         : int  , s2_lmax     : int
-            , s2_tmin         : float, s2_tmax     : float
-            , s2_rebin_stride : int  , s2_stride   : int
-            , thr_csum_s2     : float, thr_sipm_s2 : float
-            , pmt_samp_wid    : float
-            , sipm_samp_wid   : float
+def hypathia( files_in         : OneOrManyFiles
+            , file_out         : str
+            , compression      : str
+            , event_range      : EventRangeType
+            , print_mod        : int
+            , detector_db      : str
+            , run_number       : int
+            , sipm_noise_cut   : float
+            , filter_padding   : int
+            , thr_sipm         : float
+            , thr_sipm_type    : SiPMThreshold
+            , pmt_wfs_rebin    : int
+            , pmt_pe_rms       : float
+            , s1_lmin          : int  , s1_lmax     : int
+            , s1_tmin          : float, s1_tmax     : float
+            , s1_rebin_stride  : int  , s1_stride   : int
+            , thr_csum_s1      : float
+            , s2_lmin          : int  , s2_lmax     : int
+            , s2_tmin          : float, s2_tmax     : float
+            , s2_rebin_stride  : int  , s2_stride   : int
+            , thr_csum_s2      : float, thr_sipm_s2 : float
+            , pmt_samp_wid     : float
+            , sipm_samp_wid    : float
+            , cutting_function : CutAlgo
+            , cutting_params   : Dict[str, Any]
             ):
 
     sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number)
@@ -112,9 +119,13 @@ def hypathia( files_in        : OneOrManyFiles
                               item="sipm")
 
     # SiPMs calibration
-    sipm_rwf_to_cal  = fl.map(calibrate_sipms(detector_db, run_number, sipm_thr),
+    sipm_rwf_to_cal  = fl.map(calibrate_sipms(detector_db, run_number),
                               item = "sipm")
 
+    # apply function depending on user input, from provided list of functions
+    apply_cut        = apply_cutting_function(cutting_function, **cutting_params)
+
+
     event_count_in  = fl.spy_count()
     event_count_out = fl.spy_count()
 
@@ -134,7 +145,7 @@ def hypathia( files_in        : OneOrManyFiles
                                              detector_db, run_number, pmt_samp_wid, sipm_samp_wid,
                                              s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin,
                                              s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2,
-                                             h5out, sipm_rwf_to_cal)
+                                             h5out, apply_cut, sipm_rwf_to_cal)
 
         result = push(source = wf_from_files(files_in, WfType.mcrd),
                       pipe   = pipe(fl.slice(*event_range, close_all=True),
diff --git a/invisible_cities/config/hypathia.conf b/invisible_cities/config/hypathia.conf
index 3aafcd8e4..d3029f25f 100644
--- a/invisible_cities/config/hypathia.conf
+++ b/invisible_cities/config/hypathia.conf
@@ -52,3 +52,10 @@ thr_sipm_s2 = 5 * pes  # Threshold for the full sipm waveform
 
 pmt_samp_wid  = 25 * ns
 sipm_samp_wid = 1 * mus
+
+cutting_function = threshold 
+cutting_params   = dict(  thr_sipm_type = thr_sipm_type 
+                        , thr_sipm      = thr_sipm 
+                        , thr_sipm_s2   = thr_sipm_s2
+                        , detector_db   = detector_db 
+                        , run_number    = run_number)

From 7cfd74369e62f7410e3b9d21366890c656d81059 Mon Sep 17 00:00:00 2001
From: jwaiton <john.waiton@postgrad.manchester.ac.uk>
Date: Fri, 27 Feb 2026 11:47:05 +0000
Subject: [PATCH 23/49] Update docstrings, type annotation

---
 invisible_cities/reco/wfm_functions.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 82015d0b7..b9130311b 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -5,6 +5,7 @@
 import numpy as np
 from typing  import Optional
 from typing  import Callable
+from typing  import Tuple
 
 from .. core.core_functions import define_window, to_col_vector
 from .. calib               import calib_sensors_functions as csf
@@ -158,9 +159,11 @@ def median_std_method(wfs    : np.ndarray,
 def charge_threshold_method(wfs             : np.ndarray,
                             indices         : np.ndarray,
                             zeroing_thr     : Optional[float] = 2.,
-                            integration_thr : Optional[float] = 5.) -> np.ndarray:
+                            integration_thr : Optional[float] = 5.) -> Tuple[np.ndarray, np.ndarray]:
     """
-    Selects the SiPMs whose time summed waveforms within the s2 windows are above a threshold.
+    Selects the SiPMs whose time summed waveforms within the s2 windows are above two thresholds:
+        - initial zero suprresion threshold (setting values in each waveform below a value to 0)
+        - threshold over each selected integrated slice of the waveforms
 
     Parameters
     ----------
@@ -169,9 +172,7 @@ def charge_threshold_method(wfs             : np.ndarray,
 
     Returns
     -------
-    2D array of shape (n_sipms, n_time_bins) containing the waveforms of
-    each SiPM with all values below threshold set to zero, with a mask applied depending
-    on summed waveforms over certain indices.
+    Tuple of np arrays including all passing sipm ids and the corresponding waveforms
     """
     thr = to_col_vector(np.full(wfs.shape[0], zeroing_thr))
 

From 502559d26eb2a05f5749d0e0f6b92e9cea327f0a Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Tue, 3 Mar 2026 15:55:12 +0000
Subject: [PATCH 24/49] add pyrrha cutting method as an option in `CutAlgo` and
 `apply_cutting_function`

---
 invisible_cities/cities/components.py | 3 ++-
 invisible_cities/types/symbols.py     | 1 +
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index 3b9447b14..584855370 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -802,6 +802,8 @@ def apply_cutting_function(algo, **cutting_params):
 
     if algo is CutAlgo.threshold:
         func = threshold_sipm_selection(**cutting_params)
+    elif algo is CutAlgo.pyrrha:
+        func = pyrrha_sipm_selection(**cutting_params)
     else:
         # temporary solution, think of a nicer method
         func = threshold_sipm_selection(**cutting_params)
@@ -834,7 +836,6 @@ def threshold_sipm_selection(wfs, indices):
     return threshold_sipm_selection
 
 
-
 def calibrate_with_mean(dbfile, run_number):
     DataSiPM   = load_db.DataSiPM(dbfile, run_number)
     adc_to_pes = np.abs(DataSiPM.adc_to_pes.values)
diff --git a/invisible_cities/types/symbols.py b/invisible_cities/types/symbols.py
index 123679ff1..c8c43b3da 100644
--- a/invisible_cities/types/symbols.py
+++ b/invisible_cities/types/symbols.py
@@ -132,6 +132,7 @@ class XYReco(AutoNameEnumBase):
 
 class CutAlgo(AutoNameEnumBase):
     threshold  = auto()
+    pyrrha     = auto()
 
 class WfType(AutoNameEnumBase):
     rwf  = auto()

From 5d5b78670ca5f01a0b5bd53e44c591417ed8c4c8 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Tue, 3 Mar 2026 16:30:24 +0000
Subject: [PATCH 25/49] pass `load_db` parameters into `make_sipm_selection`
 function

---
 invisible_cities/reco/wfm_functions.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index b9130311b..9ad067549 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -273,7 +273,9 @@ def make_sipm_selection(wfs                 : np.ndarray,
                         selection_func      : Callable,
                         selection_kwargs    : dict,
                         proximity_threshold : float,
-                        padding_radius      : float) -> np.ndarray:
+                        padding_radius      : float,
+                        run_number          : int,
+                        detector_db         : str) -> np.ndarray:
     """
     SiPM selection pipeline, applies SiPM cuts based on user input.
     A first selection of SiPMs is made, isolated SiPMs are removed
@@ -291,7 +293,7 @@ def make_sipm_selection(wfs                 : np.ndarray,
     -------
     sipm_ids_with_signal : Array of shape (n_sipms,) with boolean values indicating which SiPMs are selected.
     """
-    detector_info = load_db.DataSiPM('next100', 0)
+    detector_info = load_db.DataSiPM(detector_db, run_number)
     sipm_x = np.array(detector_info.X)
     sipm_y = np.array(detector_info.Y)
 

From 04a0867838d6168c634736a223e4fab661acec11 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Tue, 3 Mar 2026 16:33:13 +0000
Subject: [PATCH 26/49] add documentation to `make_sipm_selection`

---
 invisible_cities/reco/wfm_functions.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 9ad067549..7bf3eaef7 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -288,6 +288,8 @@ def make_sipm_selection(wfs                 : np.ndarray,
     selection_kwargs    : Dictionary of arguments passed to selection_func.
     proximity_threshold : Threshold used to identify isolated SiPMs.
     padding_radius      : Radial padding added to each SiPM that passes the selections.
+    run_number          : Run number used to load the detector database.
+    detector_db         : Database used to load the detector geometry.
 
     Returns
     -------

From 4f243a3245cb3ffc5ba69742f03f542f15058493 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Tue, 3 Mar 2026 16:58:33 +0000
Subject: [PATCH 27/49] implement `spatial_selection_method` (previously named
 `make_sipm_selection`) into `pyrrha_sipm_selection`

---
 invisible_cities/cities/components.py  | 22 ++++++++++++++++++++++
 invisible_cities/reco/wfm_functions.py | 16 ++++++++--------
 2 files changed, 30 insertions(+), 8 deletions(-)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index 584855370..3b3c1c3a3 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -836,6 +836,28 @@ def threshold_sipm_selection(wfs, indices):
     return threshold_sipm_selection
 
 
+def pyrrha_sipm_selection(selection_function : Callable
+                         , selection_kwargs    : dict
+                         , proximity_threshold : float
+                         , padding_radius      : float
+                         , run_number : int
+                         , detector_db : str):
+    '''
+    Function that applies a generic selection function to the sipms, which can be used to 
+    implement a spatial SiPM selection method (called Pyrrha).
+    '''
+    def pyrrha_sipm_selection(wfs, indices):
+        return wfm.spatial_selection_method(wfs, 
+                                            selection_function, 
+                                            selection_kwargs, 
+                                            proximity_threshold,
+                                            padding_radius,
+                                            run_number,
+                                            detector_db)
+
+    return pyrrha_sipm_selection
+
+
 def calibrate_with_mean(dbfile, run_number):
     DataSiPM   = load_db.DataSiPM(dbfile, run_number)
     adc_to_pes = np.abs(DataSiPM.adc_to_pes.values)
diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 7bf3eaef7..a272b3b08 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -269,15 +269,15 @@ def apply_circular_padding(selected_ids_no_isolated : np.ndarray,
     return sipm_ids_with_signal
 
 
-def make_sipm_selection(wfs                 : np.ndarray,
-                        selection_func      : Callable,
-                        selection_kwargs    : dict,
-                        proximity_threshold : float,
-                        padding_radius      : float,
-                        run_number          : int,
-                        detector_db         : str) -> np.ndarray:
+def spatial_selection_method(wfs                 : np.ndarray,
+                             selection_func      : Callable,
+                             selection_kwargs    : dict,
+                             proximity_threshold : float,
+                             padding_radius      : float,
+                             run_number          : int,
+                             detector_db         : str) -> np.ndarray:
     """
-    SiPM selection pipeline, applies SiPM cuts based on user input.
+    SiPM selection function, applies SiPM cuts based on user input.
     A first selection of SiPMs is made, isolated SiPMs are removed
     and padding is added around the SiPMs that are left.
 

From dfd3ba0f1e46136fd8604e70b5e8b4c13d5e9f49 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 4 Mar 2026 15:35:19 +0000
Subject: [PATCH 28/49] modify `spatial_selection_method` to return the ids of
 the selected SiPMs and their waveforms

---
 invisible_cities/reco/wfm_functions.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index a272b3b08..7e5f75ba9 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -293,7 +293,8 @@ def spatial_selection_method(wfs                 : np.ndarray,
 
     Returns
     -------
-    sipm_ids_with_signal : Array of shape (n_sipms,) with boolean values indicating which SiPMs are selected.
+    selected_ids : Array of shape (n_sipms,) containing the indices of the selected SiPMs.
+    selected_wfs : 2D array of shape (n_selected_sipms, n_time_bins) with the waveforms of the selected SiPMs.
     """
     detector_info = load_db.DataSiPM(detector_db, run_number)
     sipm_x = np.array(detector_info.X)
@@ -314,5 +315,9 @@ def spatial_selection_method(wfs                 : np.ndarray,
         sipm_y,
         padding_radius
     )
-    return sipm_ids_with_signal
+
+    selected_ids = np.where(sipm_ids_with_signal)[0]
+    selected_wfs = wfs[selected_ids]
+
+    return selected_ids, selected_wfs
 

From f901531024723c5a04bae224dd0d1c8dbe9c7877 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 4 Mar 2026 17:34:19 +0000
Subject: [PATCH 29/49] implement waveform index slicing in
 `spatial_selection_method`

---
 invisible_cities/cities/components.py  |  1 +
 invisible_cities/reco/wfm_functions.py | 18 +++++++++++-------
 2 files changed, 12 insertions(+), 7 deletions(-)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index 3b3c1c3a3..1d0771135 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -848,6 +848,7 @@ def pyrrha_sipm_selection(selection_function : Callable
     '''
     def pyrrha_sipm_selection(wfs, indices):
         return wfm.spatial_selection_method(wfs, 
+                                            indices,
                                             selection_function, 
                                             selection_kwargs, 
                                             proximity_threshold,
diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 7e5f75ba9..22e8de9c1 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -270,6 +270,7 @@ def apply_circular_padding(selected_ids_no_isolated : np.ndarray,
 
 
 def spatial_selection_method(wfs                 : np.ndarray,
+                             indices         : np.ndarray,
                              selection_func      : Callable,
                              selection_kwargs    : dict,
                              proximity_threshold : float,
@@ -300,24 +301,27 @@ def spatial_selection_method(wfs                 : np.ndarray,
     sipm_x = np.array(detector_info.X)
     sipm_y = np.array(detector_info.Y)
 
-    selected_ids = selection_func(wfs, **selection_kwargs)
+    slice_ = slice(indices[0], indices[-1] + 1)
+    wfs_   = wfs[:, slice_]
 
-    selected_ids_no_isolated = kill_isolated_sipms(
-        selected_ids,
+    starting_ids_ = selection_func(wfs_, **selection_kwargs)
+
+    selected_ids_no_isolated_ = kill_isolated_sipms(
+        starting_ids_,
         sipm_x,
         sipm_y,
         proximity_threshold
     )
 
-    sipm_ids_with_signal = apply_circular_padding(
-        selected_ids_no_isolated,
+    sipm_ids_with_signal_ = apply_circular_padding(
+        selected_ids_no_isolated_,
         sipm_x,
         sipm_y,
         padding_radius
     )
 
-    selected_ids = np.where(sipm_ids_with_signal)[0]
-    selected_wfs = wfs[selected_ids]
+    selected_ids = np.where(sipm_ids_with_signal_)[0]
+    selected_wfs = wfs_[selected_ids]
 
     return selected_ids, selected_wfs
 

From 6eb7b88443067a82783d1380d6d534ed29444773 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 4 Mar 2026 17:34:52 +0000
Subject: [PATCH 30/49] ammend function to return sliced waveforms

---
 invisible_cities/reco/peak_functions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/invisible_cities/reco/peak_functions.py b/invisible_cities/reco/peak_functions.py
index 1a0d829ef..7b64d81e8 100644
--- a/invisible_cities/reco/peak_functions.py
+++ b/invisible_cities/reco/peak_functions.py
@@ -34,7 +34,7 @@ def select_wf_slices_above_time_integrated_thr(wfs, indices, thr):
     wfs_   = wfs[:, slice_]
 
     selected_ids = np.where(np.sum(wfs_, axis = 1) >= thr)[0]
-    selected_wfs = wfs[selected_ids]
+    selected_wfs = wfs_[selected_ids]
 
     return selected_ids, selected_wfs
 

From 16147f5cc85507f5e79eaf3ab4cdfc7a319be4c5 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Wed, 4 Mar 2026 17:41:00 +0000
Subject: [PATCH 31/49] revert back to full waveform

---
 invisible_cities/reco/peak_functions.py | 2 +-
 invisible_cities/reco/wfm_functions.py  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/invisible_cities/reco/peak_functions.py b/invisible_cities/reco/peak_functions.py
index 7b64d81e8..1a0d829ef 100644
--- a/invisible_cities/reco/peak_functions.py
+++ b/invisible_cities/reco/peak_functions.py
@@ -34,7 +34,7 @@ def select_wf_slices_above_time_integrated_thr(wfs, indices, thr):
     wfs_   = wfs[:, slice_]
 
     selected_ids = np.where(np.sum(wfs_, axis = 1) >= thr)[0]
-    selected_wfs = wfs_[selected_ids]
+    selected_wfs = wfs[selected_ids]
 
     return selected_ids, selected_wfs
 
diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 22e8de9c1..0e9e879a6 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -321,7 +321,7 @@ def spatial_selection_method(wfs                 : np.ndarray,
     )
 
     selected_ids = np.where(sipm_ids_with_signal_)[0]
-    selected_wfs = wfs_[selected_ids]
+    selected_wfs = wfs[selected_ids]
 
     return selected_ids, selected_wfs
 

From 146013ec75558a254be4ed5f9043560263b7e8b9 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Thu, 5 Mar 2026 15:53:31 +0000
Subject: [PATCH 32/49] implement `SiPMSelectionMethod` class to call specific
 selection methods for pyrrha

---
 invisible_cities/types/symbols.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/invisible_cities/types/symbols.py b/invisible_cities/types/symbols.py
index c8c43b3da..43a650493 100644
--- a/invisible_cities/types/symbols.py
+++ b/invisible_cities/types/symbols.py
@@ -134,6 +134,9 @@ class CutAlgo(AutoNameEnumBase):
     threshold  = auto()
     pyrrha     = auto()
 
+class SiPMSelectionMethod(AutoNameEnumBase):
+    median_std_method = auto()
+
 class WfType(AutoNameEnumBase):
     rwf  = auto()
     mcrd = auto()

From 3fd5f14ebc58ffc9faa0aeef99bfa66a6bcd7a21 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Thu, 5 Mar 2026 15:54:17 +0000
Subject: [PATCH 33/49] modify `pyrrha_sipm_selection` and
 `spatial_selection_method` to use SiPMSelectionMethod class

---
 invisible_cities/cities/components.py  |  9 +++++----
 invisible_cities/reco/wfm_functions.py | 16 +++++++++++-----
 2 files changed, 16 insertions(+), 9 deletions(-)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index 1d0771135..578101d8c 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -84,6 +84,7 @@
 from .. types  .ic_types          import         types_dict_tracks
 from .. types  .symbols           import                    WfType
 from .. types  .symbols           import                   CutAlgo
+from .. types  .symbols           import       SiPMSelectionMethod
 from .. types  .symbols           import               RebinMethod
 from .. types  .symbols           import                SiPMCharge
 from .. types  .symbols           import                   BlsMode
@@ -836,12 +837,12 @@ def threshold_sipm_selection(wfs, indices):
     return threshold_sipm_selection
 
 
-def pyrrha_sipm_selection(selection_function : Callable
+def pyrrha_sipm_selection(selection_method     : SiPMSelectionMethod
                          , selection_kwargs    : dict
                          , proximity_threshold : float
                          , padding_radius      : float
-                         , run_number : int
-                         , detector_db : str):
+                         , run_number          : int
+                         , detector_db         : str):
     '''
     Function that applies a generic selection function to the sipms, which can be used to 
     implement a spatial SiPM selection method (called Pyrrha).
@@ -849,7 +850,7 @@ def pyrrha_sipm_selection(selection_function : Callable
     def pyrrha_sipm_selection(wfs, indices):
         return wfm.spatial_selection_method(wfs, 
                                             indices,
-                                            selection_function, 
+                                            selection_method, 
                                             selection_kwargs, 
                                             proximity_threshold,
                                             padding_radius,
diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 0e9e879a6..123aa8854 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -13,6 +13,8 @@
 from .. database            import load_db
 from .. reco.peak_functions import select_wf_slices_above_time_integrated_thr
 
+from .. types  .symbols     import SiPMSelectionMethod
+
 def to_adc(wfs, adc_to_pes):
     """
     Convert waveform in pes to adc.
@@ -270,8 +272,8 @@ def apply_circular_padding(selected_ids_no_isolated : np.ndarray,
 
 
 def spatial_selection_method(wfs                 : np.ndarray,
-                             indices         : np.ndarray,
-                             selection_func      : Callable,
+                             indices             : np.ndarray,
+                             selection_method    : SiPMSelectionMethod,
                              selection_kwargs    : dict,
                              proximity_threshold : float,
                              padding_radius      : float,
@@ -285,8 +287,8 @@ def spatial_selection_method(wfs                 : np.ndarray,
     Parameters
     ----------
     wfs                 : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
-    selection_func      : Function returning a boolean mask of energetic SiPMs.
-    selection_kwargs    : Dictionary of arguments passed to selection_func.
+    selection_method    : Method used to select SiPMs.
+    selection_kwargs    : Dictionary of arguments passed to the selection function.
     proximity_threshold : Threshold used to identify isolated SiPMs.
     padding_radius      : Radial padding added to each SiPM that passes the selections.
     run_number          : Run number used to load the detector database.
@@ -304,7 +306,11 @@ def spatial_selection_method(wfs                 : np.ndarray,
     slice_ = slice(indices[0], indices[-1] + 1)
     wfs_   = wfs[:, slice_]
 
-    starting_ids_ = selection_func(wfs_, **selection_kwargs)
+    if selection_method is SiPMSelectionMethod.median_std_method:
+        starting_ids_ = median_std_method(wfs_, **selection_kwargs)
+    else:
+        # temporary solution, think of a nicer method
+        starting_ids_ = median_std_method(wfs_, **selection_kwargs)
 
     selected_ids_no_isolated_ = kill_isolated_sipms(
         starting_ids_,

From 83fe91471fe1767558d03af8b09023ef77187a17 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Thu, 5 Mar 2026 15:58:12 +0000
Subject: [PATCH 34/49] include another method for initial SiPM selection and
 raise an error for incorrect input

---
 invisible_cities/reco/wfm_functions.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 123aa8854..4c9cb70d2 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -308,9 +308,10 @@ def spatial_selection_method(wfs                 : np.ndarray,
 
     if selection_method is SiPMSelectionMethod.median_std_method:
         starting_ids_ = median_std_method(wfs_, **selection_kwargs)
+    elif selection_method is SiPMSelectionMethod.top_n_method:
+        starting_ids_ = top_n_method(wfs_, **selection_kwargs)
     else:
-        # temporary solution, think of a nicer method
-        starting_ids_ = median_std_method(wfs_, **selection_kwargs)
+        raise ValueError(f"Selection method {selection_method} not recognized.")
 
     selected_ids_no_isolated_ = kill_isolated_sipms(
         starting_ids_,

From 300dfb118fbe6dc13ba0f3f1b2f86cca6bccce35 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Tue, 10 Mar 2026 17:33:52 +0000
Subject: [PATCH 35/49] remove time index dependency of SiPM cutting functions

---
 invisible_cities/cities/components.py   | 11 +++++------
 invisible_cities/reco/peak_functions.py | 10 ++++------
 invisible_cities/reco/wfm_functions.py  | 13 ++++---------
 3 files changed, 13 insertions(+), 21 deletions(-)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index 578101d8c..7c08f5eff 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -809,8 +809,8 @@ def apply_cutting_function(algo, **cutting_params):
         # temporary solution, think of a nicer method
         func = threshold_sipm_selection(**cutting_params)
 
-    def apply_cutting_function(wfs, indices):
-        return func(wfs, indices)
+    def apply_cutting_function(wfs):
+        return func(wfs)
 
     return apply_cutting_function
 
@@ -831,8 +831,8 @@ def threshold_sipm_selection(thr_sipm_type
         # extract sipm threshold
         sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number)
 
-    def threshold_sipm_selection(wfs, indices):
-        return wfm.charge_threshold_method(wfs, indices, zeroing_thr = sipm_thr, integration_thr=thr_sipm_s2)
+    def threshold_sipm_selection(wfs):
+        return wfm.charge_threshold_method(wfs, zeroing_thr = sipm_thr, integration_thr=thr_sipm_s2)
 
     return threshold_sipm_selection
 
@@ -847,9 +847,8 @@ def pyrrha_sipm_selection(selection_method     : SiPMSelectionMethod
     Function that applies a generic selection function to the sipms, which can be used to 
     implement a spatial SiPM selection method (called Pyrrha).
     '''
-    def pyrrha_sipm_selection(wfs, indices):
+    def pyrrha_sipm_selection(wfs):
         return wfm.spatial_selection_method(wfs, 
-                                            indices,
                                             selection_method, 
                                             selection_kwargs, 
                                             proximity_threshold,
diff --git a/invisible_cities/reco/peak_functions.py b/invisible_cities/reco/peak_functions.py
index 1a0d829ef..d945fc19c 100644
--- a/invisible_cities/reco/peak_functions.py
+++ b/invisible_cities/reco/peak_functions.py
@@ -92,18 +92,16 @@ def build_pmt_responses(indices, times, widths, ccwf,
 
 def build_sipm_responses(indices, times, widths,
                          sipm_wfs, rebin_stride, apply_cut):
-
+    _, _, sipm_wfs = pick_slice_and_rebin(indices , times, widths,
+                                           sipm_wfs, rebin_stride,
+                                           pad_zeros = False)
     if apply_cut is not None:
         # apply cut before slicing and rebinning
         (sipm_ids,
-         sipm_wfs)   = apply_cut(sipm_wfs, indices)
+         sipm_wfs)   = apply_cut(sipm_wfs)
     else:
         # give all sipm ids as index if no cut is applied
         sipm_ids     = np.arange(sipm_wfs.shape[0])
-    # rebin
-    _, _, sipm_wfs = pick_slice_and_rebin(indices , times, widths,
-                                           sipm_wfs, rebin_stride,
-                                           pad_zeros = False)
 
     return SiPMResponses(sipm_ids, sipm_wfs)
 
diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 4c9cb70d2..53e8725af 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -11,7 +11,7 @@
 from .. calib               import calib_sensors_functions as csf
 from .. sierpe              import blr
 from .. database            import load_db
-from .. reco.peak_functions import select_wf_slices_above_time_integrated_thr
+from .. reco.peak_functions import select_wfs_above_time_integrated_thr
 
 from .. types  .symbols     import SiPMSelectionMethod
 
@@ -159,7 +159,6 @@ def median_std_method(wfs    : np.ndarray,
 
 
 def charge_threshold_method(wfs             : np.ndarray,
-                            indices         : np.ndarray,
                             zeroing_thr     : Optional[float] = 2.,
                             integration_thr : Optional[float] = 5.) -> Tuple[np.ndarray, np.ndarray]:
     """
@@ -182,7 +181,7 @@ def charge_threshold_method(wfs             : np.ndarray,
     zwfs = np.where(wfs > thr, wfs, 0)
 
     # returns selected ids and waveforms above integral
-    return select_wf_slices_above_time_integrated_thr(zwfs, indices, integration_thr)
+    return select_wfs_above_time_integrated_thr(zwfs, integration_thr)
 
 
 def top_n_method(wfs : np.ndarray,
@@ -272,7 +271,6 @@ def apply_circular_padding(selected_ids_no_isolated : np.ndarray,
 
 
 def spatial_selection_method(wfs                 : np.ndarray,
-                             indices             : np.ndarray,
                              selection_method    : SiPMSelectionMethod,
                              selection_kwargs    : dict,
                              proximity_threshold : float,
@@ -303,13 +301,10 @@ def spatial_selection_method(wfs                 : np.ndarray,
     sipm_x = np.array(detector_info.X)
     sipm_y = np.array(detector_info.Y)
 
-    slice_ = slice(indices[0], indices[-1] + 1)
-    wfs_   = wfs[:, slice_]
-
     if selection_method is SiPMSelectionMethod.median_std_method:
-        starting_ids_ = median_std_method(wfs_, **selection_kwargs)
+        starting_ids_ = median_std_method(wfs, **selection_kwargs)
     elif selection_method is SiPMSelectionMethod.top_n_method:
-        starting_ids_ = top_n_method(wfs_, **selection_kwargs)
+        starting_ids_ = top_n_method(wfs, **selection_kwargs)
     else:
         raise ValueError(f"Selection method {selection_method} not recognized.")
 

From ab88fdbda26e60020102e04ebef3927702213228 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Tue, 10 Mar 2026 17:34:34 +0000
Subject: [PATCH 36/49] remove `select_wf_slices_above_time_integrated_thr`

---
 invisible_cities/reco/peak_functions.py | 14 --------------
 1 file changed, 14 deletions(-)

diff --git a/invisible_cities/reco/peak_functions.py b/invisible_cities/reco/peak_functions.py
index d945fc19c..c440cda96 100644
--- a/invisible_cities/reco/peak_functions.py
+++ b/invisible_cities/reco/peak_functions.py
@@ -25,20 +25,6 @@ def indices_and_wf_above_threshold(wf, thr):
     return ZsWf(indices_above_thr, wf_above_thr)
 
 
-def select_wf_slices_above_time_integrated_thr(wfs, indices, thr):
-    '''
-    function that integrates over certain time slices, and passes the waveform
-    based on the aforementioned slices passing the threshold
-    '''
-    slice_ = slice(indices[0], indices[-1] + 1)
-    wfs_   = wfs[:, slice_]
-
-    selected_ids = np.where(np.sum(wfs_, axis = 1) >= thr)[0]
-    selected_wfs = wfs[selected_ids]
-
-    return selected_ids, selected_wfs
-
-
 def select_wfs_above_time_integrated_thr(wfs, thr):
     selected_ids = np.where(np.sum(wfs, axis=1) >= thr)[0]
     selected_wfs = wfs[selected_ids]

From 76b961c2b2bb7e7626dc83ffc9958d6f6bb73327 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Tue, 10 Mar 2026 17:35:14 +0000
Subject: [PATCH 37/49] update docstrings

---
 invisible_cities/reco/wfm_functions.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index 53e8725af..a5977aeab 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -168,8 +168,9 @@ def charge_threshold_method(wfs             : np.ndarray,
 
     Parameters
     ----------
-    wfs       : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
-    threshold : Charge threshold in PE, default 5.
+    wfs             : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
+    zeroing_thr     : Charge threshold for zero suppression in PE, default 2.
+    integration_thr : Charge threshold for total SiPM waveform in PE, default 5.
 
     Returns
     -------

From b7f45d4ac7da6a173c7d7aecc56809c1852eecdf Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Thu, 7 May 2026 16:22:42 +0100
Subject: [PATCH 38/49] remove unused variable from hypathia

---
 invisible_cities/cities/hypathia.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/invisible_cities/cities/hypathia.py b/invisible_cities/cities/hypathia.py
index 08bbd3be8..68f038c6c 100644
--- a/invisible_cities/cities/hypathia.py
+++ b/invisible_cities/cities/hypathia.py
@@ -84,8 +84,6 @@ def hypathia( files_in         : OneOrManyFiles
             , cutting_params   : Dict[str, Any]
             ):
 
-    sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number)
-
     #### Define data transformations
     sd = sensor_data(files_in[0], WfType.mcrd)
 

From be156fe871fc5cecc02c01dda42bf2dac0034be7 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Thu, 7 May 2026 16:53:57 +0100
Subject: [PATCH 39/49] raise error for incorrect cutting function input

---
 invisible_cities/cities/components.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index 7c08f5eff..6a69741e1 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -806,8 +806,7 @@ def apply_cutting_function(algo, **cutting_params):
     elif algo is CutAlgo.pyrrha:
         func = pyrrha_sipm_selection(**cutting_params)
     else:
-        # temporary solution, think of a nicer method
-        func = threshold_sipm_selection(**cutting_params)
+        raise ValueError(f"Unsupported cutting algorithm: {algo!r}. Expected one of {list(CutAlgo)}")
 
     def apply_cutting_function(wfs):
         return func(wfs)

From 6dd6bab5e801c3c123d67868d0449b6a38e57d9c Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Thu, 7 May 2026 18:26:43 +0100
Subject: [PATCH 40/49] add fixture and test for `charge_threshold_method()`

---
 invisible_cities/reco/wfm_functions_test.py | 33 +++++++++++++++++++++
 1 file changed, 33 insertions(+)

diff --git a/invisible_cities/reco/wfm_functions_test.py b/invisible_cities/reco/wfm_functions_test.py
index 35ab09a88..6b5c2c468 100644
--- a/invisible_cities/reco/wfm_functions_test.py
+++ b/invisible_cities/reco/wfm_functions_test.py
@@ -5,6 +5,7 @@
 import tables as tb
 
 from pytest import mark
+from pytest import fixture
 
 from .. database import load_db
 
@@ -49,3 +50,35 @@ def test_compare_cwf_blr(dbnew, ICDATADIR):
                                    event_list=range(NEVT), window_size=300)
 
     assert max(diff) < 0.15
+
+@fixture
+def sipm_wfs_for_sipm_energy_selection_testing():
+    """
+    2D array of SiPM waveforms. All but two SiPMs have 100 PEs integrated charge.
+    Two outliers:
+    - SiPM 2: 150 PEs, more than 3 sigma above the median
+    - SiPM 7: 120 PEs, between 1-2 sigma above the median
+    """
+    n_sipms     = 10
+    n_time_bins = 100
+
+    wfs = np.ones((n_sipms, n_time_bins), dtype=np.float32)
+
+    wfs[2, :] = 1.5
+    wfs[7, :] = 1.2
+    # median ~ 100, std ~ 15.5
+    return wfs, [2, 7], [2]
+
+
+def test_median_std_method(sipm_wfs_for_sipm_energy_selection_testing):
+    """
+    Test function median_std_method(). The test asserts that the function correctly 
+    identifies the outliers based on different standard deviation thresholds. 
+    """
+    wfs, expected_outliers_1sigma, expected_outliers_3sigma = sipm_wfs_for_sipm_energy_selection_testing
+
+    passing_sipms_1sigma = np.where(wfm.median_std_method(wfs, nsigma=1))[0].tolist()
+    passing_sipms_3sigma = np.where(wfm.median_std_method(wfs, nsigma=3))[0].tolist()
+
+    assert passing_sipms_1sigma == expected_outliers_1sigma
+    assert passing_sipms_3sigma == expected_outliers_3sigma
\ No newline at end of file

From b920832db019da8606f86fe56aa71bb6f15df414 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Thu, 7 May 2026 18:32:27 +0100
Subject: [PATCH 41/49] add test for `charge_threshold_method()`

---
 invisible_cities/reco/wfm_functions_test.py | 16 +++++++++++++++-
 1 file changed, 15 insertions(+), 1 deletion(-)

diff --git a/invisible_cities/reco/wfm_functions_test.py b/invisible_cities/reco/wfm_functions_test.py
index 6b5c2c468..2eb0d20a3 100644
--- a/invisible_cities/reco/wfm_functions_test.py
+++ b/invisible_cities/reco/wfm_functions_test.py
@@ -81,4 +81,18 @@ def test_median_std_method(sipm_wfs_for_sipm_energy_selection_testing):
     passing_sipms_3sigma = np.where(wfm.median_std_method(wfs, nsigma=3))[0].tolist()
 
     assert passing_sipms_1sigma == expected_outliers_1sigma
-    assert passing_sipms_3sigma == expected_outliers_3sigma
\ No newline at end of file
+    assert passing_sipms_3sigma == expected_outliers_3sigma
+
+
+def test_threshold_method(sipm_wfs_for_sipm_energy_selection_testing):
+    """
+    Test function threshold_method(). The test asserts that the function correctly
+    kills the SiPMs below a certain charge threshold. 
+    """
+    wfs, expected_outliers_110pes, expected_outliers_130pes = sipm_wfs_for_sipm_energy_selection_testing
+
+    passing_sipms_110pes, _ = wfm.charge_threshold_method(wfs, zeroing_thr=0, integration_thr=110)
+    passing_sipms_130pes, _ = wfm.charge_threshold_method(wfs, zeroing_thr=0, integration_thr=130)
+
+    assert passing_sipms_110pes.tolist() == expected_outliers_110pes
+    assert passing_sipms_130pes.tolist() == expected_outliers_130pes
\ No newline at end of file

From 75aa5792ec5e96f95656194bb01fda8d2be484ba Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Sat, 9 May 2026 02:31:18 +0100
Subject: [PATCH 42/49] simplify naming

---
 invisible_cities/reco/wfm_functions_test.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions_test.py b/invisible_cities/reco/wfm_functions_test.py
index 2eb0d20a3..e6dce88bf 100644
--- a/invisible_cities/reco/wfm_functions_test.py
+++ b/invisible_cities/reco/wfm_functions_test.py
@@ -52,7 +52,7 @@ def test_compare_cwf_blr(dbnew, ICDATADIR):
     assert max(diff) < 0.15
 
 @fixture
-def sipm_wfs_for_sipm_energy_selection_testing():
+def sipm_wfs_for_sipm_selection_testing():
     """
     2D array of SiPM waveforms. All but two SiPMs have 100 PEs integrated charge.
     Two outliers:
@@ -70,12 +70,12 @@ def sipm_wfs_for_sipm_energy_selection_testing():
     return wfs, [2, 7], [2]
 
 
-def test_median_std_method(sipm_wfs_for_sipm_energy_selection_testing):
+def test_median_std_method(sipm_wfs_for_sipm_selection_testing):
     """
     Test function median_std_method(). The test asserts that the function correctly 
     identifies the outliers based on different standard deviation thresholds. 
     """
-    wfs, expected_outliers_1sigma, expected_outliers_3sigma = sipm_wfs_for_sipm_energy_selection_testing
+    wfs, expected_outliers_1sigma, expected_outliers_3sigma = sipm_wfs_for_sipm_selection_testing
 
     passing_sipms_1sigma = np.where(wfm.median_std_method(wfs, nsigma=1))[0].tolist()
     passing_sipms_3sigma = np.where(wfm.median_std_method(wfs, nsigma=3))[0].tolist()
@@ -84,12 +84,12 @@ def test_median_std_method(sipm_wfs_for_sipm_energy_selection_testing):
     assert passing_sipms_3sigma == expected_outliers_3sigma
 
 
-def test_threshold_method(sipm_wfs_for_sipm_energy_selection_testing):
+def test_threshold_method(sipm_wfs_for_sipm_selection_testing):
     """
     Test function threshold_method(). The test asserts that the function correctly
     kills the SiPMs below a certain charge threshold. 
     """
-    wfs, expected_outliers_110pes, expected_outliers_130pes = sipm_wfs_for_sipm_energy_selection_testing
+    wfs, expected_outliers_110pes, expected_outliers_130pes = sipm_wfs_for_sipm_selection_testing
 
     passing_sipms_110pes, _ = wfm.charge_threshold_method(wfs, zeroing_thr=0, integration_thr=110)
     passing_sipms_130pes, _ = wfm.charge_threshold_method(wfs, zeroing_thr=0, integration_thr=130)

From 5ce0864727803a9cf02124dd9b2a950deef63e90 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Sat, 9 May 2026 02:37:39 +0100
Subject: [PATCH 43/49] add test for `top_n_method`

---
 invisible_cities/reco/wfm_functions_test.py | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/invisible_cities/reco/wfm_functions_test.py b/invisible_cities/reco/wfm_functions_test.py
index e6dce88bf..a4884224f 100644
--- a/invisible_cities/reco/wfm_functions_test.py
+++ b/invisible_cities/reco/wfm_functions_test.py
@@ -95,4 +95,21 @@ def test_threshold_method(sipm_wfs_for_sipm_selection_testing):
     passing_sipms_130pes, _ = wfm.charge_threshold_method(wfs, zeroing_thr=0, integration_thr=130)
 
     assert passing_sipms_110pes.tolist() == expected_outliers_110pes
-    assert passing_sipms_130pes.tolist() == expected_outliers_130pes
\ No newline at end of file
+    assert passing_sipms_130pes.tolist() == expected_outliers_130pes
+
+
+def test_top_n_method(sipm_wfs_for_sipm_selection_testing):
+    """
+    Test function top_n_method(). The test asserts that the function  selects the correct
+    number of SiPMs andcorrectly identifies the top N SiPMs based on their integrated charge. 
+    """
+    wfs, expected_outliers_top2, expected_outliers_top1 = sipm_wfs_for_sipm_selection_testing
+
+    passing_sipms_top2 = np.where(wfm.top_n_method(wfs, n=2))[0].tolist()
+    passing_sipms_top1 = np.where(wfm.top_n_method(wfs, n=1))[0].tolist()
+
+    assert len(passing_sipms_top2) == 2
+    assert len(passing_sipms_top1) == 1
+
+    assert passing_sipms_top2 == expected_outliers_top2
+    assert passing_sipms_top1 == expected_outliers_top1
\ No newline at end of file

From ebbe61ad8020241de2b3f9d8cb62920d13ee0379 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Sun, 10 May 2026 01:09:37 +0100
Subject: [PATCH 44/49] add fixture and test for `kill_isolated_sipms()`

---
 invisible_cities/reco/wfm_functions_test.py | 60 ++++++++++++++++++++-
 1 file changed, 59 insertions(+), 1 deletion(-)

diff --git a/invisible_cities/reco/wfm_functions_test.py b/invisible_cities/reco/wfm_functions_test.py
index a4884224f..3a9acaa0e 100644
--- a/invisible_cities/reco/wfm_functions_test.py
+++ b/invisible_cities/reco/wfm_functions_test.py
@@ -112,4 +112,62 @@ def test_top_n_method(sipm_wfs_for_sipm_selection_testing):
     assert len(passing_sipms_top1) == 1
 
     assert passing_sipms_top2 == expected_outliers_top2
-    assert passing_sipms_top1 == expected_outliers_top1
\ No newline at end of file
+    assert passing_sipms_top1 == expected_outliers_top1
+
+
+@fixture
+def sipm_grid_for_isolation_testing():
+    """
+    5x5 grid of with 10mm spacing containing 4 SiPMs:
+    - SiPM  7 (x=20, y=10) has nearest neighbours 12 and 13
+    - SiPM 12 (x=20, y=20) has nearest neighbours 7 and 13
+    - SiPM 13 (x=30, y=20) has nearestneighbours 12 and 7
+    - SiPM 24 (x=40, y=40) has no nearest neighbours
+    With proximity_threshold=15mm SiPMs 7, 12, 13 survive and SiPM 24 is killed.
+    With proximity_threshold=5mm, all SiPMs are killed.
+    """
+    spacing = 10
+    xs = np.arange(5) * spacing  # [0, 10, 20, 30, 40]
+    ys = np.arange(5) * spacing
+
+    grid_x, grid_y = np.meshgrid(xs, ys)
+    sipm_x = grid_x.flatten().astype(np.float32)  # shape (25,)
+    sipm_y = grid_y.flatten().astype(np.float32)
+
+    # visual representation of the grid with SiPMs plotted as their IDs:
+    # X  X  X  X  24
+    # X  X  X  X  X
+    # X  X  12 13 X
+    # X  X  7  X  X
+    # X  X  X  X  X
+    selected_ids = np.zeros(25, dtype=bool)
+    cluster_ids  = [7, 12, 13]
+    isolated_id  = [24]
+    for idx in cluster_ids + isolated_id:
+        selected_ids[idx] = True
+
+    # diagonal SiPM distance ~ 14mm
+    # only SiPMs with a nearest neighbour pass assuming proximity_threshold=15mm
+    # SiPM 24 should be killed
+    expected_survivors_15mm = np.zeros(25, dtype=bool)
+    for idx in cluster_ids:
+        expected_survivors_15mm[idx] = True
+
+    # given a SiPM distance of 10mm, there should be no survivors with proximity_threshold=5mm
+    expected_survivors_5mm = np.zeros(25, dtype=bool)
+
+    return sipm_x, sipm_y, selected_ids, expected_survivors_15mm, expected_survivors_5mm
+
+
+def test_kill_isolated_sipms(sipm_grid_for_isolation_testing):
+    """"
+    Test function kill_isolated_sipms(). The test asserts that the function correctly
+    identifies and removes isolated SiPMs based on their proximity to other SiPMs.
+    """
+    sipm_x, sipm_y, selected_ids, expected_survivors_15mm, expected_survivors_5mm = sipm_grid_for_isolation_testing
+    
+    surviving_sipms_15mm = wfm.kill_isolated_sipms(selected_ids, sipm_x, sipm_y, proximity_threshold=15.0)
+    surviving_sipms_5mm = wfm.kill_isolated_sipms(selected_ids, sipm_x, sipm_y, proximity_threshold=5.0)
+
+    assert surviving_sipms_15mm.tolist() == expected_survivors_15mm.tolist()
+    assert surviving_sipms_5mm.tolist() == expected_survivors_5mm.tolist()
\ No newline at end of file

From 51c2e630a3471da3164a030c6ca9e53eebc9cbb9 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Sun, 10 May 2026 11:49:38 +0100
Subject: [PATCH 45/49] amend `apply_circular_padding()` to not kill SiPMs when
 padding is set to 0

---
 invisible_cities/reco/wfm_functions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index a5977aeab..e0576c6e1 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -266,7 +266,7 @@ def apply_circular_padding(selected_ids_no_isolated : np.ndarray,
     for i in np.where(selected_ids_no_isolated)[0]:
         x, y = sipm_x[i], sipm_y[i]
         distances = np.sqrt((sipm_x - x)**2 + (sipm_y - y)**2)
-        sipm_ids_with_signal |= distances < padding_radius
+        sipm_ids_with_signal |= distances <= padding_radius
 
     return sipm_ids_with_signal
 

From 021e6efd967378450b47ca6d51b51f36ae01b1a1 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Sun, 10 May 2026 11:53:14 +0100
Subject: [PATCH 46/49] add test for `apply_circular_padding()` and modify
 fixture to accomodate for it

---
 invisible_cities/reco/wfm_functions_test.py | 45 ++++++++++++++++++---
 1 file changed, 40 insertions(+), 5 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions_test.py b/invisible_cities/reco/wfm_functions_test.py
index 3a9acaa0e..e2c029148 100644
--- a/invisible_cities/reco/wfm_functions_test.py
+++ b/invisible_cities/reco/wfm_functions_test.py
@@ -116,7 +116,7 @@ def test_top_n_method(sipm_wfs_for_sipm_selection_testing):
 
 
 @fixture
-def sipm_grid_for_isolation_testing():
+def sipm_grid_for_isolation_and_padding_testing():
     """
     5x5 grid of with 10mm spacing containing 4 SiPMs:
     - SiPM  7 (x=20, y=10) has nearest neighbours 12 and 13
@@ -125,6 +125,9 @@ def sipm_grid_for_isolation_testing():
     - SiPM 24 (x=40, y=40) has no nearest neighbours
     With proximity_threshold=15mm SiPMs 7, 12, 13 survive and SiPM 24 is killed.
     With proximity_threshold=5mm, all SiPMs are killed.
+
+    With a padding of 12mm around the cluster of SiPMs 7, 12, 13 you also
+    include SiPMs 2, 6, 8, 11, 14, 17, 18.
     """
     spacing = 10
     xs = np.arange(5) * spacing  # [0, 10, 20, 30, 40]
@@ -156,18 +159,50 @@ def sipm_grid_for_isolation_testing():
     # given a SiPM distance of 10mm, there should be no survivors with proximity_threshold=5mm
     expected_survivors_5mm = np.zeros(25, dtype=bool)
 
-    return sipm_x, sipm_y, selected_ids, expected_survivors_15mm, expected_survivors_5mm
+    # adding a padding of 12mm around the cluster of SiPMs 7, 12, 13 would include 
+    # SiPMs 2, 6, 8, 11, 14, 17, 18
+    # X  X  X  X  X
+    # X  X  17 18 X
+    # X  11 12 13 14
+    # X  6  7  8  X
+    # X  X  2  X  X
+    padded_cluster_ids = cluster_ids + [2, 6, 8, 11, 14, 17, 18]
+    expected_survivors_padding12 = np.zeros(25, dtype=bool)
+    for idx in padded_cluster_ids:
+        expected_survivors_padding12[idx] = True
+
+    # adding a padding of 0mm around the cluster of SiPMs 7, 12, 13 would include only the cluster itself
+    expected_survivors_padding0 = np.zeros(25, dtype=bool)
+    for idx in cluster_ids:
+        expected_survivors_padding0[idx] = True
 
+    return (sipm_x, sipm_y, selected_ids, expected_survivors_15mm, expected_survivors_5mm, 
+            expected_survivors_padding12, expected_survivors_padding0)
 
-def test_kill_isolated_sipms(sipm_grid_for_isolation_testing):
+
+def test_kill_isolated_sipms(sipm_grid_for_isolation_and_padding_testing):
     """"
     Test function kill_isolated_sipms(). The test asserts that the function correctly
     identifies and removes isolated SiPMs based on their proximity to other SiPMs.
     """
-    sipm_x, sipm_y, selected_ids, expected_survivors_15mm, expected_survivors_5mm = sipm_grid_for_isolation_testing
+    sipm_x, sipm_y, selected_ids, expected_survivors_15mm, expected_survivors_5mm, _, _ = sipm_grid_for_isolation_and_padding_testing
     
     surviving_sipms_15mm = wfm.kill_isolated_sipms(selected_ids, sipm_x, sipm_y, proximity_threshold=15.0)
     surviving_sipms_5mm = wfm.kill_isolated_sipms(selected_ids, sipm_x, sipm_y, proximity_threshold=5.0)
 
     assert surviving_sipms_15mm.tolist() == expected_survivors_15mm.tolist()
-    assert surviving_sipms_5mm.tolist() == expected_survivors_5mm.tolist()
\ No newline at end of file
+    assert surviving_sipms_5mm.tolist() == expected_survivors_5mm.tolist()
+
+
+def test_apply_circular_padding(sipm_grid_for_isolation_and_padding_testing):
+    """
+    Test function apply_circular_padding(). The test asserts that the function correctly applies 
+    a circular padding around selected SiPMs to include neighboring SiPMs within the specified radius.
+    """
+    sipm_x, sipm_y, _, selected_ids, _, expected_survivors_padding12, expected_survivors_padding0 = sipm_grid_for_isolation_and_padding_testing
+
+    padded_sipms_12mm = wfm.apply_circular_padding(selected_ids, sipm_x, sipm_y, padding_radius=12.0)
+    padded_sipms_0mm = wfm.apply_circular_padding(selected_ids, sipm_x, sipm_y, padding_radius=0.0)
+
+    assert padded_sipms_12mm.tolist() == expected_survivors_padding12.tolist()
+    assert padded_sipms_0mm.tolist() == expected_survivors_padding0.tolist()
\ No newline at end of file

From 58b88aaac82c425e634aa77d554c7984a62089e8 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Mon, 11 May 2026 16:14:57 +0100
Subject: [PATCH 47/49] add `zero_wfs_below_threshold()` to zero wavefrom
 entries below a threshold outside of `charge_threshold_method()`

---
 invisible_cities/reco/wfm_functions.py | 22 +++++++++++++++++++---
 1 file changed, 19 insertions(+), 3 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py
index e0576c6e1..0f2feba26 100644
--- a/invisible_cities/reco/wfm_functions.py
+++ b/invisible_cities/reco/wfm_functions.py
@@ -138,6 +138,24 @@ def compare_cwf_blr(cwf, pmtblr, event_list, window_size=500):
     return np.array(DIFF)
 
 
+def zero_wfs_below_threshold(wfs         : np.ndarray,
+                             zeroing_thr : Optional[float] = 2.) -> np.ndarray:
+    """
+    Zeroes the entries of the input waveforms that are below a given threshold.
+
+    Parameters
+    ----------
+    wfs         : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM.
+    zeroing_thr : Charge threshold for zero suppression in PE, default 2.
+
+    Returns
+    -------
+    2D array of shape (n_sipms, n_time_bins) containing the input waveforms with entries below threshold set to zero.
+    """
+    thr = to_col_vector(np.full(wfs.shape[0], zeroing_thr))
+    return np.where(wfs > thr, wfs, 0)
+
+
 def median_std_method(wfs    : np.ndarray,
                       nsigma : Optional[float] = 3.) -> np.ndarray:
     """
@@ -176,10 +194,8 @@ def charge_threshold_method(wfs             : np.ndarray,
     -------
     Tuple of np arrays including all passing sipm ids and the corresponding waveforms
     """
-    thr = to_col_vector(np.full(wfs.shape[0], zeroing_thr))
-
     # zero entries below threshold
-    zwfs = np.where(wfs > thr, wfs, 0)
+    zwfs = zero_wfs_below_threshold(wfs, zeroing_thr)
 
     # returns selected ids and waveforms above integral
     return select_wfs_above_time_integrated_thr(zwfs, integration_thr)

From bd1b9d95f9d08720f7641a8472a418b8cc6bf4f6 Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Mon, 11 May 2026 16:33:47 +0100
Subject: [PATCH 48/49] update `calibrate_sipms()` and corresponding test it to
 rely on `zero_wfs_below_threshold()` for zeroing waveform entries

---
 .../calib/calib_sensors_functions.py          |  9 ++----
 .../calib/calib_sensors_functions_test.py     | 28 ++++++++++---------
 2 files changed, 17 insertions(+), 20 deletions(-)

diff --git a/invisible_cities/calib/calib_sensors_functions.py b/invisible_cities/calib/calib_sensors_functions.py
index ba2c8854e..38d522f20 100644
--- a/invisible_cities/calib/calib_sensors_functions.py
+++ b/invisible_cities/calib/calib_sensors_functions.py
@@ -136,18 +136,13 @@ def pmt_subtract_maw(cwfs, n_maw=100):
     return cwfs - maw
 
 
-def calibrate_sipms(sipm_wfs, adc_to_pes, *, thr = None, bls_mode=BlsMode.mode):
+def calibrate_sipms(sipm_wfs, adc_to_pes, *, bls_mode=BlsMode.mode):
     """
     Subtract baseline, and calibrates waveforms to pes.
     """
     bls  = subtract_baseline(sipm_wfs, bls_mode=bls_mode)
     cwfs = calibrate_wfs(bls, adc_to_pes)
-    if thr is None:
-        return cwfs
-    else:
-        # if you apply a threshold here, apply as usual
-        thr  = to_col_vector(np.full(sipm_wfs.shape[0], thr))
-        return np.where(cwfs > thr, cwfs, 0)
+    return cwfs
 
 
 def subtract_mean  (wfs): return subtract_baseline(wfs, bls_mode=BlsMode.mean  )
diff --git a/invisible_cities/calib/calib_sensors_functions_test.py b/invisible_cities/calib/calib_sensors_functions_test.py
index bdbfc673e..1b3948561 100644
--- a/invisible_cities/calib/calib_sensors_functions_test.py
+++ b/invisible_cities/calib/calib_sensors_functions_test.py
@@ -13,6 +13,7 @@
 
 from .          import calib_sensors_functions as csf
 from .. core    import core_functions          as cf
+from .. reco    import wfm_functions           as wfm
 from .. sierpe  import fee                     as FE
 from .. sierpe  import waveform_generator      as wfg
 
@@ -207,10 +208,11 @@ def test_calibrate_sipms_stat(oscillating_waveform_with_baseline,
      baseline )              = oscillating_waveform_with_baseline
     #n_maw                    = n_samples // 500
 
-    ccwfs = csf.calibrate_sipms(wfs, adc_to_pes, thr=nsigma * noise_sigma, bls_mode=BlsMode.mode)
+    ccwfs = csf.calibrate_sipms(wfs, adc_to_pes, bls_mode=BlsMode.mode)
+    zeroed_ccwfs = wfm.zero_wfs_below_threshold(ccwfs, zeroing_thr=nsigma * noise_sigma)
 
-    number_of_zeros = np.count_nonzero(ccwfs == 0)
-    assert number_of_zeros > fraction * ccwfs.size
+    number_of_zeros = np.count_nonzero(zeroed_ccwfs == 0)
+    assert number_of_zeros > fraction * zeroed_ccwfs.size
 
 
 def test_calibrate_sipms_common_threshold(toy_sipm_signal):
@@ -218,8 +220,8 @@ def test_calibrate_sipms_common_threshold(toy_sipm_signal):
      signal_zs_common_threshold, _,
      common_threshold, _) = toy_sipm_signal
 
-    zs_wf = csf.calibrate_sipms(signal_adc, adc_to_pes,
-                                thr=common_threshold, bls_mode=BlsMode.mode)
+    ccwf = csf.calibrate_sipms(signal_adc, adc_to_pes, bls_mode=BlsMode.mode)
+    zs_wf = wfm.zero_wfs_below_threshold(ccwf, zeroing_thr=common_threshold)
 
     for actual, expected in zip(zs_wf, signal_zs_common_threshold):
         assert actual == approx(expected)
@@ -230,10 +232,9 @@ def test_calibrate_sipms_individual_thresholds(toy_sipm_signal):
      _, signal_zs_individual_thresholds,
      _, individual_thresholds) = toy_sipm_signal
 
+    ccwf = csf.calibrate_sipms(signal_adc, adc_to_pes, bls_mode=BlsMode.mode)
+    zs_wf = wfm.zero_wfs_below_threshold(ccwf, zeroing_thr=individual_thresholds)
 
-    zs_wf = csf.calibrate_sipms(signal_adc, adc_to_pes,
-                                thr=individual_thresholds,
-                                bls_mode=BlsMode.mode)
     for actual, expected in zip(zs_wf, signal_zs_individual_thresholds):
         assert actual == approx(expected)
 
@@ -309,11 +310,12 @@ def test_area_of_sum_equals_sum_of_areas_pmts(square_pmt_and_sipm_waveforms):
 
 def test_area_of_sum_equals_sum_of_areas_sipms(square_pmt_and_sipm_waveforms):
     _, nsensors, _, _, _, sipms_wfm, _ = square_pmt_and_sipm_waveforms
-    adc_to_pes = np.full(nsensors, 100, dtype=float)
-    cwfs       = csf.calibrate_sipms(sipms_wfm, adc_to_pes, thr=10, bls_mode=BlsMode.mode)
-    stot       = np.sum(cwfs[0]) * nsensors
-    sums       = np.sum(cwfs, axis=1)
-    stot2      = reduce(add, sums)
+    adc_to_pes  = np.full(nsensors, 100, dtype=float)
+    cwfs        = csf.calibrate_sipms(sipms_wfm, adc_to_pes, bls_mode=BlsMode.mode)
+    zeroed_cwfs = wfm.zero_wfs_below_threshold(cwfs, zeroing_thr=10)
+    stot        = np.sum(zeroed_cwfs[0]) * nsensors
+    sums        = np.sum(zeroed_cwfs, axis=1)
+    stot2       = reduce(add, sums)
     assert stot == approx(stot2, rel=1e-3)
 
 

From 66e94a7e34ea9d1328b1d74b9643fc1c1045aa9d Mon Sep 17 00:00:00 2001
From: Ian0sborne <ian.osborne@student.manchester.ac.uk>
Date: Mon, 11 May 2026 16:49:45 +0100
Subject: [PATCH 49/49] add test for `zero_wfs_below_threshold()` with small
 tweak in fixture

---
 invisible_cities/reco/wfm_functions_test.py | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/invisible_cities/reco/wfm_functions_test.py b/invisible_cities/reco/wfm_functions_test.py
index e2c029148..4ef8209b3 100644
--- a/invisible_cities/reco/wfm_functions_test.py
+++ b/invisible_cities/reco/wfm_functions_test.py
@@ -67,7 +67,20 @@ def sipm_wfs_for_sipm_selection_testing():
     wfs[2, :] = 1.5
     wfs[7, :] = 1.2
     # median ~ 100, std ~ 15.5
-    return wfs, [2, 7], [2]
+    return wfs, [2, 7], [2], [0, 1, 3, 4, 5, 6, 8, 9]
+
+
+def test_zero_wfs_below_threshold(sipm_wfs_for_sipm_selection_testing):
+    """
+    Test function zero_wfs_below_threshold(). The test asserts that the function correctly
+    sets to zero the entries of the waveforms that are below a specified threshold.
+    """
+    wfs, passing_wfs_11_ids, _, zeroed_wfs_11_ids  = sipm_wfs_for_sipm_selection_testing
+
+    zeroed_wfs_11 = wfm.zero_wfs_below_threshold(wfs, zeroing_thr=1.1)
+
+    assert np.all(zeroed_wfs_11[zeroed_wfs_11_ids] == 0)
+    assert np.all(zeroed_wfs_11[passing_wfs_11_ids] == wfs[passing_wfs_11_ids])
 
 
 def test_median_std_method(sipm_wfs_for_sipm_selection_testing):
@@ -75,7 +88,7 @@ def test_median_std_method(sipm_wfs_for_sipm_selection_testing):
     Test function median_std_method(). The test asserts that the function correctly 
     identifies the outliers based on different standard deviation thresholds. 
     """
-    wfs, expected_outliers_1sigma, expected_outliers_3sigma = sipm_wfs_for_sipm_selection_testing
+    wfs, expected_outliers_1sigma, expected_outliers_3sigma, _ = sipm_wfs_for_sipm_selection_testing
 
     passing_sipms_1sigma = np.where(wfm.median_std_method(wfs, nsigma=1))[0].tolist()
     passing_sipms_3sigma = np.where(wfm.median_std_method(wfs, nsigma=3))[0].tolist()
@@ -89,7 +102,7 @@ def test_threshold_method(sipm_wfs_for_sipm_selection_testing):
     Test function threshold_method(). The test asserts that the function correctly
     kills the SiPMs below a certain charge threshold. 
     """
-    wfs, expected_outliers_110pes, expected_outliers_130pes = sipm_wfs_for_sipm_selection_testing
+    wfs, expected_outliers_110pes, expected_outliers_130pes, _ = sipm_wfs_for_sipm_selection_testing
 
     passing_sipms_110pes, _ = wfm.charge_threshold_method(wfs, zeroing_thr=0, integration_thr=110)
     passing_sipms_130pes, _ = wfm.charge_threshold_method(wfs, zeroing_thr=0, integration_thr=130)
@@ -103,7 +116,7 @@ def test_top_n_method(sipm_wfs_for_sipm_selection_testing):
     Test function top_n_method(). The test asserts that the function  selects the correct
     number of SiPMs andcorrectly identifies the top N SiPMs based on their integrated charge. 
     """
-    wfs, expected_outliers_top2, expected_outliers_top1 = sipm_wfs_for_sipm_selection_testing
+    wfs, expected_outliers_top2, expected_outliers_top1, _ = sipm_wfs_for_sipm_selection_testing
 
     passing_sipms_top2 = np.where(wfm.top_n_method(wfs, n=2))[0].tolist()
     passing_sipms_top1 = np.where(wfm.top_n_method(wfs, n=1))[0].tolist()