S2 merging exclude sparse peaklets

straxen/plugins/events/events.py

@@ -121,6 +121,7 @@ def get_window_size(self):
     def _is_triggering(self, peaks):
         _is_triggering = peaks["area"] > self.trigger_min_area
         _is_triggering &= peaks["n_competing"] <= self.trigger_max_competing
+        # have to consider the peak with type 20
         if self.exclude_s1_as_triggering_peaks:
             _is_triggering &= peaks["type"] == 2
         else:

straxen/plugins/merged_s2s/merged_s2s.py

@@ -14,9 +14,15 @@ class MergedS2s(strax.OverlapWindowPlugin):
     """Merge together peaklets if peak finding favours that they would form a single peak
     instead."""
 
-    __version__ = "1.1.0"
-
-    depends_on: Tuple[str, ...] = ("peaklets", "peaklet_classification", "lone_hits")
+    __version__ = "2.0.0"
+
+    # Use DEFAULT_POSREC_ALGO here?
+    depends_on: Tuple[str, ...] = (
+        "peaklets",
+        "peaklet_positions_cnf",
+        "peaklet_classification",
+        "lone_hits",
+    )
     data_kind = "merged_s2s"
     provides = "merged_s2s"
 
@@ -52,6 +58,18 @@ class MergedS2s(strax.OverlapWindowPlugin):
         ),
     )
 
+    position_density_s2_area_threshold = straxen.URLConfig(
+        default=1e4,
+        infer_type=False,
+        help="Threshold of S2 area for using the position reconstruction density",
+    )
+
+    dr_cnf_avg_threshold = straxen.URLConfig(
+        default=14,
+        infer_type=False,
+        help="Threshold for the average distance to the center of the top array allowed for S2s",
+    )
+
     n_top_pmts = straxen.URLConfig(type=int, help="Number of top TPC array PMTs")
 
     n_tpc_pmts = straxen.URLConfig(type=int, help="Number of TPC PMTs")
@@ -131,7 +149,6 @@ def compute(self, peaklets, lone_hits):
             end_merge_at,
             max_buffer=int(self.s2_merge_max_duration // np.gcd.reduce(peaklets["dt"])),
         )
-        merged_s2s["type"] = 2
 
         # Updated time and length of lone_hits and sort again:
         lh = np.copy(lone_hits)
@@ -141,20 +158,29 @@ def compute(self, peaklets, lone_hits):
         lh["length"] = lh["right_integration"] - lh["left_integration"]
         lh = strax.sort_by_time(lh)
 
-        _n_top_pmts = self.n_top_pmts if "data_top" in self.dtype.names else -1
+        _store_data_top = "data_top" in self.dtype.names
         _store_data_start = "data_start" in self.dtype.names
         strax.add_lone_hits(
             merged_s2s,
             lh,
             self.to_pe,
-            n_top_channels=_n_top_pmts,
+            n_top_channels=self.n_top_pmts if _store_data_top else -1,
             store_data_start=_store_data_start,
         )
 
         strax.compute_widths(merged_s2s)
 
-        if (_n_top_pmts <= 0) or (not _store_data_start):
-            merged_s2s = drop_data_field(merged_s2s, self.dtype, "_drop_data_field_merged_s2s")
+        merged_s2s["type"] = 2
+        assert self.merge_without_s1
+        area_top = peaklets["area_per_channel"][:, : self.n_top_pmts].sum(axis=1)
+        # by using pos-rec density, merge_without_s1 must be true
+        dr_cnf_avg = self.weighted_average_dr(area_top, peaklets, merged_s2s)
+        mask = merged_s2s["area"] < self.position_density_s2_area_threshold
+        mask &= dr_cnf_avg > self.dr_cnf_avg_threshold
+        merged_s2s["type"][mask] = 20
+
+        # if (not _store_data_top) or (not _store_data_start):
+        merged_s2s = drop_data_field(merged_s2s, self.dtype, "_drop_data_field_merged_s2s")
 
         return merged_s2s
 
@@ -226,6 +252,26 @@ def get_merge_instructions(
 
         return merge_start, merge_stop_exclusive
 
+    def weighted_average_dr(self, area_top, peaklets, merged_s2s):
+        """Weighted average of the distance to the center of the top array for the peaklets."""
+        mask = area_top > 0
+        mask &= ~np.isnan(peaklets["x_cnf"])
+        mask &= ~np.isnan(peaklets["y_cnf"])
+        dr_cnf_avg = np.full(len(merged_s2s), np.nan)
+        windows = strax.touching_windows(peaklets, merged_s2s)
+        for i, window in enumerate(windows):
+            ps = slice(*window)
+            if mask[ps].sum() == 0:
+                continue
+            weights = area_top[ps][mask[ps]]
+            x_cnf = peaklets["x_cnf"][ps][mask[ps]]
+            y_cnf = peaklets["y_cnf"][ps][mask[ps]]
+            x_cnf_avg = np.average(x_cnf, weights=weights)
+            y_cnf_avg = np.average(y_cnf, weights=weights)
+            dr_cnf = np.sqrt((x_cnf - x_cnf_avg) ** 2 + (y_cnf - y_cnf_avg) ** 2)
+            dr_cnf_avg[i] = np.average(dr_cnf, weights=weights)
+        return dr_cnf_avg
+
 
 @numba.njit(cache=True, nogil=True)
 def _filter_s1_starts(start_merge_at, types, end_merge_at):

straxen/plugins/peaks/peak_proximity.py

@@ -12,7 +12,7 @@ class PeakProximity(strax.OverlapWindowPlugin):
     """Look for peaks around a peak to determine how many peaks are in proximity (in time) of a
     peak."""
 
-    __version__ = "0.4.0"
+    __version__ = "0.5.0"
 
     depends_on = "peak_basics"
     dtype = [
@@ -56,8 +56,12 @@ def get_window_size(self):
         return self.peak_max_proximity_time
 
     def compute(self, peaks):
-        windows = strax.touching_windows(peaks, peaks, window=self.nearby_window)
-        n_left, n_tot = self.find_n_competing(peaks, windows, fraction=self.min_area_fraction)
+        # later we can even remove type 0
+        mask = np.isin(peaks["type"], [0, 1, 2])
+        windows = strax.touching_windows(peaks[mask], peaks, window=self.nearby_window)
+        n_left, n_tot = self.find_n_competing(
+            peaks[mask], peaks, windows, fraction=self.min_area_fraction
+        )
 
         t_to_prev_peak = np.ones(len(peaks), dtype=np.int64) * self.peak_max_proximity_time
         t_to_prev_peak[1:] = peaks["time"][1:] - peaks["endtime"][:-1]
@@ -77,15 +81,17 @@ def compute(self, peaks):
 
     @staticmethod
     @numba.jit(nopython=True, nogil=True, cache=True)
-    def find_n_competing(peaks, windows, fraction):
+    def find_n_competing(peaks_in_roi, peaks, windows, fraction):
         n_left = np.zeros(len(peaks), dtype=np.int32)
         n_tot = n_left.copy()
+        areas_in_roi = peaks_in_roi["area"]
         areas = peaks["area"]
 
+        dig = np.searchsorted(peaks_in_roi["center_time"], peaks["center_time"])
         for i, peak in enumerate(peaks):
             left_i, right_i = windows[i]
             threshold = areas[i] * fraction
-            n_left[i] = np.sum(areas[left_i:i] > threshold)
-            n_tot[i] = n_left[i] + np.sum(areas[i + 1 : right_i] > threshold)
+            n_left[i] = np.sum(areas_in_roi[left_i : dig[i]] > threshold)
+            n_tot[i] = n_left[i] + np.sum(areas_in_roi[dig[i] : right_i] > threshold)
 
         return n_left, n_tot