Veto plugins

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 *
 

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]

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

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)

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