Resolution function fitting

It was shown in Ref. [17] that the raw resolution function (obtained by the optical code) can not be used for a reliable resonance analysis. The reason is that the residual statistical fluctuations from the computationally intensive FLUKA+MNCP simulations of the neutron production and transport through the spallation target are not negligible relative to the fluctuations in the experimental data. Thus, using the raw resolution function in the analysis of the experimental data artificially and unnecessarily increases the involved statistical uncertainties. (Examples of smearing the initially smooth spectra by the raw resolution function may be found in later Figs. 5 and 6.) Clearly, a smoothed form of the resolution function is necessary, so as to avoid this adverse effect. It is also highly desirable that the smoothed form be efficiently parameterized, that is, that it is more compact than just a densely interpolated resolution function matrix filled with the values of a smoothed function (a matrix such as those from Figs. 1 and 2).

One way of proceeding would be to identify an analytical parametrization of the entire resolution function matrix. An example of such parametrization for a resolution function of EAR1 from Phase-1 of the n_TOF operation can be found in Ref. [2]. However, such an analytical form is difficult to identify and may no longer be appropriate when the alterations are introduced to the neutron production process. For example, the replacement of a spallation target between Phase-1 and Phase-2 of the n_TOF operation [4] notably affected the shape of a resolution function, rendering previous parametrization invalid. Furthermore, a resolution function for EAR2 differs from the one for EAR1 and requires its own dedicated parametrization. To avoid a tedious manual identification of new analytical forms, we take advantage of the machine learning techniques, in particular of the deep feedforward neural networks. The idea stems from a fact that the multilayer feedforward neural networks act as the universal approximators, capable of approximating any sufficiently well behaved function to any desired degree of accuracy [22, 23]. In other words, such networks can be thought of as “black box” fitting functions capable of modeling any function of practical importance. The application of neural networks to this task is a part of ongoing efforts to introduce the machine learning techniques into a widespread practice at n_TOF [24-26]. The possibility of applying the convolutional neural networks in unfolding the effects of the resolution function is also being investigated, with very promising results on the horizon [27].

We demonstrate the proof-of-concept by fitting a resolution function of EAR1, from Phase-3 of the n_TOF operation. To this end, we used the neural network training capabilities of the TMultiLayerPerceptron class [28] from root. Using root allows for a seamless integration of the end result (a trained neural network) within a vast majority of the data analysis codes from n_TOF. We provide a basic example of the code usage among the openly available data files [29].

Our goal is to fit a single form of a resolution function (either RT, $R_{\mathcal{E}}$ or Rλ) and reconstruct all other forms from this single fit by applying the appropriate transformations. This will ensure a perfect consistency between all forms of a resolution function, which would not necessarily be satisfied by fitting each form separately. Because of the described advantages of the Rλ representation (a uniformity of relevant λ values and insensitivity to a nominal flight path L), it is an obvious choice for fitting.

There is another consideration to be taken into account, that will allow for a greater flexibility in reconstructing particular forms of a resolution function. The resolution functions of the n_TOF facility (for different experimental areas) are affected by two separable contributions: (1) a neutron production and transport process inside a spallation target, as well as a neutron transport outside of it; and (2) a time distribution (a finite time width) of the primary proton beam from the CERN Proton Synchrotron irradiating the spallation target. The proton beam time distribution is Gaussian in shape with a standard deviation of σT=7 ns. Let ${\mathcal{R}}_T$ designate a resolution function in time of flight, without the effects of the proton beam width (as if σT=0). ${\mathcal{R}}_T$ is easily extracted from the raw results of the FLUKA+MNCP simulations processed by an optical transport code. A resolution function RT affected by the proton beam width is then obtained by a simple convolution with a temporal proton beam profile (a normalized Gaussian): $R_T\left(E\mathrm{,}{T}^{\prime }\right)=\frac{1}{{\sigma }_T\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }}{\mathcal{R}}_T\left(E\mathrm{,}\tau \right)\exp \left[-\frac{{\left(T^{\prime }-\tau \right)}^2}{2{\sigma }_T^2}\right]\text{d}\tau \mathrm{.}$ (8) On account of a linear relationship between T and λ from Eq. (4), a representation Rλ of a resolution function affected by a proton beam width may still be expressed as a convolution of a resolution function ${\mathcal{R}}_{\lambda }$ with an instantaneous proton beam, and a λ-equivalent of a temporal beam profile. From Eq. (4), it follows that the λ-width of this profile is dependent on a true neutron energy, as ${\sigma }_{\lambda }\left(E\right)=v_E{\sigma }_T$. Using Eq. (5) this dependence can be expressed as: ${\sigma }_{\lambda }\left(E\right)=c{\beta }_E{\sigma }_T=\frac{\sqrt{E\left(E+2m{c}^2\right)}}{E+m{c}^2}c{\sigma }_T\mathrm{.}$ (9) Thus, the required convolution takes the following form: $R_{\lambda }\left(E\mathrm{,}{\lambda }^{\prime }\right)=\frac{1}{{\sigma }_{\lambda }\left(E\right)\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }}{\mathcal{R}}_{\lambda }\left(E\mathrm{,}\text{Λ}\right)\exp \left[-\frac{{\left({\lambda }^{\prime }-\text{Λ}\right)}^2}{2{\sigma }_{\lambda }^2\left(E\right)}\right]\text{d}\text{Λ}\mathrm{.}$ (10) Equations (6) and (10) constitute all the transformation rules required for reconstructing any desired form of a resolution function from a proton-beam-width-free form ${\mathcal{R}}_{\lambda }$. Therefore, we apply a neural network fitting precisely to this, which is the most convenient representation ${\mathcal{R}}_{\lambda }$. The arguments E and ${\lambda }^{\prime }$ of a two-dimensional resolution function ${\mathcal{R}}_{\lambda }\left(E\mathrm{,}{\lambda }^{\prime }\right)$ play a role of the network inputs. A single scalar-function value for each pair of arguments plays a role of a single network output. However, since the true neutron energies E and the resolution function values ${\mathcal{R}}_{\lambda }$ span multiple orders of magnitude, for reasons of numerical stability we take the logarithms of these quantities for the input and output neutrons.

Here, we describe a neural network training and optimization procedure. First, the 450×60 numerical resolution function matrix was constructed, spanning 450 uniformly distributed λ values between -2 m and 7 m (2 cm steps), together with 60 isolethargically distributed ${{\displaystyle \log }}_{10}E$ values between 10^-3 eV and 10⁹ eV (five points per decade), creating a dataset of 27 000 points. It should be noted that the trained neural network is not be used for any kind of extrapolation outside of the range. It serves only as a smooth and compact parametrization throughout the physically meaningful range, to be used only for the resolution function reconstruction within a range entirely covered by the numerical matrix being fitted. Therefore, no cross-validation was required during a training procedure. In other words, all 27 000 entries could be used as training counts, without needing to reserve any of them for separate testing. Thus the quality of the trained network could be entirely assessed against the training data. We used a few simple and effective steps: we observed the saturation rate of a loss function during the training, and visually compared the agreement between the final fit and the raw data (see Fig. 4, soon to be discussed). We also closely monitored a distribution of fitting residuals.

We tested several training methods available in the TMultiLayerPerceptron class, using a few preliminary network structures. We did not observe any significant improvement– either in the computational efficiency or in the quality of the final results – over the default Broyden-Fletcher-Goldfarb-Shanno (BFGS) method with the default hyperparameter values [28]. Hence, we opted for the default training parameters, employing the sigmoid activation function. In TMultiLayerPerceptron, the binary cross-entropy is the only loss function associated with the sigmoid activation function. By testing a range of network structures, we identified an optimal structure consisting of three hidden layers, each composed of 15 neurons. The optimal structure was easily identified. Less than three hidden layers do not seem to have enough flexibility to recover the quick variations in the resolution function, even with highly increased number of neurons. Aside from the visually obvious overfitting, more than three hidden layers do not bring any further improvement in the resolution function reconstruction. Once the optimal three-layer structure was identified, the number of neurons was varied in steps of five, in different combinations throughout the layers (e.g. 15-10-20). By searching for the simplest structure providing a satisfactory resolution function reconstruction – without further improvement with an increasing number of neutrons – we quickly converged upon the optimal 15-15-15 structure, which is shown in Fig. 3. The network inputs are $x={\lambda }^{\prime }/\text{(1 cm)}$ and $y={{\displaystyle \log }}_{10}[E/\text{(1eV)}]$, with a single output $z={{\displaystyle \log }}_{10}\left[{\mathcal{R}}_{\lambda }/\left({\text{1cm}}^{-1}\right)\right]$. In general, the final state of a trained neural network depends on a particular training run owing to a random initialization of its weights and biases. With the optimal network structure and after a sufficient number of training epoch, these variations are essentially negligible. For the final training we used 10⁴ training epochs. A single-threaded training on AMD Ryzen 7 7735HS (3.2 GHz) CPU under Linux Mint 21.2 Cinnamon takes 1 min per 100 epochs, making a total of 100 min for 10⁴ epochs.

Fig. 3Full-screen View

Neural network structure used for modeling a resolution function of EAR1 from Phase-3 of n_TOF operation. Each fully connected hidden layer consists of 15 neurons. The inputs correspond to $x={\lambda }^{\prime }/\text{(1 cm)}$ and $y={{\displaystyle \log }}_{10}\left[E/\text{(1eV)}\right]$. A single output is $z={{\displaystyle \log }}_{10}\left[{\mathcal{R}}_{\lambda }/\left({\text{1 cm}}^{-1}\right)\right]$

Figure 4 compares a raw resolution function with a fitted one, for several values of a true neutron energy E. The effect of a proton beam width is negligible below σλ=1 cm. As per Eq. (9), a time width of σT=7 ns implies an equivalent energy limit of E=10 keV (corresponding, for EAR1 flight path of L=180 m, to the times of flight below 0.13 ms). Hence, at lower neutron energies the proton beam width of 7 ns does not have any significant effect upon the resolution function. For this reason a resolution function smeared by means of a numerical convolution from Eq. (10) is shown only at energies higher than 10 keV.

Fig. 4Full-screen View

(Color online) Raw resolution function fitted by a single neural network, at neutron energies of (a) 40 meV, (b) 10 keV and (c) 2 MeV. A raw and fitted function correspond to an instantaneous proton beam irradiating a spallation target. A beam width of 7 ns does not have a notable effect below 10 keV

Numerical resolution function reconstruction

To reconstruct the resolution function forms RT and $R_{\mathcal{E}}$ smeared by a proton beam width, one relies on the transformation rules from Eq. (6). In the general case, a smeared resolution function Rλ should first be obtained from an unsmeared fit ${\mathcal{R}}_{\lambda }$ – by means of Eq. (10) – and only then should a required transformation from Eq. (6) be applied. The reason is this. Combining Eqs. (1) and (4) yields a nonlinear relation between $\mathcal{E}$ and λ: $\mathcal{E}=m{c}^2\left\{{\left[1-{\left(\frac{v_E}{c}\frac{L}{L+\lambda }\right)}^2\right]}^{-1/2}-1\right\}\mathrm{.}$ (11) Because of this nonlinearity, the transformation between an unsmeared resolution function ${\mathcal{R}}_{\mathcal{E}}$ and a smeared $R_{\mathcal{E}}$ can no longer be expressed as a formal convolution, which bears upon the questions of computational complexity in applying a smearing transformation. We address these questions in Sect. 3.3.

Because of a simple relation between λ and T from Eq. (4), a transformation rule for RT from Eq. (6) involves a very simple derivative Dλ/DT=vE, yielding: $R_T\left(E\mathrm{,}{T}^{\prime }\right)=c{\beta }_E\times R_{\lambda }\left(E\mathrm{,}c{\beta }_E{T}^{\prime }-L\right)\mathrm{,}$ (12) with βE being defined by Eq. (5). A transformation for $R_{\mathcal{E}}$ is slightly more complex. By first obtaining $\lambda (\mathcal{E})$ and $\text{d}\lambda /\text{d}\mathcal{E}$ from Eq. (11), we obtain: $R_{\mathcal{E}}\left(E\mathrm{,}{\mathcal{E}}^{\prime }\right)=\frac{L{m}^2{c}^4}{{\left({\mathcal{E}}^{\prime }+m{c}^2\right)}^3}\frac{{\beta }_E}{{\beta }_{{\mathcal{E}}^{\prime }}^3}\times R_{\lambda }\left[E\mathrm{,}L\left(\frac{{\beta }_E}{{\beta }_{{\mathcal{E}}^{\prime }}}-1\right)\right]\mathrm{,}$ (13) where $\lambda (\mathcal{E})=L({\beta }_E/{\beta }_{\mathcal{E}}-1)$. By construction, these transformations preserve a norm of a resolution function.

We use a transformed resolution function from Eq. (13) to demonstrate the effect upon and an agreement with the n_TOF experimental data. Figure 5 shows two selected resonances from a recent measurement of the ⁵³Cr(n,γ) reaction in EAR1 [30]. Although the measurement was performed during n_TOF Phase-4, for presentation purposes, we use here a resolution function from Phase-3 as a first approximation of the resolution function from Phase-4. The plot shows a resolution-function-free reaction yield – manually constructed by appropriately scaling a neutron capture cross section from ENDF/B-VIII.0 database [31] – together with two resolution resolution-function-smeared yields, compared with the experimental data. (The preliminary experimental data are shown, as their analysis is not yet complete. The region between the resonances is still affected by the residual background contributions, not all of which have yet been subtracted. The relative background contribution inside the resonances is negligible for visual purposes, so that a meaningful visual comparison with the ENDF resonances can still be made.) Smeared yields have been obtained either by applying the raw, unparameterized resolution function or the one fitted by the neural network. Reaction yields (Y) transform precisely as the differential spectra from Eq. (7) and have been obtained by a transformation: $Y_{\mathcal{E}}\left({\mathcal{E}}^{\prime }\right)={\displaystyle {\int }_0^{\infty }}{Y}_E\left(E^{\prime }\right)R_{\mathcal{E}}\left(E^{\prime }\mathrm{,}{\mathcal{E}}^{\prime }\right)\text{d}{E}^{\prime }\mathrm{,}$ (14) wherein YE is a yield constructed from ENDF/B-VIII.0 data. Resolution function $R_{\mathcal{E}}$ was calculated according to Eq. (13), starting either from the raw or the fitted resolution function Rλ. The difference between two smeared yields illustrates our initial claim that the raw resolution function should not be used in the accurate data analysis, as it unnecessarily introduces artificial fluctuations into otherwise smooth data, or enhances the existing ones.

Fig. 5Full-screen View

(Color online) Selected resonances from n_TOF measurement of the ⁵³Cr(n,γ) reaction, compared to the reaction yields based on ENDF/B-VIII.0 data, with of without the resolution function having been applied

A simple and comprehensive insight into effects of applying either the raw or the smoothed resolution function may be obtained by applying them to a constant spectrum of unit height. Figure 6 shows these results for the $\mathcal{E}$-spectrum spanning 10 orders of magnitude in neutron energy, completely analogously to a resonant yield from Fig. 5. A detrimental effect of applying the raw resolution function is immediately evident, as opposed to a smooth end-result from the fitted resolution function. In that, Figure 6 shows the result of applying a resolution function unsmeared by a proton beam width (as if σT=0). The reason is that smearing the raw resolution function smooths a high energy part, thus obscuring the statistical fluctuations inherent in the simulated data, while these fluctuations are precisely what we are trying to exhibit here. Other than smoothing a high energy part of the spectrum obtained with the raw resolution function, a realistic value of σT=7 ns does not affect the shape of the folded spectra in a visible way (as they are displayed in Fig. 6). It should be taken into account that the data from Fig. 6 are displayed in 100 bins per decade. With finer binning, the fluctuations from the raw resolution function become even more pronounced. We note in passing that the n_TOF data are often analyzed in thousands of bins per decade (see, for example, Refs.[32, 33] for 5000 bins per decade applied to the neutron capture data or Ref. [34] for 2000 bins per decade applied to the fission data). While the resonance plots like those from Fig. 5 show how the underlying data are locally deformed by a resolution function, a plot from Fig. 6 clearly shows that the global data trend may also be affected.

Fig. 6Full-screen View

(Color online) Effect of applying either the raw or the fitted resolution function to a unit spectrum spanning 10 orders of magnitude in neutron energy. The data are shown in 100 bins per decade. A level of fluctuations introduced by the raw resolution function increases with denser binning

The resolution function affected spectra from Fig. 6 seem to be systematically higher than unity. This might suggest that the application of a resolution function violates a norm preservation, that is, it violates a conservation of the total number of underlying counts. This violation is only apparent and is addressed in the Appendix.

Resolution function class

To facilitate the use of a newly fitted resolution function at n_TOF, we have written a self-contained C++ class [35], serving as a simple interface for the evaluation of any desired resolution function form (Rλ, RT or $R_{\mathcal{E}}$), starting from a fitted proton-beam-width-free form ${\mathcal{R}}_{\lambda }^{*}$ (not necessarily normalized, as denoted by asterisk). The class serves the following functions:

(1) initialization of a trained neural network ($\to {\mathcal{R}}_{\lambda }^{*}$);

(2) numerical normalization ($\to {\mathcal{R}}_{\lambda }$);

(3) smearing due to a proton beam width ($\to R_{\lambda }$);

(4) transformation to alternative forms ($\to R_T$ or $R_{\mathcal{E}}$);

(5) computationally efficient implementation of the above operations.

The class houses the parameters of a trained neural network and initializes a TMultiLayerPerceptron object, allowing for a direct evaluation of a fitted, proton-beam-width-free (unnormalized) resolution function ${\mathcal{R}}_{\lambda }^{*}$.

In evaluating any form of a resolution function, the first operation to be performed is a numerical normalization such that ${\displaystyle {\int }_{-\infty }^{\infty }}{\mathcal{R}}_{\lambda }\left(E\mathrm{,}{\lambda }^{\prime }\right)\text{d}{\lambda }^{\prime }=1$, for any required value of E. The reason for performing a normalization at all is the fact that a fitting function without a priori imposed norm does not necessarily preserve a norm of the fitted data. Therefore, even if already-normalized data were fitted, the end-result needs to be a posteriori normalized. This also allows for the unnormalized raw data to be fitted, which may improve the quality of the fit, since the raw data may be unnormalized not only by a constant factor, but by an arbitrary E-dependent function, affecting a global twodimensional trend that the trained neural network needs to reproduce. The reason for a numerical normalization to be performed first is a computational efficiency, owing to the fact that both the convolutional smearing from Eq. (10) and the kinematic parameter transformations from Eqs. (12) and (13) preserve the norm. To clarify the point, let ${\mathbb{T}}_X$ denote a total transformation (a composition of smearing and kinematic-parameter-transformation operations) acting on unnormalized ${\mathcal{R}}_{\lambda }^{*}$ so as to produce an unnormalized, smeared and transformed resolution function $R_X^{*}$, in a sense: $R_X^{*}={\mathbb{T}}_X\left\{{\mathcal{R}}_{\lambda }^{*}\right\}$. Owing to a norm preserving property of all operations from ${\mathbb{T}}_X$, a normalized resolution function can be obtained in two ways: $R_X\left(E\mathrm{,}{X}^{\prime }\right)=={\mathbb{T}}_X\left\{\frac{{\mathcal{R}}_{\lambda }^{*}\left[E\mathrm{,}{\lambda }^{\prime }(X^{\prime })\right]}{{\displaystyle {\sum }_{\text{Λ}}{\mathcal{R}}_{\lambda }^{*}}\left(E\mathrm{,}\text{Λ}\right)\text{Δ}\text{Λ}}\right\}=\frac{{\mathbb{T}}_X\left\{{\mathcal{R}}_{\lambda }^{*}\left[E\mathrm{,}{\lambda }^{\prime }(X^{\prime })\right]\right\}}{{\displaystyle \sum {}_x}{\mathbb{T}}_X\left\{{\mathcal{R}}_{\lambda }^{*}\left[E\mathrm{,}\lambda (x)\right]\right\}\text{Δ}x}\mathrm{,}$ (15) where the sum from either denominator represents a discrete numerical integration. When the normalization is performed first, as in the first expression, the computationally expensive operations from ${\mathbb{T}}_X$ need to be performed only once in order to obtain a single point from $R_X$. If, on the other hand, the normalization was performed last (the second expression), then the same computationally expensive operations need to be performed unnecessarily many times during the norm calculation by means of a discrete summation. Furthermore, evaluating a trained network with many neural links (Fig. 3) is also moderately expensive. For this reason, we trigger a fresh normalization only when a new value of E is registered, which is different from a previous one requested by the user. This avoids the unnecessary repetitions of the same summation from Eq. (15).

For a user-supplied value of a temporal proton beam width σT, a normalized resolution function Rλ is smeared by performing a discrete version of a convolution from Eq. (10). This is the most significant bottle-neck of the numerical calculations, since a naive convolution algorithm is computationally expensive. When users request single points of a resolution function (one point at a time), little can be done to speed up a convolution algorithm itself. In this case, our class offers a possibility of a bilinear interpolation, which will soon be described. However, significant algorithmic improvements are available when a convolution needs to be calculated (i.e. smearing needs to be performed) over a grid of uniformly spaced λ-points. While a naive convolution algorithm is of excessive $\mathcal{O}\left(N^2\right)$ computational complexity – N being a number of required points – the famous Fast Fourier Transform algorithm manages a job in only $\mathcal{O}(N{{\displaystyle \log }}_{2}N)$ steps [36]. In this case, we employ an efficient Fast Fourier Transform implementation from Ref. [37], keeping the overall computational workload at a manageable level.

Having thus obtained a smeared Rλ, only a kinematic parameter transformation from Eq. (12) or (13) remains to be performed whenever the resolution function forms RT or $R_{\mathcal{E}}$ is requested. If the evaluation of RT over a grid of uniformly spaced T-points is needed, the computational advantages of the Fast Fourier Transform may again be relied on. This is because the smeared RT may still be expressed as a formal convolution. In fact, it is the original convolution from Eq. (8). At the same time, due to a linearity between λ and T from Eq. (4), a uniform grid of T-points corresponds to a uniform grid of λ-points, allowing for the Fast Fourier Transform to be applied. However, due to the n_TOF beam spanning more than 10 orders of magnitude in neutron energy, the measured time of flight spectra are spread over the multiple orders of magnitude. Therefore, the isolethargic spacing of T-points is commonly used, corresponding to an isolethargic spacing of λ-points (at a given E), for which a regular Fast Fourier Transform algorithm can no longer be used. On the other hand, whether a uniform of isolethargic spacing of $\mathcal{E}$-points is used (as is common practice), neither corresponds to a uniform spacing of λ-points, owing to a nonlinear relationship from Eq. (11). To extend a computational efficiency to all cases, our class allows users to activate the interpolation mode, allowing any form of the resolution function to be evaluated by means of a bilinear interpolation between pre-calculated resolution function points. Thus, a smeared and normalized resolution function Rλ is evaluated on a dense (E, λ) grid, allowing all the normalization and smearing operations to be performed in a single go, without later repetitions. This grid is isolethargically spaced over E and, more importantly, uniformly spaced over λ, so that the Fast Fourier Transform can be taken full advantage of during a smearing stage. The user-requested values of a resolution function are calculated by first obtaining a required $R_{\lambda }\left(E\mathrm{,}{\lambda }^{\prime }\right)$ point by a bilinear interpolation between four closest points from a pre-calculated grid (a simple bilinear interpolation algorithm may also be found in Ref. [37]). If necessary, a kinematic parameter transformation from Eq. (12) or (13) is then performed. For a sufficiently dense grid of pre-evaluated points it makes little numerical difference whether the kinematic parameter transformation is performed before or after the interpolation. However, prior to interpolation it has to be executed at four grid-points; after the interpolation it needs to be applied only once. We have selected a computationally efficient procedure.