Selection of Experimental data

The first step in the nuclear data evaluation process is usually to carefully select experimental data since using all the experiments from the EXFOR database [26] without any selection normally leads to the computation of very large chi squares between model calculations and experiments. The large chi squares can be attributed partly to the presence of discrepant and outlier experiments as well as the presence of experiments with unreported or under-reported uncertainties. In this work, outliers were treated using a binary accept/reject approach. For example (as carried out also in Ref. [15]), experiments that were observed to be inconsistent with all or most of the data sets available for a particular channel and energy range were assigned a binary value of zero and therefore were not considered. In addition, experiments that deviate from the trend of the evaluations from the major nuclear data libraries (when available), were not considered. Similarly, experiments without reported uncertainties were penalized by assigning them a binary value of zero except in the cases where the considered experimental data set(s) under consideration were the only experiment(s) available in the energy range of interest. In this case, a 10% relative uncertainty is assumed for each data point of that particular experimental data set. Also, in a situation where these experimental data sets with unreported uncertainties have been considered and reported to be of reasonable quality in Ref. [27], a 10% relative uncertainty was assumed. In some cases, as carried out also in Ref. [27], experiments that were found to be close to the threshold energies and hence, usually difficult to measure, were not considered in our optimization procedure.

We note that the selection and rejection of experiments introduces bias into the evaluation process as is normally the case with all other evaluated nuclear data libraries. Ref. [4] however argues that, instead of rejecting discrepant and outlier experiments as carried out in this work, subjective weights which take into account the quality of each experimental data set, could instead, be assigned to each data set. In this way, ‘bad’ experiments would be assigned with smaller weights and hence contribute less to the optimization and therefore, no experimental information would be discarded. In Ref. [28], the use of Marginal Likelihood Optimization (MLO) for the automatic correction of the uncertainties of inconsistent experiments was proposed. Another approach as presented in Ref. [29], has been to identify Unrecognized Sources of Uncertainties (USU) in experiments and try to include them in evaluations. These approaches were however, not utilized in this work.

In this work, three experimental data types were utilized in the evaluation procedure: (1) reaction cross sections, (2) residual production cross sections, and (3) the elastic angular distributions. In the case of the reaction cross sections, the following eight reaction channels were considered: (p,non-el), (p,n), (p,3n), (p,4n), (p,2np)g, (p,2np)m, (p,γ) and (p,xn) cross sections. These channels were selected because (1) experimental data were available within the considered energy range for the considered cross sections of ⁵⁹Co and (2) because we desire a general purpose file which is optimized to many reaction channels with experimental data as much as possible. A total of 169 reaction cross section experimental data points were used in the optimization. For the residual production cross sections, a total of 141 experimental data points were considered for the following cross sections:

• ⁵⁹Co(p,x)⁴⁶Sc,

• ⁵⁹Co(p,x)⁴⁸V,

• ⁵⁹Co(p,x)⁵²Mn,

• ⁵⁹Co(p,x)⁵⁵Fe,

• ⁵⁹Co(p,x)⁵⁵Co,

• ⁵⁹Co(p,x)⁵⁶Co,

• ⁵⁹Co(p,x)⁵⁷Co,

• ⁵⁹Co(p,x)⁵⁸Co,

• ⁵⁹Co(p,x)⁵⁷Ni.

In the case of the elastic angular distributions, the angles were considered from 1 to 180 degrees while the incident energies considered were from 5 to 40 MeV and a total of 185 experimental data points were considered.

Model Calculations

Model calculations were performed using the TALYS version 1.9 [30] code. TALYS is a state of the art nuclear reactions code used for the analysis of nuclear reactions for a number of incident particles. These particles include neutrons, protons, deuterons, tritons, ³He- and α particles within the 1 keV to 200 MeV energy range [30].

In Table 1, similar to Ref. [20], the models considered in this work are listed. A total of 52 different physical models were randomly varied. A model as used here represents either a complete nuclear reaction model such as the Jeukenne-Lejeune-Mahaux optical model, or a sub-model, and in some cases, components of a model or sub-model. On the other hand, a model set or combination, represents a vector of these models or sub-models, coupled together in the TALYS code for nuclear reaction calculations. For example, statepot is a flag used to invoke the optical model parameterisation for each excited state in a Distorted Wave Born Approximation (DWBA) or coupled-channels calculation. As can be seen in Table 1, the TALYS code contains 4 pre-equilibrium models, 6 level density models, 8 gamma-strength function models, 4 mass models, and 4 Jeukenne-Lejeune-Mahaux (JLM) optical models, among others. Several of these models are linked together for nuclear reaction calculation with the outputs of some models used as input to other models. Even though all these models are available in TALYS, the default models have been used in proton evaluations for the proton sub-library of the TENDL library [20].

Since we assumed that all models were equally important a priori, our algorithm sampled each model type within it’s bounds. To achieved a non-informative prior as much as possible, each model type was sampled from a uniform distribution within it’s bounds. Models were then randomly drawn from each model distribution to create a total of about 200 model combinations. In the case of the level density (ld) model for example, there are 6 different models available in TALYS. Hence, these models were each assigned a unique ID and then sampled from model 1 to 6 in such a way that each ld model was equally likely to be selected. In addition, parameters to each model combination were also randomly generated to create multiply TALYS input files each with unique parameters, which were then run with the TALYS code to produce random cross section curves as well as angular distributions.

It should be noted that each model combination together with their unique parameter set gives different cross section curves and angular distributions. In Fig. 1 for example, ⁵⁹Co(p,n) cross section curves computed with three different model combinations are presented. Also, the calculated cross sections are compared with the default TALYS models (TENDL-2019) and the JENDL/He-2007 and JENDL-4.0/HE libraries as well as with experimental data. In the case of the models labelled (A) in the figure, the following models: ldmodel 6: Microscopic level densities (temperature dependent Hartree-Fock-Bogolyubov (HFB), Gogny force) from Hilaire’s combinatorial tables; preeqmode 4: Multi-step direct/compound model; strength 4: Hartree-Fock-Bogolyubov tables and massmodel 1: Möller table + other models, were used. For the model combination labelled (B), the following models were utilized: ldmodel 2: Back-shifted Fermi gas model; preeqmode 2: Exciton model:- Numerical transition rates with energy-dependent matrix element; strength 1: Kopecky-Uhl generalized Lorentzian, + other models. No mass model was invoked in this case of (B). In the case of models labelled (C), the cross section curves were produced using ldmodel 3: Generalized superfluid model; preeqmode 2: Exciton model:- Numerical transition rates with energy-dependent matrix element; strength 4: Hartree-Fock-Bogolyubov (HFB) tables; massmodel 1: Möller table, + other models. The other models as stated here are default TALYS models that were included in the model calculations. For each model combination, two cross section curves are shown: (1) with default parameters and the other with perturb parameters. These initial calculations give an idea of the shape and possible spread of the cross section curves for each model combination. It was observed in this work that the shape of the cross section curves produced with each model combination was not affected with variation of the parameters and hence, a smaller number of parameters variations can be carried out around each model combination in the case of the initial or parent generation. It can be seen from the figure that model combination (B) gives a better fit to the available experimental data than model combinations (A) and (C). However, care must be taken as there is no guarantee that a model set which produces a good fit to experimental data for a particular channel would produce good fits for other channels. The ‘best’ file in Fig. 1 represents the optimal file selected from the parent generation (Gen. 0) taking into consideration both reaction and residual production cross sections as well as the elastic angular distributions. The parent generation here refers to the initial random nuclear data (ND) files produced from the variation of both models and their parameters and signifies the initial generation from which all subsequent generations were produced. Detailed description of the optimization procedure is presented in Sects. 2.3 and 2.4.

Fig. 1Full-screen View

(Color online) ⁵⁹Co(p,n) cross section computed with three different model combinations compared with cross sections computed with the default TALYS models (TENDL-2019), the JENDL/He-2007 and JENDL-4.0/HE libraries, as well as experimental data. The ‘best’ file here is the optimal file selected from the parent generation taking into consideration, reaction and residual production cross sections as well as the elastic angular distributions. For each model combination A, B and C, calculations were carried out with default and then perturbed model parameters.

By randomly varying parameters of each model combination using the TALYS code package [30], a set of random nuclear data files in the Evaluated Nuclear Data File (ENDF) format were produced. Each TALYS input file contains a set of these models as presented in Table 1 as well as the parameters to these models (see Table 2). In a situation where a model type is not explicitly listed in the TALYS input file, a default TALYS model was used.

Similar to Table 1, a list of selected model parameters with their corresponding parameter widths (uncertainties) is presented in Table 2. A similar table is presented in Ref. [20]. The parameter uncertainties were obtained from the TENDL project [20, 30]. We note here that these parameter widths or uncertainties were obtained using default models in TALYS. We however think that this is a good enough starting point for our iBMC methodology. These parameter widths were used as prior uncertainties for the sampling of the parameters from a uniform distribution for each model combination used in this work. From the table, the parameter widths (uncertainties) are given as a fraction (%) of their absolute values except in the case of level density parameter a, and the gπ and gν parameters of the pre-equlibrium model, where the uncertainties are given in terms of the mass number A. As can be seen from the table, the parameters have been grouped under the following model types: the optical model, pre-equilibrium, gamma-ray strength function and level density models. In the case of the optical model, the phenomenological optical model parameters are used for a potential of the Koning-Delaroche form [32] while the parameters for the semi-microscopic optical model are used for the Jeukenne-Lejeune-Mahaux (JLM) microscopic optical model potential [33]. Since the fission cross section was not considered in this work, fission models and their parameters are not presented in the table.

From Table 2, the geometrical parameters $r_V^p$ and $a_V^p$ denotes the radius and diffuseness parameters of the real central potential, $a_D^p$ is the surface diffusivity and $r_c^p$ is the Coulomb radius. $r_{SO}^p$, $a_{SO}^p$, $v_{SO1}^p$, $v_{SO2}^p$, $w_{SO1}^p$, $w_{SO2}^p$ are the spin-orbit potential parameters and v₁ - v₄ are adjustable parameters used in the computation of the depth of the real central potential. w₁ - w₂ and d₁ - d₄ are the adjustable parameters used in the computation of the volume and surface absorption optical model potential respectively. The superscript p denotes proton induced reactions while the subscripts V, D, SO, and W denotes the volume-central, the surface-central, spin-orbit potentials and the imaginary depth of the optical model respectively. A more detailed description of these parameters can be found in Refs. [30-32]. In the case of the JLM model parameters, λV, λW, λ_V1, and λ_W1 denotes the overall real, imaginary, real isovector, and imaginary isovector potential depth normalization factors of the semi-microscopic optical model potential respectively [33]. ${\lambda }_{V_{SO}}$ and ${\lambda }_{W_{SO}}$ are the real and imaginary spin orbit (SO) potential depth normalization factors respectively. It should be noted that all parameters were assumed to be energy independent in our calculations.

In the case of the pre-equilibrium model, gπ and gν are the single-particle state densities, M² is the average squared matrix element used in the exciton model and Rγ is the parameter for the pre-equilibrium γ-decay. C_strip, C_break and C_knock are the stripping, break-up and knock-out contributions used to scale the complex-particle pre-equilibrium cross section per outgoing particle respectively [31]. E_surf is the effective well depth for surface effects in the exciton model. Rνν, Rπν, Rππ, and Rνπ are respectively, the neutron-neutron ratio, the proton-neutron ratio, the proton-proton ratio and the neutron-proton ratio of the matrix element used in the two component exciton model. In the case of the level density model parameters, k_rot is the rotational enhancement factor of the level density model, σ² is the spin cut-off parameter which represents the width of the angular momentum distribution of the level density, a is the level-density parameter which is related to the density of single-particle states near the Fermi energy, and E₀ and T are the the back-shift energy and the temperature used in the constant temperature model respectively. Rσ is a global adjustable constant for the spin cut-off parameter. In the case of the gamma ray strength function, Γ_γ is the average radiative capture width and ${\sigma }_{E\mathcal{l}}$, $E_{E\mathcal{l}}$ and ${\Gamma }_{E\mathcal{l}}$ are the strength, energy and width of the giant resonance, respectively. Where $E\mathcal{l}$ represents the electric (E) multipole type and $\mathcal{l}$ the multipolarity [30]. These gamma-ray strength function parameters are used in the description of the gamma emission channel [30]. It should be noted here that not all the TALYS model parameters are listed in Table 2, a more complete list of all the model parameters can be found in Ref. [31] while the parameter uncertainties can be found in Refs. [4, 30]. As can be seen from the table, rather large parameter uncertainties were assigned to the k_rot, C_strip, C_break and C_knock in order to take into account the shortcomings of the corresponding pre-equilibrium models [20]. As mentioned earlier, these parameters were varied simultaneously for each model combination to create multiple TALYS input files.

The incident proton energy grid for model calculations was selected taking into consideration the incident energies of the available experimental data as well as the availability of computational resources. It must be stated here that in order to observe the potential coverage of the parameter space for each model combination, two initial TALYS calculations are normally executed using two parameter vectors with extreme values, e.g. with: (1) the lower and (2) the upper bound values of each parameter. The lower and upper bounds corresponds to the minimum and maximum values of each parameter from the associated uniform prior parameter distributions. From these calculations, we are able to observe (visually) the trend and potential spread of the cross sections and the angular distributions of interest with respect to the scattered experimental data available. This was done because as mentioned earlier, the parameter widths as presented in Table 2 were obtained by comparing the cross section curves produced with default TALYS model combination with varying parameters and hence, might not be applicable to other model combinations under consideration. Additionally, by observing the lower and upper bound curves, we are able to determine if any changes in the parameters within the uncertainties as given in Table 2, would significantly affect the spread of the considered cross sections. We note that since no parameter correlations were included in these initial calculations, this approach might not describe the entire coverage of the parameter space. This however, is helpful in determining the potential uncertainty band for each cross section. Also, it must be stated here that we normally start with calculations with relatively small number of excitation energy bins of 20. The bin size is then increased to 60 for subsequent iterations in order to improve the accuracy of the TALYS results but at a higher computational cost. We note also that the total number of 200 model combinations used in this work does not cover the entire model space. We are however of the opinion that the number is adequate for demonstrating the applicability of the proposed method.

The 200 model combinations were each run with the TALYS code while varying all the model parameters to produced a large set of random ENDF nuclear data files. The random ENDF files produced were converted into x–y tables for comparison with experimental data. The reason for going through the TALYS → ENDF → x–y tables route instead of using the x–y tables directly produced by TALYS is because, during processing of TALYS results into the ENDF format, renormalizations are sometimes applied to the results using auxiliary codes such as the autonorm utility within the T6 code package [30] which is used to smoothen the results of TALYS at the incident energies where different model results are joined together [30]. These normalizations should normally not affect the results significantly, however, in order to enable a fair comparison between our evaluations and that of other nuclear data libraries which were only available in the ENDF format, the TALYS → ENDF → x–y route was used.

Bayesian calibration and selecting the winning model

The ultimate goal of a Bayesian calibration is to maximize the likelihood that model outputs are statistically consistent with experimental data [34]. The first step in this work involves the pre-selection of models and their parameters, followed by the quantification of the uncertainties of model parameters and the determination of their distributions. As mentioned earlier, the parameter uncertainties were adopted from the TENDL library project. In this work, it was assumed that we have no prior information on the models as well as on their parameters and hence, the models and their parameters were both sampled from uniform distributions. This assumption is entirely not true since the model types used (see Table 1) were pre-selected using expert judgment and to some extent, model sensitivity analysis. The model sensitivity analysis involved holding constant all other models as the default TALYS models while changing the model of interest, one-at-a-time. The spread in the cross sections give an indication of the sensitivity of each model to the cross section of interest. For example, in Fig. 2, excitation functions for the ⁵⁹Co(p,3n) cross section computed with the four mass models implemented in the TALYS code, is presented. No changes in the ⁵⁹Co(p,3n) cross section were observed after changing and running the different mass models one-at-a-time for the energy range under consideration. This gives an indication that the mass models implemented in TALYS are not sensitive to the ⁵⁹Co(p,3n) cross section. The ‘Best file (5th Gen.)’ in the figure represents the final evaluation obtained in this work while the TENDL-2019 represents the evaluation from the TENDL-2019 library.

Fig. 2Full-screen View

(Color online) ⁵⁹Co(p,3n) cross section computed using the different mass models implemented in TALYS. All other models were maintained as the default TALYS models while changing each mass models one-at-a-time. The ‘Best file (5th Gen.)’ represents our evaluation. It can be observed that changing the mass models did not have any significant impact on the (p,3n) cross section.

The pre-selection of models was carried out in order to limit the model space because of computational resource constraints. We do however understand that by using expert judgement to select a subset of the models available in the TALYS code system [30], we introduce some user bias into the initial model selection process. We however note here that since we are limited by the models available in TALYS in this case, selection bias cannot be entirely excluded from the process. A possible solution would have been to use all the models available in TALYS as well as additional models from other nuclear reaction codes such as EMPIRE [35]. This would increase the model space tremendously but would ultimately lead to a higher computation cost. A more detailed study on the use of Bayesian model selection and the Occam’s Razor in the selection of model combinations within a Total Monte Carlo framework is proposed and presented in Ref. [36].

Furthermore, it was observed in this study that some model combinations produced model outputs with non-smooth curves which appear unphysical and therefore were excluded from subsequent model runs. For example, the ⁵⁹Co(p,γ) cross section computed with three different model combinations (A, B and C) are compared with experimental data and the TENDL-2017 evaluation. Resonance-like structures which are difficult to explain from nuclear reaction theory were observed for ⁵⁹Co(p,γ) between about 4 to 10 MeV for model combination (A) as can be seen in Fig. 3, and therefore, this model combination was excluded from subsequent calculations. Note that model combinations (A), (B) and (C) are the same model combinations presented earlier in Fig. 1. As can be seen from Fig. 3, a relatively ‘bad’ model set such as model combination (A) can be difficult to identify by using only the χ² as the goodness of fit estimator since it can be observed that the cross section curves produced by this model combination was able to reproduce quite favorably, the experimental data from Butler (1957). It can also be observed from the figure that the cross section curves passes through some experimental data points from Drake (1973) which results in a relatively low χ² making it difficult to identify this model combination as a ‘bad’ model set. Because of this, visual inspection is sometimes needed to identify these ‘bad’ model sets. This should however, be carried out on an isotope-by-isotope and channel-by-channel basis.

Fig. 3Full-screen View

(Color online) ⁵⁹Co(p,γ) cross section computed with three different model combinations (A, B and C), are compared with cross sections computed with the default TALYS models (TENDL-2017), the JENDL/He-2007 and JENDL-4.0/HE libraries, as well as experimental data. The ‘best’ file is the file selected from the parent generation taking into consideration the reaction and residual production cross sections, and the elastic angular distributions. It can be observed that the curves from model combination (A) have non-smooth curves between 4 to 12 MeV which are difficult to explain from nuclear reaction theory.

In the case of the input parameters to the models, all the parameters available within the TALYS code were varied. In this way, we were able to largely exclude parameter selection bias from our analyses. We however note that since only the TALYS code was used, we are constrained by the parameters available in TALYS. We also note that not all the model parameters are sensitive to the considered cross sections or the elastic angular distributions and therefore, these parameters could have been identified through a parameter sensitivity analyses and then excluded from the parameter variation process. This approach was however not carried out in this work.

Now, suppose that we have a set of J models, ${\overrightarrow{M}}_j$, where j=1,2,...,J, given a set of experimental data (${\overrightarrow{\sigma }}_E$) with corresponding uncertainties ($\Delta {\overrightarrow{\sigma }}_E$). Each model combination also consists of a vector of K model parameters, ${\overrightarrow{p}}_k$, where k=1,2,...,K. As mentioned earlier, a model in this case refers to a vector of different models and sub-models as presented in Table 1 while the parameter set denotes a vector of parameters to these models (see Table 2). As previously mentioned, we assume that all models are of equal importance a priori and therefore each model is characterized by a uniform prior distribution ($P({\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)$). As carried out also in Ref. [15], we assume that we have no prior knowledge on the model parameters and hence, the parameters ${\overrightarrow{p}}_k$ were also drawn from a uniform distribution. Now, if $L({\overrightarrow{\sigma }}_E|{\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)$ is our likelihood function for model ${\overrightarrow{M}}_j$ and parameter set ${\overrightarrow{p}}_k$, the likelihood function can be given within the Bayesian Monte Carlo [4] and Unified Monte Carlo (UMCB-G and UMC-B) [8, 9] approaches as: $L({\overrightarrow{\sigma }}_E|{\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)\propto \exp \left(-\frac{{\chi }_{G(k\mathrm{,}j)}^2}{2}\right)\mathrm{,}$ (1) where ${\chi }_{G(k\mathrm{,}j)}$ is the global chi square given in Eq. 5. Given that $P({\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)$ is our combined prior distribution of models and their parameters, and $L({\overrightarrow{\sigma }}_E|{\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)$ is the likelihood function given in Eq. 1, we can now compute our posterior distribution ($P({\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k|{\overrightarrow{\sigma }}_E)$) as follows: $P({\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k|{\overrightarrow{\sigma }}_E)=\frac{L({\overrightarrow{\sigma }}_E|{\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)P({\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)}{P({\overrightarrow{\sigma }}_E)}\mathrm{,}$ (2) where $P({\overrightarrow{\sigma }}_E)$, which is the marginal likelihood also referred to as the model evidence, is simply a normalisation constant and therefore not considered in the optimization. Eq. 2 becomes: $P({\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k|{\overrightarrow{\sigma }}_E)=L({\overrightarrow{\sigma }}_E|{\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)P({\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)\mathrm{.}$ (3) Based on Eq. 1, we assign each random nuclear data file with a weight equal to the likelihood function, also called Bayesian Monte Carlo (BMC) weights (see Ref. [4]): $w_{k\mathrm{,}j}=\exp \left(-\frac{{\chi }_{G(k\mathrm{,}j)}^2}{2}\right)\mathrm{.}$ (4) From Eqs. 1 and 4, ${\overrightarrow{\sigma }}_E$ is our experimental data, and ${\chi }_{G(k\mathrm{,}j)}^2$ and wk,j respectively are the global reduced chi square and BMC weights for parameter set k and model combination j. ${\chi }_{G(k\mathrm{,}j)}^2$ is our single but multiple criteria objective function obtained as a linear combination of the individual reduced χ² computed for each considered experimental data type (β) and given by: ${\chi }_{G(k\mathrm{,}j)}^2=\frac{1}{N_{\beta }}{\displaystyle \sum _{\beta =1}^{N_{\beta }}}{\chi }_{G(k\mathrm{,}j)}^2(\beta )\mathrm{,}$ (5) where ${\chi }_{(k\mathrm{,}j)}^2(\beta )$ is a vector of reduced chi square values computed using the experimental data type β and Nβ is the total number of considered experimental data types. For Nβ experimental data types, Eq. 5 can be given as: ${\chi }_{G(\text{k}\mathrm{,}j)}^2=\frac{{\chi }_k^2(\text{x}s)+{\chi }_k^2(\text{r}p)+{\chi }_k^2(\text{D}A)}{N_{\beta }}\mathrm{,}$ (6) where Nβ is equal to three in this case, ${\chi }_k^2(\text{x}s)$ and ${\chi }_k^2(\text{r}p)$ are the reduced chi squares computed for the reaction and the residual production cross sections respectively, and ${\chi }_k^2(\text{D}A)$ is the reduced chi square computed for the elastic angular distributions. Eq. 6 can be written as: ${\chi }_{\text{G}(k\mathrm{,}j)}^2=w_{{\beta }_1}^E{\chi }_k^2(\text{x}s)+w_{{\beta }_2}^E{\chi }_k^2(\text{r}p)+w_{{\beta }_3}^E{\chi }_k^2(\text{D}A)\mathrm{,}$ (7) where the $w_{{\beta }_1}^E$, $w_{{\beta }_2}^E$ and $w_{{\beta }_3}^E$ are the coefficients of the linear combination presented in Eq. 7 and represents the weighting factors for the experimental data types: (1) reaction cross sections, (2) residual production cross sections, and (3) elastic angular distribution respectively. In Eq. 7, the experimental data types are equally weighted, implying that equal weighting factors were assigned to each: $w_{{\beta }_1}^E=w_{{\beta }_2}^E=w_{{\beta }_3}^E=1/N_{\beta }$. This is normally the case for a general purpose nuclear data library evaluation where all experimental data types are assigned equal importance. Alternatively, the weighting factors for each experimental data type can be determined by the evaluator depending on the needs of the evaluation. Since the weighting factors are normalized, they must sum up to 1: ${\displaystyle \sum _{\beta =1}^{N_{\beta }}}{w}_{\beta }^E=1.$ (8) Eq. 5 was computed assuming that there were no correlations between the different experimental data types considered. We note that this assumption is simplistic since in some cases, similar or the same instruments, methods, and authors were involved in the experiments and measurements of more than one experimental data type, which could introduce cross correlations between these experimental data types. For example, F. Ditroi was involved in the measurement and analysis of the reaction cross section, ⁵⁹Co(p,3n), between 16.0–69.8 MeV as well as in the measurement of the residual production cross section, ⁵⁹Co(p,x)⁵⁶Co, between 15.0–69.8 MeV. However, these cross correlations are not readily available and therefore not used in this work.

Similarly, in the computation of the individual reduced chi squares in Eq. 6 such as the ${\chi }_k^2(\text{x}s)$, the experimental data sets were assumed to be uncorrelated. To include experimental correlations, the generalized chi square as presented in Ref. [4] should have been used. However, experimental correlations are scarce especially with respect to proton induced reactions and when available, are usually incomplete. Therefore, the reduced chi square with respect to the reaction cross sections (${\chi }_k^2(\text{x}s)$) for example, was computed using Eq. 9 (similar expressions have been presented also in Refs. [4] and [17]): ${\chi }_k^2(xs)=\frac{1}{N_c}{\displaystyle \sum _{c=1}^{N_c}}\frac{1}{N_m}{\displaystyle \sum _{m=1}^{N_m}}\frac{1}{N_{pt}}{\displaystyle \sum _{i=1}^{N_{p}t}}{\left(\frac{{\sigma }_{T(k)}^{ci}-{\sigma }_E^{cmi}}{\Delta {\sigma }_E^{cmi}}\right)}^2\mathrm{,}$ (9) Where Nc is the total number of considered channels c, Nm is the total number of experimental data sets m, and Npt is the total number of considered data points for each experimental data set; σ_Ecmi and ${\sigma }_{T(k)}^{ci}$ are respectively, the vectors of the experimental and TALYS calculated cross sections at the energy i, for the data set m, and channel c. Similarly, $\Delta {\sigma }_E^{cmi}$ is the corresponding experimental uncertainty at energy i, data set m and channel c. In cases where there were no matches in energy (or in angle in the case of the elastic angular distributions) between the TALYS calculations and the considered experiments, similar to what was carried out in Refs. [15, 17], we linearly interpolate to fill in the missing TALYS values. The same approach as presented in Eq. 9 was applied for the computation of ${\chi }_k^2(\text{r}p)$ and ${\chi }_k^2(\text{D}A)$, however, in the case of the ${\chi }_k^2(\text{D}A)$, we match TALYS calculations with that of the experiments in both angle and in energy but we only interpolate on the angle.

In Refs. [15, 17], the reduced chi square presented in Eq. 9, was computed by averaging the chi square values for each channel over all the considered experimental data points. This approach, as stated also in Ref. [4], assigns equal weights (aside their uncertainties) to all the experimental data points which may lead to a situation where an experimental data set with a large number of measurements completely dominates the goodness of fit estimations. Also, channels with many different experimental data sets but fewer measurements would be assigned smaller weights compared to those with fewer experimental data sets but with many measurements. In an attempt to assign channels with many different experimental data sets with larger weights (aside their uncertainties) in line with what was carried out in Ref. [4], we averaged the chi square over each experimental data set by dividing by Npt as presented in Eq. 9. In this way, each experimental data set contributes equally to the goodness of fit estimation. Further, since the measurements from each data set are known to be highly correlated, by averaging over each experimental data set, we effectively combine the information from each experimental set into a single goodness of fit estimate. The correlations come about as a result of the fact that usually, the same equipment as well as methods (and authors) were used or involved in these measurements. It is also known that the addition of correlated experiments is not as effective in reducing the uncertainty in our calibration as the uncorrelated or independent experiments [37]. Experimental data set as used here refers to one or more measurements carried out at a specific energy or energy range, for a particular channel and isotope, with a unique EXFOR ID in the EXFOR database.

Statistical information of the posterior distribution in Bayesian estimations as presented in Eqs. 2 and 3, can normally be summarized by computing central tendency statistics [38], where the posterior mean is used as the best estimate (with its corresponding variance). However, given a large sample size and assuming that our prior was sampled from a uniform distribution, as stated in Ref. [39], the posterior Probability density function (PDF) can be asymptotically approximated by a Gaussian PDF centered on the Maximum a Posteriori (MAP) estimate. Therefore, in this work, as used also in Bayesian Model Selection (BMS) [40], the best or winning model becomes the model (and parameter) set that maximizes the posterior probability which is also known as the Maximum a Posteriori (MAP) estimate (i.e the mode of the posterior distribution), given as: $L_{\text{MAP}}=*\underset{m}{argmax}\left[L({\overrightarrow{\sigma }}_E|{\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)P({\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)\right]$ (10) where L_MAP is the Maximum a Posteriori (MAP) estimate and the index m denotes the considered models with their parameters. However, in Bayesian statistics, as stated in Ref. [38], the maximum likelihood estimation (MLE) "is a special case of Bayesian Estimation (BE) in which (1) the estimate is based on the mode of the posterior distribution, and (2) all the parameters values are equally likely". The MLE is therefore viewed as a special case of the MAP estimate in the case where a uniform distribution is assumed for the prior distribution of the parameters [38]. Since in this work, we assume a uniform distribution for the models as well as for their parameters, the MLE was used. Therefore, given a uniform distribution of models and their parameters, Eq. 10 reduces to the maximum likelihood estimate denoted by L_MLE and given as: $L_{\text{MLE}}=*\underset{m}{argmax}\left[L({\overrightarrow{\sigma }}_E|{\overrightarrow{M}}_j\mathrm{,}{\overrightarrow{p}}_k)\right]$ (11) The L_MAP and L_MLE estimates were computed and compared for selected model parameters in the case of p+⁵⁹Co and found to be the same (see Table 8. Therefore, from Eq. 11, the model combination (with its parameter set) which maximizes the likelihood function was considered as the winning model set and the file that makes the experimental data most probable.

Iterative Bayes Procedure

The algorithm for the Iterative Bayesian procedure proposed in this work is presented in Table 3. The idea of the iterative procedure is to minimize the bias between our experimental observables and the corresponding model outputs in an iterative fashion. As can be seen from the table, we start with the selection of model combinations and their parameters (including parameter uncertainties and distributions) - see Table 1 for more information on the considered models. Next, we select the energy grid for the TALYS calculations. The energy grid was chosen such that there was a large number of matches in incident energy between the TALYS calculations and that of the corresponding experimental data. As mentioned earlier, we linearly interpolate in energy in the case of the reaction and residual production cross sections, and in angle, in the case of the elastic angular distributions, for the purpose of filling in the missing TALYS values.

Next, a set of random combinations of the models were generated from uniform distributions. In addition, the model parameters to each model combination were also sampled from uniform distributions to create a large set of TALYS input files. These input files were then run with the TALYS code system to produce a set of random nuclear data files in the ENDF format. These initial random nuclear data files constitute the parent generation (Gen. 0). The ENDF files produced were processed into x–y tables for the purpose of comparison with experimental data from the EXFOR database as well as with other nuclear data libraries, which were obtained in the ENDF format. The selection of experimental data for the considered channels has been presented earlier in Sect. 2.1.

Next, we select the ‘best’ model combination (with the parameter set) by identifying the random nuclear data file with the largest likelihood function value. Using this best file as our new central file or ‘best’ estimate, we re-sample around this file, however, this time only varying the model parameters of the selected model to produce another set of random nuclear data files referred to as the 1st generation (Gen. 1). The idea here is that, once we are able to identify the ‘best’ model combination, we proceed with subsequent iterations as though the models selected were the true models. To proceed with the many different models for each iteration would imply that no model combination is good enough and hence, an average should instead, be taken over all or a selection of the models. This approach has been proposed and presented in a dedicated paper [22].

Next, is to update the model parameters and their uncertainties. However, before updating the parameter uncertainties, we first compute the Effective Sample Size (ESS). This was done in order to ensure that a sufficient number of effective samples are available for the computation of the posterior uncertainties and covariances. The effective sample size (ESS) is given as the squared sum of the file weights divided by the summed squares of these weights: $ESS=\frac{{\left({\displaystyle {\sum }_{k=1}^n{w}_k}\right)}^2}{{\displaystyle {\sum }_{k=1}^n{w}_k^2}}$ (12) The ESS is a useful metric for understanding the information value of each random nuclear data file and therefore gives an indication on how many random samples have significant impact on the posterior distribution. It has been noticed in this work that with an ESS ≤ 10%, very few random nuclear data files have significant impact on the posterior distribution and hence, as a rule of the thumb, we only update the parameter uncertainties for an ESS ≥ 10%. We then repeat our model parameter variation step (steps 3 through to 14 in Table 3) in an iterative fashion until we reach convergence.

Convergence was determined by monitoring the evolution of the maximum likelihood estimate computed for each iteration. If $L_{\text{MLE}}^{(n)}$ represents the maximum likelihood estimate (MLE) for the last (n) iteration and $L_{\text{MLE}}^{(n-1)}$, the MLE value of the last but one iteration, then ΔL_MLE, the relative change in the likelihood function can be given as: $\Delta L_{\text{(MLE)}}=\frac{L_{\text{MLE}}^{(n)}-L_{\text{MLE}}^{(n-1)}}{L_{\text{LME}}^{(n)}}⩽{\epsilon }_G$ (13) The iteration is said to have converged when the relative change in the maximum likelihood estimates (ΔL_MLE) for the last 2 iterations is within a target global tolerance denoted by εG. In this work, we chose to set the value of εG to: εG ≤ 0.05 (5%). It should be noted however that this value was chosen arbitrarily. Ideally, a smaller tolerance value (εG) is preferred however, small values of εG can be difficult to achieve for a number of reasons: (1) because of the limitation of the models used, (2) the condition of Pareto optimality in multi-objective optimization, (3) computational resource constraint, and (4) when target accuracy is reached. In the case of (1), since the fit can only be improved with parameter variations within the limits of the selected models, we are constrained by the deficiencies of the winning model set. For example, if the winning model set is unable to reproduce the shape of experimental data for a particular channel, only varying the parameters would not be sufficient for improving the fits to experimental data. As mentioned earlier, the effects of model defects would then have to be taken into account in these cases.

With (2), a point worthy of note is that within the limits of our models in a multi-objective optimization procedure (as in this case), a point is reached in the iteration process, where further variation of the parameters is not able to improve the fits to the reaction cross sections for example, without making the fits to the residual production cross sections or elastic angular distributions worse-off. This condition is referred to as Pareto optimality [41]. This is because, multi-objective optimization problems as in our case, involves the simultaneous optimization of multiple competing objectives, which most often, leads to trade-off solutions largely known as Pareto-optimal solutions. A possible solution would be to assign subjective weights to each criteria or the chi square computed for each experimental data type as presented in Eq. 7 depending on the needs of the evaluation. For example, if the target of the evaluation is the production of radio-isotopes, relatively larger weights could be assigned to the residual production cross sections in the computation of the global chi square. This would place more importance on the residual cross sections in the optimization procedure. Since the goal of this work was however to provide a general purpose evaluation which should compare reasonably well with experimental data for all or a large number of the cross sections, as well as the residual production cross section and angular distributions, equal weights were assigned to each individual experimental data type. As a general rule of the thumb, when no further improvements are possible with the global likelihood function, the iteration should be stopped since the Pareto optimum might have been reached.

For (3), the creation of random nuclear data files can be computationally expensive depending on the energy grid and the TALYS excitation bin size used. For example, more than 6,000 random nuclear data files were produced and used in this work. This translated into several months of computational time. Therefore, the value of εG should be chosen taken into consideration the accuracy gained as well as the computational resources available. In the case of (4), the value of εG can be chosen with a particular objective or target accuracy in mind. For example, one underlying objective of this work has been to create an evaluation that globally outperforms and significantly improves the TENDL-2019 evaluation for the p+⁵⁹Co. Hence, the iteration can be stopped after this objective is achieved. To check if the final results were robust to a change in εG, the adjusted results were compared with experimental data for different changes in εG: 15%, 10% and 5% and it was observed in the case of p+⁵⁹Co and p+¹¹¹Cd that the pareto optimum was reached at about 5% for the reaction cross sections (see Table 5 for example). We note here that only static or energy-independent model parameters were used for the entire energy range under consideration, hence, this does not offer enough flexibility in the adjustment or fitting procedure. The use of energy-dependent parameters are proposed for future work.

Next, we return the final solution for Bayesian calibration i.e. the final ‘best’ file. This file was then compared against available microscopic experimental data from the EXFOR database as well as with evaluations from other nuclear data libraries (if available). For testing and validation of the evaluation, we also compare our evaluation with experimental data from the channels that were not used in the optimization. This was done so as not to use the same experimental data for both adjustment and validation. Where available, the ‘best’ file should also be taken through integral validation where the evaluation is tested against a large set of integral benchmark experiments. These benchmarks are however, not readily available for proton induced reactions and therefore, not used in this work. Finally, we update also, the final parameter uncertainties (see Sect. 2.5).

Updating model parameter uncertainties

Once we are satisfied with the ‘best’ file, we move to the next step where we infer parameter uncertainties and covariance information to this file. We note here that in Bayesian inference, each parameter is described by a probability distribution instead of a point estimate. Therefore, by varying the model parameters around this file, we ensure that the prior parameter distribution is centered around the ‘best’ file. The covariance information associated with this file is then contained in the distribution of the random nuclear data files produced. Alternatively, as done in most nuclear data libraries, the covariance information associated with the evaluation can be stored in MF31-40 (in ENDF terminology).

Under appropriate regularity conditions, the posterior can be approximated to a normal distribution with a mean and standard deviation [42]. Therefore, if we assume that the parameter posterior distribution is Gaussian, a weighted variance $(va{r}_w({\overrightarrow{p}}_{\iota }))$ for parameter ($p_{\iota }$) can be computed as follows: $va{r}_w({\overrightarrow{p}}_{\iota })=\frac{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k\mathrm{.}{\overrightarrow{p}}_{\iota k}{}^2}}{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}}-{\overline{p}}_{\iota }{}_w^2\mathrm{,}$ (14) where wk and pιk are the file weight and value of the parameter ($p_{\iota }$) of the random file k respectively, ${\overline{p}}_{\iota }{}_w$ is the weighted average value of the parameter, $p_{\iota }$, which can given by: ${\overline{p}}_{\iota }{}_w=\frac{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k\mathrm{.}{\overrightarrow{p}}_{\iota k}}}{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}}\mathrm{.}$ (15) As mentioned earlier, given a large sample size and also that, the parameters are sampled from uniform distributions, the MLE is approximately equal or close to the mean values of each parameter. It is shown later in Table 8 that the mean values of the posterior parameter distributions were closed to the MLE estimates with respect to the optical model parameters as expected. Similarly, we can compute a weighted parameter covariance matrix between parameters $p_{\iota }$ and pτ as follows: $co{v}_w(p_{\iota }\mathrm{,}{p}_{\tau })=\frac{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}\left({\overrightarrow{p}}_{\iota k}-{\overline{p}}_{\iota }\right)\left({\overrightarrow{p}}_{\tau k}-{\overline{p}}_{\tau }\right)}{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}}\mathrm{,}$ (16) where ${\overline{p}}_{\iota }$ and ${\overline{p}}_{\tau }$ are the mean values of $p_{\iota }$ and pτ parameters respectively, pτ k is the value of the parameter (pτ) in the random file k and covw(pι,pτ) is the weighted covariance matrix between parameters $p_{\iota }$ and pτ.

The use of the global likelihood function used in this work for updating the parameter uncertainties as presented in Eq. 4 raises a number of concerns:

1. Parameter sensitivities were not explicitly taken into account in the computation of the the global reduced chi square as presented in Eq. 6. This is particularly important since each parameter would normally have a different sensitivity to the different cross sections or the angular distributions. This could lead to a situation where the weighted averaged parameters computed from the posterior distribution would not necessarily provide the best central values. In these cases, the parameter mean and the MLE values would be different. A possible solution would be to include parameter sensitivities as channel weights in the computation of χ²(xs), χ²(rp) and χ²(DA). In this way, if a particular channel is only sensitive to 2 or 3 parameters, this information would be included in the computation of the global reduced chi square. The inclusion of parameter sensitivities is however beyond the scope of this work.

2. The use of energy - independent parameters as used in this work does not give enough flexibility to the adjustment procedure. It should be noted here, as mentioned also in Ref. [4], that the predictive power of TALYS is energy dependent and therefore, the use of energy-dependent parameters in the computation of weights would give more flexibility to the adjustment procedure as parameters can be mapped to the weights in energy and/or in angle, thereby improving the ability to better constrain the model parameters to experimental data. An example on the use of energy-dependent parameters in the treatment of model defects for nuclear data evaluation is presented in Ref. [43].

3. Experimental correlations were not included in the optimization. Since some of the experimental errors are known to be correlated, these correlations should have been included by using the generalized chi square presented in Ref. [4]. However, since these correlations were not readily available, no correlations were used.

Updating the cross sections

Similar to section 2.5, the weighted mean of the cross sections which corresponds to (or must be close to the values of the ‘best’ file), can be given as: ${\overline{{\sigma }_{\text{T}(k)}^c}}_w=\frac{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k\mathrm{.}{\sigma }_{\text{T}(k)}^c}}{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}}\mathrm{,}$ (17) where ${\overline{{\sigma }_{\text{T}(k)}^c}}_w$ is the weighted mean of the TALYS (T) calculated cross sections taken from K random files and ${\sigma }_{\text{T}(k)}^c$ is a vector of TALYS calculated cross sections of interest. w denotes weighted and wk represents the weights as a function of random files. The corresponding weighted variance ($va{r}_w({\overline{{\sigma }_{\text{T}(k)}^c}}_w)$), similar to Eq. 14 can be expressed as: $va{r}_w({\overline{{\sigma }_{\text{T}(k)}^c}}_w)=\frac{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k\mathrm{.}\overrightarrow{{\sigma }_{\text{T}(k)}^c{}^2}}}{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}}-{\overline{{\sigma }_{\text{T}(k)}^c}}^2\mathrm{.}$ (18) Similarly, a weighted or the posterior covariance matrix between cross sections at energy a (${\sigma }_{{\text{T}}_a}^c$) and b (${\sigma }_{{\text{T}}_b}^c$) can be given as: $co{v}_w({\sigma }_{{\text{T}}_a}^c\mathrm{,}{\sigma }_{{\text{T}}_b}^c)=\frac{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}\left(\overrightarrow{{\sigma }_{{\text{T}}_{a(k)}}^c}-\overline{{\sigma }_{{\text{T}}_a}^c}\right)\left(\overrightarrow{{\sigma }_{{\text{T}}_{b(k)}}^c}-\overline{{\sigma }_{{\text{T}}_b}^c}\right)}{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}}\mathrm{.}$ (19) The sample weighted correlation coefficient (rw) can be given as [44]: $r_w=\frac{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}\left(\overrightarrow{{\sigma }_{{\text{T}}_{a(k)}}^c}-\overline{{\sigma }_{{\text{T}}_a}^c}\right)\left(\overrightarrow{{\sigma }_{{\text{T}}_{b(k)}}^c}-\overline{{\sigma }_{{\text{T}}_b}^c}\right)}{\sqrt{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}{\left(\overrightarrow{{\sigma }_{{\text{T}}_{a(k)}}^c}-\overline{{\sigma }_{{\text{T}}_a}^c}\right)}^2}{\sqrt{{\displaystyle {\sum }_{k=1}^K{\overrightarrow{w}}_k}\left(\overrightarrow{{\sigma }_{{\text{T}}_{b(k)}}^c}-\overline{{\sigma }_{{\text{T}}_b}^c}\right)}}^2}\mathrm{.}$ (20) It must be stated here that even though these covariances and correlation matrices were produced in this work, the same information is contained in the large set of random cross sections produced. These random nuclear data files can be processed and used in a Total Monte Carlo approach for uncertainty propagation to applications [3, 15, 13].

Final parameter uncertainties

To show the convergence of the posterior parameter distributions of the 5th generation, we present a plot example of the convergence of the posterior mean and standard deviation for the vso2adjust ($v_{so2}^p$) parameter in Fig. 10. It can be seen from the figure that even though over a thousand samples were produced, the mean and the standard deviation converges well after about 600 runs or samples. The fluctuations observed after 600 samples were found to be below 1%. Similar plots (not shown) have been produced for the other model parameters and it was observed that a large number of the parameters converged after about 800 samples. For uncertainty propagation to applications purposes, within the Total Monte Carlo approach, it has been observed earlier that convergence could be reached after 300 samples [60].

Fig. 10Full-screen View

(Color online) Plot example showing the convergence of the posterior mean and the standard deviation of the vso2adjust optical model parameter distribution (for p + ⁵⁹Co). The parameter distribution used here was obtained by re-sampling the parameters around the ‘best’ file values obtained from the 5th generation.

In Fig. 11, a distribution of the file weights computed for p + ⁵⁹Co using Eq. 4, is presented. As expected, a large number of files were assigned with low or insignificant file weights and therefore contributed less to the final (or posterior) uncertainties. This is consistent with earlier results presented with respect to neutron induced reactions in Refs. [15, 10]. The low and insignificant weights obtained could be attributed to the following: (1) the presence of outlier experiments as well as experiments with unreported or under-reported uncertainties, (2) sampling from a region of the parameter space where the likelihood is probably low, (3) the difficulty in optimizing our fits simultaneously to three different experimental data types considered, (4) the presence of model defects, and (5) not considering experimental correlations in the computation of the χ². In the case of (1), data considered as outlier were assigned with a binary value of zero and therefore not considered in the optimization. Even though outlier experiments were discarded, it was still observed that our models had difficulty in reproducing experimental data especially at the threshold energies. It is observed in Ref. [4] that, the nuclear model calculations can easily be several factors off experimental data at the threshold energies as well as at the high-energy tail of the excitation function. Our solution here has been to exclude a few data points at the threshold regions for cross sections. Furthermore, instead of rejecting outlier experiments, it is proposed that approaches for the automatic correction of the uncertainties of inconsistent experiments and the identification of Unrecognized Sources of Uncertainties (USU) in experiments proposed in Refs. [28] and [29] respectively, be incorporated into the iBMC methodology in future work.

Fig. 11Full-screen View

(Color online) Distribution of global file weights for p + ⁵⁹Co computed using Eq. 4. The weights have been normalized with the maximum weight in order to relate them to 1.

To address (2), an accept/reject method could be used to reject samples with insignificant file weights as carried out in Ref. [15]. This approach could however introduce some bias depending on the cut-off parameter. Therefore, instead of rejecting samples, we instead, compute an Effective Sample Size (ESS) using Eq. 12 in order to ensure that a sufficient number of effective samples were available for the computation of the final posterior moments. In the case of (3), it was observed that optimising our model calculations to three experimental data types simultaneously can be difficult. A solution would be to optimise each experimental data type or cross section individually and then combine them into an evaluation. This would however, lead to inconsistent evaluations as the sum rules must be obeyed. With respect to (5), experimental correlations and covariances are often not readily available and therefore were not considered in this work. The inclusion of these correlations is however recommended for future work.

Table 8 presents a comparison between the posterior mean and the Maximum Likelihood Estimate (MLE) values for a selected number of model parameters in the case of the 5th generation (p + ⁵⁹Co). The corresponding updated uncertainties obtained from each posterior parameter distribution are also presented. The mean values were obtained by taking the weighted average of the posterior distribution of each parameter using the computed file weights while the MLE values are the parameter values of the TALYS input file with the maximum likelihood estimate (i.e. the final ‘best’ file). It was observed in this work that the MAP estimate was the same as the MLE values as expected and therefore, not included in Table 8.

Similar to the prior parameter uncertainty values given in Table 2, the updated (posterior) uncertainties are given as a fraction (%) of the posterior mean values. It can be seen from the table that the MLE and the posterior mean values are close for most of the optical model parameters presented as expected - a relative difference of less than 1% was obtained for most of the optical parameters. This is not surprising, as mentioned earlier, under appropriate regularity conditions, the posterior distribution should be centered at or close to the MLE values. However, large relative change between the MLE and the mean values of more than 10% were recorded for example, in the case of the spin cut-off parameter (e.g. σ²(⁵⁶Fe)) of the level density model and the single-particle state densities (g_π(⁵⁶Fe)). For these parameters, it was observed that the prior parameter distributions were not flat but were rather skewed towards the right. This could be attributed to the inability of our algorithm to sample the entire space of these parameters and would therefore require a larger number of samples.

From the table, a relatively significant reduction in the prior uncertainty is observed for a number of parameters. For example, a prior uncertainty of 5% was reduced to 3.5% for the avdadjust parameter and a 20% prior uncertainty was reduced to 8.2% and 7.4% for the wso1adjust and wso2adjust parameters respectively. Also, a prior uncertainty of 10% and 2% were reduced to 4.5% and 1.3% respectively for the rcadjust and the rvadjust parameters. Small or negligible reductions were however observed for parameters such as v1adjust: (from 2% to 1.8%), and in the case of v2adjust and v3adjust parameters, from 3% to 2.7%. For the d1adjust, d2adjust and d3adjust parameters, the prior uncertainties were reduced from 10% for each parameter, to 7.6%, 7.7% and 8.4% respectively. A large reduction in posterior spread (uncertainties) gives an indication that the parameter under consideration has a large sensitivity to the channels used in computing the global likelihood function. On the other hand, small or no reduction in the posterior uncertainty implies that the file weights do not have impact on the posterior parameter distribution of interest, signifying no sensitivity of the parameter under consideration to the channels used in the computation of the weights. This is because, by using a uniform or relatively flat priors (as was in this case), the posterior distribution is determined solely (or almost) by the experimental data through the likelihood function.

In Fig. 12, prior and posterior distributions for selected optical model parameters in the case of the 5th generation: avadjust ($a_V^p$), avdadjust ($a_D^p$), v1adjust ($v_1^p$), w1adjust ($w_1^p$), are presented. As can be seen from the figure, the prior samples were drawn from uniform parameter distributions. The posterior distributions describes the updated information for the considered parameters after experimental data were taken into account. In general, it can be seen that the priors are relatively flat as expected. However, the posterior distribution peaks sharply at the mean value which incidentally is very close to the MLE and the MAP values as expected. The posterior (or updated) uncertainties for the parameters in Fig. 12 are as follows: $a_V^p$ - 1.8% (reduced from 2%), $a_D^p$ - 3.5% (reduced from 4%), $v_1^p$ - 1.8% (reduced from 2%), $r_V^p$ - 1.3% (reduced from 2%). An Effective Sample Size (ESS) of 110 which represents 10.8% of the total sample size, was obtained. Our inability to significantly reduced the posterior uncertainties for some of the parameters can be attributed to the insensitivity of these parameters to the cross sections and the elastic angular distributions used in computing the file weights.

Fig. 12Full-screen View

(Color online) Prior and posterior parameter distributions for selected optical model parameters: avadjust ($a_V^p$), avdadjust ($a_D^p$), v1adjust ($v_1^p$), rvadjust ($r_V^p$) in the case of the p + ⁵⁹Co (Gen. 5). An Effective Sample Size (ESS) of 110 (representing 10.8% of the total sample size) was obtained.

As a natural consequence of the Monte Carlo method used, correlations and other statistical information such as mean and variances can be extracted from both the prior and posterior distributions. In Fig. 13, prior and posterior (weighted) correlation matrices for a selected number of optical model parameters are presented. The correlation matrices presented are symmetric and give a measure of the strength of the linear dependence between two parameters and vary from -1 (perfect negative correlation) to +1 (perfect positive correlation). It can be seen from the figure that the prior correlation values are generally relatively low. This is expected since the parameters were sampled in an uncorrelated manner. Additional correlations were then introduced through the use of the same experimental data used in the computation of the file weights as can be seen in the weighted (posterior) correlation matrix in Fig.13. A negative correlation was observed between the rvadjust ($r_V^p$) and v1adjust ($v_1^p$) parameters for example. This is expected as there is a well known inverse relationship between the $r_V^p$ and $v_1^p$ parameters which can be given as follows: $v_1^p*r_V^p=constant$ [61]. It must also be noted that most of the optical models showed similar behavior as can be seen in Fig. 13. Similar negative correlations were also observed in Ref. [61]. We note here that these parameter correlations were however, automatically taken into account in the simultaneous variations of the model parameters. Similarly, energy-energy correlations for each cross section can be obtained. As proof of principle, evaluated posterior correlation matrices for selected proton induced reaction cross sections ((p,nonel), (p,n), Co(p,γ) and (p,4n)) of p + ⁵⁹Co are presented in Fig. 14. Similar matrices (not shown) were obtained for the residual production cross sections and the elastic angular distributions. These matrices were obtained by combining the global file weights with the cross sections. The high correlations as seen in the figure can be attributed largely to the energy-energy correlations that come from the use of the same theoretical models in the production of the random cross sections. A similar observation was made in Ref. [60]. Additional correlations come from the use of the same experimental data in the computation of the file weights. This however, has less impact since the correlations were introduced through the global likelihood function which also took into account, experimental residual cross section data as well as the elastic angular distributions. These weighted correlations and covariance information can be used for uncertainty propagation to applications. It must be stated here that even though these parameter correlations and covariances are available, the final cross sections were produced from the direct variation of model parameters and not from perturbing the cross sections.

Fig. 13Full-screen View

(Color online) Prior and posterior (weighted) correlation matrices for a selected number of optical model parameters for p + ⁵⁹Co (Gen. 5).

Fig. 14Full-screen View

(Color online) Evaluated posterior correlation matrix for selected proton induced cross sections of p + ⁵⁹Co (Gen. 5). The incident energy used corresponds to the incident energies of the experimental data used in the computation of the file weights.