optmap: add experimental photon_threshold_per_hit

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

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
     )