From 08cd212173637de5ea80e326419b40f4ae932a83 Mon Sep 17 00:00:00 2001 From: dachengx Date: Fri, 6 Dec 2024 08:35:24 -0600 Subject: [PATCH 1/6] S2 merging exclude sparse peaklets --- straxen/plugins/events/events.py | 1 + straxen/plugins/merged_s2s/merged_s2s.py | 44 ++++++++++++++++++++++-- straxen/plugins/peaks/peak_proximity.py | 1 + 3 files changed, 43 insertions(+), 3 deletions(-) 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..98850167c 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,12 @@ class MergedS2s(strax.OverlapWindowPlugin): ), ) + 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 +143,13 @@ def compute(self, peaklets, lone_hits): end_merge_at, max_buffer=int(self.s2_merge_max_duration // np.gcd.reduce(peaklets["dt"])), ) + + 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) merged_s2s["type"] = 2 + merged_s2s["type"][dr_cnf_avg > self.dr_cnf_avg_threshold] = 20 # Updated time and length of lone_hits and sort again: lh = np.copy(lone_hits) @@ -226,6 +244,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(peaklets), 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..4e97823df 100644 --- a/straxen/plugins/peaks/peak_proximity.py +++ b/straxen/plugins/peaks/peak_proximity.py @@ -57,6 +57,7 @@ def get_window_size(self): def compute(self, peaks): windows = strax.touching_windows(peaks, peaks, window=self.nearby_window) + # have to consider the peak with type 20 n_left, n_tot = self.find_n_competing(peaks, windows, fraction=self.min_area_fraction) t_to_prev_peak = np.ones(len(peaks), dtype=np.int64) * self.peak_max_proximity_time From 63c89e852a4cbdc1807ebbc0bf33f652f7dac69d Mon Sep 17 00:00:00 2001 From: dachengx Date: Fri, 6 Dec 2024 08:56:38 -0600 Subject: [PATCH 2/6] Set 1E4 PE threshold of whether apply the algorithm --- straxen/plugins/merged_s2s/merged_s2s.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/straxen/plugins/merged_s2s/merged_s2s.py b/straxen/plugins/merged_s2s/merged_s2s.py index 98850167c..9f3dd4c41 100644 --- a/straxen/plugins/merged_s2s/merged_s2s.py +++ b/straxen/plugins/merged_s2s/merged_s2s.py @@ -58,6 +58,12 @@ 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, @@ -144,13 +150,6 @@ def compute(self, peaklets, lone_hits): max_buffer=int(self.s2_merge_max_duration // np.gcd.reduce(peaklets["dt"])), ) - 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) - merged_s2s["type"] = 2 - merged_s2s["type"][dr_cnf_avg > self.dr_cnf_avg_threshold] = 20 - # Updated time and length of lone_hits and sort again: lh = np.copy(lone_hits) del lone_hits @@ -174,6 +173,15 @@ def compute(self, peaklets, lone_hits): if (_n_top_pmts <= 0) or (not _store_data_start): merged_s2s = drop_data_field(merged_s2s, self.dtype, "_drop_data_field_merged_s2s") + 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) + merged_s2s["type"] = 2 + mask = merged_s2s["area"] < self.position_density_s2_area_threshold + mask &= dr_cnf_avg > self.dr_cnf_avg_threshold + merged_s2s["type"][mask] = 20 + return merged_s2s @staticmethod From f7067437db6905e8002f2dfa297bf8b5a0c771c1 Mon Sep 17 00:00:00 2001 From: dachengx Date: Fri, 6 Dec 2024 09:19:30 -0600 Subject: [PATCH 3/6] Fix bug --- straxen/plugins/merged_s2s/merged_s2s.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/straxen/plugins/merged_s2s/merged_s2s.py b/straxen/plugins/merged_s2s/merged_s2s.py index 9f3dd4c41..8b948a8e5 100644 --- a/straxen/plugins/merged_s2s/merged_s2s.py +++ b/straxen/plugins/merged_s2s/merged_s2s.py @@ -257,7 +257,7 @@ def weighted_average_dr(self, area_top, peaklets, merged_s2s): mask = area_top > 0 mask &= ~np.isnan(peaklets["x_cnf"]) mask &= ~np.isnan(peaklets["y_cnf"]) - dr_cnf_avg = np.full(len(peaklets), np.nan) + 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) From 1c752b7dabe93e523a955db3b695a9207d62f0f6 Mon Sep 17 00:00:00 2001 From: dachengx Date: Fri, 6 Dec 2024 10:39:10 -0600 Subject: [PATCH 4/6] Fix bug --- straxen/plugins/merged_s2s/merged_s2s.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/straxen/plugins/merged_s2s/merged_s2s.py b/straxen/plugins/merged_s2s/merged_s2s.py index 8b948a8e5..1e7dc57fb 100644 --- a/straxen/plugins/merged_s2s/merged_s2s.py +++ b/straxen/plugins/merged_s2s/merged_s2s.py @@ -158,21 +158,18 @@ 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") - 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 @@ -182,6 +179,9 @@ def compute(self, peaklets, lone_hits): 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 @staticmethod From 715d888fb82335c303615612fb695b4ac7ba3ee8 Mon Sep 17 00:00:00 2001 From: dachengx Date: Sat, 7 Dec 2024 08:53:59 -0600 Subject: [PATCH 5/6] Minor change --- straxen/plugins/merged_s2s/merged_s2s.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/straxen/plugins/merged_s2s/merged_s2s.py b/straxen/plugins/merged_s2s/merged_s2s.py index 1e7dc57fb..560c1c78d 100644 --- a/straxen/plugins/merged_s2s/merged_s2s.py +++ b/straxen/plugins/merged_s2s/merged_s2s.py @@ -170,11 +170,11 @@ def compute(self, peaklets, lone_hits): strax.compute_widths(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) - merged_s2s["type"] = 2 mask = merged_s2s["area"] < self.position_density_s2_area_threshold mask &= dr_cnf_avg > self.dr_cnf_avg_threshold merged_s2s["type"][mask] = 20 From c954cfff12d66fc1a51ca7ea24cc4a2c1fd7c049 Mon Sep 17 00:00:00 2001 From: dachengx Date: Tue, 10 Dec 2024 09:25:11 -0600 Subject: [PATCH 6/6] Calculate `n_competing` only of type 0, 1, 2 --- straxen/plugins/peaks/peak_proximity.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/straxen/plugins/peaks/peak_proximity.py b/straxen/plugins/peaks/peak_proximity.py index 4e97823df..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,9 +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) - # have to consider the peak with type 20 - 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] @@ -78,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