11867531211

Committed 16 Nov 2024 05:23AM UTC coverage: 89.785% (+0.02%) from 89.77%

Build # 11867531211

Build Type

push

github

Committed by

web-flow

Commit Message

Merge branch 'sr1_leftovers' into master (#1478)

* Update pytest.yml (#1431)

* Update pytest.yml

Specify the strax to be v1.6.5

* add install base_env

* Add force reinstall

* Update definition of the SE Score (previously the SE density) for SR1 WIMP (#1430)

* Update score definition. Modify file names.

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Modify file names in the initialization file.

* Rearrangenames. Move sr phase assignment elsewhere.

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Add url configs and modify code style.

* Modify the parameter names.

* Fix data type in url config.

* Add docstring for the eps used to prevent divide by zero.

* Reformmated with precommit.

* Add docstrings. Remove redundant code.

* Add docstring for the 2D Gaussian.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Dacheng Xu <dx2227@columbia.edu>

* Copy https://github.com/XENONnT/straxen/pull/1417 (#1438)

* Bump to v2.2.6 (#1441)

* Bump version: 2.2.4 → 2.2.6

* Update HISTORY.md

* Constraint strax version

---------

Co-authored-by: Kexin Liu <lkx21@mails.tsinghua.edu.cn>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Run Details

72 of 72 new or added lines in 4 files covered. (100.0%)

7 existing lines in 3 files now uncovered.

8614 of 9594 relevant lines covered (89.79%)

1.8 hits per line

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

77.47

/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
import scipy.stats as sps

# 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.1"

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

    run_doc = straxen.URLConfig(
        default="run_doc://comments?run_id=plugin.run_id",
        infer_type=False,
        help=(
            "Comments dictionary from the run metadata. "
            "It is used for discriminate between LED on and LED "
            "off runs."
        ),
    )

    defualt_run_comments = straxen.URLConfig(
        default=["auto, SC user: ", "pulserScript: "],
        type=list,
        help=("List of default comments the automatic script for runs in PMT calibration. "),
    )

    noise_run_comments = straxen.URLConfig(
        default=["SPE_calibration_step0", "Gain_calibration_step3"],
        type=list,
        help=("List of comments for noise runs in PMT calibration. "),
    )

    led_cal_record_length = straxen.URLConfig(
        default=160, infer_type=False, help="Length (samples) of one record without 0 padding."
    )

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

    minimum_led_position = straxen.URLConfig(
        default=60,
        infer_type=False,
        help=(
            "The minimum sample index to consider for LED hits. Hits before this sample are "
            "ignored."
        ),
    )

    fixed_position = straxen.URLConfig(
        default=88,
        infer_type=False,
        help=(
            "Fixed ADC sample upon which the integration window is defined. "
            "This is used as default when no hits or less than a certain amount "
            "are identified."
        ),
    )

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

    area_averaging_length = straxen.URLConfig(
        default=7,
        infer_type=False,
        help=(
            "The total length of the averaging window for the area calculation. "
            "To mitigate a possiple bias from noise, the area is integrated multiple times with "
            "sligntly different window lengths and then averaged. area_averaging_length should "
            "be divisible by area_averaging_step."
        ),
    )

    area_averaging_step = straxen.URLConfig(
        default=1,
        infer_type=False,
        help=(
            "The step size used for the different windows, averaged for the area calculation. "
            "To mitigate a possiple bias from noise, the area is integrated multiple times with "
            "sligntly different window lengths and then averaged. area_averaging_length should "
            "be divisible by area_averaging_step."
        ),
    )

    noise_window = straxen.URLConfig(
        default=(10, 50), 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",
    )

    led_cal_hit_min_height_over_noise = straxen.URLConfig(
        default=6,
        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 averaged in integration windows", "area"), np.float32),
        (("Amplitude in LED window", "amplitude_led"), np.float32),
        (("Amplitude in off LED window", "amplitude_noise"), np.float32),
        (("Channel", "channel"), np.int16),
        (("Start time of the interval (ns since unix epoch)", "time"), np.int64),
        (("Time resolution in ns", "dt"), np.int16),
        (("Length of the interval in samples", "length"), np.int32),
        (("Whether there was a hit found in the record", "triggered"), bool),
        (("Sample index of the hit that defines the window position", "hit_position"), np.uint8),
        (("Window used for integration", "integration_window"), np.uint8, (2,)),
        (("Baseline from the record", "baseline"), np.float32),
    ]

    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.

        """

        self.is_led_on = is_the_led_on(
            self.run_doc, self.defualt_run_comments, self.noise_run_comments
        )

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

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

        led_windows, triggered = get_led_windows(
            records,
            self.minimum_led_position,
            self.fixed_position,
            self.is_led_on,
            self.led_hit_extension,
            self.led_cal_hit_min_height_over_noise,
            self.led_cal_record_length,
            self.area_averaging_length,
        )

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

        area = get_area(records, led_windows, self.area_averaging_length, self.area_averaging_step)
        temp["area"] = area["area"]

        temp["triggered"] = triggered
        temp["hit_position"] = led_windows[:, 0] - self.led_hit_extension[0]
        temp["integration_window"] = led_windows
        temp["baseline"] = records["baseline"]
        return temp


def is_the_led_on(run_doc, defualt_run_comments, noise_run_comments):
    """Utilizing the run database metadata to determine whether the run ID corresponds to LED on or
    LED off runs.

    The LED off, or noise runs, are identified by having 'Gain_calibration_step3' or
    'SPE_calibration_step0' in the comment.

    """
    # Check if run_doc is a list with a dictionary
    if isinstance(run_doc, list) and isinstance(run_doc[0], dict):
        # Extract the dictionary
        doc = run_doc[0]

        # Check if the required keys are present
        required_keys = {"user", "date", "comment"}
        if all(key in doc for key in required_keys):
            # Check if 'comment' contains any of the noise run comments
            comment = doc["comment"]

            # Check if the comment matches the expected pattern
            if not all(x in comment for x in defualt_run_comments):
                raise ValueError("The comment does not match the expected pattern.")

            if any(noise_comment in comment for noise_comment in noise_run_comments):
                return False
            else:
                return True
        else:
            raise ValueError("The dictionary does not contain the required keys.")
    else:
        raise ValueError("The input is not a list with a single dictionary.")


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

    Subtract the pulse float(data) from baseline.

    """

    record_length_padded = 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 with 0 padding", "data"), "f4", (record_length_padded,)),
    ]

    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"] == led_cal_record_length))[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"][:, :led_cal_record_length] = (
        -1.0 * (records["data"][:, :led_cal_record_length].transpose() - bl[:]).transpose()
    )
    records["baseline"] = bl
    records["baseline_rms"] = rms
    return records


def get_led_windows(
    records,
    minimum_led_position,
    fixed_position,
    is_led_on,
    led_hit_extension,
    hit_min_height_over_noise,
    record_length,
    area_averaging_length,
):
    """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 minimum_led_position: The minimum simple index of the LED hits. Hits before this sample
        are ignored.
    :param is_led_on: Fetch from the run database. It is used for discriminate between LED on and
        LED off runs.
    :param fixed_position: Fixed ADC sample upon which the integration window is defined. Used if no
        hits are identified
    :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
    :param record_length: The length of one led_calibration record
    :param area_averaging_length: The length (samples) of the window to do the averaging on.

    """
    if len(records) == 0:  # If input is empty, return empty arrays of correct shape
        return np.empty((0, 2), dtype=np.int64), np.empty(0, dtype=bool)

    hits = strax.find_hits(
        records,
        min_amplitude=0,  # Always use the height over noise threshold.
        min_height_over_noise=hit_min_height_over_noise,
    )

    maximum_led_position = record_length - area_averaging_length - led_hit_extension[1]
    hits = hits[hits["left"] >= minimum_led_position]
    # Check if the records are sorted properly by 'record_i' first and 'time' second and sort them
    # if they are not
    record_i = hits["record_i"]
    time = hits["time"]
    if not (
        np.all(
            (record_i[:-1] < record_i[1:])
            | ((record_i[:-1] == record_i[1:]) & (time[:-1] <= time[1:]))
        )
    ):
        hits.sort(order=["record_i", "time"])

    # If there are not hits in the chunk or if
    # the run is a nosie run, with LED off,
    # the integration window is defined beased on a
    # hard-coded ADC sample
    if (not is_led_on) or (len(hits) == 0):
        default_hit_position = fixed_position
    else:
        default_hit_position = sps.mode(hits["left"])[0]

        if isinstance(default_hit_position, np.ndarray):
            default_hit_position = default_hit_position[0]

        if default_hit_position > maximum_led_position:
            default_hit_position = maximum_led_position

    triggered = np.zeros(len(records), dtype=bool)

    default_windows = np.tile(default_hit_position + np.array(led_hit_extension), (len(records), 1))
    return _get_led_windows(
        hits, default_windows, led_hit_extension, maximum_led_position, triggered
    )


@numba.jit(nopython=True)
def _get_led_windows(hits, default_windows, led_hit_extension, maximum_led_position, triggered):
    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

        triggered[hit["record_i"]] = True

        hit_left = hit["left"]
        # Limit the position of the window so it stays inside the record.
        if hit_left > maximum_led_position:
            hit_left = maximum_led_position

        left = hit_left + led_hit_extension[0]
        right = hit_left + led_hit_extension[1]

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

    return windows, triggered


_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, area_averaging_length, area_averaging_step):
    """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.
    :param area_averaging_length: The total length in records of the window over which to do the
        averaging of the areas.
    :param area_averaging_step: The increase in length for each step of the averaging.
    :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 = np.arange(0, area_averaging_length + area_averaging_step, area_averaging_step)

    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	import scipy.stats as sps	2✔
16
17	# This makes sure shorthands for only the necessary functions
18	# are made available under straxen.[...]
19	export, __all__ = strax.exporter()	2✔
20
21	channel_list = [i for i in range(494)]	2✔
22
23
24	@export	2✔
25	class LEDCalibration(strax.Plugin):	2✔
26	"""
27	Preliminary version, several parameters to set during commissioning.
28	LEDCalibration returns: channel, time, dt, length, Area,
29	amplitudeLED and amplitudeNOISE.
30
31	The new variables are:
32	- Area: Area computed in the given window, averaged over 6
33	windows that have the same starting sample and different end
34	samples.
35	- amplitudeLED: peak amplitude of the LED on run in the given
36	window.
37	- amplitudeNOISE: amplitude of the LED on run in a window far
38	from the signal one.
39	"""
40
41	__version__ = "0.3.1"	2✔
42
43	depends_on = "raw_records"	2✔
44	data_kind = "led_cal"	2✔
45	compressor = "zstd"	2✔
46	parallel = "process"	2✔
47	rechunk_on_save = False	2✔
48
49	run_doc = straxen.URLConfig(	2✔
50	default="run_doc://comments?run_id=plugin.run_id",
51	infer_type=False,
52	help=(
53	"Comments dictionary from the run metadata. "
54	"It is used for discriminate between LED on and LED "
55	"off runs."
56	),
57	)
58
59	defualt_run_comments = straxen.URLConfig(	2✔
60	default=["auto, SC user: ", "pulserScript: "],
61	type=list,
62	help=("List of default comments the automatic script for runs in PMT calibration. "),
63	)
64
65	noise_run_comments = straxen.URLConfig(	2✔
66	default=["SPE_calibration_step0", "Gain_calibration_step3"],
67	type=list,
68	help=("List of comments for noise runs in PMT calibration. "),
69	)
70
71	led_cal_record_length = straxen.URLConfig(	2✔
72	default=160, infer_type=False, help="Length (samples) of one record without 0 padding."
73	)
74
75	baseline_window = straxen.URLConfig(	2✔
76	default=(0, 40),
77	infer_type=False,
78	help="Window (samples) for baseline calculation.",
79	)
80
81	minimum_led_position = straxen.URLConfig(	2✔
82	default=60,
83	infer_type=False,
84	help=(
85	"The minimum sample index to consider for LED hits. Hits before this sample are "
86	"ignored."
87	),
88	)
89
90	fixed_position = straxen.URLConfig(	2✔
91	default=88,
92	infer_type=False,
93	help=(
94	"Fixed ADC sample upon which the integration window is defined. "
95	"This is used as default when no hits or less than a certain amount "
96	"are identified."
97	),
98	)
99
100	led_hit_extension = straxen.URLConfig(	2✔
101	default=(-8, 32),
102	infer_type=False,
103	help="The extension around the LED hit to integrate.",
104	)
105
106	area_averaging_length = straxen.URLConfig(	2✔
107	default=7,
108	infer_type=False,
109	help=(
110	"The total length of the averaging window for the area calculation. "
111	"To mitigate a possiple bias from noise, the area is integrated multiple times with "
112	"sligntly different window lengths and then averaged. area_averaging_length should "
113	"be divisible by area_averaging_step."
114	),
115	)
116
117	area_averaging_step = straxen.URLConfig(	2✔
118	default=1,
119	infer_type=False,
120	help=(
121	"The step size used for the different windows, averaged for the area calculation. "
122	"To mitigate a possiple bias from noise, the area is integrated multiple times with "
123	"sligntly different window lengths and then averaged. area_averaging_length should "
124	"be divisible by area_averaging_step."
125	),
126	)
127
128	noise_window = straxen.URLConfig(	2✔
129	default=(10, 50), infer_type=False, help="Window (samples) to analyse the noise"
130	)
131
132	channel_list = straxen.URLConfig(	2✔
133	default=(tuple(channel_list)),
134	infer_type=False,
135	help="List of PMTs. Defalt value: all the PMTs",
136	)
137
138	led_cal_hit_min_height_over_noise = straxen.URLConfig(	2✔
139	default=6,
140	infer_type=False,
141	help=(
142	"Minimum hit amplitude in numbers of baseline_rms above baseline. "
143	"Actual threshold used is max(hit_min_amplitude, hit_min_"
144	"height_over_noise * baseline_rms)."
145	),
146	)
147
148	dtype = [	2✔
149	(("Area averaged in integration windows", "area"), np.float32),
150	(("Amplitude in LED window", "amplitude_led"), np.float32),
151	(("Amplitude in off LED window", "amplitude_noise"), np.float32),
152	(("Channel", "channel"), np.int16),
153	(("Start time of the interval (ns since unix epoch)", "time"), np.int64),
154	(("Time resolution in ns", "dt"), np.int16),
155	(("Length of the interval in samples", "length"), np.int32),
156	(("Whether there was a hit found in the record", "triggered"), bool),
157	(("Sample index of the hit that defines the window position", "hit_position"), np.uint8),
158	(("Window used for integration", "integration_window"), np.uint8, (2,)),
159	(("Baseline from the record", "baseline"), np.float32),
160	]
161
162	def compute(self, raw_records):	2✔
163	"""The data for LED calibration are build for those PMT which belongs to channel list.
164
165	This is used for the different ligh levels. As default value all the PMTs are considered.
166
167	"""
168
169	self.is_led_on = is_the_led_on(	2✔
170	self.run_doc, self.defualt_run_comments, self.noise_run_comments
171	)
172
173	mask = np.where(np.in1d(raw_records["channel"], self.channel_list))[0]	2✔
174	raw_records_active_channels = raw_records[mask]	2✔
175	records = get_records(	2✔
176	raw_records_active_channels, self.baseline_window, self.led_cal_record_length
177	)
178	del raw_records_active_channels, raw_records	2✔
179
180	temp = np.zeros(len(records), dtype=self.dtype)	2✔
181	strax.copy_to_buffer(records, temp, "_recs_to_temp_led")	2✔
182
183	led_windows, triggered = get_led_windows(	2✔
184	records,
185	self.minimum_led_position,
186	self.fixed_position,
187	self.is_led_on,
188	self.led_hit_extension,
189	self.led_cal_hit_min_height_over_noise,
190	self.led_cal_record_length,
191	self.area_averaging_length,
192	)
193
194	on, off = get_amplitude(records, led_windows, self.noise_window)	2✔
195	temp["amplitude_led"] = on["amplitude"]	2✔
196	temp["amplitude_noise"] = off["amplitude"]	2✔
197
198	area = get_area(records, led_windows, self.area_averaging_length, self.area_averaging_step)	2✔
199	temp["area"] = area["area"]	2✔
200
201	temp["triggered"] = triggered	2✔
202	temp["hit_position"] = led_windows[:, 0] - self.led_hit_extension[0]	2✔
203	temp["integration_window"] = led_windows	2✔
204	temp["baseline"] = records["baseline"]	2✔
205	return temp	2✔
206
207
208	def is_the_led_on(run_doc, defualt_run_comments, noise_run_comments):	2✔
209	"""Utilizing the run database metadata to determine whether the run ID corresponds to LED on or
210	LED off runs.
211
212	The LED off, or noise runs, are identified by having 'Gain_calibration_step3' or
213	'SPE_calibration_step0' in the comment.
214
215	"""
216	# Check if run_doc is a list with a dictionary
217	if isinstance(run_doc, list) and isinstance(run_doc[0], dict):	2✔
218	# Extract the dictionary
219	doc = run_doc[0]	2✔
220
221	# Check if the required keys are present
222	required_keys = {"user", "date", "comment"}	2✔
223	if all(key in doc for key in required_keys):	2✔
224	# Check if 'comment' contains any of the noise run comments
225	comment = doc["comment"]	2✔
226
227	# Check if the comment matches the expected pattern
228	if not all(x in comment for x in defualt_run_comments):	2✔
229	raise ValueError("The comment does not match the expected pattern.")
230
231	if any(noise_comment in comment for noise_comment in noise_run_comments):	2✔
232	return False	×
233	else:
234	return True	2✔
235	else:
236	raise ValueError("The dictionary does not contain the required keys.")
237	else:
238	raise ValueError("The input is not a list with a single dictionary.")
239
240
241	def get_records(raw_records, baseline_window, led_cal_record_length):	2✔
242	"""Determine baseline as the average of the first baseline_samples of each pulse.
243
244	Subtract the pulse float(data) from baseline.
245
246	"""
247
248	record_length_padded = np.shape(raw_records.dtype["data"])[0]	2✔
249
250	_dtype = [	2✔
251	(("Start time since unix epoch [ns]", "time"), "<i8"),
252	(("Length of the interval in samples", "length"), "<i4"),
253	(("Width of one sample [ns]", "dt"), "<i2"),
254	(("Channel/PMT number", "channel"), "<i2"),
255	(
256	(
257	"Length of pulse to which the record belongs (without zero-padding)",
258	"pulse_length",
259	),
260	"<i4",
261	),
262	(("Fragment number in the pulse", "record_i"), "<i2"),
263	(
264	("Baseline in ADC counts. data = int(baseline) - data_orig", "baseline"),
265	"f4",
266	),
267	(
268	("Baseline RMS in ADC counts. data = baseline - data_orig", "baseline_rms"),
269	"f4",
270	),
271	(("Waveform data in raw ADC counts with 0 padding", "data"), "f4", (record_length_padded,)),
272	]
273
274	records = np.zeros(len(raw_records), dtype=_dtype)	2✔
275	strax.copy_to_buffer(raw_records, records, "_rr_to_r_led")	2✔
276
277	mask = np.where((records["record_i"] == 0) & (records["length"] == led_cal_record_length))[0]	2✔
278	records = records[mask]	2✔
279	bl = records["data"][:, baseline_window[0] : baseline_window[1]].mean(axis=1)	2✔
280	rms = records["data"][:, baseline_window[0] : baseline_window[1]].std(axis=1)	2✔
281	records["data"][:, :led_cal_record_length] = (	2✔
282	-1.0 * (records["data"][:, :led_cal_record_length].transpose() - bl[:]).transpose()
283	)
284	records["baseline"] = bl	2✔
285	records["baseline_rms"] = rms	2✔
286	return records	2✔
287
288
289	def get_led_windows(	2✔
290	records,
291	minimum_led_position,
292	fixed_position,
293	is_led_on,
294	led_hit_extension,
295	hit_min_height_over_noise,
296	record_length,
297	area_averaging_length,
298	):
299	"""Search for hits in the records, if a hit is found, return an interval around the hit given by
300	led_hit_extension. If no hit is found in the record, return the default window.
301
302	:param records: Array of the records to search for LED hits.
303	:param minimum_led_position: The minimum simple index of the LED hits. Hits before this sample
304	are ignored.
305	:param is_led_on: Fetch from the run database. It is used for discriminate between LED on and
306	LED off runs.
307	:param fixed_position: Fixed ADC sample upon which the integration window is defined. Used if no
308	hits are identified
309	:param led_hit_extension: The integration window around the first hit found to use. A tuple of
310	form (samples_before, samples_after) the first LED hit.
311	:param hit_min_amplitude: Minimum amplitude of the signal to be considered a hit.
312	:param hit_min_height_over_noise: Minimum height of the signal over noise to be considered a
313	hit. :return (len(records), 2) array: Integration window for each record
314	:param record_length: The length of one led_calibration record
315	:param area_averaging_length: The length (samples) of the window to do the averaging on.
316
317	"""
318	if len(records) == 0: # If input is empty, return empty arrays of correct shape	2✔
319	return np.empty((0, 2), dtype=np.int64), np.empty(0, dtype=bool)	2✔
320
321	hits = strax.find_hits(	×
322	records,
323	min_amplitude=0, # Always use the height over noise threshold.
324	min_height_over_noise=hit_min_height_over_noise,
325	)
326
327	maximum_led_position = record_length - area_averaging_length - led_hit_extension[1]	×
328	hits = hits[hits["left"] >= minimum_led_position]	×
329	# Check if the records are sorted properly by 'record_i' first and 'time' second and sort them
330	# if they are not
331	record_i = hits["record_i"]	×
332	time = hits["time"]	×
333	if not (	×
334	np.all(
335	(record_i[:-1] < record_i[1:])
336	\| ((record_i[:-1] == record_i[1:]) & (time[:-1] <= time[1:]))
337	)
338	):
339	hits.sort(order=["record_i", "time"])	×
340
341	# If there are not hits in the chunk or if
342	# the run is a nosie run, with LED off,
343	# the integration window is defined beased on a
344	# hard-coded ADC sample
345	if (not is_led_on) or (len(hits) == 0):	×
346	default_hit_position = fixed_position	×
347	else:
348	default_hit_position = sps.mode(hits["left"])[0]	×
349
350	if isinstance(default_hit_position, np.ndarray):	×
351	default_hit_position = default_hit_position[0]	×
352
353	if default_hit_position > maximum_led_position:	×
354	default_hit_position = maximum_led_position	×
355
356	triggered = np.zeros(len(records), dtype=bool)	×
357
358	default_windows = np.tile(default_hit_position + np.array(led_hit_extension), (len(records), 1))	×
359	return _get_led_windows(	×
360	hits, default_windows, led_hit_extension, maximum_led_position, triggered
361	)
362
363
364	@numba.jit(nopython=True)	2✔
365	def _get_led_windows(hits, default_windows, led_hit_extension, maximum_led_position, triggered):	2✔
366	windows = default_windows	×
367	last = -1	×
368
369	for hit in hits:	×
370	if hit["record_i"] == last:	×
371	continue # If there are multiple hits in one record, ignore after the first	×
372
373	triggered[hit["record_i"]] = True	×
374
375	hit_left = hit["left"]	×
376	# Limit the position of the window so it stays inside the record.
377	if hit_left > maximum_led_position:	×
378	hit_left = maximum_led_position	×
379
380	left = hit_left + led_hit_extension[0]	×
381	right = hit_left + led_hit_extension[1]	×
382
383	windows[hit["record_i"]] = np.array([left, right])	×
384	last = hit["record_i"]	×
385
386	return windows, triggered	×
387
388
389	_on_off_dtype = np.dtype([("channel", "int16"), ("amplitude", "float32")])	2✔
390
391
392	@numba.jit(nopython=True)	2✔
393	def get_amplitude(records, led_windows, noise_window):	2✔
394	"""Needed for the SPE computation.
395
396	Get the maximum of the signal in two different regions, one where there is no signal, and one
397	where there is.
398
399	:param records: Array of records
400	:param ndarray led_windows : 2d array of shape (len(records), 2) with the window to use as the
401	signal on area for each record. Inclusive left boundary and exclusive right boundary.
402	:param tuple noise_window: Tuple with the window, used for the signal off area for all records.
403	:return ndarray ons: 1d array of length len(records). The maximum amplitude in the led window
404	area for each record.
405	:return ndarray offs: 1d array of length len(records). The maximum amplitude in the noise area
406	for each record.
407
408	"""
409	ons = np.zeros(len(records), dtype=_on_off_dtype)	2✔
410	offs = np.zeros(len(records), dtype=_on_off_dtype)	2✔
411
412	for i, record in enumerate(records):	2✔
413	ons[i]["channel"] = record["channel"]	×
414	offs[i]["channel"] = record["channel"]	×
415
416	ons[i]["amplitude"] = np.max(record["data"][led_windows[i, 0] : led_windows[i, 1]])	×
417	offs[i]["amplitude"] = np.max(record["data"][noise_window[0] : noise_window[1]])	×
418
419	return ons, offs	2✔
420
421
422	_area_dtype = np.dtype([("channel", "int16"), ("area", "float32")])	2✔
423
424
425	@numba.jit(nopython=True)	2✔
426	def get_area(records, led_windows, area_averaging_length, area_averaging_step):	2✔
427	"""Needed for the gain computation.
428
429	Integrate the record in the defined window area. To reduce the effects of the noise, this is
430	done with 6 different window lengths, which are then averaged.
431
432	:param records: Array of records
433	:param ndarray led_windows : 2d array of shape (len(records), 2) with the window to use as the
434	integration boundaries.
435	:param area_averaging_length: The total length in records of the window over which to do the
436	averaging of the areas.
437	:param area_averaging_step: The increase in length for each step of the averaging.
438	:return ndarray area: 1d array of length len(records) with the averaged integrated areas for
439	each record.
440
441	"""
442	area = np.zeros(len(records), dtype=_area_dtype)	2✔
443	end_pos = np.arange(0, area_averaging_length + area_averaging_step, area_averaging_step)	2✔
444
445	for i, record in enumerate(records):	2✔
446	area[i]["channel"] = record["channel"]	×
447	for right in end_pos:	×
448	area[i]["area"] += np.sum(record["data"][led_windows[i, 0] : led_windows[i, 1] + right])	×
449
450	area[i]["area"] /= float(len(end_pos))	×
451
452	return area	2✔
453
454
455	@export	2✔
456	class nVetoExtTimings(strax.Plugin):	2✔
457	"""Plugin which computes the time difference `delta_time` from pulse timing of `hitlets_nv` to
458	start time of `raw_records` which belong the `hitlets_nv`.
459
460	They are used as the external trigger timings.
461
462	"""
463
464	__version__ = "0.0.1"	2✔
465
466	depends_on = ("raw_records_nv", "hitlets_nv")	2✔
467	provides = "ext_timings_nv"	2✔
468	data_kind = "hitlets_nv"	2✔
469
470	compressor = "zstd"	2✔
471
472	channel_map = straxen.URLConfig(	2✔
473	track=False,
474	type=immutabledict,
475	help="immutabledict mapping subdetector to (min, max) channel number.",
476	)
477
478	def infer_dtype(self):	2✔
479	dtype = []	2✔
480	dtype += strax.time_dt_fields	2✔
481	dtype += [	2✔
482	(("Delta time from trigger timing [ns]", "delta_time"), np.int16),
483	(
484	("Index to which pulse (not record) the hitlet belongs to.", "pulse_i"),
485	np.int32,
486	),
487	]
488	return dtype	2✔
489
490	def setup(self):	2✔
491	self.nv_pmt_start = self.channel_map["nveto"][0]	2✔
492	self.nv_pmt_stop = self.channel_map["nveto"][1] + 1	2✔
493
494	def compute(self, hitlets_nv, raw_records_nv):	2✔
495	rr_nv = raw_records_nv[raw_records_nv["record_i"] == 0]	2✔
496	pulses = np.zeros(len(rr_nv), dtype=self.pulse_dtype())	2✔
497	pulses["time"] = rr_nv["time"]	2✔
498	pulses["endtime"] = rr_nv["time"] + rr_nv["pulse_length"] * rr_nv["dt"]	2✔
499	pulses["channel"] = rr_nv["channel"]	2✔
500
501	ext_timings_nv = np.zeros_like(hitlets_nv, dtype=self.dtype)	2✔
502	ext_timings_nv["time"] = hitlets_nv["time"]	2✔
503	ext_timings_nv["length"] = hitlets_nv["length"]	2✔
504	ext_timings_nv["dt"] = hitlets_nv["dt"]	2✔
505	self.calc_delta_time(	2✔
506	ext_timings_nv, pulses, hitlets_nv, self.nv_pmt_start, self.nv_pmt_stop
507	)
508
509	return ext_timings_nv	2✔
510
511	@staticmethod	2✔
512	def pulse_dtype():	2✔
513	pulse_dtype = []	2✔
514	pulse_dtype += strax.time_fields	2✔
515	pulse_dtype += [(("PMT channel", "channel"), np.int16)]	2✔
516	return pulse_dtype	2✔
517
518	@staticmethod	2✔
519	def calc_delta_time(ext_timings_nv_delta_time, pulses, hitlets_nv, nv_pmt_start, nv_pmt_stop):	2✔
520	"""Numpy access with fancy index returns copy, not view This for-loop is required to
521	substitute in one by one."""
522	hitlet_index = np.arange(len(hitlets_nv))	2✔
523	pulse_index = np.arange(len(pulses))	2✔
524	for ch in range(nv_pmt_start, nv_pmt_stop):	2✔
525	mask_hitlets_in_channel = hitlets_nv["channel"] == ch	2✔
526	hitlet_in_channel_index = hitlet_index[mask_hitlets_in_channel]	2✔
527
528	mask_pulse_in_channel = pulses["channel"] == ch	2✔
529	pulse_in_channel_index = pulse_index[mask_pulse_in_channel]	2✔
530
531	hitlets_in_channel = hitlets_nv[hitlet_in_channel_index]	2✔
532	pulses_in_channel = pulses[pulse_in_channel_index]	2✔
533	hit_in_pulse_index = strax.fully_contained_in(hitlets_in_channel, pulses_in_channel)	2✔
534	for h_i, p_i in zip(hitlet_in_channel_index, hit_in_pulse_index):	2✔
535	if p_i == -1:	2✔
536	continue	×
537	res = ext_timings_nv_delta_time[h_i]	2✔
538
539	res["delta_time"] = (	2✔
540	hitlets_nv[h_i]["time"]
541	+ hitlets_nv[h_i]["time_amplitude"]
542	- pulses_in_channel[p_i]["time"]
543	)
544	res["pulse_i"] = pulse_in_channel_index[p_i]	2✔

XENONnT / straxen / 11867531211

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