diff --git a/straxen/plugins/events/events.py b/straxen/plugins/events/events.py index 3696ea31c..fa0da5a11 100644 --- a/straxen/plugins/events/events.py +++ b/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: diff --git a/straxen/plugins/merged_s2s/merged_s2s.py b/straxen/plugins/merged_s2s/merged_s2s.py index 8cdd1a43e..560c1c78d 100644 --- a/straxen/plugins/merged_s2s/merged_s2s.py +++ b/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): diff --git a/straxen/plugins/peaks/peak_proximity.py b/straxen/plugins/peaks/peak_proximity.py index ddcefa440..9181eebb0 100644 --- a/straxen/plugins/peaks/peak_proximity.py +++ b/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