Precision Neutron Counting with Tail Pulse Pileup Rejection

From: B. Fornal, “Neutron Dark Decay,” Universe 9, 449 (2023). https://doi.org/10.3390/universe9100449

Overview

A free neutron undergoes beta decay into a proton, an electron, and an electron antineutrino, with a lifetime of about 880 s. The most precise measurements of the neutron lifetime come from two different methods, known as beam and bottle measurements.

• In the beam method, a beam of cold neutrons passes through a Penning trap, which captures and counts the protons from neutrons that decay inside it. The lifetime follows from the ratio of the number of neutrons in the trap volume, determined with a calibrated neutron flux monitor, to the rate of decay protons.

• In the bottle method, ultracold neutrons are stored in a trap (e.g. the magneto-gravitational trap of UCNτ), and the neutrons that survive after different storage times are counted.

The two methods give lifetimes that differ by about 9–10 s, a discrepancy of about 4σ (e.g. 887.7 ± 2.2 s from the beam method [3] versus 877.75 ± 0.28 (stat.) +0.22/−0.16 (syst.) s from UCNτ [4]). This neutron lifetime puzzle is what the next generation of beam and bottle experiments aims to resolve. The figure below illustrates the beam method studied in this project.

From BL3 NSF presentation

The bottle method has reached an uncertainty of a few tenths of a second, whereas the current beam measurement has an uncertainty of a couple of seconds. The next-generation beam experiment, BL3, aims to reduce the beam uncertainty to the sub-second level. A key ingredient is the absolute calibration of the neutron flux monitor. A new version of the “alpha–gamma” calibration system [2], Alpha-Gamma 2.0, is being developed at Illinois for this purpose. It counts the alpha particles and the 478 keV gamma rays produced when neutrons from the beam are captured in a boron (¹⁰B) or lithium (⁶Li) target, and uses these counts to calibrate the neutron counters.

The alphas and gammas are detected by silicon and high-purity germanium detectors, respectively. In these semiconductor detectors, the incident radiation creates electron–hole pairs, whose total charge is proportional to the deposited energy. This short charge pulse is integrated by a charge-sensitive preamplifier. Its output is a tail pulse, with a fast rise (nanoseconds) and a slow exponential decay (tens to hundreds of microseconds), whose height is proportional to the energy. A digitizer samples this signal, and a digital trapezoidal filter, often running in the digitizer’s firmware, converts each tail pulse into a trapezoid whose flat-top height measures the energy, with reduced noise.

We know which energies to expect from the alphas and gammas, so the task is to count precisely the number of events in each detector within a given energy window. One complication is pile-up: if two particles deposit energy within a short time of each other, their tail pulses, and hence their trapezoids, overlap, which makes the combined shape difficult to analyze. This project focuses on detecting and resolving the energies of piled-up pulses, using data from test radiation sources such as cobalt-60, and simulated data.

Part 1: Understanding BL3

Question 1 - Why do we care about the neutron lifetime?

From this brief article, or from other articles you find (good search terms are “neutron lifetime” and “neutron lifetime puzzle”; see also Ref. [6]), summarize in your own words what other physics depends on the free-neutron lifetime. Consider, for example, Big Bang nucleosynthesis and tests of the Standard Model through the CKM matrix element \(V_{ud}\). Which other fields or researchers might be interested in this quantity? Cite your sources.

Question 2 - Why is energy resolution important?

For both the alpha and the gamma detectors, we count only events within a narrow window of the energy spectrum, because the events we want are monoenergetic. Neutron capture on ¹⁰B produces, in 94% of captures, a 1.47 MeV alpha particle together with a 478 keV gamma ray from the excited ⁷Li nucleus. (The total kinetic energy released in this branch is 2.31 MeV, shared between the alpha and the ⁷Li recoil.)

• Assume the detectors are inside a stainless-steel vacuum chamber. Which detector will see more background events, and why? Refer to the Wikipedia article on radiation or another source.

• What physical feature of the background spectrum could leak into our narrow energy window? (Hint: this effect is common to all gamma-ray spectra of this type. Look up the features of a typical gamma-ray spectrum, such as the Compton continuum.) This is a precision measurement, so it is important to minimize the background and understand its origin.

• After reading the introduction to pile-up below, explain how pile-up can move the events we are interested in out of the energy window if it is not corrected. Unlike the background, this happens in both detectors.

Part 2: Understanding the Trapezoidal Filter

To extract a precise lifetime, BL3 needs to measure several particle rates with precision better than 0.01%. Real signals are distorted by noise and pile-up, which makes it hard to resolve pulses accurately. Instead of reading the energy directly from the peak height of the tail pulse, a trapezoidal filter is very useful for this purpose.

The trapezoidal filter is a widely used digital signal-processing technique in nuclear and particle physics, particularly for shaping and analyzing signals from charge-sensitive preamplifiers. Its main purpose is to extract energy information from pulses while minimizing noise and preserving timing resolution. For a tail pulse \(v(n)\) sampled at times \(n\) (in units of the sampling period), with an exponential decay time \(\tau\) (in the same units), the filter can be written recursively as [1]:

\[ \Large d^{k,l}(n) = v(n) - v(n-k) - v(n-l) + v(n-k-l) \]

\[ \Large s(n) = \sum_{i=0}^{n} \left[ \sum_{j=0}^{i} d^{k,l}(j) + M\, d^{k,l}(i) \right], \qquad M = \frac{1}{e^{1/\tau} - 1} \approx \tau \]

Here \(s(n)\) is the filter output, \(k\) is the rise (peaking) time, \(l - k\) is the length of the flat top (gap time), and \(M\) corrects for the exponential decay of the preamplifier pulse (“pole-zero correction”). With the correct \(M\), a single tail pulse is turned into a trapezoid with a linear rise of length \(k\), a flat top of length \(l - k\), and a linear fall. Its flat-top height is proportional to the pulse amplitude. Without the \(M\) term, the output does not have a flat top and does not return to the baseline.

The filter’s output resembles a trapezoid, characterized by a linear rise, a flat top, and a linear fall. Key parameters influencing its performance include the peaking time (or rise time), which determines how quickly the signal reaches the flat top and affects noise suppression and resolution; the gap time, which defines the width of the flat top and allows for ballistic deficit correction; and the decay time correction, or energy correction, which compensates for exponential decay in the preamplifier output, ensuring accurate amplitude measurement. Together, these parameters must be carefully tuned to balance energy resolution, noise rejection, and signal pile-up recovery based on the specific characteristics of the detector and electronics.

While the trapezoidal filter is effective at shaping pulses and improving energy resolution, it does not inherently resolve pulse pile-up, a common issue in high-rate environments where two or more events occur within a short time window. When pile-up occurs, the resulting waveform becomes a superposition of overlapping pulses, which the trapezoidal filter may interpret as a single distorted pulse, leading to inaccurate timing and energy estimation. To address this limitation, machine learning (ML) models, particularly recurrent neural networks (RNNs) and convolutional neural networks (CNNs), offer powerful alternatives. CNNs excel at recognizing local patterns in waveforms, such as distinguishing overlapping peaks, while RNNs can capture temporal dependencies to model the sequence of pulses over time. By training these models on labeled waveform data, they can learn to accurately identify the number of peaks, estimate each pulse’s timing, and recover the true energy of individual events—even in complex pile-up scenarios. This approach can significantly enhance signal processing performance beyond the capabilities of traditional trapezoidal filtering.

Shown below is the effect of pile-up where the left pulse is alone, but the right pulse is piled up with a secondary pulse spaced close together.

From: Yang et al., “Efficient pile-up correction based on pulse-tail prediction for high count rates,” Nucl. Instrum. Meth. A 1029, 166376 (2022). https://doi.org/10.1016/j.nima.2022.166376

Question 3 - Advantage of trapezoids

Briefly explain in your own words what a trapezoidal filter is. Implement the filter above in a few lines of NumPy, and apply it to a simulated tail pulse \(v(n) = A\,e^{-(n - n_0)/\tau}\) for \(n \ge n_0\) (e.g. \(\tau = 500\), \(k = 50\), \(l = 80\)), with and without the \(M\) term. Plot the results.

What is the advantage of trapezoidal filtering for resolving the energy spectrum of a detector, compared with using only the peak height of the tail pulse? (Hint: these filters are designed for real detectors that make precision measurements with analog real-world components. Which real-world effects complicate the measurement, e.g. electronic noise or variations in the rise time (“ballistic deficit”)?)

Question 4 - Parameter estimation

Shown below is a simulated pulse before convolution in blue and the resulting trapezoid produced from the convolution in the equation above.

From this image, estimate the energy (amplitude), the peaking time and the gap (flat-top) time of the trapezoidal filter. Briefly describe each parameter.

Part 3: A Pre-ML Approach to Data Analysis

One way to analyze the tail-pulse data without a trapezoidal filter and without machine learning is pulse fitting. We use the gradient of the tail pulse to “guess” the number of pulses, then fit the waveform with a sum of overlapping asymmetric pulses. Each pulse is described by a function with a rise time (tr), a decay time (td), an arrival time (x0), and an amplitude (A). td and tr are shared by all pulses within a fit. All times are in units of the sampling period (“samples”).

After fitting, we can histogram the amplitudes to obtain the energy distribution. However, there are some challenges:

• The choice of pulse model matters, because it affects the accuracy and resolution of the fitted amplitudes.

• Detecting the pulses is difficult, especially with pile-up, because the overlapping gradients hide the information about how many pulses are present.

Before dealing with those challenges, let’s implement this algorithm with the pulse model

\[ \Large y(x)=A\,\frac{e^{-(x-x_0)/t_d}}{1+e^{-(x-x_0)/t_r}} \]

Below we analyze a single pulse from real data acquired with the detector. We first subtract the baseline (the digitizer output in the absence of a pulse).

To estimate the number of pulses more reliably in the presence of noise, we apply a low-pass (Savitzky–Golay) filter to the waveform and store it in data_grad. We use this filtered waveform to compute the gradient, while the fit is performed on the original waveform after light Gaussian smoothing.

Plot the waveform and its gradient, and print the fitted energy (amplitude).

Note: for this part we use real waveforms recorded in the experiment (trial_9.npy: 177,436 waveforms of 1000 samples each, in ADC units).

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import glob
from scipy.optimize import curve_fit, least_squares
from scipy.signal import savgol_filter, find_peaks
from scipy.ndimage import gaussian_filter1d
from tqdm.notebook import tqdm

plt.style.use('seaborn-v0_8')

The data files are hosted on Google Drive and downloaded with gdown (pip install gdown if needed). The file is 338 MB.

# Import the data
import gdown

file_id = '1OF9F9_JQ4YygE3FcPWBtNBSLpl1XxuFg'
file_name = 'trial_9.npy'
gdown.download(f'https://drive.google.com/uc?id={file_id}', file_name, quiet=False)

data = np.load(file_name)
data = data[:10000] # We take first 10000 data for faster runtime
data_grad = savgol_filter(data, window_length=50, polyorder=3, axis=1)
data = gaussian_filter1d(data, sigma=1.8, axis=1)

def pulse(x, x0, td, tr, A): # single pulse
    y = np.zeros_like(x)
    y = A * np.exp(-(x-x0) / td) / (1 + np.exp(-(x-x0) / tr))
    return y

def detect_pulses_grad_fast(waveform, grad_threshold=1.7, window=10):
    grad = np.gradient(waveform)
    kernel = np.ones(window)
    conv = np.convolve(grad > grad_threshold, kernel, mode='same')
    peaks, _ = find_peaks(conv, height=window)
    return peaks

def multi_pulse_shared_tr_td(x, *params):
    td = params[0]
    tr = params[1]
    N = (len(params) - 2) // 2
    x0s = params[2:2 + N]
    As = params[2 + N:]
    y = np.zeros_like(x)
    for x0, A in zip(x0s, As):
        y += A * np.exp(-(x-x0) / td) / (1 + np.exp(-(x-x0) / tr))
    return y

def fit_waveform_multi(waveform, smooth_grad, index=None):
    x = np.arange(len(waveform), dtype=float)
    residual = waveform.copy()
    baseline = int(np.mean(residual[:200]))
    residual = waveform - baseline
    peak_indices = detect_pulses_grad_fast(smooth_grad)

    delta_x0 = 50
    filtered = [(idx, residual[idx] - residual[max(0, idx - delta_x0)]) for idx in peak_indices]
    filtered = [(idx, A0) for idx, A0 in filtered if A0 > 0]

    if len(filtered) == 0:
        return [], [], residual, -100000

    peak_indices, A0s = zip(*filtered)
    peak_indices = list(peak_indices)
    A0s = list(A0s)
    p0 = [11508, 15] + list(peak_indices) + list(A0s)
    N = len(peak_indices)

    x0_lowers = [max(0, idx - delta_x0) for idx in peak_indices]
    x0_uppers = [min(len(waveform), idx + delta_x0) for idx in peak_indices]
    bounds_lower = [1e3, 1] + x0_lowers + [0] * N
    bounds_upper = [1e5, 200] + x0_uppers + [5e3] * N

    try:
        popt, _ = curve_fit(multi_pulse_shared_tr_td, x, residual,
                            p0=p0, maxfev=5000, bounds=(bounds_lower, bounds_upper))

    except Exception as e:
        print("Fitting failed:", e)
        return [], [], residual, -100000

    fitted_waveform = multi_pulse_shared_tr_td(x, *popt)
    residual = residual - fitted_waveform

    td, tr = popt[0], popt[1]
    x0s, As = popt[2:2+N], popt[2+N:]
    pulses = [(x0, td, tr, A) for x0, A in zip(x0s, As)]

    return pulses, popt, residual, baseline

index = 0

pulses, popt, residual, baseline = fit_waveform_multi(data[index], data_grad[index])
print(pulses)
print(popt)
print(residual)
print(baseline)

# Plot out the waveform, its gradient, and print the energy

Now, the code below can find the energy distribution from looping over all the waveforms. Because the gradient depends on the energy, some of the events on the lower energy side of the spectrum would fail to fit. Ignore those.

# Gamma
def find_energy():
    all_pulses = []
    all_tds = []
    all_trs = []
    all_energies = []

    for i, (raw, grad) in tqdm(enumerate(zip(data, data_grad)), total=len(data)):
        pulses, popt, residual, baseline = fit_waveform_multi(raw, grad, i)
        if baseline == -100000:
            continue
        all_pulses.append(pulses)
        if popt is not None and len(popt) >= 2:
            all_tds.append(popt[0])
            all_trs.append(popt[1])
        for p in pulses:
            all_energies.append(p[-1]) # Assuming p[-1] is the amplitude

    # Histogram 1: Full energy range
    plt.figure(figsize=(6, 4))
    energy_bins = np.arange(min(all_energies), max(all_energies))
    energy_hist, energy_edges = np.histogram(all_energies, bins=energy_bins)
    plt.hist(all_energies, bins=energy_bins, edgecolor='k')
    plt.xlabel('Pulse Amplitude (Energy)')
    plt.ylabel('Counts')
    plt.title('Energy Spectrum')
    plt.yscale('log')
    plt.tight_layout()
    plt.show()

    # Histogram 2: Zoomed to mean (-1.5 std, +5 std)
    plt.figure(figsize=(6, 4))
    zoom_bins = np.arange(int(np.mean(all_energies) - 1.5*np.std(all_energies)),
                          int(np.mean(all_energies) + 5*np.std(all_energies)))
    zoom_hist, zoom_edges = np.histogram(all_energies, bins=zoom_bins)
    plt.hist(all_energies, bins=zoom_bins, edgecolor='k')
    plt.xlabel('Pulse Amplitude (Energy)')
    plt.ylabel('Counts')
    plt.title('Zoomed Energy Spectrum')
    plt.yscale('log')
    plt.tight_layout()
    plt.show()

    # Histogram: td
    plt.figure(figsize=(6, 4))
    td_bins = np.linspace(min(all_tds), max(all_tds), 50)
    td_hist, td_edges = np.histogram(all_tds, bins=td_bins)
    plt.hist(all_tds, bins=td_bins, edgecolor='k')
    plt.xlabel('Decay Time (td)')
    plt.ylabel('Counts')
    plt.title('Distribution of td')
    plt.yscale('log')
    plt.tight_layout()
    plt.show()

    # Histogram: tr
    plt.figure(figsize=(6, 4))
    tr_bins = np.linspace(min(all_trs), max(all_trs), 50)
    tr_hist, tr_edges = np.histogram(all_trs, bins=tr_bins)
    plt.hist(all_trs, bins=tr_bins, edgecolor='k')
    plt.xlabel('Rise Time (tr)')
    plt.ylabel('Counts')
    plt.title('Distribution of tr')
    plt.yscale('log')
    plt.tight_layout()
    plt.show()

find_energy()

Now change the fitted pulse shape to

\[ \Large y(x)=A\left(1-e^{-(x-x_0)/t_r}\right)e^{-(x-x_0)/t_d}\,\Theta(x-x_0) \]

where \(\Theta\) is the Heaviside step function, and find the energy distribution again. Which model describes the data better? Compare the fit residuals and the widths of the peaks in the energy spectra.

def pulse(x, x0, td, tr, A):
    y = np.zeros_like(x)

    # Implement this
    raise NotImplementedError()

def multi_pulse_shared_tr_td(x, *params):
    td = params[0]
    tr = params[1]
    N = (len(params) - 2) // 2
    x0s = params[2:2 + N]
    As = params[2 + N:]
    y = np.zeros_like(x)
    for x0, A in zip(x0s, As):

        # Implement this
        raise NotImplementedError()

index = 0

pulses, popt, residual, baseline = fit_waveform_multi(data[index], data_grad[index])
print(pulses)
print(popt)
print(residual)
print(baseline)

# Plot out the waveform, its gradient, and print the energy

find_energy()

Part 4: Understanding the Dataset with ML

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier, GradientBoostingClassifier, GradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, accuracy_score
import pandas as pd
import ast
from collections import Counter
from matplotlib.colors import LogNorm
from sklearn.multioutput import MultiOutputRegressor

plt.style.use('seaborn-v0_8')

For the ML part, we split the problem into two tasks: finding the number of pulses in a waveform, and finding their energies.

For this part we use simulated waveforms (data_set2.parquet, 233 MB, 99,996 events with 1–6 pulses each, 16,666 events per class). Each row contains:

• raw: the simulated tail-pulse signal (1000 samples; decay time \(\tau\) = 500 samples);

• waveform: the same signal after trapezoidal filtering;

• num_pulses, times, heights: the true number of pulses, their arrival times (sample index), and their heights (arbitrary units, between 10 and 25). These are the labels for supervised learning.

Note that the simulated waveforms contain no noise, unlike the real data of Part 3. Keep this in mind when you interpret your results.

Question 5 - An ML approach

For this step, we train a model to find the number of pulses in a waveform. First, download and read the data:

# Import the data
import gdown

file_id = '1_0vTiUCjNfjSs9GaBBn0pQFWk0DEkpn6'
file_name = 'data_set2.parquet'
gdown.download(f'https://drive.google.com/uc?id={file_id}', file_name, quiet=False)

data = pd.read_parquet(file_name)
print(data.head())
for i in range(1, 10):
    print(f"Events with {i} pulse(s): ", len(data[data['num_pulses'] == i]))

The dataframe is a bit tricky to work with, so here is code that splits it into a training set and a test set (80%/20% of each class):

train_data = pd.DataFrame()
test_data = pd.DataFrame()
for i in range(1, 7):
    train_data = pd.concat([
        train_data,
        data[data['num_pulses'] == i].iloc[0:int(len(data[data['num_pulses'] == i])*0.8)]
    ])
    test_data = pd.concat([
        test_data,
        data[data['num_pulses'] == i].iloc[int(len(data[data['num_pulses'] == i])*0.8):]
    ])

X_train = np.array([row["waveform"] for _, row in train_data.iterrows()])
y_train = np.array(train_data["num_pulses"])

y_train_cls = y_train.astype(int)

X_test = np.array([row["waveform"] for _, row in test_data.iterrows()])
y_test = np.array(test_data["num_pulses"])
y_test_cls = y_test.astype(int)

Train a classifier on the raw waveform arrays of these sets (e.g. a random forest, or a 1D CNN), and plot the confusion matrix.

A classifier applied directly to the 1000 samples has to discover for itself which aspects of the waveform matter. We can help it by extracting a few physically motivated features that summarize the waveform:

def extract_features(waveform):
    slope = np.diff(waveform)

    pos_slope = slope[slope > 0]
    neg_slope = slope[slope < 0]

    features = {
        "mean_val": np.mean(waveform),
        "std_val": np.std(waveform),
        "max_val": np.max(waveform),

        "mean_pos_slope": np.mean(pos_slope) if len(pos_slope) > 0 else 0,
        "std_pos_slope": np.std(pos_slope) if len(pos_slope) > 0 else 0,
        "max_pos_slope": np.max(pos_slope) if len(pos_slope) > 0 else 0,

        "mean_neg_slope": np.mean(neg_slope) if len(neg_slope) > 0 else 0,
        "std_neg_slope": np.std(neg_slope) if len(neg_slope) > 0 else 0,
        "max_neg_slope": np.min(neg_slope) if len(neg_slope) > 0 else 0,  # most negative
    }

    return features

features_list = []
labels = []

for waveform, label in zip(X_train, y_train_cls):
    feats = extract_features(waveform)
    feats["num_pulses"] = label
    features_list.append(feats)

df_features = pd.DataFrame(features_list)

Make a pairplot of the features. Which features do you think will be the most useful?

import seaborn as sns
import matplotlib.pyplot as plt

sns.set(style='whitegrid', context='notebook')
sns.pairplot(df_features, hue='num_pulses', corner=True, plot_kws={'alpha': 0.5})

plt.suptitle("Waveform Feature Clustering by Pulse Count", y=1.02)
plt.show()

Now train a classifier on these features (df_features, with the column num_pulses as the label) and call it clf. Then run the evaluation cell that follows. How does it compare with the classifier trained on the raw waveform?

# YOUR CODE HERE: define and train clf on the training features

from sklearn.metrics import classification_report, confusion_matrix, ConfusionMatrixDisplay, f1_score, roc_auc_score, roc_curve, auc
from sklearn.preprocessing import label_binarize

features_list = []
labels = []

for waveform, label in zip(X_test, y_test_cls):
    feats = extract_features(waveform)
    feats["num_pulses"] = label
    features_list.append(feats)

df_test_features = pd.DataFrame(features_list)
X_test_feats = df_test_features.drop(columns='num_pulses')
y_test_feats = df_test_features['num_pulses']

y_pred = clf.predict(X_test_feats)

print(classification_report(y_test_feats, y_pred))
ConfusionMatrixDisplay(confusion_matrix(y_test, y_pred)).plot()
plt.title("confusion Matrix")
plt.show()

f1_scores = f1_score(y_test_feats, y_pred, average=None)
classes = np.unique(y_test_feats)

plt.figure()
plt.bar(classes, f1_scores)
plt.xlabel('Class (num_pulses)')
plt.ylabel('F1 Score')
plt.title('F1 Score per Class')
plt.xticks(classes)
plt.grid(True)
plt.show()

y_score = clf.predict_proba(X_test_feats)
y_test_bin = label_binarize(y_test_feats, classes=classes)

plt.figure()
for i, cls in enumerate(classes):
    fpr, tpr, _ = roc_curve(y_test_bin[:, i], y_score[:, i])
    roc_auc = auc(fpr, tpr)
    plt.plot(fpr, tpr, label=f'Class {cls} (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], 'k--')
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title("ROC AUC per Class")
plt.legend(loc="lower right")
plt.grid(True)
plt.show()

This works much better, but it does not give us the energies (heights) of the pulses. Let’s move on to the last part.

Things to think about: the simulated waveforms are noise-free. Add Gaussian noise at the level of the real data (estimate it from the baseline region of the trial_9.npy waveforms, after scaling to the same units) and test how your classifier’s performance changes.

Question 6 - Advanced approach

Train a model that predicts the energies (heights) of the pulses in a piled-up waveform, and use it to plot an energy spectrum. Compare it with the true spectrum of the heights in the test set, and with a naive estimate, e.g. the maximum of each waveform.

The number of pulses varies from event to event, so think about how to represent the output. Some options:

• predict a fixed-length vector of 6 heights, sorted by arrival time and padded with zeros, possibly using your Question 5 classifier to decide how many to keep;

• predict, for every time sample, the probability that a pulse starts there and its height (a 1D CNN or U-Net “segmentation” approach). Then read off the pulses as peaks of the probability.

Evaluate the energy resolution (e.g. the distribution of predicted minus true height) as a function of the number of pulses. Optional: apply your model to the real waveforms of Part 3, and compare the energy spectrum with the one from pulse fitting.

References

[1] V. T. Jordanov, G. F. Knoll, “Digital synthesis of pulse shapes in real time for high resolution radiation spectroscopy,” Nucl. Instrum. Meth. A 345, 337 (1994). https://doi.org/10.1016/0168-9002(94)91011-1

[2] A. T. Yue et al., “Precision determination of absolute neutron flux,” Metrologia 55, 460 (2018). https://doi.org/10.1088/1681-7575/aac283

[3] A. T. Yue et al., “Improved Determination of the Neutron Lifetime,” Phys. Rev. Lett. 111, 222501 (2013). https://doi.org/10.1103/PhysRevLett.111.222501

[4] F. M. Gonzalez et al. (UCNτ Collaboration), “Improved Neutron Lifetime Measurement with UCNτ,” Phys. Rev. Lett. 127, 162501 (2021). https://doi.org/10.1103/PhysRevLett.127.162501

[5] H. Yang et al., “Efficient pile-up correction based on pulse-tail prediction for high count rates,” Nucl. Instrum. Meth. A 1029, 166376 (2022). https://doi.org/10.1016/j.nima.2022.166376

[6] B. Fornal, “Neutron Dark Decay,” Universe 9, 449 (2023). https://doi.org/10.3390/universe9100449

Acknowledgments

• Initial version: Calvin Nettelhorst and Ryan Huang with some guidance from Mark Neubauer

• Revised: 2026

© Copyright 2026