Merge pull request #261 from AxFoundation/lone_hit_integration

Routines for lone hit integration

strax/dtypes.py

@@ -88,6 +88,10 @@ def record_dtype(samples_per_record=DEFAULT_RECORD_LENGTH):
         'left'), np.int16),
     (('Index of first sample in record just beyond hit (exclusive bound)',
         'right'), np.int16),
+    (('For lone hits, index of sample in record where integration starts',
+      'left_integration'), np.int16),
+    (('For lone hits, index of first sample beyond integration region',
+      'right_integration'), np.int16),
     (('Internal (temporary) index of fragment in which hit was found',
         'record_i'), np.int32),
     (('ADC threshold applied in order to find hits',

strax/processing/peak_building.py

@@ -253,3 +253,98 @@ def find_peak_groups(peaks, gap_threshold,
         left_extension=left_extension, right_extension=right_extension,
         min_channels=1, min_area=0)
     return fake_peaks['time'], strax.endtime(fake_peaks)
+
+
+##
+# Lone hit integration
+##
+
+@numba.njit(nogil=True, cache=True)
+def _find_hit_integration_bounds(
+        lone_hits, peaks, records, save_outside_hits, n_channels):
+    """"Update lone hits to include integration bounds
+
+    save_outside_hits: in ns!!
+    """
+    result = np.zeros((len(lone_hits), 2), dtype=np.int64)
+    if not len(lone_hits):
+        return result
+
+    # By default, use save_outside_hits to determine bounds
+    result[:, 0] = lone_hits['time'] - save_outside_hits[0]
+    result[:, 1] = strax.endtime(lone_hits) + save_outside_hits[1]
+
+    NO_EARLIER_HIT = -1
+    last_hit_index = np.ones(n_channels, dtype=np.int32) * NO_EARLIER_HIT
+
+    n_peaks = len(peaks)
+    FAR_AWAY = 9223372036_854775807   # np.iinfo(np.int64).max, April 2262
+    peak_i = 0
+
+    for hit_i, h in enumerate(lone_hits):
+        ch = h['channel']
+
+        # Find end of previous peak and start of next peak
+        # (note peaks are disjoint from any lone hit, even though
+        # lone hits may not be disjoint from each other)
+        while peak_i < n_peaks and peaks[peak_i]['time'] < h['time']:
+            peak_i += 1
+        prev_p_end = strax.endtime(peaks[peak_i - 1]) if peak_i != 0 else 0
+        next_p_start = peaks[peak_i]['time'] if peak_i != n_peaks else FAR_AWAY
+
+
+        # Ensure we do not integrate parts of peaks
+        # or (at least for now) beyond the record in which the hit was found
+        r = records[h['record_i']]
+        result[hit_i][0] = max(prev_p_end,
+                               r['time'],
+                               result[hit_i][0])
+        result[hit_i][1] = min(next_p_start,
+                               strax.endtime(r),
+                               result[hit_i][1])
+
+        if last_hit_index[ch] != NO_EARLIER_HIT:
+            # Ensure previous hit does not integrate the over-threshold region
+            # of this hit
+            result[last_hit_index[ch]][1] = min(result[last_hit_index[ch]][1],
+                                                h['time'])
+            # Ensure this hit doesn't integrate anything the previous hit
+            # already integrated
+            result[hit_i][0] = max(result[last_hit_index[ch]][1],
+                                   result[hit_i][0])
+
+        last_hit_index[ch] = hit_i
+
+    # Convert to index in record and store
+    t0 = records[lone_hits['record_i']]['time']
+    dt = records[lone_hits['record_i']]['dt']
+    for hit_i, h in enumerate(lone_hits):
+        h['left_integration'] = (result[hit_i, 0] - t0[hit_i]) // dt[hit_i]
+        h['right_integration'] = (result[hit_i, 1] - t0[hit_i]) // dt[hit_i]
+
+
+@export
+@numba.njit(nogil=True, cache=True)
+def integrate_lone_hits(
+        lone_hits, records, peaks, save_outside_hits, n_channels):
+    """Update the area of lone_hits to the integral in ADCcounts x samples
+
+    :param lone_hits: Hits outside of peaks
+    :param records: Records in which hits and peaks were found
+    :param peaks: Peaks
+    :param save_outside_hits: (left, right) *TIME* with wich we should extend
+    the integration window of hits
+    the integration region
+    :param n_channels: number of channels
+
+    TODO: this doesn't extend the integration range beyond record boundaries
+    """
+    _find_hit_integration_bounds(
+        lone_hits, peaks, records, save_outside_hits, n_channels)
+    for hit_i, h in enumerate(lone_hits):
+        r = records[h['record_i']]
+        start, end = h['left_integration'], h['right_integration']
+        # TODO: when we add amplitude multiplier, adjust this too!
+        h['area'] = (
+                r['data'][start:end].sum()
+                + (r['baseline'] % 1) * (end - start))

strax/utils.py

@@ -221,7 +221,7 @@ def profile_threaded(filename):
 @export
 def to_str_tuple(x) -> ty.Tuple[str]:
     if isinstance(x, str):
-        return x,
+        return (x,)
     elif isinstance(x, list):
         return tuple(x)
     elif isinstance(x, tuple):