IBD yield prediction

The Huber-Mueller model is a theoretical framework for predicting the antineutrino spectra produced by the fission reactions of four isotopes in reactors. Each of these isotopic antineutrino spectrum can be parameterized using the exponent of a fifth-order polynomial as follows: $s_l\left(E_{\nu }\right)=\exp \left({\displaystyle \sum _{p=1}^6}{\alpha }_{lp}{E}_{\nu }^{p-1}\right)\mathrm{,}$ (1) where $l=\left\{{}^{235}\text{U}{\mathrm{,}}^{238}\text{U}{\mathrm{,}}^{239}\text{Pu}{\mathrm{,}}^{241}\text{Pu}\right\}$, Eν is the ${\overline{\nu }}_{\text{e}}$ energy, and the αlps are polynomial coefficients for the isotope l. The αlp coefficients for ²³⁵U, ²³⁹Pu, and ²⁴¹Pu were derived using the conversion method by Huber [21], whereas the αlp coefficients for ²³⁸U were obtained using the summation method of Mueller et al. [16]. To incorporate the RAA in this study, we modified the isotopic antineutrino spectrum in Eq. (1) as follows: $S_l\left(E_{\nu }\right)=s_l\left(E_{\nu }\right)r_{\text{RAA}}\left(E_{\nu }\right)\mathrm{,}$ (2) where r_RAA(Eν) is the ratio of the RAA between the spectra measured in the Daya Bay experiment [32] and the Huber-Mueller model prediction. To evaluate r_RAA(Eν), we performed cubic spline interpolation within the provided energy range of 1.8-8 MeV and set it uniformly to 1 for energy values above 8 MeV.

The antineutrino yield per fission can be expressed as $\varphi \left(E_{\nu }\mathrm{,}t\right)={\displaystyle \sum _l}{f}_l\left(t\right)S_l\left(E_{\nu }\right)\mathrm{,}$ (3) where the fission fraction $f_l\left(t\right)$ represents the relative contribution of the isotope l to the fission reaction at time t. The event rate of antineutrinos emitted from the reactor core can be calculated as $\frac{\text{d}N}{\text{d}{E}_{\nu }}=\frac{W(t)}{{\displaystyle \sum _l}{f}_l\left(t\right){ϵ}_l}\varphi \left(E_{\nu }\mathrm{,}t\right)\mathrm{,}$ (4) where W(t) is the thermal power of the reactor at time t, ϵl is the mean energy released per fission of the isotope l, and the values for ϵl were obtained from Ref. [36].

In the standard three-flavor neutrino oscillation framework, the survival probability P_ee of ${\overline{\nu }}_{\text{e}}$ after propagating a distance L is given by [37] $\begin{array}{ll}{P}_{\text{ee}}\left(L\mathrm{,}{E}_{\nu }\right)\hfill & =P\left({\overline{\nu }}_{\text{e}}\to {\overline{\nu }}_{\text{e}};L\mathrm{,}{E}_{\nu }\right)\hfill \\ \hfill & =1-{{\displaystyle \cos }}^4{\theta }_{13}{{\displaystyle \sin }}^2\left(2{\theta }_{12}\right){{\displaystyle \sin }}^2\left({\text{Δ}}_{21}\right)\hfill \\ \hfill & -{{\displaystyle \cos }}^2{\theta }_{12}{{\displaystyle \sin }}^2\left(2{\theta }_{13}\right){{\displaystyle \sin }}^2\left({\text{Δ}}_{31}\right)\hfill \\ \hfill & -{{\displaystyle \sin }}^2{\theta }_{12}{{\displaystyle \sin }}^2\left(2{\theta }_{13}\right){{\displaystyle \sin }}^2\left({\text{Δ}}_{32}\right)\mathrm{,}\hfill \end{array}$ (5) where the θijs represent the neutrino mixing angles. The oscillation phases Δ_ij are given by $\begin{array}{l}{\text{Δ}}_{ij}=\frac{\text{Δ}{m}_{ij}^{2}L}{4E_{\nu }}\simeq \frac{1.267\text{Δ}{m}_{ij}^2\left[{\text{eV}}^{\text{2}}\right]L[\text{m}]}{E_{\nu }\left[\text{MeV}\right]}\mathrm{,}\hfill \end{array}$ (6) where $\text{Δ}{m}_{ij}^2$ denotes the mass-squared difference between the two mass eigenstates mi and mj, i.e., $\text{Δ}{m}_{ij}^2\equiv m_i^2-m_j^2$.

For short-baseline reactor neutrino experiments, considering that the term involving Δ₂₁ is negligible and $\text{Δ}{m}_{31}^2{m}_{32}^2$, Eq. (5) can be simplified to $\begin{array}{l}{P}_{\text{ee}}\left(L\mathrm{,}{E}_{\nu }\right)\approx 1-{{\displaystyle \sin }}^2\left(2{\theta }_{13}\right){{\displaystyle \sin }}^2\left(\frac{\text{Δ}{m}_{31}^{2}L}{4E_{\nu }}\right)\mathrm{.}\hfill \end{array}$ (7) Unless otherwise specified, ${{\displaystyle \sin }}^2{\theta }_{13}=(2.20\pm 0.07)\times {10}^{-2}$, $\text{Δ}{m}_{32}^2=(2.437\pm 0.033)\times {10}^{-3}{\text{eV}}^2$, and $\text{Δ}{m}_{21}^2=(7.53\pm 0.18)\times {10}^{-5}{\text{eV}}^2$ in this study based on the values from the Particle Data Group (PDG) 2022 [37].

As the ${\overline{\nu }}_{\text{e}}$ emitted by the reactor propagate to the LS detector, some of them engage in IBD reactions with the free target protons in the LS, which are denoted as ${\overline{\nu }}_{\text{e}}+p\to e^{+}+n$. In this process, the positron e⁺ rapidly deposits its energy and annihilates the surrounding electron e^- to form two 0.511 MeV gammas, generating a prompt signal. The neutron n scatters within the detector until it is thermalized and subsequently captured by hydrogen (99%) or carbon (1%) within ~200 μs, thereby releasing a 2.22 or 4.95 MeV gamma, respectively, and yielding a delayed signal [38]. An IBD event is identified by the prompt-delayed signal pair during such a brief interval. The measured IBD event number Mk in the k-th reconstructed prompt energy E_rec bin observed at a detector within the data acquisition time T_DAQ is therefore given by [5] $\begin{array}{c}{\text{M}}_k=\frac{N_{\text{p}}\epsilon }{4\pi L^2}{\displaystyle \underset{E_{\text{rec}}^k}{\overset{E_{\text{rec}}^{k+1}}{\int }}}\text{d}{E}_{\text{rec}}{\displaystyle \underset{T_{\text{DAQ}}}{\int }}\text{d}t{\displaystyle \underset{E_{\text{thr}}}{\int }}\text{d}{E}_{\nu }\\ \times \frac{\text{d}N}{\text{d}{E}_{\nu }}{P}_{\text{ee}}\left(L\mathrm{,}{E}_{\nu }\right){\sigma }_{\text{IBD}}\left(E_{\nu }\right)G\left(E_{\nu }\mathrm{,}{E}_{\text{rec}}\right)\end{array}$ (8) where N_p is the number of free target protons in the LS, ε is the detection efficiency of the detector, the IBD threshold energy E_thr~m_n-m_p+m_n~1.8 MeV, σ_IBD(Eν) is the cross-section of the IBD taken from Ref. [39], and $G\left(E_{\text{rec}}\mathrm{,}{E}_{\nu }\right)$ is a normalized Gaussian smearing function, which includes the energy resolution effect.

To simplify the calculation, we assumed that the detector has no energy leakage or LS nonlinearity [32]. Thus, the prompt energy $E_{\text{pro}}\sim E_{\nu }-0.78$ MeV, and E_rec is expected to obey the distribution $G\left(E_{\text{rec}}\mathrm{,}{E}_{\nu }\right)$ defined as follows [5]: $\begin{array}{l}G\left(E_{\nu }\mathrm{,}{E}_{\text{rec}}\right)\simeq \frac{1}{\sqrt{2\pi }{\delta }_{E_{\text{pro}}}}\exp \left\{-\frac{{\left(E_{\text{pro}}-E_{\text{rec}}\right)}^2}{2{\left({\delta }_{E_{\text{pro}}}\right)}^2}\right\}\mathrm{,}\hfill \end{array}\mathrm{.}$ (9) The energy resolution ${\delta }_{E_{\text{pro}}}$ is parameterized as $\begin{array}{l}\frac{{\delta }_{E_{\text{pro}}}}{E_{\text{pro}}}=\sqrt{{\left(\frac{p_0}{\sqrt{E_{\text{pro}}}}\right)}^2+p_1^2+{\left(\frac{p_2}{E_{\text{pro}}}\right)}^2}\mathrm{,}\hfill \end{array}$ (10) where p₀ quantifies the statistical fluctuations in the photons detected by the detector, p₁ is predominantly influenced by residual effects resulting from the spatial nonuniformity and temporal instability correction of the detector, and p₂ quantifies the effects associated with the photomultiplier tube (PMT), notably the PMT dark noise [5, 38].

For simplicity, we set p₀ = 0.08, p₁ = 0, and p₂ = 0 in this study for an energy resolution of 8% at 1 MeV in the detector. Therefore, under full reactor power and classical fission fractions conditions [32], the detector observes the energy spectrum of the IBD events (i.e., the reconstructed prompt energy spectrum) distorted by the RAA in one day, as shown in Fig. 1, and approximately 7473 IBD events are recorded.

Fig. 1Full-screen View

The virtual detector observes the reconstructed prompt energy spectrum distorted by the RAA with an exposure of $(2.9\times 5\times 1){\text{GW}}_{\text{th}}\cdot \text{tonnes}\cdot \text{day}$. The insert shows the cross-section of the IBD reaction [39] and the four isotope antineutrino spectra obtained by modifying the Huber-Muller model according to Eq. (2)

Simulated samples and targets in dataset

Considering the significant computational resources and time required for the integral terms in Eq. (8), Eq. (8) is typically converted to a discrete summation or matrix multiplication equivalent form in practical computations. In this study, the integral form of the reconstructed prompt energy spectrum is rewritten as an element of the row matrix $M_{1\times N_{E_{\text{rec}}}}$, which is given in Eq. (11). Each element of the matrix represents the measured IBD event number in the corresponding energy bin. $\begin{array}{ll}{M}_{1\times N_{E_{rec}}}\hfill & =X_{1\times 4}\cdot S_{4\times N_{E_{\nu }}}\cdot P_{N_{E_{\nu }}\times N_{E_{\nu }}}\cdot {\sigma }_{N_{E_{\nu }}\times N_{E_{\nu }}}\cdot R_{N_{E_{\nu }}\times N_{E_{\text{rec}}}}\hfill \\ \hfill & =X_{1\times 4}\cdot S_{4\times N_{E_{\nu }}}\cdot P\sigma R_{N_{E_{\nu }}\times N_{E_{\text{rec}}}}\mathrm{,}\hfill \end{array}$ (11) where the subscripts denote the dimensions of the corresponding matrices and the number 4 indicates the four isotopes. $N_{E_{\text{rec}}}$ is the number of bins in the reconstructed prompt energy spectrum while $N_{E_{\nu }}$ is the number of terms in the discretized sum for integration over Eν, which is also the number of bins in the extracted isotopic antineutrino spectrum. The ranges of Eν and E_rec are 1.8-10 MeV and 0.8-10 MeV, respectively. In this study, $N_{E_{\text{rec}}}$ was set to 80 based on the limit of the virtual detector energy resolution whereas $N_{E_{\nu }}$ was optimized to 401 after balancing model performance and computational cost. The element Xl in X_1×4 can be expressed as $\begin{array}{l}{X}_l={\displaystyle \sum _u^{N_t}}\frac{N_p\epsilon W\left(t^u\right)f_l\left(t^u\right)}{4\pi L^2{\displaystyle \sum {}_l}{f}_l\left(t^u\right){ϵ}_l}\text{Δ}t\text{Δ}{E}_{\nu }\text{Δ}{E}_{\text{rec}}\mathrm{,}\hfill \end{array}$ (12) where T_DAQ in Eq. (8) is divided into Nt time units of Δt, u is the time unit index, and $\text{Δ}{E}_{\nu }$ and $\text{Δ}{E}_{\text{rec}}$ are the bin widths of the extracted isotopic antineutrino spectrum and reconstructed prompt energy spectrum, respectively. In Eq. (12), W, fl, and ϵl are reactor-related parameters. W and fl vary as the reactor evolves while ϵl and the remaining parameters are constants. Xl is therefore referred to as the reactor dynamic evolution information.

Each row of $S_{4\times N_{E_{\nu }}}$ represents the binned antineutrino spectrum for isotope l, as given by $S_l\left(E_{\nu }\right)$. Both $P_{N_{E_{\nu }}\times N_{E_{\nu }}}$ and ${\sigma }_{N_{E_{\nu }}\times N_{E_{\nu }}}$ are diagonal matrices whose diagonal elements are given by $P_{\text{ee}}\left(L\mathrm{,}{E}_{\nu }\right)$ and σ_IBD(Eν), respectively. The role of $R_{N_{E_{\nu }}\times N_{E_{\text{rec}}}}$ is to map each Eν to a spectrum of E_rec. $R_{N_{E_{\nu }}\times N_{E_{\text{rec}}}}$ is therefore also referred to as the detector response matrix [Rqk], which is defined as follows: $\begin{array}{c}{R}_{qk}=R\left(E_{\nu }^q\mathrm{,}{E}_{\text{rec}}^k\right)=G\left(E_{\nu }^q\mathrm{,}\frac{E_{\text{rec}}^k+E_{\text{rec}}^{k+1}}{2}\right)\\ =G\left(E_{\nu }^q\mathrm{,}{E}_{\text{rec}}^k+\frac{\text{Δ}{E}_{\text{rec}}}{2}\right)\mathrm{,}\end{array}$ (13) where q is the index for binning Eν, $q\in \left[\mathrm{1,}\mathrm{2,}\cdots \mathrm{,}{N}_{E_{\nu }}\right]$, and $k\in \left[\mathrm{1,}\mathrm{2,}\cdots \mathrm{,}{N}_{E_{\text{rec}}}\right]$. In contexts that do not involve the oscillation parameters or unfolding, $P_{N_{E_{\nu }}\times N_{E_{\nu }}}$, ${\sigma }_{N_{E_{\nu }}\times N_{E_{\nu }}}$, and $R_{N_{E_{\nu }}\times N_{E_{\text{rec}}}}$ can be pre-multiplied to obtain the matrix $P\sigma R_{N_{E_{\nu }}\times N_{E_{\text{rec}}}}$.

The matrix multiplication relation in Eq. (11) provides the mathematical foundation for constructing the FNN architecture presented in Sect. 3. Furthermore, X_1×4 and $M_{1\times N_{E_{\text{rec}}}}$ respectively constitute a sample and its associated target in our dataset, which serve as a feature-label pair for supervised learning in the FNN model implemented in this study.

As described in Eq. (12), the fission fraction varies dynamically with burn-up as the reactor operates. In each reactor core refueling cycle, the cycle burn-up can be calculated as [32] $\text{Burn}-\text{up}=\frac{W\cdot D}{M_{{\text{U}}^{\text{ini}}}}\mathrm{,}$ (14) where W, D, and $M_{{\text{U}}^{\text{ini}}}$ represent the total thermal power of the reactor, the number of days since the refueling cycle started, and the mass of the initial uranium fuel loaded into the reactor, respectively. The unit for burn-up is ${\text{GW}}_{\text{th}}\cdot \text{day}\cdot {\text{tonne}}^{-1}$. Given that the real-time power output of the reactor is dynamic and cannot exceed its maximum capacity of 2.9 GW_th for safe operation, we used a random number generator for a normal distribution with a mean of 2.9 GW_th and downward fluctuation of 0.5% to determine the daily average power output of the virtual reactor [33]. By incorporating the fission fraction evolution data of the isotopes during a complete burn-up cycle from Ref. [32], we obtained the evolution of the fission fractions for the four main isotopes as a function of the operation day, as shown in Fig. 2.

Fig. 2Full-screen View

The evolution of fissile fractions for the four main isotopes in the reactor core as a function of operation day, which includes one complete refueling cycle [32]. The cumulative fission fraction of the four main isotopes used in this experiment is normalized to unity, and other isotopes contributing less than 0.3% are excluded from our analysis

Under the assumption that the thermal power and fissile fractions for the four main isotopes of the reactor are constant within each day, we accumulated the exposure over each 3-day interval as a sample to create a dataset of 600 simulated samples and their corresponding targets for subsequent analysis.

Mathematical foundations of FNN model

The NN is a powerful machine learning model that has been widely explored and applied across various fields. The universal approximation theorem [45, 46] implies that any continuous function can be approximated with arbitrary precision using an appropriate NN, even if the NN is an FNN with only one hidden layer containing a sufficient number of neurons. However, the internal structure and parameters of the NN in such scenarios often lack physical meaning or interpretability. This results in black-box models, which are not fully trusted by high-energy physicists. Therefore, we designed and implemented a white-box NN model in this study for converting the mathematical mapping function in Eq. (11) to a FNN model.

An FNN is typically composed of one to several single-layer perceptrons, which are considered the fundamental building units of the FNN and play a vital role in its overall functionality [47]. Each perceptron in the FNN follows the computational flow shown in Fig. 3 to process data. Forward and backward propagation are two phases in the NN training process that interact to optimize network performance.

Fig. 3Full-screen View

An example illustration of the structure of a single-layer perceptron along with forward (black flow arrows) and backward (red flow arrows) propagation

During the forward propagation phase, the perceptron performs computation by computing the dot product of the input vector $\overleftarrow{x}={\left[x_1\mathrm{,}{x}_2\mathrm{,...,}{x}_N\right]}^{\text{T}}$ with the weight coefficient vector $\overleftarrow{w}={\left[w_1\mathrm{,}{w}_2\mathrm{,...,}{w}_N\right]}^{\text{T}}$, adding the bias b, and applying the activation function ψ to yield the activation result y as the output. The discrepancy between the output y and target $\widehat{y}$ is then calculated using the loss function $L(y\mathrm{,}\widehat{y})$. Forward propagation provides the foundation for evaluating network performance. Backward propagation in turn determines how the network parameters (weights and bias) are updated to reduce loss. It can be described as $\begin{array}{l}{{\omega }^{\prime }}_{\mu }={\omega }_{\mu }-\eta \times \left[{\nabla }_{{\omega }_{\mu }}L(y\mathrm{,}\widehat{y})+\lambda {\omega }_{\mu }\right]\mathrm{,}\hfill \end{array}$ (15) $\begin{array}{l}{b}^{\prime }=b-\eta \times {\nabla }_{b}L(y\mathrm{,}\widehat{y})\mathrm{,}\hfill \end{array}$ (16) where ωμ and ${{\omega }^{\prime }}_{\mu }$ represent the μ-th weight coefficient of the current and subsequent steps, respectively; b and b’ the biases of the current and subsequent steps, respectively; and η and λ are the learning rate and weight decay rate, respectively. This iterative update process of the parameters based on the computed gradients allows the NN to learn and improve its predictions over time.

To allow matrix multiplication in the perceptrons, the bias b must be eliminated, i.e., set to zero. The absence of negative values in our data flow justifies the use of the default rectified linear unit (ReLU) activation function, which is defined as $\psi (z)=\max \left\{\mathrm{0,}z\right\}$. This setup also permits the perceptrons to be chained to perform successive matrix dot product operations, which is integral to the development of our FNN model.

As shown in Fig. 4, the architecture of the FNN model consists of three layers comprising, from left to right, the input, hidden, and output layers with four, $N_{E_{\nu }}$, and $N_{E_{\text{rec}}}$ neurons, respectively. The neurons in adjacent layers are connected using a fully connected approach; that is, each neuron in one layer is connected to every neuron in the subsequent layer with no connections between neurons within the same layer. The training process of the model starts from the input layer, at which each neuron receives the reactor dynamic evolution information corresponding to its fission isotope. The output of the hidden layer is the scaled total spectrum [hq] of antineutrinos emitted by the reactor. The output layer then provides a predicted reconstructed prompt energy spectrum [yk]. The two weight coefficient matrices W⁽¹⁾ and W⁽²⁾ correspond to the transposes of the matrices $S_{4\times N_{E_{\nu }}}$ and $P\sigma R_{N_{E_{\nu }}\times N_{E_{\text{rec}}}}$, respectively. The matrix W⁽¹⁾ contains the fission isotope antineutrino spectra to be extracted, which are learned during training. In contrast, the matrix W⁽²⁾ is fixed as $P\sigma R_{N_{E_{\nu }}\times N_{E_{\text{rec}}}}^T$ because it is assumed to be a constant matrix without uncertainties in this study. The FNN is therefore a supervised learning model that iteratively refines W⁽¹⁾ to minimize the discrepancies between its outputs and corresponding targets.

Fig. 4Full-screen View

The FNN is a white-box model that describes the mapping relation between the reactor dynamic evolution information and reconstructed prompt energy spectrum. The architecture of the FNN model includes an input layer, hidden layer, output layer, and two sets of weight coefficient matrices W⁽¹⁾ and W⁽²⁾. The weight values between neurons associated with connections of the same color form the rows of the weight coefficient matrix

Training strategy

All the samples generated in Sect. 2.2 were utilized solely to train the FNN model. The validation and testing processes were omitted. This approach was chosen because our aim is to minimize the loss function during the training process to determine the optimal W⁽¹⁾ for extracting the four main isotopic antineutrino spectra. Our focus is on optimizing spectra extraction performance rather than evaluating model performance across various datasets, as well as on simplifying the process and aligning with our primary research objective.

The loss function is a fundamental component in deep-learning models. It serves as the criterion for evaluating how well the model predictions match the actual outcomes and provides a numerical indicator of model accuracy. The Combined Neyman-Pearson (CNP) chi-square model is a statistical model frequently employed in HEP experiments to quantify the error between predicted and measured values [48]. Based on this model, we define the loss function for the FNN model as $\begin{array}{l}{\chi }_{\text{CNP}}^2={\displaystyle \sum _{k=1}^{N_{E_{\text{rec}}}}}\frac{{\left[M_k-y_k\left(W^{\left(1\right)}\right)\right]}^2}{3/\left[\frac{1}{M_k}+\frac{2}{y_k\left(W^{\left(1\right)}\right)}\right]}\mathrm{,}\hfill \end{array}$ (17) where Mk is the IBD event number in the k-th bin for the measured reconstructed prompt neutrino energy spectrum given by Eq. (8) and yk is the corresponding predicted value output of the model. We used this loss function to guide the optimization process of W⁽¹⁾ during the training process so that the FNN was driven towards increasingly precise predictions.

After defining the loss function, it is essential to select a suitable optimizer, learning rate schedule, batch size, and epoch, among other hyperparameters. Following hyperparameter tuning using the Optuna framework [49] and extensive testing, we developed two training strategies denoted as the short- and long-epoch strategies to investigate the performance of the FNN model in extracting the antineutrino spectra of the four fission isotopes from the reconstructed prompt energy spectrum [50]. As shown in Table 1, a critical commonality between these two strategies is the segmentation of the hidden layer in the FNN model into multiple partitions or parallel hidden layers. This setup allows distinct learning and weight decay rates to be assigned to each partition to facilitate differential performance outcomes. Because the focus in this study is not on the isotope antineutrino spectra above 8 MeV, i.e., in the (303, 401] partition or the matrix $P\sigma R_{N_{E_{\nu }}\times N_{E_{\text{rec}}}}$, we fixed their learning and weight decay rates to zero and disabled the gradient calculations for the corresponding weight coefficients. Additionally, we set the initial value of W⁽¹⁾ based on the Huber-Mueller model.

As indicated by their names, the main distinction between the short- and long-epoch strategies lies in the epochs. The short-epoch strategy leverages the AdamW [51] optimizer with non-zero weight decay rates for faster loss reduction. In contrast, in the long-epoch strategy, the Adam [52] optimizer is applied without weight decay, i.e., the weight decay rates are set to zero. Superior convergence results were obtained using the long-epoch strategy. The results are presented and discussed in Sect. 4. As illustrated in Table 1, these circumstances also led to minor differences in the configurations of the learning rate schedulers. Nonetheless, the same metric, i.e., the sum of the losses for all samples denoted as ${\chi }_{{\displaystyle \sum }\text{CNP}}^2$, was monitored in both schedulers.

We also extracted the antineutrino spectra of the four fission isotopes using the χ² minimization method to provide a comparison and benchmark for the FNN model. We employed the Minuit2 minimization library from ROOT [53] to implement this method. ${\chi }_{{\displaystyle \sum }\text{CNP}}^2$ was used as the objective function to be minimized to find the best fit. The same dataset as that for the FNN model was used as the measured value in this fitting process. In contrast, the predicted value was derived from Eq. (11) where the $S_{4\times N_{E_{\nu }}}$ matrix elements corresponding to $\le 8$ MeV are the parameters to be fitted and the remaining elements considered as fixed parameters in the fitting procedure. We adopted the “Combined” minimizer algorithm to minimize the objective function with initial fitting values from the Huber-Mueller model and fitting step sizes of 1% of the order of magnitude of these values. We set the tolerance for the fitting procedure to 1×10^-30. The fitting stopped automatically only when the improvement in the ${\chi }_{{\displaystyle \sum }\text{CNP}}^2$ value between consecutive iterations fell below this threshold.

The FNN model was implemented using PyTorch [54], a Python-based deep learning library that supports both CPU and GPU platforms and is one of the mainstream tools for developing and training NN models. A NVIDIA GeForce RTX 3060 Ti GPU platform was used to deploy the FNN model, whereas tasks involving Optuna and ROOT were performed on two identical servers, each of which was equipped with two 28-core Intel(R) Xeon(R) Gold 6330 CPUs @ 2.00 GHz.