10168286066

Committed 30 Jul 2024 07:07PM UTC coverage: 90.957% (-0.2%) from 91.135%

Build # 10168286066

Build Type

Pull #1401

github

Committed by

web-flow

Commit Message

Merge 49c9fb3c8 into 4d6c2b7ba

Pull Request Pull Request #1401: Dynamic led window

Run Details

48 of 71 new or added lines in 2 files covered. (67.61%)

1 existing line in 1 file now uncovered.

8992 of 9886 relevant lines covered (90.96%)

1.82 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

99.24

/straxen/plugins/afterpulses/afterpulse_processing.py

import numba
import numpy as np
import strax
import straxen

export, __all__ = strax.exporter()


@export
class LEDAfterpulseProcessing(strax.Plugin):
    """Plugin for processing LED afterpulses.

    Detect LED pulses and afterpulses (APs) in raw_records waveforms. Compute the AP datatype.

    """

    __version__ = "0.6.0"
    depends_on = "raw_records"
    data_kind = "afterpulses"
    provides = "afterpulses"
    compressor = "zstd"
    parallel = "process"
    rechunk_on_save = True

    gain_model = straxen.URLConfig(
        infer_type=False,
        help="PMT gain model. Specify as (model_type, model_config)",
    )

    n_tpc_pmts = straxen.URLConfig(
        type=int,
        help="Number of PMTs in TPC",
    )

    LED_hit_left_boundary = straxen.URLConfig(
        default=50,
        infer_type=False,
        help="Left boundary after which the first hit marks the position of the LED window",
    )

    LED_hit_right_boundary = straxen.URLConfig(
        default=110,
        infer_type=False,
        help="Right boundary after which the LED hit will no longer be searched",
    )

    LED_window_width = straxen.URLConfig(
        default=20,
        infer_type=False,
        help="Width of the window in which hits are merged into the LED hit after the first hit",
    )

    baseline_samples = straxen.URLConfig(
        default=40,
        infer_type=False,
        help="Number of samples to use at start of WF to determine the baseline",
    )

    hit_min_amplitude = straxen.URLConfig(
        track=True,
        infer_type=False,
        default="cmt://hit_thresholds_tpc?version=ONLINE&run_id=plugin.run_id",
        help=(
            "Minimum hit amplitude in ADC counts above baseline. "
            "Specify as a tuple of length n_tpc_pmts, or a number,"
            'or a string like "legacy-thresholds://pmt_commissioning_initial" which means calling'
            "hitfinder_thresholds.py"
            'or url string like "cmt://hit_thresholds_tpc?version=ONLINE" which means'
            "calling cmt."
        ),
    )

    hit_min_height_over_noise = straxen.URLConfig(
        default=4,
        infer_type=False,
        help=(
            "Minimum hit amplitude in numbers of baseline_rms above baseline."
            "Actual threshold used is max(hit_min_amplitude, hit_min_"
            "height_over_noise * baseline_rms)."
        ),
    )

    save_outside_hits = straxen.URLConfig(
        default=(3, 20),
        infer_type=False,
        help="Save (left, right) samples besides hits; cut the rest",
    )

    def infer_dtype(self):
        dtype = dtype_afterpulses()
        return dtype

    def setup(self):
        self.to_pe = self.gain_model
        self.hit_thresholds = self.hit_min_amplitude
        self.hit_left_extension, self.hit_right_extension = self.save_outside_hits

    def compute(self, raw_records):
        # Convert everything to the records data type -- adds extra fields.
        records = strax.raw_to_records(raw_records)
        del raw_records

        # calculate baseline and baseline rms
        strax.baseline(records, baseline_samples=self.baseline_samples, flip=True)

        # find all hits
        hits = strax.find_hits(
            records,
            min_amplitude=self.hit_thresholds,
            min_height_over_noise=self.hit_min_height_over_noise,
        )

        # sort hits by record_i and time, then find LED hit and afterpulse
        # hits within the same record
        hits_ap = find_ap(
            hits,
            records,
            LED_hit_left_boundary=self.LED_hit_left_boundary,
            LED_hit_right_boundary=self.LED_hit_right_boundary,
            LED_window_width=self.LED_window_width,
            hit_left_extension=self.hit_left_extension,
            hit_right_extension=self.hit_right_extension,
        )

        hits_ap["area_pe"] = hits_ap["area"] * self.to_pe[hits_ap["channel"]]
        hits_ap["height_pe"] = hits_ap["height"] * self.to_pe[hits_ap["channel"]]

        return hits_ap


@export
def find_ap(
    hits,
    records,
    LED_hit_left_boundary,
    LED_hit_right_boundary,
    LED_window_width,
    hit_left_extension,
    hit_right_extension,
):
    """Find afterpulses (APs) in the given hits data within specified LED hit boundaries and
    extensions.

    Parameters
    ----------
    hits :
        Array containing hit data.
    records :
        Array containing record data.
    LED_hit_left_boundary :
        Left boundary of the LED hit window.
    LED_hit_right_boundary :
        Right boundary of the LED hit window.
    LED_window_width :
        Extension to the right of the first hit found in the LED hit window
        within which hits are merged into the LED hit.
    hit_left_extension :
        Extension to the left of the hit window.
    hit_right_extension :
        Extension to the right of the hit window.

    Returns
    -------
    Array containing afterpulse data.

    Notes
    -----
    - Hits to the left of the LED_hit_left_boundary are ignored.
    - If no hit is found between LED_hit_left_boundary and LED_hit_right_boundary
    the record is skipped.
    - The merged LED hits are also saved and can be selected for by having
    t_delay = 0 by definition.

    """

    buffer = np.zeros(len(hits), dtype=dtype_afterpulses())

    if not len(hits):
        return buffer

    # sort hits first by record_i, then by time
    hits_sorted = np.sort(hits, order=("record_i", "time"))
    res = _find_ap(
        hits_sorted,
        records,
        LED_hit_left_boundary,
        LED_hit_right_boundary,
        LED_window_width,
        hit_left_extension,
        hit_right_extension,
        buffer=buffer,
    )
    return res


@numba.jit(nopython=True, nogil=True, cache=True)
def _find_ap(
    hits,
    records,
    LED_hit_left_boundary,
    LED_hit_right_boundary,
    LED_window_width,
    hit_left_extension,
    hit_right_extension,
    buffer=None,
):
    # hits need to be sorted by record_i, then time!
    offset = 0

    is_LED = False
    t_LED = None
    t_LED_hit = None

    prev_record_i = hits[0]["record_i"]
    record_data = records[prev_record_i]["data"]
    record_len = records[prev_record_i]["length"]
    baseline_fpart = records[prev_record_i]["baseline"] % 1

    for h_i, h in enumerate(hits):
        if h["record_i"] > prev_record_i:
            # start of a new record
            is_LED = False
            # only increment buffer if the old one is not empty! this happens
            # when no (LED) hit is found in the previous record
            if not buffer[offset]["time"] == 0:
                offset += 1
            prev_record_i = h["record_i"]
            record_data = records[prev_record_i]["data"]
            baseline_fpart = records[prev_record_i]["baseline"] % 1

        res = buffer[offset]

        if h["left"] < LED_hit_left_boundary:
            # if hit is before LED hit window: discard
            continue

        if (h["left"] < LED_hit_right_boundary) and not is_LED:
            # this is the first hit in the LED hit window
            fill_hitpars(
                res,
                h,
                hit_left_extension,
                hit_right_extension,
                record_data,
                record_len,
                baseline_fpart,
            )

            t_LED_hit = res["sample_10pc_area"]
            t_LED = res["sample_10pc_area"]
            is_LED = True

            continue

        if not is_LED:
            # No hit found in the LED hit window: skip to next record
            continue

        if h["left"] < (t_LED_hit + LED_window_width):
            # This hit is still inside the LED window: extend the LED hit
            fill_hitpars(
                res,
                h,
                hit_left_extension,
                hit_right_extension,
                record_data,
                record_len,
                baseline_fpart,
                extend=True,
            )

            # set the LED time in the current WF
            t_LED = res["sample_10pc_area"]

            continue

        # Here begins a new hit after the LED window

        # if a hit is completely inside the previous hit's right_extension,
        # then skip it (because it's already included in the previous hit)
        if h["right"] <= res["right_integration"]:
            continue

        # if a hit only partly overlaps with the previous hit's right_
        # extension, merge them (extend previous hit by this one)
        if h["left"] <= res["right_integration"]:
            fill_hitpars(
                res,
                h,
                hit_left_extension,
                hit_right_extension,
                record_data,
                record_len,
                baseline_fpart,
                extend=True,
            )

            res["tdelay"] = res["sample_10pc_area"] - t_LED

            continue

        # an actual new hit increases the buffer index
        offset += 1
        res = buffer[offset]

        fill_hitpars(
            res,
            h,
            hit_left_extension,
            hit_right_extension,
            record_data,
            record_len,
            baseline_fpart,
        )

        res["tdelay"] = res["sample_10pc_area"] - t_LED

    return buffer[:offset]


@export
@numba.jit(nopython=True, nogil=True, cache=True)
def get_sample_area_quantile(data, quantile, baseline_fpart):
    """Return the index of the first sample in the hit where the integrated area of the hit to that
    index is above the specified quantile of the total area.

    Parameters:
    - data :
        Array containing the baselined waveform data of the hit.
    - quantile :
        The quantile (0 to 1) representing the threshold for the area.
    - baseline_fpart :
        Fractional part of the baseline (baseline % 1) of the record

    Return:
    - int: The index of the first sample where the area exceeds the quantile of the total area.
           If no such sample is found, returns 0.

    Notes:
    - If no quantile is found where the area exceeds the threshold, it returns 0. This is
    usually caused by real events in the baseline window, which can result in a negative
    area.

    """

    area = 0
    area_tot = data.sum() + len(data) * baseline_fpart

    for d_i, d in enumerate(data):
        area += d + baseline_fpart
        if area > (quantile * area_tot):
            return d_i
        if d_i == len(data) - 1:
            # if no quantile was found, something went wrong
            # (negative area due to wrong baseline, caused by real events that
            # by coincidence fall in the first samples of the trigger window)
            # print('no quantile found: set to 0')
            return 0

    # What happened here?!
    return 0


@numba.jit(nopython=True, nogil=True, cache=True)
def fill_hitpars(
    result,
    hit,
    hit_left_extension,
    hit_right_extension,
    record_data,
    record_len,
    baseline_fpart,
    extend=False,
):
    if not extend:  # fill first time only
        result["time"] = hit["time"] - hit_left_extension * hit["dt"]
        result["dt"] = hit["dt"]
        result["channel"] = hit["channel"]
        result["left"] = hit["left"]
        result["record_i"] = hit["record_i"]
        result["threshold"] = hit["threshold"]
        result["left_integration"] = hit["left"] - hit_left_extension
        result["height"] = hit["height"]

    # fill always (if hits are merged, only these will be updated)
    result["right"] = hit["right"]
    result["right_integration"] = hit["right"] + hit_right_extension
    if result["right_integration"] > record_len:
        result["right_integration"] = record_len  # cap right_integration at end of record
    result["length"] = result["right_integration"] - result["left_integration"]

    hit_data = record_data[result["left_integration"] : result["right_integration"]]
    result["area"] = hit_data.sum() + result["length"] * baseline_fpart
    result["sample_10pc_area"] = result["left_integration"] + get_sample_area_quantile(
        hit_data, 0.1, baseline_fpart
    )
    result["sample_50pc_area"] = result["left_integration"] + get_sample_area_quantile(
        hit_data, 0.5, baseline_fpart
    )
    if len(hit_data):
        result["max"] = result["left_integration"] + hit_data.argmax()

    if extend:  # only when merging hits
        result["height"] = max(result["height"], hit["height"])


@export
def dtype_afterpulses():
    """The afterpulse datatype.

    Return:
    - The afterpulse datatype

    """
    dtype_ap = [
        (("Channel/PMT number", "channel"), "<i2"),
        (("Time resolution in ns", "dt"), "<i2"),
        (("Start time of the interval (ns since unix epoch)", "time"), "<i8"),
        (("Length of the interval in samples", "length"), "<i4"),
        (("Integral in ADC x samples", "area"), "<i4"),
        (("Pulse area in PE", "area_pe"), "<f4"),
        (("Sample index in which hit starts", "left"), "<i2"),
        (
            (
                "Sample index in which hit area succeeds 10% of total area",
                "sample_10pc_area",
            ),
            "<i2",
        ),
        (
            (
                "Sample index in which hit area succeeds 50% of total area",
                "sample_50pc_area",
            ),
            "<i2",
        ),
        (("Sample index of hit maximum", "max"), "<i2"),
        (
            (
                "Index of first sample in record just beyond hit (exclusive bound)",
                "right",
            ),
            "<i2",
        ),
        (("Height of hit in ADC counts", "height"), "<i4"),
        (("Height of hit in PE", "height_pe"), "<f4"),
        (("Delay of hit w.r.t. LED hit in same WF, in samples", "tdelay"), "<i2"),
        (
            (
                "Internal (temporary) index of fragment in which hit was found",
                "record_i",
            ),
            "<i4",
        ),
        (
            ("Index of sample in record where integration starts", "left_integration"),
            np.int16,
        ),
        (
            ("Index of first sample beyond integration region", "right_integration"),
            np.int16,
        ),
        (("ADC threshold applied in order to find hits", "threshold"), np.float32),
    ]
    return dtype_ap

1	import numba	2✔
2	import numpy as np	2✔
3	import strax	2✔
4	import straxen	2✔
5
6	export, __all__ = strax.exporter()	2✔
7
8
9	@export	2✔
10	class LEDAfterpulseProcessing(strax.Plugin):	2✔
11	"""Plugin for processing LED afterpulses.
12
13	Detect LED pulses and afterpulses (APs) in raw_records waveforms. Compute the AP datatype.
14
15	"""
16
17	__version__ = "0.6.0"	2✔
18	depends_on = "raw_records"	2✔
19	data_kind = "afterpulses"	2✔
20	provides = "afterpulses"	2✔
21	compressor = "zstd"	2✔
22	parallel = "process"	2✔
23	rechunk_on_save = True	2✔
24
25	gain_model = straxen.URLConfig(	2✔
26	infer_type=False,
27	help="PMT gain model. Specify as (model_type, model_config)",
28	)
29
30	n_tpc_pmts = straxen.URLConfig(	2✔
31	type=int,
32	help="Number of PMTs in TPC",
33	)
34
35	LED_hit_left_boundary = straxen.URLConfig(	2✔
36	default=50,
37	infer_type=False,
38	help="Left boundary after which the first hit marks the position of the LED window",
39	)
40
41	LED_hit_right_boundary = straxen.URLConfig(	2✔
42	default=110,
43	infer_type=False,
44	help="Right boundary after which the LED hit will no longer be searched",
45	)
46
47	LED_window_width = straxen.URLConfig(	2✔
48	default=20,
49	infer_type=False,
50	help="Width of the window in which hits are merged into the LED hit after the first hit",
51	)
52
53	baseline_samples = straxen.URLConfig(	2✔
54	default=40,
55	infer_type=False,
56	help="Number of samples to use at start of WF to determine the baseline",
57	)
58
59	hit_min_amplitude = straxen.URLConfig(	2✔
60	track=True,
61	infer_type=False,
62	default="cmt://hit_thresholds_tpc?version=ONLINE&run_id=plugin.run_id",
63	help=(
64	"Minimum hit amplitude in ADC counts above baseline. "
65	"Specify as a tuple of length n_tpc_pmts, or a number,"
66	'or a string like "legacy-thresholds://pmt_commissioning_initial" which means calling'
67	"hitfinder_thresholds.py"
68	'or url string like "cmt://hit_thresholds_tpc?version=ONLINE" which means'
69	"calling cmt."
70	),
71	)
72
73	hit_min_height_over_noise = straxen.URLConfig(	2✔
74	default=4,
75	infer_type=False,
76	help=(
77	"Minimum hit amplitude in numbers of baseline_rms above baseline."
78	"Actual threshold used is max(hit_min_amplitude, hit_min_"
79	"height_over_noise * baseline_rms)."
80	),
81	)
82
83	save_outside_hits = straxen.URLConfig(	2✔
84	default=(3, 20),
85	infer_type=False,
86	help="Save (left, right) samples besides hits; cut the rest",
87	)
88
89	def infer_dtype(self):	2✔
90	dtype = dtype_afterpulses()	2✔
91	return dtype	2✔
92
93	def setup(self):	2✔
94	self.to_pe = self.gain_model	2✔
95	self.hit_thresholds = self.hit_min_amplitude	2✔
96	self.hit_left_extension, self.hit_right_extension = self.save_outside_hits	2✔
97
98	def compute(self, raw_records):	2✔
99	# Convert everything to the records data type -- adds extra fields.
100	records = strax.raw_to_records(raw_records)	2✔
101	del raw_records	2✔
102
103	# calculate baseline and baseline rms
104	strax.baseline(records, baseline_samples=self.baseline_samples, flip=True)	2✔
105
106	# find all hits
107	hits = strax.find_hits(	2✔
108	records,
109	min_amplitude=self.hit_thresholds,
110	min_height_over_noise=self.hit_min_height_over_noise,
111	)
112
113	# sort hits by record_i and time, then find LED hit and afterpulse
114	# hits within the same record
115	hits_ap = find_ap(	2✔
116	hits,
117	records,
118	LED_hit_left_boundary=self.LED_hit_left_boundary,
119	LED_hit_right_boundary=self.LED_hit_right_boundary,
120	LED_window_width=self.LED_window_width,
121	hit_left_extension=self.hit_left_extension,
122	hit_right_extension=self.hit_right_extension,
123	)
124
125	hits_ap["area_pe"] = hits_ap["area"] * self.to_pe[hits_ap["channel"]]	2✔
126	hits_ap["height_pe"] = hits_ap["height"] * self.to_pe[hits_ap["channel"]]	2✔
127
128	return hits_ap	2✔
129
130
131	@export	2✔
132	def find_ap(	2✔
133	hits,
134	records,
135	LED_hit_left_boundary,
136	LED_hit_right_boundary,
137	LED_window_width,
138	hit_left_extension,
139	hit_right_extension,
140	):
141	"""Find afterpulses (APs) in the given hits data within specified LED hit boundaries and
142	extensions.
143
144	Parameters
145	----------
146	hits :
147	Array containing hit data.
148	records :
149	Array containing record data.
150	LED_hit_left_boundary :
151	Left boundary of the LED hit window.
152	LED_hit_right_boundary :
153	Right boundary of the LED hit window.
154	LED_window_width :
155	Extension to the right of the first hit found in the LED hit window
156	within which hits are merged into the LED hit.
157	hit_left_extension :
158	Extension to the left of the hit window.
159	hit_right_extension :
160	Extension to the right of the hit window.
161
162	Returns
163	-------
164	Array containing afterpulse data.
165
166	Notes
167	-----
168	- Hits to the left of the LED_hit_left_boundary are ignored.
169	- If no hit is found between LED_hit_left_boundary and LED_hit_right_boundary
170	the record is skipped.
171	- The merged LED hits are also saved and can be selected for by having
172	t_delay = 0 by definition.
173
174	"""
175
176	buffer = np.zeros(len(hits), dtype=dtype_afterpulses())	2✔
177
178	if not len(hits):	2✔
179	return buffer	2✔
180
181	# sort hits first by record_i, then by time
182	hits_sorted = np.sort(hits, order=("record_i", "time"))	2✔
183	res = _find_ap(	2✔
184	hits_sorted,
185	records,
186	LED_hit_left_boundary,
187	LED_hit_right_boundary,
188	LED_window_width,
189	hit_left_extension,
190	hit_right_extension,
191	buffer=buffer,
192	)
193	return res	2✔
194
195
196	@numba.jit(nopython=True, nogil=True, cache=True)	2✔
197	def _find_ap(	2✔
198	hits,
199	records,
200	LED_hit_left_boundary,
201	LED_hit_right_boundary,
202	LED_window_width,
203	hit_left_extension,
204	hit_right_extension,
205	buffer=None,
206	):
207	# hits need to be sorted by record_i, then time!
208	offset = 0	2✔
209
210	is_LED = False	2✔
211	t_LED = None	2✔
212	t_LED_hit = None	2✔
213
214	prev_record_i = hits[0]["record_i"]	2✔
215	record_data = records[prev_record_i]["data"]	2✔
216	record_len = records[prev_record_i]["length"]	2✔
217	baseline_fpart = records[prev_record_i]["baseline"] % 1	2✔
218
219	for h_i, h in enumerate(hits):	2✔
220	if h["record_i"] > prev_record_i:	2✔
221	# start of a new record
222	is_LED = False	2✔
223	# only increment buffer if the old one is not empty! this happens
224	# when no (LED) hit is found in the previous record
225	if not buffer[offset]["time"] == 0:	2✔
226	offset += 1	2✔
227	prev_record_i = h["record_i"]	2✔
228	record_data = records[prev_record_i]["data"]	2✔
229	baseline_fpart = records[prev_record_i]["baseline"] % 1	2✔
230
231	res = buffer[offset]	2✔
232
233	if h["left"] < LED_hit_left_boundary:	2✔
234	# if hit is before LED hit window: discard
235	continue	2✔
236
237	if (h["left"] < LED_hit_right_boundary) and not is_LED:	2✔
238	# this is the first hit in the LED hit window
239	fill_hitpars(	2✔
240	res,
241	h,
242	hit_left_extension,
243	hit_right_extension,
244	record_data,
245	record_len,
246	baseline_fpart,
247	)
248
249	t_LED_hit = res["sample_10pc_area"]	2✔
250	t_LED = res["sample_10pc_area"]	2✔
251	is_LED = True	2✔
252
253	continue	2✔
254
255	if not is_LED:	2✔
256	# No hit found in the LED hit window: skip to next record
NEW 257	continue	×
258
259	if h["left"] < (t_LED_hit + LED_window_width):	2✔
260	# This hit is still inside the LED window: extend the LED hit
261	fill_hitpars(	2✔
262	res,
263	h,
264	hit_left_extension,
265	hit_right_extension,
266	record_data,
267	record_len,
268	baseline_fpart,
269	extend=True,
270	)
271
272	# set the LED time in the current WF
273	t_LED = res["sample_10pc_area"]	2✔
274
275	continue	2✔
276
277	# Here begins a new hit after the LED window
278
279	# if a hit is completely inside the previous hit's right_extension,
280	# then skip it (because it's already included in the previous hit)
281	if h["right"] <= res["right_integration"]:	2✔
282	continue	2✔
283
284	# if a hit only partly overlaps with the previous hit's right_
285	# extension, merge them (extend previous hit by this one)
286	if h["left"] <= res["right_integration"]:	2✔
287	fill_hitpars(	2✔
288	res,
289	h,
290	hit_left_extension,
291	hit_right_extension,
292	record_data,
293	record_len,
294	baseline_fpart,
295	extend=True,
296	)
297
298	res["tdelay"] = res["sample_10pc_area"] - t_LED	2✔
299
300	continue	2✔
301
302	# an actual new hit increases the buffer index
303	offset += 1	2✔
304	res = buffer[offset]	2✔
305
306	fill_hitpars(	2✔
307	res,
308	h,
309	hit_left_extension,
310	hit_right_extension,
311	record_data,
312	record_len,
313	baseline_fpart,
314	)
315
316	res["tdelay"] = res["sample_10pc_area"] - t_LED	2✔
317
318	return buffer[:offset]	2✔
319
320
321	@export	2✔
322	@numba.jit(nopython=True, nogil=True, cache=True)	2✔
323	def get_sample_area_quantile(data, quantile, baseline_fpart):	2✔
324	"""Return the index of the first sample in the hit where the integrated area of the hit to that
325	index is above the specified quantile of the total area.
326
327	Parameters:
328	- data :
329	Array containing the baselined waveform data of the hit.
330	- quantile :
331	The quantile (0 to 1) representing the threshold for the area.
332	- baseline_fpart :
333	Fractional part of the baseline (baseline % 1) of the record
334
335	Return:
336	- int: The index of the first sample where the area exceeds the quantile of the total area.
337	If no such sample is found, returns 0.
338
339	Notes:
340	- If no quantile is found where the area exceeds the threshold, it returns 0. This is
341	usually caused by real events in the baseline window, which can result in a negative
342	area.
343
344	"""
345
346	area = 0	2✔
347	area_tot = data.sum() + len(data) * baseline_fpart	2✔
348
349	for d_i, d in enumerate(data):	2✔
350	area += d + baseline_fpart	2✔
351	if area > (quantile * area_tot):	2✔
352	return d_i	2✔
353	if d_i == len(data) - 1:	2✔
354	# if no quantile was found, something went wrong
355	# (negative area due to wrong baseline, caused by real events that
356	# by coincidence fall in the first samples of the trigger window)
357	# print('no quantile found: set to 0')
358	return 0	2✔
359
360	# What happened here?!
361	return 0	2✔
362
363
364	@numba.jit(nopython=True, nogil=True, cache=True)	2✔
365	def fill_hitpars(	2✔
366	result,
367	hit,
368	hit_left_extension,
369	hit_right_extension,
370	record_data,
371	record_len,
372	baseline_fpart,
373	extend=False,
374	):
375	if not extend: # fill first time only	2✔
376	result["time"] = hit["time"] - hit_left_extension * hit["dt"]	2✔
377	result["dt"] = hit["dt"]	2✔
378	result["channel"] = hit["channel"]	2✔
379	result["left"] = hit["left"]	2✔
380	result["record_i"] = hit["record_i"]	2✔
381	result["threshold"] = hit["threshold"]	2✔
382	result["left_integration"] = hit["left"] - hit_left_extension	2✔
383	result["height"] = hit["height"]	2✔
384
385	# fill always (if hits are merged, only these will be updated)
386	result["right"] = hit["right"]	2✔
387	result["right_integration"] = hit["right"] + hit_right_extension	2✔
388	if result["right_integration"] > record_len:	2✔
389	result["right_integration"] = record_len # cap right_integration at end of record	2✔
390	result["length"] = result["right_integration"] - result["left_integration"]	2✔
391
392	hit_data = record_data[result["left_integration"] : result["right_integration"]]	2✔
393	result["area"] = hit_data.sum() + result["length"] * baseline_fpart	2✔
394	result["sample_10pc_area"] = result["left_integration"] + get_sample_area_quantile(	2✔
395	hit_data, 0.1, baseline_fpart
396	)
397	result["sample_50pc_area"] = result["left_integration"] + get_sample_area_quantile(	2✔
398	hit_data, 0.5, baseline_fpart
399	)
400	if len(hit_data):	2✔
401	result["max"] = result["left_integration"] + hit_data.argmax()	2✔
402
403	if extend: # only when merging hits	2✔
404	result["height"] = max(result["height"], hit["height"])	2✔
405
406
407	@export	2✔
408	def dtype_afterpulses():	2✔
409	"""The afterpulse datatype.
410
411	Return:
412	- The afterpulse datatype
413
414	"""
415	dtype_ap = [	2✔
416	(("Channel/PMT number", "channel"), "<i2"),
417	(("Time resolution in ns", "dt"), "<i2"),
418	(("Start time of the interval (ns since unix epoch)", "time"), "<i8"),
419	(("Length of the interval in samples", "length"), "<i4"),
420	(("Integral in ADC x samples", "area"), "<i4"),
421	(("Pulse area in PE", "area_pe"), "<f4"),
422	(("Sample index in which hit starts", "left"), "<i2"),
423	(
424	(
425	"Sample index in which hit area succeeds 10% of total area",
426	"sample_10pc_area",
427	),
428	"<i2",
429	),
430	(
431	(
432	"Sample index in which hit area succeeds 50% of total area",
433	"sample_50pc_area",
434	),
435	"<i2",
436	),
437	(("Sample index of hit maximum", "max"), "<i2"),
438	(
439	(
440	"Index of first sample in record just beyond hit (exclusive bound)",
441	"right",
442	),
443	"<i2",
444	),
445	(("Height of hit in ADC counts", "height"), "<i4"),
446	(("Height of hit in PE", "height_pe"), "<f4"),
447	(("Delay of hit w.r.t. LED hit in same WF, in samples", "tdelay"), "<i2"),
448	(
449	(
450	"Internal (temporary) index of fragment in which hit was found",
451	"record_i",
452	),
453	"<i4",
454	),
455	(
456	("Index of sample in record where integration starts", "left_integration"),
457	np.int16,
458	),
459	(
460	("Index of first sample beyond integration region", "right_integration"),
461	np.int16,
462	),
463	(("ADC threshold applied in order to find hits", "threshold"), np.float32),
464	]
465	return dtype_ap	2✔

XENONnT / straxen / 10168286066

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous