Function-fitting method

Saturation-induced signal truncation caused by the dynamic range limitation of the DAQ system leads to an irreversible loss of complete waveform morphology and charge integral data. Because the charge integral is a critical parameter directly proportional to the energy deposition and the Dis factor of detection events, the precise reconstruction of saturated signals is essential for comprehensive waveform recovery.

The hypothesis is that signal waveforms with different amplitudes exhibit analogous shape characteristics. To test this hypothesis, an artificial neural network (ANN) methodology was employed [40]. This study utilized a generalized radial basis function (GRBF) neural network, a topology validated for robust signal reconstruction performance [40, 41]. Four amplitude ranges were designated as training cohorts: 50–100, 100–150, 150–200, and 200–saturation (244.5 ADC units), with training sets containing 10,000, 10,000, 10,000, and 88 events, respectively. The scarcity of signals in the 200–244.5 range arises from the limited number of high-energy gamma-ray events emitted after neutron capture by Gd nuclei. All training samples were randomly selected from the detection datasets to ensure statistical representativeness, with a reduced population in the final cohort reflecting the natural amplitude distribution sparsity. The top panel of Fig. 2 displays the ANN outputs as characteristic scatter distributions: black inverted triangles (50–100), pink squares (100–150), purple circles (150–200), and blue pentagrams (200–saturation). These clustered distributions validate the waveform similarity hypothesis across the amplitude ranges.

Fig. 2Full-screen View

Comparison of event waveform morphology between ANN and function-fitting outputs. The upper panel illustrates the waveform shapes of four events with different amplitudes reconstructed by ANN and function fitting, with the fitting function applied over the 0–60 ns range. The lower panel shows the reconstruction performance of the function-fitting method for the saturated events. The X-axis represents the signal generation time within a 100-ns time window. For Dis factor calculation, relative timing (e.g., 15 ns post-peak) is typically adopted

Parallel reconstruction efforts apply a function-fitting method to the saturated signals. Red curves in Fig. 2 (top) represents the fits to all ANN output points using Eq. (1), where p₁, p₂, and p₃ denote the peak amplitude, peak time, and half-height width, respectively. Table 1 compares the ANN and function-fitting reconstruction results, showing an integration average deviation within a few thousandths across all amplitude ranges. The bottom panel of Fig. 2 demonstrates the reconstruction of a saturated signal: black scatters represent raw data, the red curve is the fitted function, and 244.5 ADC units mark the DAQ dynamic range limit. Both panels in Fig. 2 confirm the ability of the function to approximate the waveform morphology. However, the function aligns closely with low-amplitude waveforms but deviates slightly at higher amplitudes than the others.(1)

Analysis of the bias value of the function fitting

To quantify the systematic deviations introduced by the fitting function during the waveform reconstruction, we defined two distinct bias categories for systematic error characterization. This analysis specifically focuses on the reconstructed charge integral bias of physical events, which is calculated by statistical averaging over multiple events processed using the fitting procedure. Parameter estimation errors from the fitting process were not analyzed independently, as these errors were inherently embedded within the statistical ensemble treatment and thus contributed intrinsically to the observed charge integral bias.

The first type of bias arises from the imperfect alignment between Eq. (1) and the actual waveform, leading to systematic deviations when fitting the high-amplitude waveforms. The analysis framework employs unsaturated signals as baseline references across four amplitude ranges: 50–100, 100–150, 150–200, and 200–saturation (244.5 ADC units). The quantitative evaluation of the charge integral bias magnitudes within these ranges enables the systematic estimation of the uncertainty for saturated signal reconstruction. Equation (2) is further defined to characterize the first bias type:(2)where Q_f is the integral under the model of the fitting function, and Q_r is the area integral of the real waveform. Bias denotes the degree of bias of the modeled integral of the fitted function from the actual integral of the real waveform.

Figure 3(a) displays the amplitude-dependent bias distributions via color-coded histograms for the four amplitude ranges. The lower inset shows the expectation values and standard deviations derived from the Gaussian fits. As shown in Fig. 3(a), waveform integral biases using Eq. (1) remained below 0.006 (0.6%), with decreasing trends as the amplitudes increased. For the highest amplitude range (200–saturation), the bias measures 0.0023 and follows a decreasing pattern. This confirms the first-type bias for the saturated signal reconstruction using Eq. (1) is bounded at 0.0023. Systematic analysis validates that first-type biases remain within acceptable thresholds for practical applications, justifying their exclusion from subsequent analyses.

Fig. 3Full-screen View

Bias results of two types under function fitting. The left inset (a) illustrates the first type of bias, with the histogram below displaying the statistical information of the mean values. The right inset (b) illustrates the second type of bias, with the histogram below showing statistical information of the mean bias

The second type of bias stems from reconstruction errors in the saturated signals caused by waveform truncation owing to information loss. Larger saturated signal amplitudes result in more severe truncation, exacerbating information loss, and thereby increasing the reconstruction bias.

Prior to analyzing this bias, the Cut parameter was defined to quantify the truncation extent using Eq. (3).(3)where 244.5 ADC units denotes the DAQ dynamic range limit, and A represents the fitted peak amplitude. A higher Cut value indicates greater waveform truncation, whereas a lower Cut corresponds to minimal truncation.

To quantify the integral bias between the actual and reconstructed waveforms, an Artificial Cut (AC) was applied to waveforms within the 100–150 amplitude range. Note that the Cut defined in Eq. (3) represents the percentage of the actual saturated signal truncation level relative to the full waveform amplitude, whereas AC denotes the same percentage-based Cut applied to unsaturated waveforms. The AC mimics truncation-induced information loss, which is analogous to the natural Cut. This amplitude range provides sufficient event statistics (10,000 events per Cut value from 10% to 80%) to reliably characterize the potential truncation levels in the saturated signals. Integral bias is calculated similarly to the first bias type (Eq. (2)), with the results shown in Fig. 3(b). Figure 3(b) displays the integral bias distributions for varying Cut values via color-coded histograms, with red dashed lines denoting Gaussian fits. The lower inset presents the Gaussian-fitted mean values. The results demonstrate that higher Cut values yield larger biases and broader distributions, validating the hypothesis that reconstructing heavily truncated waveforms introduces greater deviations. For 80% Cut, the bias reaches 10.7%, whereas 10% Cut results in 0.2% bias. This confirms that the reconstruction function significantly affects the impacts waveform fidelity, necessitating bias correction in both energy and Dis factor calculations for saturated signals.

Results of reconstruction

Following the evaluation of the reconstruction biases, saturated signal events were reconstructed. The correction factor (CF) is defined as CF=1-Bias, where Bias quantifies the systematic deviations described in the previous section. The reconstructed charge integral was calculated using Eq. (4), which incorporates CF to correct for the bias effects.(4)The Gd-LS detector employs a dual PMT readout system, with PMTs positioned at the front and back ends to collect the light signals. Consequently, each event generated two distinct waveforms. The total charge integral Q_total is defined as the mean of the charge integrals from both PMTs (Q₁ and Q₂), as expressed in Eq. (5):(5)The Dis factor is computed using Eq. (6), where Q_1part and Q_2part represent the charge integrals within a 15 ns window following the signal peak (see Fig. 2, bottom). Typically, the Dis factor for gamma events is below 0.15, while that for neutron recoil events exceeds this threshold [16]. Error propagation methods were applied to estimate the uncertainties in Q_total and the Dis factor, accounting for statistical variances in Q_1part, Q_2part, and temporal measurements.(6)Figure 4 presents the reconstruction results for the saturated events. The left panel shows a two-dimensional scatter plot of the Q-value versus Dis for the reconstructed events, with all Dis values strictly below 0.15, aligning with the waveform morphology of gamma events. Because the energy of gamma events in the environment is typically below 10 MeV, corresponding to a Q-value of less than 4000, while the reconstructed events exhibit Q-values concentrated around 10,000, this indicates that these saturated events likely originate from cosmic-ray muons traversing the detector. The right panel presents a similar Q-Dis scatter plot, including low-energy background signals, where the red curve marks the gamma-like event trend. This distribution confirmed the classification of the reconstructed saturated events as gamma-like signals induced by cosmic-ray muon interactions. Consequently, the 26 reconstruction events were designated as muon event candidates.

Fig. 4Full-screen View

(Color online) Reconstructed results of saturated signals using the function-fitting method. Left panel: Two-dimensional distribution of charge-integrated Q-values and Dis for 26 saturated events. Right panel: Two-dimensional scatter plot of the Q value and Dis events containing low-energy background signals. The red curve was obtained by performing a Gaussian fit on the Dis distribution within different Q-value ranges. The value of the color ruler represents the density of points within a grid range

Event selection

Based on the detector geometry, cosmic-ray muons deposit approximately 50 MeV of energy over a mean trajectory length of 30 cm within the detector. To suppress background contamination from environmental gamma events and peripheral muon interactions, muon candidates were selected with energy thresholds exceeding 10 MeV. Under this selection criterion, PMT flasher events and electronic noise are the dominant background sources in the high-energy region. PMT flasher events originate from spontaneous light emission caused by base discharge within PMTs, whereas electronic noise stems from baseline fluctuations in the readout electronics, which are induced by external interference. Both background types exhibited anomalously high waveform amplitudes localized to specific PMT channels. To discriminate these backgrounds from genuine muon signals, a parameter r related to the waveform amplitude is defined, which equals the ratio of the difference between the right and left waveform amplitudes to their sum (), where L and R denote the signal amplitudes from the left and right PMTs, respectively.

We compared the reconstruction events with a large sample of pre-identified low-energy gamma events. These gamma events originate from neutron capture by Gd nuclei, which emit characteristic gamma rays. Characterized by low energy (typically below 5 MeV) and selected using time correlation criteria (described in Sect. 1), these events are uncontaminated by PMT flasher artifacts or electronic noise. Figure 5 compares the amplitude asymmetry parameter r between the low-energy gamma events and the reconstructed events. The left panel shows the distribution of r, revealing a higher density of reconstruction events than gamma events in the 0.2–0.4 range, suggesting the presence of background events. To reject background events, the right panel presents a two-dimensional scatter plot of the Q-value versus r, where the Q-values of the reconstruction events are scaled by a factor of 0.1 for enhanced visual discrimination, with the gamma event distribution defining the candidate region. Within the 0.2–0.4 range, five low-energy events were located outside the candidate region, leading to their classification as background. Furthermore, all reconstruction events exhibited large waveform amplitudes (with the smallest being approximately 244.5 ADC units), and the probability of simultaneous flasher event occurrence in the dual-PMT system was extremely low. Therefore, we attributed these background events to electronic noise. After stringent selection, 21 events remained as validated muon candidates.

Fig. 5Full-screen View

(Color online) Additional event distribution predictions and post-cut two-dimensional scatter plots of r versus Q-values. Left panel: Comparative distributions of r between gamma events and reconstructed saturated events; normalization was performed using the leftmost value of the distribution as the baseline. Right panel: Two-dimensional distribution of r and Q-values for gamma events versus the reconstructed saturated events. The value of the color ruler represents the density of points within the grid range. To facilitate a direct comparison, the Q-values of the reconstructed events were scaled by a factor of 0.1. This analysis enables robust removal of electronic noises

Simulation setting

In the particle source configuration, the GeneralParticleSource in Geant4 was used to generate cosmic-ray muons, with the cosmic-ray muon energy and incidence angle inputs derived from measurements conducted at CJPL by Guo et al. [14]. The total muon flux was 3.53 × 10^-10 cm^-2 s^-1 with an average energy of 340 GeV. The distributions of muon energy, azimuthal and zenithal angles are illustrated schematically in Fig. 6(a), (b), and (c). In addition, the impact of the mountainous structure above the laboratory on muons was considered in recent measurements. Based on this consideration, more muon events originated from the southwest and southeast directions. The location of the particle source is shown in Fig. 6(d). Muons were generated in a circular plane 5 m directly above the detector, with a radius of 5 m.

Fig. 6Full-screen View

(Color online) Particle source and detector simulation modeling used in Geant4 simulation. (a) Muon energy spectrum characterization. (b) Azimuthal distribution of incident muons. (c) Zenith angle distribution of incident muons. (d): Geometric configuration of detector modeling in Geant4

The detector model in the detector geometry configuration is illustrated in Fig. 6(d). A PE room with dimensions 8 m (length) × 5.6 m (width) × 4.1 m (height) was simulated using 1 m thick PE sheets (green squares). The detector was positioned in one corner of the PE room and was fully enclosed by PE sheets. A 0.6 m × 0.6 m × 0.05 m lead plate (pink square) was installed at the front of the detector, whereas 0.005 m thick aluminum plates (yellow squares) shielded the rear and sides. The detector body comprised a cylindrical structure with a diameter of 0.3 m and a length of 0.4 m (brown cylinder) surrounded by a PTFE reflective layer (gray cylinder). In the experimental setup, glass light guides and PMTs were mounted at both ends of the detector body.

The active area of detector

The detector’s active area is estimated by(7)where S represents the active area of the detector, T denotes the operational live time (412 days in this experiment, or 3.56 × 10⁷ s), and corresponds to the simulated muon flux (3.53 × 10^-10 cm^-2 s^-1). The parameter Nμ is defined as the mean event count with energy deposition exceeding 10 MeV, derived from the simulation.

An exposure unit is defined as 412 days of continuous detector operation, which is equivalent to 9,685 muons traversing a circular region with a 5-meter radius above the detector. Through 3,000 simulated exposure units, the event count distribution shown in Fig. 7. The simulations yielded a mean of 19.2 events per exposure unit. From this result, the mean event count was calculated as Nμ = 19.2, and the calculated active area S = 1528 cm².

Fig. 7Full-screen View

Distribution of the number of events with energy deposition exceeding 10 MeV in detectors

The energy deposition spectrum of the detector and the source of muon-induced neutrons

To estimate the composition of the reconstructed energy spectrum, correct the detector energy response, and investigate potential neutron background contributions, simulations were conducted with muon events equivalent to 1,000 times the experimental exposure. The simulation results are presented in Fig. 8. The upper panel of Fig. 8 displays the electron-equivalent energy (MeVee) spectrum below 100 MeVee in the detector, normalized to events per day per MeVee. This spectrum demonstrates that muon events dominate the visible events above 10 MeVee. In the 10–20 MeVee range, a small fraction of muon-induced electron events contributed to the spectrum. Additionally, recoil protons produced by the elastic scattering between muon-induced neutrons and hydrogen nuclei generate visible energy after the quenching effects. As shown by the pink band in the upper panel, muon-induced neutron events are approximately two orders of magnitude less frequent than muon events, with most exhibiting energies below 5 MeVee. The lower panel of Fig. 8 illustrates the spatial distribution of muon-induced neutron production within the PE room under these simulation conditions. The results revealed that while most neutrons originated from interactions at the PE room walls, detectable neutron-induced events primarily occured inside the detector. This spatial discrepancy arises because the 1-meter-thick polyethylene shielding effectively attenuates external neutrons, whereas internal neutron events result from muons directly penetrating the detector.

Fig. 8Full-screen View

The electron-equivalent energy spectrum (upper panel) and distribution of muon-induced neutron production (lower panel) under simulated conditions with 1,000 times the experimental muon exposure. The upper panel shows the total visible energy events (black solid line), muon events (red solid line), muon-induced electron events (blue solid line), and recoil proton events after the quenching effect (pink band). In the lower panel, the X and Y coordinates represent the horizontal and vertical lengths of the PE room, where black points mark all muon-induced neutron production positions and red points denote visible neutron events

Muon-induced neutrons detection efficiency

The neutron detection efficiency is defined as the ratio of neutrons undergoing reactions in the LS to the total number of neutrons entering the detector, expressed as:(8)where N₁ represents the number of neutrons that undergo elastic scattering with hydrogen nuclei in the LS, and N_Total is the total number of neutrons entering the detector. These values were derived from simulations, and 10,000 muon-induced neutrons were generated (sampled according to the neutron energy spectrum), and 5312 events satisfied the selection criteria. Thus, the efficiency was calculated as ϵ = 53.1%, with a 0.7% uncertainty originating from the MC simulation.

Additionally, the geometric efficiency ϵ_g is defined as the ratio of neutrons entering from outside the sensitive volume of the detector to the total number of neutrons. This definition accounts for the fact that only external neutrons can be effectively detected, as the internal neutrons produced by muons traversing the detector are obscured by higher-amplitude muon signals and cannot be distinguished. The lower inset of Fig. 8 illustrates the simulated distribution of muon-induced neutrons, showing that external neutrons constitute 20.0% of the visible neutron population. Consequently, ϵ_g is determined to be 20.0%, where the 2.1% uncertainty stems from the MC statistic.

Correction of detector energy response

When the Hamamatsu R5912-02 PMT receives a large number of photons, the anode output current increases sharply. In the region of the last dynode, electron accumulation forms a space charge, generating a reverse electric field that weakens the accelerating field strength and reduces the secondary electron emission efficiency. Consequently, the detector exhibited nonlinear energy responses above 10 MeV. Correcting this nonlinearity is essential for accurately reconstructing the true energy of events.

The energy response below 10 MeV was calibrated using gamma sources ⁶⁰Co and ¹³⁷Cs. For energies above 10 MeV, the calibration was achieved by comparing the Q-values of the reconstructed saturated events with the simulated muon energy deposition peaks. These three calibration points are shown in the upper inset of Fig. 9. The small inset within the upper panel provides a magnified view of the low-energy region, where the red and green lines represent cubic polynomial and linear fits, respectively. The error of the calibration point was determined by the error of the mean value obtained from a Gaussian fit to the simulated peak and the uncertainty of the peak position of the integral spectrum. The minimal discrepancy between the two fits below 10 MeV deviates by only a few percent near 10 MeV, justifying the approximation of the response as linear in this range. The lower inset of Fig. 9 displays the reconstructed energy of muon candidate events using the corrected energy response curve.

Fig. 9Full-screen View

The correction results of the detector’s energy response. The upper inset illustrates the differences between the cubic polynomial and linear fitting calibration points, where the zero-point calibration was performed using periodic pulse signals generated by a signal generator. The lower inset shows the two-dimensional distribution of energy versus Dis for reconstructed muon candidate events using the corrected calibration curve. Cuts applied indicates that electronic noise events have been deducted

Systematic uncertainties

Table 2 summarizes the total systematic uncertainties. The dominant systematic uncertainties arise from the reconstructed energy, energy scale, radius of the quartz glass vessel, wall thickness of the PE room, PE shielding layer, and Pb shielding layer.

The uncertainty in the energy scale originates from the uncertainty in the calibration point. The uncertainty of the calibration points arises from two aspects: the mean value error obtained through Gaussian fitting of the simulated energy spectrum peak and the uncertainty of the peak position of the integral spectrum. The uncertainty in the peak position in the integral spectrum is given by the difference between the distributions of the upper and lower errors in the reconstructed integral. Consequently, the uncertainty in the energy scale was determined by performing calibration curve fits using both the upper and lower error bounds of the calibration points and comparing the resulting flux differences between these two cases. The corresponding flux uncertainty is estimated to be 4.6%.

The uncertainty in the reconstructed saturated event energy affects the count in different energy bins. Its impact on the flux uncertainty is realized by affecting the best-fit value (k in the next section). To estimate this conservatively, we constructed two energy spectra: the first based on the reconstructed energy plus the upper error, and the second based on the reconstructed energy minus the lower error. By comparing the flux differences between these two distributions, we estimated the uncertainty in the reconstructed energy, and the corresponding flux uncertainty was estimated to be 5.1%.

The systematic uncertainties from the material dimensions were quantified through simulations. For example, a 0.2 cm uncertainty in the quartz glass vessel radius contributes 1.2% to the muon flux uncertainty, calculated as the relative change in reconstructed muon events when the radius varies by ± 0.2 cm. Finally, the total uncertainties of muon flux, muon-induced electron flux, and muon-induced neutron flux, derived through error propagation, were calculated as 7.0%, 7.1% and 7.2%, respectively.

The flux measurements of muons, muon-induced electrons, and muon-induced neutrons

It is assumed that the energy spectrum of cosmic-ray muons in the detector comprises contributions from three distinct components: muons, muon-induced electrons, and muon-induced neutron recoils. By employing the simulated energy distributions of the detector, we can fit the energy distribution of the experimental data points to derive both the spectral composition and flux in the measured spectrum. This fitting procedure is mathematically expressed as:(9)where k denotes the best-fit value, and f(E_sim) represents a step function approximation of the simulated energy distribution, constructed using 20 equal-width bins.

The left panel of Fig. 10 shows the measured energy spectrum of the detector. The χ²/n.d.f. of this fit is calculated to be 8.25/11 for the binning shown in the figure. k is determined to be 1.03 ± 0.19, where the uncertainty is dominated by the statistics of data. The right panel of Fig. 10 displays the neutron and alpha recoil events. Data in the energy region below 10 MeVee were derived from prior studies, and spontaneous fission of ²³⁸U was considered [16, 19]. Data above 10 MeVee originate from neutron recoil events selected via n-γ discrimination applied to the reconstructed saturated events. The uncertainties were calculated using Poisson statistics at the 90% confidence level. In addition, the pink area in the figure represents the energy spectrum of the recoil protons after calculating the quenching effect in the simulation, which originates from neutrons entering from outside the detector.

Fig. 10Full-screen View

(Color online) Experimental measurement results. Left panel: Measured event distribution in the energy region above 10 MeVee. The data points represent the contributions from muons and muon-induced electron energy in the detector. The solid line corresponds to the simulated energy spectrum, and the red dashed line represents the curve of the data fitting. Right panel: Measured energy spectrum of nuclear recoil events below 100 MeVee. The error bars above 10 MeVee represent the statistical uncertainties of the measured data in this study, with upper limits calculated at the 90% Poisson confidence level. Data points below 10 MeVee correspond to recoil proton events that are measured in ambient environments. The pink band shows the simulated spectrum of the recoil protons after applying quenching corrections. The inset in the right panel is a magnified view of the region below 5 MeVee, where the blue and brown regions denote the environmental gamma-ray and alpha background, respectively

Furthermore, because muon-induced electrons and neutrons originate externally to the detector, the fluxes of muon, muon-induced electrons, and muon-induced neutrons can be expressed as follows:(10)(11)where bins represents the i-th energy bin, represents the flux value of muons in the i-th bin, Wi is their width, and Wi = 5 MeV. Therefore, the total flux can be calculated as the sum of the flux values within all bins. The fluxes were calculated as follows:(12)(13)(14)As expected, the measured cosmic-ray muon flux was consistent with the results reported by Guo et al. [14]. The calculated muon-induced electron flux exceeds the muon event flux, as a single muon traversing the material generates multiple electrons along its path, with 93.2% of these electron events exhibiting energy depositions below 10 MeVee. For proton recoil events caused by muon-induced neutrons, after considering the quenching effect, 7.9% of the events were distributed above 10 MeVee. This shows that the flux results above 10 MeVee are below the measurement threshold of the detector (1 count per 412 days). Additionally, measurements at CJPL indicate that the neutron yield from muons passing through LS is approximately 0.028 per muon [42]. Therefore, no muon-induced neutron events above 10 MeVee were detected in the 412 day measurement results. At the 90% Poisson confidence level, the upper limit for such events above 10 MeVee is 2.3 counts, corresponding to a flux upper limit of 3.52 × 10^-9 cm^-2 s^-1.