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

84.46

/straxen/plugins/led_cal/led_calibration.py

"""
Dear nT analyser,
if you want to complain please contact:
    chiara@physik.uzh.ch,
    gvolta@physik.uzh.ch,
    kazama@isee.nagoya-u.ac.jp
    torben.flehmke@fysik.su.se
"""

from immutabledict import immutabledict
import strax
import straxen
import numba
import numpy as np

# This makes sure shorthands for only the necessary functions
# are made available under straxen.[...]
export, __all__ = strax.exporter()

channel_list = [i for i in range(494)]


@export
class LEDCalibration(strax.Plugin):
    """
    Preliminary version, several parameters to set during commissioning.
    LEDCalibration returns: channel, time, dt, length, Area,
    amplitudeLED and amplitudeNOISE.

    The new variables are:
        - Area: Area computed in the given window, averaged over 6
          windows that have the same starting sample and different end
          samples.
        - amplitudeLED: peak amplitude of the LED on run in the given
          window.
        - amplitudeNOISE: amplitude of the LED on run in a window far
          from the signal one.
    """

    __version__ = "0.3.0"

    depends_on = "raw_records"
    data_kind = "led_cal"
    compressor = "zstd"
    parallel = "process"
    rechunk_on_save = False

    baseline_window = straxen.URLConfig(
        default=(0, 40),
        infer_type=False,
        help="Window (samples) for baseline calculation.",
    )

    default_led_window = straxen.URLConfig(
        default=(78, 107),
        infer_type=False,
        help="Default window (samples) to integrate and get the maximum amplitude \
            if no hit was found in the record.",
    )

    led_hit_extension = straxen.URLConfig(
        default=(-5, 24),
        infer_type=False,
        help="The extension around the LED hit to integrate.",
    )

    noise_window = straxen.URLConfig(
        default=(10, 48), infer_type=False, help="Window (samples) to analyse the noise"
    )

    channel_list = straxen.URLConfig(
        default=(tuple(channel_list)),
        infer_type=False,
        help="List of PMTs. Defalt value: all the PMTs",
    )

    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)."
        ),
    )

    dtype = [
        ("area", np.float32, "Area averaged in integration windows"),
        ("amplitude_led", np.float32, "Amplitude in LED window"),
        ("amplitude_noise", np.float32, "Amplitude in off LED window"),
        ("channel", np.int16, "Channel"),
        ("time", np.int64, "Start time of the interval (ns since unix epoch)"),
        ("dt", np.int16, "Time resolution in ns"),
        ("length", np.int32, "Length of the interval in samples"),
    ]

    def compute(self, raw_records):
        """The data for LED calibration are build for those PMT which belongs to channel list.

        This is used for the different ligh levels. As default value all the PMTs are considered.

        """
        mask = np.where(np.in1d(raw_records["channel"], self.channel_list))[0]
        rr = raw_records[mask]
        r = get_records(rr, baseline_window=self.baseline_window)
        del rr, raw_records

        temp = np.zeros(len(r), dtype=self.dtype)
        strax.copy_to_buffer(r, temp, "_recs_to_temp_led")

        led_windows = get_led_windows(
            r,
            self.default_led_window,
            self.led_hit_extension,
            self.hit_min_amplitude,
            self.hit_min_height_over_noise,
        )

        on, off = get_amplitude(r, led_windows, self.noise_window)
        temp["amplitude_led"] = on["amplitude"]
        temp["amplitude_noise"] = off["amplitude"]

        area = get_area(r, led_windows)
        temp["area"] = area["area"]
        return temp


def get_records(raw_records, baseline_window):
    """Determine baseline as the average of the first baseline_samples of each pulse.

    Subtract the pulse float(data) from baseline.

    """

    record_length = np.shape(raw_records.dtype["data"])[0]

    _dtype = [
        (("Start time since unix epoch [ns]", "time"), "<i8"),
        (("Length of the interval in samples", "length"), "<i4"),
        (("Width of one sample [ns]", "dt"), "<i2"),
        (("Channel/PMT number", "channel"), "<i2"),
        (
            (
                "Length of pulse to which the record belongs (without zero-padding)",
                "pulse_length",
            ),
            "<i4",
        ),
        (("Fragment number in the pulse", "record_i"), "<i2"),
        (
            ("Baseline in ADC counts. data = int(baseline) - data_orig", "baseline"),
            "f4",
        ),
        (
            ("Baseline RMS in ADC counts. data = baseline - data_orig", "baseline_rms"),
            "f4",
        ),
        (("Waveform data in raw ADC counts", "data"), "f4", (record_length,)),
    ]

    records = np.zeros(len(raw_records), dtype=_dtype)
    strax.copy_to_buffer(raw_records, records, "_rr_to_r_led")

    mask = np.where((records["record_i"] == 0) & (records["length"] == 160))[0]
    records = records[mask]
    bl = records["data"][:, baseline_window[0] : baseline_window[1]].mean(axis=1)
    rms = records["data"][:, baseline_window[0] : baseline_window[1]].std(axis=1)
    records["data"][:, :160] = -1.0 * (records["data"][:, :160].transpose() - bl[:]).transpose()
    records["baseline"] = bl
    records["baseline_rms"] = rms
    return records


def get_led_windows(
    records,
    default_window,
    led_hit_extension,
    hit_min_amplitude,
    hit_min_height_over_noise,
):
    """Search for hits in the records, if a hit is found, return an interval around the hit given by
    led_hit_extension. If no hit is found in the record, return the default window.

    :param records: Array of the records to search for LED hits.
    :param default_window: Default window to use if no LED hit is found given as a tuple (start,
        end)
    :param led_hit_extension: The integration window around the first hit found to use. A tuple of
        form (samples_before, samples_after) the first LED hit.
    :param hit_min_amplitude: Minimum amplitude of the signal to be considered a hit.
    :param hit_min_height_over_noise: Minimum height of the signal over noise to be considered a
        hit. :return (len(records), 2) array: Integration window for each record

    """
    hits = strax.find_hits(
        records,
        min_amplitude=hit_min_amplitude,
        min_height_over_noise=hit_min_height_over_noise,
    )
    default_windows = np.tile(default_window, (len(records), 1))
    return _get_led_windows(hits, default_windows, led_hit_extension)


@numba.jit(nopython=True)
def _get_led_windows(hits, default_windows, led_hit_extension):
    windows = default_windows
    last = -1
    for hit in hits:
        if hit["record_i"] == last:
            continue  # If there are multiple hits in one record, ignore after the first

        left = hit["left"] + led_hit_extension[0]
        # Limit the position of the window so it stays inside the record.
        if left < default_windows[0, 0]:
            left = default_windows[0, 0]
        elif left > 150 - (led_hit_extension[1] - led_hit_extension[0]):
            left = 150 - (led_hit_extension[1] - led_hit_extension[0])

        right = hit["left"] + led_hit_extension[1]
        if right < default_windows[0, 1]:
            right = default_windows[0, 1]
        elif right > 150:
            right = 150

        windows[hit["record_i"]] = np.array([left, right])
        last = hit["record_i"]

    return windows


_on_off_dtype = np.dtype([("channel", "int16"), ("amplitude", "float32")])


@numba.jit(nopython=True)
def get_amplitude(records, led_windows, noise_window):
    """Needed for the SPE computation.

    Get the maximum of the signal in two different regions, one where there is no signal, and one
    where there is.

    :param records: Array of records
    :param ndarray led_windows : 2d array of shape (len(records), 2) with the window to use as the
        signal on area for each record. Inclusive left boundary and exclusive right boundary.
    :param tuple noise_window: Tuple with the window, used for the signal off area for all records.
    :return ndarray ons: 1d array of length len(records). The maximum amplitude in the led window
        area for each record.
    :return ndarray offs: 1d array of length len(records). The maximum amplitude in the noise area
        for each record.

    """
    ons = np.zeros(len(records), dtype=_on_off_dtype)
    offs = np.zeros(len(records), dtype=_on_off_dtype)

    for i, record in enumerate(records):
        ons[i]["channel"] = record["channel"]
        offs[i]["channel"] = record["channel"]

        ons[i]["amplitude"] = np.max(record["data"][led_windows[i, 0] : led_windows[i, 1]])
        offs[i]["amplitude"] = np.max(record["data"][noise_window[0] : noise_window[1]])

    return ons, offs


_area_dtype = np.dtype([("channel", "int16"), ("area", "float32")])


@numba.jit(nopython=True)
def get_area(records, led_windows):
    """Needed for the gain computation.

    Integrate the record in the defined window area. To reduce the effects of the noise, this is
    done with 6 different window lengths, which are then averaged.

    :param records: Array of records
    :param ndarray led_windows : 2d array of shape (len(records), 2) with the window to use as the
        integration boundaries.
    :return ndarray area: 1d array of length len(records) with the averaged integrated areas for
        each record.

    """
    area = np.zeros(len(records), dtype=_area_dtype)
    end_pos = [2 * i for i in range(6)]

    for i, record in enumerate(records):
        area[i]["channel"] = record["channel"]
        for right in end_pos:
            area[i]["area"] += np.sum(record["data"][led_windows[i, 0] : led_windows[i, 1] + right])

        area[i]["area"] /= float(len(end_pos))

    return area


@export
class nVetoExtTimings(strax.Plugin):
    """Plugin which computes the time difference `delta_time` from pulse timing of `hitlets_nv` to
    start time of `raw_records` which belong the `hitlets_nv`.

    They are used as the external trigger timings.

    """

    __version__ = "0.0.1"

    depends_on = ("raw_records_nv", "hitlets_nv")
    provides = "ext_timings_nv"
    data_kind = "hitlets_nv"

    compressor = "zstd"

    channel_map = straxen.URLConfig(
        track=False,
        type=immutabledict,
        help="immutabledict mapping subdetector to (min, max) channel number.",
    )

    def infer_dtype(self):
        dtype = []
        dtype += strax.time_dt_fields
        dtype += [
            (("Delta time from trigger timing [ns]", "delta_time"), np.int16),
            (
                ("Index to which pulse (not record) the hitlet belongs to.", "pulse_i"),
                np.int32,
            ),
        ]
        return dtype

    def setup(self):
        self.nv_pmt_start = self.channel_map["nveto"][0]
        self.nv_pmt_stop = self.channel_map["nveto"][1] + 1

    def compute(self, hitlets_nv, raw_records_nv):
        rr_nv = raw_records_nv[raw_records_nv["record_i"] == 0]
        pulses = np.zeros(len(rr_nv), dtype=self.pulse_dtype())
        pulses["time"] = rr_nv["time"]
        pulses["endtime"] = rr_nv["time"] + rr_nv["pulse_length"] * rr_nv["dt"]
        pulses["channel"] = rr_nv["channel"]

        ext_timings_nv = np.zeros_like(hitlets_nv, dtype=self.dtype)
        ext_timings_nv["time"] = hitlets_nv["time"]
        ext_timings_nv["length"] = hitlets_nv["length"]
        ext_timings_nv["dt"] = hitlets_nv["dt"]
        self.calc_delta_time(
            ext_timings_nv, pulses, hitlets_nv, self.nv_pmt_start, self.nv_pmt_stop
        )

        return ext_timings_nv

    @staticmethod
    def pulse_dtype():
        pulse_dtype = []
        pulse_dtype += strax.time_fields
        pulse_dtype += [(("PMT channel", "channel"), np.int16)]
        return pulse_dtype

    @staticmethod
    def calc_delta_time(ext_timings_nv_delta_time, pulses, hitlets_nv, nv_pmt_start, nv_pmt_stop):
        """Numpy access with fancy index returns copy, not view This for-loop is required to
        substitute in one by one."""
        hitlet_index = np.arange(len(hitlets_nv))
        pulse_index = np.arange(len(pulses))
        for ch in range(nv_pmt_start, nv_pmt_stop):
            mask_hitlets_in_channel = hitlets_nv["channel"] == ch
            hitlet_in_channel_index = hitlet_index[mask_hitlets_in_channel]

            mask_pulse_in_channel = pulses["channel"] == ch
            pulse_in_channel_index = pulse_index[mask_pulse_in_channel]

            hitlets_in_channel = hitlets_nv[hitlet_in_channel_index]
            pulses_in_channel = pulses[pulse_in_channel_index]
            hit_in_pulse_index = strax.fully_contained_in(hitlets_in_channel, pulses_in_channel)
            for h_i, p_i in zip(hitlet_in_channel_index, hit_in_pulse_index):
                if p_i == -1:
                    continue
                res = ext_timings_nv_delta_time[h_i]

                res["delta_time"] = (
                    hitlets_nv[h_i]["time"]
                    + hitlets_nv[h_i]["time_amplitude"]
                    - pulses_in_channel[p_i]["time"]
                )
                res["pulse_i"] = pulse_in_channel_index[p_i]

1	"""
2	Dear nT analyser,
3	if you want to complain please contact:
4	chiara@physik.uzh.ch,
5	gvolta@physik.uzh.ch,
6	kazama@isee.nagoya-u.ac.jp
7	torben.flehmke@fysik.su.se
8	"""
9
10	from immutabledict import immutabledict	2✔
11	import strax	2✔
12	import straxen	2✔
13	import numba	2✔
14	import numpy as np	2✔
15
16	# This makes sure shorthands for only the necessary functions
17	# are made available under straxen.[...]
18	export, __all__ = strax.exporter()	2✔
19
20	channel_list = [i for i in range(494)]	2✔
21
22
23	@export	2✔
24	class LEDCalibration(strax.Plugin):	2✔
25	"""
26	Preliminary version, several parameters to set during commissioning.
27	LEDCalibration returns: channel, time, dt, length, Area,
28	amplitudeLED and amplitudeNOISE.
29
30	The new variables are:
31	- Area: Area computed in the given window, averaged over 6
32	windows that have the same starting sample and different end
33	samples.
34	- amplitudeLED: peak amplitude of the LED on run in the given
35	window.
36	- amplitudeNOISE: amplitude of the LED on run in a window far
37	from the signal one.
38	"""
39
40	__version__ = "0.3.0"	2✔
41
42	depends_on = "raw_records"	2✔
43	data_kind = "led_cal"	2✔
44	compressor = "zstd"	2✔
45	parallel = "process"	2✔
46	rechunk_on_save = False	2✔
47
48	baseline_window = straxen.URLConfig(	2✔
49	default=(0, 40),
50	infer_type=False,
51	help="Window (samples) for baseline calculation.",
52	)
53
54	default_led_window = straxen.URLConfig(	2✔
55	default=(78, 107),
56	infer_type=False,
57	help="Default window (samples) to integrate and get the maximum amplitude \
58	if no hit was found in the record.",
59	)
60
61	led_hit_extension = straxen.URLConfig(	2✔
62	default=(-5, 24),
63	infer_type=False,
64	help="The extension around the LED hit to integrate.",
65	)
66
67	noise_window = straxen.URLConfig(	2✔
68	default=(10, 48), infer_type=False, help="Window (samples) to analyse the noise"
69	)
70
71	channel_list = straxen.URLConfig(	2✔
72	default=(tuple(channel_list)),
73	infer_type=False,
74	help="List of PMTs. Defalt value: all the PMTs",
75	)
76
77	hit_min_amplitude = straxen.URLConfig(	2✔
78	track=True,
79	infer_type=False,
80	default="cmt://hit_thresholds_tpc?version=ONLINE&run_id=plugin.run_id",
81	help=(
82	"Minimum hit amplitude in ADC counts above baseline. "
83	"Specify as a tuple of length n_tpc_pmts, or a number,"
84	'or a string like "legacy-thresholds://pmt_commissioning_initial" which means calling'
85	"hitfinder_thresholds.py"
86	'or url string like "cmt://hit_thresholds_tpc?version=ONLINE" which means'
87	"calling cmt."
88	),
89	)
90
91	hit_min_height_over_noise = straxen.URLConfig(	2✔
92	default=4,
93	infer_type=False,
94	help=(
95	"Minimum hit amplitude in numbers of baseline_rms above baseline."
96	"Actual threshold used is max(hit_min_amplitude, hit_min_"
97	"height_over_noise * baseline_rms)."
98	),
99	)
100
101	dtype = [	2✔
102	("area", np.float32, "Area averaged in integration windows"),
103	("amplitude_led", np.float32, "Amplitude in LED window"),
104	("amplitude_noise", np.float32, "Amplitude in off LED window"),
105	("channel", np.int16, "Channel"),
106	("time", np.int64, "Start time of the interval (ns since unix epoch)"),
107	("dt", np.int16, "Time resolution in ns"),
108	("length", np.int32, "Length of the interval in samples"),
109	]
110
111	def compute(self, raw_records):	2✔
112	"""The data for LED calibration are build for those PMT which belongs to channel list.
113
114	This is used for the different ligh levels. As default value all the PMTs are considered.
115
116	"""
117	mask = np.where(np.in1d(raw_records["channel"], self.channel_list))[0]	2✔
118	rr = raw_records[mask]	2✔
119	r = get_records(rr, baseline_window=self.baseline_window)	2✔
120	del rr, raw_records	2✔
121
122	temp = np.zeros(len(r), dtype=self.dtype)	2✔
123	strax.copy_to_buffer(r, temp, "_recs_to_temp_led")	2✔
124
125	led_windows = get_led_windows(	2✔
126	r,
127	self.default_led_window,
128	self.led_hit_extension,
129	self.hit_min_amplitude,
130	self.hit_min_height_over_noise,
131	)
132
133	on, off = get_amplitude(r, led_windows, self.noise_window)	2✔
134	temp["amplitude_led"] = on["amplitude"]	2✔
135	temp["amplitude_noise"] = off["amplitude"]	2✔
136
137	area = get_area(r, led_windows)	2✔
138	temp["area"] = area["area"]	2✔
139	return temp	2✔
140
141
142	def get_records(raw_records, baseline_window):	2✔
143	"""Determine baseline as the average of the first baseline_samples of each pulse.
144
145	Subtract the pulse float(data) from baseline.
146
147	"""
148
149	record_length = np.shape(raw_records.dtype["data"])[0]	2✔
150
151	_dtype = [	2✔
152	(("Start time since unix epoch [ns]", "time"), "<i8"),
153	(("Length of the interval in samples", "length"), "<i4"),
154	(("Width of one sample [ns]", "dt"), "<i2"),
155	(("Channel/PMT number", "channel"), "<i2"),
156	(
157	(
158	"Length of pulse to which the record belongs (without zero-padding)",
159	"pulse_length",
160	),
161	"<i4",
162	),
163	(("Fragment number in the pulse", "record_i"), "<i2"),
164	(
165	("Baseline in ADC counts. data = int(baseline) - data_orig", "baseline"),
166	"f4",
167	),
168	(
169	("Baseline RMS in ADC counts. data = baseline - data_orig", "baseline_rms"),
170	"f4",
171	),
172	(("Waveform data in raw ADC counts", "data"), "f4", (record_length,)),
173	]
174
175	records = np.zeros(len(raw_records), dtype=_dtype)	2✔
176	strax.copy_to_buffer(raw_records, records, "_rr_to_r_led")	2✔
177
178	mask = np.where((records["record_i"] == 0) & (records["length"] == 160))[0]	2✔
179	records = records[mask]	2✔
180	bl = records["data"][:, baseline_window[0] : baseline_window[1]].mean(axis=1)	2✔
181	rms = records["data"][:, baseline_window[0] : baseline_window[1]].std(axis=1)	2✔
182	records["data"][:, :160] = -1.0 * (records["data"][:, :160].transpose() - bl[:]).transpose()	2✔
183	records["baseline"] = bl	2✔
184	records["baseline_rms"] = rms	2✔
185	return records	2✔
186
187
188	def get_led_windows(	2✔
189	records,
190	default_window,
191	led_hit_extension,
192	hit_min_amplitude,
193	hit_min_height_over_noise,
194	):
195	"""Search for hits in the records, if a hit is found, return an interval around the hit given by
196	led_hit_extension. If no hit is found in the record, return the default window.
197
198	:param records: Array of the records to search for LED hits.
199	:param default_window: Default window to use if no LED hit is found given as a tuple (start,
200	end)
201	:param led_hit_extension: The integration window around the first hit found to use. A tuple of
202	form (samples_before, samples_after) the first LED hit.
203	:param hit_min_amplitude: Minimum amplitude of the signal to be considered a hit.
204	:param hit_min_height_over_noise: Minimum height of the signal over noise to be considered a
205	hit. :return (len(records), 2) array: Integration window for each record
206
207	"""
208	hits = strax.find_hits(	2✔
209	records,
210	min_amplitude=hit_min_amplitude,
211	min_height_over_noise=hit_min_height_over_noise,
212	)
213	default_windows = np.tile(default_window, (len(records), 1))	2✔
214	return _get_led_windows(hits, default_windows, led_hit_extension)	2✔
215
216
217	@numba.jit(nopython=True)	2✔
218	def _get_led_windows(hits, default_windows, led_hit_extension):	2✔
219	windows = default_windows	2✔
220	last = -1	2✔
221	for hit in hits:	2✔
NEW 222	if hit["record_i"] == last:	×
NEW 223	continue # If there are multiple hits in one record, ignore after the first	×
224
NEW 225	left = hit["left"] + led_hit_extension[0]	×
226	# Limit the position of the window so it stays inside the record.
NEW 227	if left < default_windows[0, 0]:	×
NEW 228	left = default_windows[0, 0]	×
NEW 229	elif left > 150 - (led_hit_extension[1] - led_hit_extension[0]):	×
NEW 230	left = 150 - (led_hit_extension[1] - led_hit_extension[0])	×
231
NEW 232	right = hit["left"] + led_hit_extension[1]	×
NEW 233	if right < default_windows[0, 1]:	×
NEW 234	right = default_windows[0, 1]	×
NEW 235	elif right > 150:	×
NEW 236	right = 150	×
237
NEW 238	windows[hit["record_i"]] = np.array([left, right])	×
NEW 239	last = hit["record_i"]	×
240
241	return windows	2✔
242
243
244	_on_off_dtype = np.dtype([("channel", "int16"), ("amplitude", "float32")])	2✔
245
246
247	@numba.jit(nopython=True)	2✔
248	def get_amplitude(records, led_windows, noise_window):	2✔
249	"""Needed for the SPE computation.
250
251	Get the maximum of the signal in two different regions, one where there is no signal, and one
252	where there is.
253
254	:param records: Array of records
255	:param ndarray led_windows : 2d array of shape (len(records), 2) with the window to use as the
256	signal on area for each record. Inclusive left boundary and exclusive right boundary.
257	:param tuple noise_window: Tuple with the window, used for the signal off area for all records.
258	:return ndarray ons: 1d array of length len(records). The maximum amplitude in the led window
259	area for each record.
260	:return ndarray offs: 1d array of length len(records). The maximum amplitude in the noise area
261	for each record.
262
263	"""
264	ons = np.zeros(len(records), dtype=_on_off_dtype)	2✔
265	offs = np.zeros(len(records), dtype=_on_off_dtype)	2✔
266
267	for i, record in enumerate(records):	2✔
NEW 268	ons[i]["channel"] = record["channel"]	×
NEW 269	offs[i]["channel"] = record["channel"]	×
270
NEW 271	ons[i]["amplitude"] = np.max(record["data"][led_windows[i, 0] : led_windows[i, 1]])	×
NEW 272	offs[i]["amplitude"] = np.max(record["data"][noise_window[0] : noise_window[1]])	×
273
274	return ons, offs	2✔
275
276
277	_area_dtype = np.dtype([("channel", "int16"), ("area", "float32")])	2✔
278
279
280	@numba.jit(nopython=True)	2✔
281	def get_area(records, led_windows):	2✔
282	"""Needed for the gain computation.
283
284	Integrate the record in the defined window area. To reduce the effects of the noise, this is
285	done with 6 different window lengths, which are then averaged.
286
287	:param records: Array of records
288	:param ndarray led_windows : 2d array of shape (len(records), 2) with the window to use as the
289	integration boundaries.
290	:return ndarray area: 1d array of length len(records) with the averaged integrated areas for
291	each record.
292
293	"""
294	area = np.zeros(len(records), dtype=_area_dtype)	2✔
295	end_pos = [2 * i for i in range(6)]	2✔
296
297	for i, record in enumerate(records):	2✔
NEW 298	area[i]["channel"] = record["channel"]	×
NEW 299	for right in end_pos:	×
NEW 300	area[i]["area"] += np.sum(record["data"][led_windows[i, 0] : led_windows[i, 1] + right])	×
301
NEW 302	area[i]["area"] /= float(len(end_pos))	×
303
304	return area	2✔
305
306
307	@export	2✔
308	class nVetoExtTimings(strax.Plugin):	2✔
309	"""Plugin which computes the time difference `delta_time` from pulse timing of `hitlets_nv` to
310	start time of `raw_records` which belong the `hitlets_nv`.
311
312	They are used as the external trigger timings.
313
314	"""
315
316	__version__ = "0.0.1"	2✔
317
318	depends_on = ("raw_records_nv", "hitlets_nv")	2✔
319	provides = "ext_timings_nv"	2✔
320	data_kind = "hitlets_nv"	2✔
321
322	compressor = "zstd"	2✔
323
324	channel_map = straxen.URLConfig(	2✔
325	track=False,
326	type=immutabledict,
327	help="immutabledict mapping subdetector to (min, max) channel number.",
328	)
329
330	def infer_dtype(self):	2✔
331	dtype = []	2✔
332	dtype += strax.time_dt_fields	2✔
333	dtype += [	2✔
334	(("Delta time from trigger timing [ns]", "delta_time"), np.int16),
335	(
336	("Index to which pulse (not record) the hitlet belongs to.", "pulse_i"),
337	np.int32,
338	),
339	]
340	return dtype	2✔
341
342	def setup(self):	2✔
343	self.nv_pmt_start = self.channel_map["nveto"][0]	2✔
344	self.nv_pmt_stop = self.channel_map["nveto"][1] + 1	2✔
345
346	def compute(self, hitlets_nv, raw_records_nv):	2✔
347	rr_nv = raw_records_nv[raw_records_nv["record_i"] == 0]	2✔
348	pulses = np.zeros(len(rr_nv), dtype=self.pulse_dtype())	2✔
349	pulses["time"] = rr_nv["time"]	2✔
350	pulses["endtime"] = rr_nv["time"] + rr_nv["pulse_length"] * rr_nv["dt"]	2✔
351	pulses["channel"] = rr_nv["channel"]	2✔
352
353	ext_timings_nv = np.zeros_like(hitlets_nv, dtype=self.dtype)	2✔
354	ext_timings_nv["time"] = hitlets_nv["time"]	2✔
355	ext_timings_nv["length"] = hitlets_nv["length"]	2✔
356	ext_timings_nv["dt"] = hitlets_nv["dt"]	2✔
357	self.calc_delta_time(	2✔
358	ext_timings_nv, pulses, hitlets_nv, self.nv_pmt_start, self.nv_pmt_stop
359	)
360
361	return ext_timings_nv	2✔
362
363	@staticmethod	2✔
364	def pulse_dtype():	2✔
365	pulse_dtype = []	2✔
366	pulse_dtype += strax.time_fields	2✔
367	pulse_dtype += [(("PMT channel", "channel"), np.int16)]	2✔
368	return pulse_dtype	2✔
369
370	@staticmethod	2✔
371	def calc_delta_time(ext_timings_nv_delta_time, pulses, hitlets_nv, nv_pmt_start, nv_pmt_stop):	2✔
372	"""Numpy access with fancy index returns copy, not view This for-loop is required to
373	substitute in one by one."""
374	hitlet_index = np.arange(len(hitlets_nv))	2✔
375	pulse_index = np.arange(len(pulses))	2✔
376	for ch in range(nv_pmt_start, nv_pmt_stop):	2✔
377	mask_hitlets_in_channel = hitlets_nv["channel"] == ch	2✔
378	hitlet_in_channel_index = hitlet_index[mask_hitlets_in_channel]	2✔
379
380	mask_pulse_in_channel = pulses["channel"] == ch	2✔
381	pulse_in_channel_index = pulse_index[mask_pulse_in_channel]	2✔
382
383	hitlets_in_channel = hitlets_nv[hitlet_in_channel_index]	2✔
384	pulses_in_channel = pulses[pulse_in_channel_index]	2✔
385	hit_in_pulse_index = strax.fully_contained_in(hitlets_in_channel, pulses_in_channel)	2✔
386	for h_i, p_i in zip(hitlet_in_channel_index, hit_in_pulse_index):	2✔
387	if p_i == -1:	2✔
388	continue	×
389	res = ext_timings_nv_delta_time[h_i]	2✔
390
391	res["delta_time"] = (	2✔
392	hitlets_nv[h_i]["time"]
393	+ hitlets_nv[h_i]["time_amplitude"]
394	- pulses_in_channel[p_i]["time"]
395	)
396	res["pulse_i"] = pulse_in_channel_index[p_i]	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