diff --git a/straxen/plugins/__init__.py b/straxen/plugins/__init__.py index fac3e30ab..0d7061089 100644 --- a/straxen/plugins/__init__.py +++ b/straxen/plugins/__init__.py @@ -39,6 +39,9 @@ from . import veto_events from .veto_events import * +from . import veto_veto_regions +from .veto_veto_regions import * + from . import acqmon_processing from .acqmon_processing import * diff --git a/straxen/plugins/nveto_recorder.py b/straxen/plugins/nveto_recorder.py index fb2517ee2..e89815375 100644 --- a/straxen/plugins/nveto_recorder.py +++ b/straxen/plugins/nveto_recorder.py @@ -105,29 +105,37 @@ def compute(self, raw_records_nv, start, end): min_amplitude=self.config['hit_min_amplitude_nv']) del temp_records - # First we have to split rr into records and lone records: # Please note that we consider everything as a lone record which # does not satisfy the coincidence requirement - intervals = coincidence(hits, - self.config['coincidence_level_recorder_nv'], - self.config['resolving_time_recorder_nv'], - self.config['pre_trigger_time_nv']) + intervals = find_coincidence(hits, + self.config['coincidence_level_recorder_nv'], + self.config['resolving_time_recorder_nv'], + self.config['pre_trigger_time_nv']) del hits # Always save the first and last resolving_time nanoseconds (e.g. 600 ns) since we cannot guarantee the gap # size to be larger. (We cannot use an OverlapingWindow plugin either since it requires disjoint objects.) if len(intervals): - intervals_with_bounds = np.zeros((len(intervals) + 2, 2), dtype=np.int64) - intervals_with_bounds[1:-1, :] = intervals - intervals_with_bounds[0, :] = start, min(start + self.config['resolving_time_recorder_nv'], intervals[0, 0]) - intervals_with_bounds[-1, :] = max(end - self.config['resolving_time_recorder_nv'], intervals[-1, 1]), end + intervals_with_bounds = np.zeros(len(intervals) + 2, dtype=strax.time_fields) + intervals_with_bounds['time'][1:-1] = intervals['time'] + intervals_with_bounds['endtime'][1:-1] = intervals['endtime'] + intervals_with_bounds['time'][0] = start + intervals_with_bounds['endtime'][0] = min(start + self.config['resolving_time_recorder_nv'], + intervals['time'][0]) + intervals_with_bounds['time'][-1] = max(end - self.config['resolving_time_recorder_nv'], + intervals['endtime'][-1]) + intervals_with_bounds['endtime'][-1] = end del intervals else: - intervals_with_bounds = np.zeros((0, 2), dtype=np.int64) + intervals_with_bounds = np.zeros((0, 2), dtype=strax.time_fields) neighbors = strax.record_links(raw_records_nv) - mask = pulse_in_interval(raw_records_nv, neighbors, *np.transpose(intervals_with_bounds)) + mask = pulse_in_interval(raw_records_nv, + neighbors, + intervals_with_bounds['time'], + intervals_with_bounds['endtime'],) + rr, lone_records = straxen.mask_and_not(raw_records_nv, mask) # Compute some properties of the lone_records: @@ -326,7 +334,7 @@ def pulse_in_interval(raw_records, record_links, start_times, end_times): @export -def coincidence(records, nfold=4, resolving_time=300, pre_trigger=0): +def find_coincidence(records, nfold=4, resolving_time=300, pre_trigger=0): """ Checks if n-neighboring events are less apart from each other then the specified resolving time. @@ -351,10 +359,12 @@ def coincidence(records, nfold=4, resolving_time=300, pre_trigger=0): # In case of a "single-fold" coincidence every thing gives # the start of a new interval: start_times = records['time'] - intervals = _merge_intervals(start_times-pre_trigger, - resolving_time+pre_trigger) + intervals = np.zeros(len(start_times), dtype=strax.time_fields) + intervals['time'] = start_times - pre_trigger + intervals['endtime'] = start_times + resolving_time + intervals = merge_intervals(intervals) else: - intervals = np.zeros((0, 2), np.int64) + intervals = np.zeros(0, dtype=strax.time_fields) return intervals @@ -377,7 +387,7 @@ def _coincidence(rr, nfold=4, resolving_time=300): t_diff = np.diff(start_times, prepend=start_times[0]) # 2. Now we have to check if n-events are within resolving time: - # -> use moving average with size n to accumulate time between n-pulses + # -> use moving sum with size n to accumulate time between n-pulses # -> check if accumulated time is below resolving time # generate kernel: @@ -393,34 +403,48 @@ def _coincidence(rr, nfold=4, resolving_time=300): return start_times[mask] -@numba.njit(nogil=True, cache=True) -def _merge_intervals(start_time, resolving_time): +@export +def merge_intervals(intervals): """ - Function which merges overlapping time intervals into a single one. + Function which merges overlapping intervals into a single one. - Note: - If start times of two intervals are exactly resolving_time apart - from each other they will be merged into a single interval. + :param intervals: Any numpy array with strax time fields. + :returns: New intervals with time and endtime according to the + overlapping intervals. """ - if not len(start_time): - # If the input is empty return empty array: - return np.zeros((0, 2), dtype=np.int64) - - # check for gaps larger than resolving_time: - # The gaps will indicate the starts of new intervals - gaps = np.diff(start_time) > resolving_time - - last_element = np.argwhere(gaps).flatten() - first_element = 0 - # Creating output - # There is one more interval than gaps - intervals = np.zeros((np.sum(gaps) + 1, 2), dtype=np.int64) - - # Looping over all intervals, except for the last one: - for ind, le in enumerate(last_element): - intervals[ind] = (start_time[first_element], start_time[le]+resolving_time) - first_element = le + 1 - - # Now we have to deal with the last gap: - intervals[-1] = (start_time[first_element], start_time[first_element:][-1] + resolving_time) - return intervals + res = np.zeros(len(intervals), dtype=strax.time_fields) + + if len(intervals): + res = _merge_intervals(intervals['time'], + strax.endtime(intervals), + res) + return res + + +@numba.njit(cache=True, nogil=True) +def _merge_intervals(start, end, res): + offset = 0 + + interval_start = start[0] + interval_end = end[0] + for next_interval_start, next_interval_end in zip(start[1:], end[1:]): + if interval_end >= next_interval_start: + # Interval overlaps, updated only end: + interval_end = next_interval_end + continue + + # Intervals do not overlap, save interval: + res[offset]['time'] = interval_start + res[offset]['endtime'] = interval_end + offset += 1 + + # New interval: + interval_start = next_interval_start + interval_end = next_interval_end + + # Save last interval: + res[offset]['time'] = interval_start + res[offset]['endtime'] = interval_end + offset += 1 + + return res[:offset] diff --git a/straxen/plugins/veto_events.py b/straxen/plugins/veto_events.py index 6de801c36..e3a2ecbed 100644 --- a/straxen/plugins/veto_events.py +++ b/straxen/plugins/veto_events.py @@ -149,19 +149,10 @@ def find_veto_events(hitlets: np.ndarray, :returns: events, hitelt_ids_per_event """ # Find intervals which satisfy requirement: - intervals = straxen.plugins.nveto_recorder.coincidence(hitlets, - coincidence_level, - resolving_time, - left_extension, - - ) - - # Create some preliminary events: - event_intervals = np.zeros(len(intervals), - dtype=strax.time_fields - ) - event_intervals['time'] = intervals[:, 0] - event_intervals['endtime'] = intervals[:, 1] + event_intervals = straxen.plugins.nveto_recorder.find_coincidence(hitlets, + coincidence_level, + resolving_time, + left_extension,) # Find all hitlets which touch the coincidence windows: # (we cannot use fully_contained in here since some muon signals diff --git a/straxen/plugins/veto_veto_regions.py b/straxen/plugins/veto_veto_regions.py new file mode 100644 index 000000000..9e4f336ee --- /dev/null +++ b/straxen/plugins/veto_veto_regions.py @@ -0,0 +1,159 @@ +import strax +import straxen +import numpy as np +import numba + +export, __all__ = strax.exporter() +MV_PREAMBLE = 'Muno-Veto Plugin: Same as the corresponding nVETO-Plugin.\n' + + +@export +@strax.takes_config( + strax.Option( + 'min_veto_area_nv', default=5, type=float, track=True, + help='Minimal area required in pe to trigger veto.'), + strax.Option( + 'min_veto_hits_nv', default=10, type=int, track=True, + help='Minimal number of hitlets in n/mveto_event to trigger a veto.'), + strax.Option( + 'min_veto_channel_nv', default=0, type=int, track=True, + help='Minimal number PMT channel contributing to the n/mveto_event.'), + strax.Option( + 'veto_left_extension_nv', default=500_000, type=int, track=True, + help='Veto time in ns left t the start of a vetoing event.'), + strax.Option( + 'veto_right_extension_nv', default=0, type=int, track=True, + help='Veto time in ns right to the end of a vetoing event.'), +) +class nVETOVetoRegions(strax.OverlapWindowPlugin): + """ + Plugin which defines the time intervals in which peaks should be + tagged as vetoed. An event must surpass all three criteria to trigger + a veto. + """ + __version__ = '0.0.1' + + depends_on = 'events_nv' + provides = 'veto_regions_nv' + data_kind = 'veto_regions_nv' + save_when = strax.SaveWhen.NEVER + + dtype = strax.time_fields + ends_with = '_nv' + + def get_window_size(self): + return 10 * (self.config['veto_left_extension_nv'] + self.config['veto_right_extension_nv']) + + def compute(self, events_nv, start, end): + vetos = create_veto_intervals(events_nv, + self.config['min_veto_area_nv'], + self.config['min_veto_hits_nv'], + self.config['min_veto_channel_nv'], + self.config['veto_left_extension_nv'], + self.config['veto_right_extension_nv']) + + # Now we have to do clip all times and end endtimes in case + # they go beyond a chunk boundary: + vetos['time'] = np.clip(vetos['time'], start, end) + vetos['endtime'] = np.clip(vetos['endtime'], start, end) + + return vetos + + +def create_veto_intervals(events, + min_area, + min_hits, + min_contributing_channels, + left_extension, + right_extension, ): + """ + Function which creates veto regions. + + :param events: nveto events + :param min_area: min area in pe required to create a veto region. + :param min_hits: same but with hitlets + :param min_contributing_channels: Minimal number of contributing + channels. + :param left_extension: Left extension of the event to define veto + region in ns. + :param right_extension: Same but right hand side after. + + :returns: numpy.structured.array containing the veto regions. + """ + res = np.zeros(len(events), + dtype=strax.time_fields) + + res = _create_veto_intervals(events, + min_area, + min_hits, + min_contributing_channels, + left_extension, + right_extension, + res, ) + res = straxen.merge_intervals(res) + + return res + + +@numba.njit(cache=True) +def _create_veto_intervals(events, + min_area, + min_hits, + min_contributing_channels, + left_extension, + right_extension, + res, ): + offset = 0 + + for ev in events: + satisfies_veto_trigger = (ev['area'] >= min_area + and ev['n_hits'] >= min_hits + and ev['n_contributing_pmt'] >= min_contributing_channels) + if not satisfies_veto_trigger: + continue + + res[offset]['time'] = ev['time'] - left_extension + res[offset]['endtime'] = ev['endtime'] + right_extension + offset += 1 + return res[:offset] + + +@strax.takes_config( + strax.Option( + 'min_veto_area_mv', default=10, type=float, track=True, + child_option=True, parent_option_name='min_veto_area_nv', + help='Minimal area required in pe to trigger veto.'), + strax.Option( + 'min_veto_hits_mv', default=0, type=int, track=True, + child_option=True, parent_option_name='min_veto_hits_nv', + help='Minimal number of hitlets in event to trigger veto.'), + strax.Option( + 'min_veto_channel_mv', default=5, type=int, track=True, + child_option=True, parent_option_name='min_veto_channel_nv', + help='Minimal number PMT channel contributing to the event.'), + strax.Option( + 'veto_left_extension_mv', default=0, type=int, track=True, + child_option=True, parent_option_name='veto_left_extension_nv', + help='Veto time in ns left t the start of a vetoing event.'), + strax.Option( + 'veto_right_extension_mv', default=1_000_000, type=int, track=True, + child_option=True, parent_option_name='veto_right_extension_nv', + help='Veto time in ns right to the end of a vetoing event.'), +) +class muVETOVetoRegions(nVETOVetoRegions): + __doc__ = MV_PREAMBLE + nVETOVetoRegions.__doc__ + __version__ = '0.0.1' + + depends_on = 'events_mv' + provides = 'veto_regions_mv' + data_kind = 'veto_regions_mv' + save_when = strax.SaveWhen.NEVER + + dtype = strax.time_fields + ends_with = '_mv' + + def get_window_size(self): + return 10 * (self.config['veto_left_mv'] + self.config['veto_right_mv']) + + def compute(self, events_mv, start, end): + return super().compue(events_mv, start, end) diff --git a/tests/test_nveto_events.py b/tests/test_nveto_events.py index 730d2a6dc..2ab59f306 100644 --- a/tests/test_nveto_events.py +++ b/tests/test_nveto_events.py @@ -114,16 +114,9 @@ def test_nveto_event_building(hitlets, """ hitlets = strax.sort_by_time(hitlets) - intervals = straxen.plugins.nveto_recorder.coincidence(hitlets, - coincidence, - 300) - - # Create event_intervals: - event_intervals = np.zeros(len(intervals), - dtype=strax.time_fields - ) - event_intervals['time'] = intervals[:, 0] - event_intervals['endtime'] = intervals[:, 1] + event_intervals = straxen.plugins.nveto_recorder.find_coincidence(hitlets, + coincidence, + 300) mes = 'Found overlapping events returned by "coincidence".' assert np.all(event_intervals['endtime'][:-1] - event_intervals['time'][1:] < 0), mes diff --git a/tests/test_nveto_recorder.py b/tests/test_nveto_recorder.py new file mode 100644 index 000000000..162aeedc4 --- /dev/null +++ b/tests/test_nveto_recorder.py @@ -0,0 +1,136 @@ +import strax +import straxen + +import numpy as np +import unittest + + +class TestMergeIntervals(unittest.TestCase): + + def setUp(self): + self.intervals = np.zeros(4, dtype=strax.time_fields) + self.intervals['time'] = [2, 3, 7, 20] + self.intervals['endtime'] = [4, 7, 8, 22] + + def test_empty_intervals(self): + intervals = np.zeros(0, dtype=strax.time_fields) + intervals = straxen.merge_intervals(intervals) + assert len(intervals) == 0, 'Empty input should return empty intervals!' + + def test_merge_overlapping_intervals(self): + intervals = straxen.merge_intervals(self.intervals) + + assert len(intervals) == 2, 'Got the wrong number of intervals!' + + time_is_correct = intervals[0]['time'] == self.intervals['time'][0] + assert time_is_correct, 'First interval has the wrong time!' + time_is_correct = intervals[0]['endtime'] == self.intervals['endtime'][-2] + assert time_is_correct, 'First interval has the wrong endtime!' + + time_is_correct = intervals[-1]['time'] == self.intervals['time'][-1] + assert time_is_correct, 'Second interval has the wrong time!' + time_is_correct = intervals[-1]['endtime'] == self.intervals['endtime'][-1] + assert time_is_correct, 'Second interval has the wrong endtime!' + + +class TestCoincidence(unittest.TestCase): + + def setUp(self): + self.intervals = np.zeros(8, dtype=strax.time_fields) + self.intervals['time'] = [3, 6, 9, 12, 15, 18, 21, 38] + self.intervals['endtime'] = [5, 8, 10, 13, 16, 19, 23, 42] + + def test_empty_inputs(self): + intervals = np.zeros(0, dtype=strax.time_fields) + intervals = straxen.find_coincidence(intervals) + assert len(intervals) == 0, 'Empty input should return empty intervals!' + + def test_without_coincidence(self): + resolving_time = 10 + truth_time = np.array([self.intervals['time'][0], + self.intervals['time'][-1]]) + truth_endtime = np.array([self.intervals['time'][-2], + self.intervals['time'][-1]]) + resolving_time + self._test_coincidence(resolving_time=resolving_time, + coincidence=1, + pre_trigger=0, + n_concidences_truth=2, + times_truth=truth_time, + endtime_truth=truth_endtime + ) + pre_trigger = 2 + self._test_coincidence(resolving_time=resolving_time, + coincidence=1, + pre_trigger=pre_trigger, + n_concidences_truth=2, + times_truth=truth_time - pre_trigger, + endtime_truth=truth_endtime + ) + + def test_even_fold(self): + resolving_time = 10 + self._test_coincidence(resolving_time=resolving_time, + coincidence=2, + pre_trigger=0, + n_concidences_truth=1, + times_truth=self.intervals['time'][0], + endtime_truth=self.intervals['time'][-3] + resolving_time, + ) + pre_trigger = 2 + self._test_coincidence(resolving_time=resolving_time, + coincidence=2, + pre_trigger=pre_trigger, + n_concidences_truth=1, + times_truth=self.intervals['time'][0] - pre_trigger, + endtime_truth=self.intervals['time'][-3] + resolving_time, + ) + + self._test_coincidence(resolving_time=resolving_time, + coincidence=4, + pre_trigger=pre_trigger, + n_concidences_truth=1, + times_truth=self.intervals['time'][0] - pre_trigger, + endtime_truth=self.intervals['time'][-5] + resolving_time, + ) + + def test_odd_fold(self): + resolving_time = 10 + self._test_coincidence(resolving_time=resolving_time, + coincidence=3, + pre_trigger=0, + n_concidences_truth=1, + times_truth=self.intervals['time'][0], + endtime_truth=self.intervals['time'][-4] + resolving_time, + ) + pre_trigger = 2 + self._test_coincidence(resolving_time=resolving_time, + coincidence=3, + pre_trigger=pre_trigger, + n_concidences_truth=1, + times_truth=self.intervals['time'][0] - pre_trigger, + endtime_truth=self.intervals['time'][-4] + resolving_time, + ) + + self._test_coincidence(resolving_time=resolving_time, + coincidence=5, + pre_trigger=pre_trigger, + n_concidences_truth=0, + times_truth=self.intervals['time'][:0], + endtime_truth=self.intervals['time'][:0], + ) + + def _test_coincidence(self, resolving_time, coincidence, pre_trigger, + n_concidences_truth, times_truth, endtime_truth): + coincidence = straxen.find_coincidence(self.intervals, + nfold=coincidence, + resolving_time=resolving_time, + pre_trigger=pre_trigger) + number_coincidence_correct = len(coincidence) == n_concidences_truth + assert number_coincidence_correct, 'Have not found the correct number of coincidences' + + time_is_correct = np.all(coincidence['time'] == times_truth) + assert time_is_correct, 'Coincidence does not have the correct time' + + endtime_is_correct = np.all(coincidence['endtime'] == endtime_truth) + print(coincidence['endtime'], endtime_truth) + assert endtime_is_correct, 'Coincidence does not have the correct endtime' diff --git a/tests/test_veto_veto_regions.py b/tests/test_veto_veto_regions.py new file mode 100644 index 000000000..51f30743c --- /dev/null +++ b/tests/test_veto_veto_regions.py @@ -0,0 +1,72 @@ +import straxen +import numpy as np +import unittest + + +class TestCreateVetoIntervals(unittest.TestCase): + def setUp(self): + dtype = straxen.plugins.veto_events.veto_event_dtype('nveto_eventbumber') + dtype += straxen.plugins.veto_events.veto_event_positions_dtype()[2:] + self.dtype = dtype + + self.events = np.zeros(4, self.dtype) + self.events['area'] = 1 + self.events['n_hits'] = 1 + self.events['n_contributing_pmt'] = 1 + self.events['time'] = [2, 5, 7, 20] + self.events['endtime'] = [3, 7, 8, 22] + + def test_empty_inputs(self): + events = np.zeros(0, self.dtype) + vetos = straxen.plugins.veto_veto_regions.create_veto_intervals(events, 0, 0, 0, 10, 10) + assert len(vetos) == 0, 'Empty input must return empty output!' + + def test_concatenate_overlapping_intervals(self): + left_extension = 1 + right_extension = 4 + vetos = straxen.plugins.veto_veto_regions.create_veto_intervals(self.events, + min_area=0, + min_hits=0, + min_contributing_channels=0, + left_extension=left_extension, + right_extension=right_extension) + assert len(vetos) == 2, 'Got the wrong number of veto intervals!' + + time_is_correct = vetos[0]['time'] == self.events['time'][0]-left_extension + assert time_is_correct, 'First veto event has the wrong time!' + time_is_correct = vetos[0]['endtime'] == self.events['endtime'][2] + right_extension + assert time_is_correct, 'First veto event has the wrong endtime!' + + time_is_correct = vetos[1]['time'] == self.events['time'][-1] - left_extension + assert time_is_correct, 'Second veto event has the wrong time!' + time_is_correct = vetos[1]['endtime'] == self.events['endtime'][-1] + right_extension + assert time_is_correct, 'Second veto event has the wrong endtime!' + + def test_thresholds(self): + events = np.copy(self.events[:1]) + + self._test_threshold_type(events, 'area', 'min_area', 2) + self._test_threshold_type(events, 'n_hits', 'min_hits', 2) + self._test_threshold_type(events, 'n_contributing_pmt', 'min_contributing_channels', 2) + + @staticmethod + def _test_threshold_type(events, field, threshold_type, threshold): + thresholds = {'min_area': 1, + 'min_hits': 1, + 'min_contributing_channels': 1} + thresholds = {key: (threshold if threshold_type == key else 1) for key in thresholds.keys()} + + vetos = straxen.plugins.veto_veto_regions.create_veto_intervals(events, + **thresholds, + left_extension=0, + right_extension=0) + print(events[field], thresholds, vetos) + assert len(vetos) == 0, f'Vetos for {threshold_type} threshold should be empty since it is below threshold!' + + events[field] = threshold + vetos = straxen.plugins.veto_veto_regions.create_veto_intervals(events, + **thresholds, + left_extension=0, + right_extension=0) + assert len(vetos) == 1, f'{threshold_type} threshold did not work, have a wrong number of vetos!' + events[field] = 1