From c39c4d9dafde9dcf2e88443bb1a3fa19b8ecc590 Mon Sep 17 00:00:00 2001
From: Manuel Huber <info@manuelhu.de>
Date: Fri, 13 Feb 2026 16:14:15 +0100
Subject: [PATCH] optmap: add experimental photon_threshold_per_hit

---
 src/reboost/optmap/convolve.py | 15 +++++++++++++--
 src/reboost/spms/pe.py         |  3 ++-
 2 files changed, 15 insertions(+), 3 deletions(-)

diff --git a/src/reboost/optmap/convolve.py b/src/reboost/optmap/convolve.py
index b6c2072d..ab74b2e6 100644
--- a/src/reboost/optmap/convolve.py
+++ b/src/reboost/optmap/convolve.py
@@ -170,6 +170,7 @@ def iterate_stepwise_depositions_numdet(
     det: str,
     map_scaling: float = 1,
     map_scaling_sigma: float = 0,
+    photon_threshold_per_hit: int = -1,
     rng: np.random.Generator | None = None,
 ):
     if edep_hits.xloc.ndim == 1:
@@ -187,6 +188,7 @@ def iterate_stepwise_depositions_numdet(
         optmap.edges,
         optmap.weights,
         ak.sum(counts),
+        photon_threshold_per_hit,
     )
 
     if res["det_no_stats"] > 0:
@@ -367,8 +369,9 @@ def _iterate_stepwise_depositions_numdet(
     optmap_edges,
     optmap_weights,
     output_length: int,
+    photon_threshold_per_hit: int = -1,
 ):
-    oob = ib = det_no_stats = 0
+    oob = ib = det_no_stats = skipped = 0
     output = np.empty(shape=output_length, dtype=np.int64)
 
     output_index = 0
@@ -380,7 +383,14 @@ def _iterate_stepwise_depositions_numdet(
             map_scaling_evt = rng.normal(loc=map_scaling, scale=map_scaling_sigma)
 
         # iterate steps inside the hit
+        photons_in_hit = 0
         for si in range(len(hit.xloc)):
+            if photon_threshold_per_hit > 0 and photons_in_hit > photon_threshold_per_hit:
+                output[output_index] = 0
+                output_index += 1
+                skipped += 1
+                continue
+
             loc = np.array([hit.xloc[si], hit.yloc[si], hit.zloc[si]], dtype=np.float64)
             # coordinates -> bins of the optical map.
             bins = np.empty(3, dtype=np.int64)
@@ -407,9 +417,10 @@ def _iterate_stepwise_depositions_numdet(
             pois_cnt = 0 if detp <= 0.0 else rng.poisson(lam=hit.num_scint_ph[si] * detp)
             output[output_index] = pois_cnt
             output_index += 1
+            photons_in_hit += pois_cnt
 
     assert output_index == output_length
-    return output, {"oob": oob, "ib": ib, "det_no_stats": det_no_stats}
+    return output, {"oob": oob, "ib": ib, "det_no_stats": det_no_stats, "skipped": skipped}
 
 
 # - run with NUMBA_FULL_TRACEBACKS=1 NUMBA_BOUNDSCHECK=1 for testing/checking
diff --git a/src/reboost/spms/pe.py b/src/reboost/spms/pe.py
index bf42e69e..4422fbbf 100644
--- a/src/reboost/spms/pe.py
+++ b/src/reboost/spms/pe.py
@@ -208,6 +208,7 @@ def number_of_detected_photoelectrons(
     spm_detector: str,
     map_scaling: float = 1,
     map_scaling_sigma: float = 0,
+    photon_threshold_per_hit: int = -1,
 ) -> ak.Array:
     """Derive the number of detected photoelectrons (p.e.) from scintillator hits using an optical map.
 
@@ -233,7 +234,7 @@ def number_of_detected_photoelectrons(
     )
 
     return convolve.iterate_stepwise_depositions_numdet(
-        hits, optmap, spm_detector, map_scaling, map_scaling_sigma
+        hits, optmap, spm_detector, map_scaling, map_scaling_sigma, photon_threshold_per_hit
     )