Bayesian updating for data adjustments and multi-level uncertainty propagation within Total Monte Carlo

TL;DR

This paper introduces a Bayesian method to combine differential and integral nuclear data within the Total Monte Carlo framework, improving data adjustments and uncertainty propagation for nuclear reactions.

Contribution

It develops a novel Bayesian approach to integrate multiple data types for nuclear data adjustment using Total Monte Carlo, enhancing accuracy and uncertainty quantification.

Findings

The method successfully adjusted $^{208}$Pb nuclear data below 20 MeV.

The adjusted data compared favorably with experimental and evaluated data.

The approach improves the reliability of nuclear data for applications.

Abstract

In this work, a method is proposed for combining differential and integral benchmark experimental data within a Bayesian framework for nuclear data adjustments and multi-level uncertainty propagation using the Total Monte Carlo method. First, input parameters to basic nuclear physics models implemented within the state of the art nuclear reactions code, TALYS, were sampled from uniform distributions and randomly varied to produce a large set of random nuclear data files. Next, a probabilistic data assimilation was carried out by computing the likelihood function for each random nuclear data file based first on only differential experimental data (1st update) and then on integral benchmark data (2nd update). The individual likelihood functions from the two updates were then combined into a global likelihood function which was used for the selection of the final 'best' file. The proposed method has been applied for the adjustment of $^{208}$Pb in the fast neutron energy region below 20 MeV. The 'best' file from the adjustments was compared with available experimental data from the EXFOR database as well as evaluations from the major nuclear data libraries and found to compare favourably.

Tables

Table 1. Table 1: Selected nuclear model parameters of the TALYS code (with their parameter widths) used for parameter variations. The parameter widths are given as a fraction (%) of their absolute values. The parameter widths of a 𝑎 a , g π subscript 𝑔 𝜋 g_{\pi} and g ν subscript 𝑔 𝜈 g_{\nu} are given in terms of the mass number A 𝐴 A . a 𝑎 a is the real central diffuseness, g π subscript 𝑔 𝜋 g_{\pi} and g ν subscript 𝑔 𝜈 g_{\nu} are the single-particle state densities as used in exciton model analyses  [ 35 ] . A complete list of all the model parameters can be found in Ref.  [ 36 ] . All parameters were varied altogether within their parameter widths in a TMC fashion.

Parameter	Uncertainty(%)	Parameter	Uncertainty(%)
$r_V^n$	1.5	$a_V^n$	2.0
$v_1^n$	1.9	$v_2^n$	3.0
$v_3^n$	3.1	$v_4^n$	5.0
$w_1^n$	9.7	$w_2^n$	10.0
$d_1^n$	9.4	$d_2^n$	10.0
$d_3^n$	9.4	$r_D^n$	3.5
$a_D^n$	4.0	$r_{SO}^n$	9.7
$a_{SO}^n$	10.0	$v_{SO1}^n$	5.0
$v_{SO2}^n$	10.0	$w_{SO1}^n$	20.0
$w_{SO2}^n$	20.0	${\Gamma }_{\gamma }$	5.0
$a{(}^{207}Pb)$	4.5	$a{(}^{206}Pb)$	6.5
$a{(}^{208}Pb)$	5.0	$a{(}^{205}Pb)$	6.5
${\sigma }^2$	19.0	$M^2$	21.0
$g_{\pi }$(${}_{}{}^{207}Pb$)	6.5	$g_{\nu }$(${}_{}{}^{207}Pb$)	6.5

Table 2. Table 2: Differential experimental data used for the adjustment of 208 Pb showing the number of data points per reaction channel used. Only data points between 5 and 20 MeV were used in the computation of file weights. Also, only the names of the first authors have been presented.

\pbox20cmCross
section	\pbox20cm MT
number	\pbox20cmData points used
5 - 20 MeV	Author	EXFOR ID	Year
(n,tot)	001	142	R.F. Carlton	13735002	1991
(n,tot)	001	116	D.G. Foster Jr	10047092	1971
(n,non-el)	003	1	G.M.Haas	11794005	1963
(n,inl)	004	46	L.C. Mihailescu	23039003	2008
(n,2n)	016	16	Frehaut	20416058	1980
(n,$\gamma$)	102	8	I. Bergqvist	10226002	1972
(n,$\gamma$)	102	7	J. Csikai	30074007	1967
(n,$\gamma$)	102	1	D. Drake	10193006	1971

Table 3. Table 3: The integral experiments used in this work showing the evaluated experimental benchmark k eff subscript 𝑘 eff k_{\rm eff} and corresponding experimental uncertainties. These benchmarks were obtained from the (ICSBEP) handbook  [ 42 ] . HMF stands for Highly Enriched Uranium (HEU) Metallic Fast, PMF for Plutonium Metallic Fast and LCT for Low Enriched Uranium (LEU) Compound Thermal benchmarks.

Benchmark category	Abbreviation	\pbox20cmExperimental
Benchmark $k_{\mathrm{eff}}$	\pbox20cmExperimental Benchmark
uncertainty (pcm)
PMF035 case 1	pmf35c1	1.0000	160
HMF027 case 1	hmf27c1	1.0000	250
HMF057 case 1	hmf57c1	1.0000	200
HMF057 case 2	hmf57c2	1.0000	230
HMF057 case 3	hmf57c3	1.0000	320
HMF057 case 4	hmf57c4	1.0000	400
HMF057 case 5	hmf57c5	1.0000	190
HMF064 case 1	hmf64c1	0.9996	80
LCT010 case 1	lct10c1	1.0000	210
LCT010 case 20	lct10c20	1.0000	280

Table 4. Table 4: Comparison of different nuclear data libraries as well as adjustments from this work using the reduced chi squared (see Eq.  3 ). Note that only experimental data from the (n,tot), (n,non-el), (n,inl), (n, γ 𝛾 \gamma ) and the (n,2n) cross sections of 208 Pb were considered. No data was available for the (n,non-el) channel in the JEFF-3.3 and JENDL-4.0 libraries and therefore, the (n,non-el) was obtained by subtracting the (n,el) from their total cross sections. As presented in the Table, the 1st update (weighted channels) represents the adjustment with differential data only where each channel was assigned a weight equal to its average cross section, while the (unweighted channels) represents adjustment where all considered channels were assigned with equal weights. In the case of the 2nd update, only the hmf57c1 was used for adjustment while in the case of ’This work (global likelihood)’, adjustments were carried out using selected experiments from EXFOR and the hmf57c1 benchmark . In the case of the global likelihood, both the weighted and unweighted channels gave the same ’best’ file. Note: the χ 2 superscript 𝜒 2 \chi^{2} results given for the ENDF/B-VIII.0, JEFF-3.3, JENDL-4.0, TENDL-2017, CENDL-3.1 and the nominal file were obtained using unweighted channels.

Libraries	(n,tot)	(n,non-el)	(n,inl)	(n,2n)	(n,$\gamma$)	Avg ${\chi }^2$
ENDF/B-VIII.0	4.50	0.29	19.34	2.34	1.13	5.52
JEFF-3.3	3.16	0.16	23.03	3.82	0.37	6.11
JENDL-4.0	3.16	0.09	25.04	3.82	0.37	6.57
TENDL-2017	3.54	0.16	3.81	8.97	8.95	5.08
CENDL-3.1	4.55	1.17	23.48	1.61	2.38	6.64
Nominal (prior) file	4.78	0.16	17.52	2.10	8.30	6.57
\pbox20cmThis work (1st update)
(unweighted channels)	5.34	0.02	10.51	4.60	3.62	4.82
\pbox20cmThis work (1st update)
(weighted channels)	3.26	0.18	6.70	3.83	69.11	16.61
\pbox20cmThis work (2nd update)
(unweighted channels)	8.02	6.40	5.52	39.45	70.66	26.01
\pbox20cmThis work (global likelihood)	10.57	4.22	7.94	0.88	0.83	4.89

Table 5. Table 5: Summary of results for the mean k eff subscript 𝑘 eff k_{\rm eff} with corresponding nuclear data uncertainty ( ± plus-or-minus \pm uncertainty on the estimated ND uncertainty) and the Effective Sampling Size (ESS) for the 1st and 2nd priors, and the posterior distributions of the 1st, 2nd and the combined updates. Weighted channels represents the case where each channel was assigned a weight equal to its average cross section over the considered energy range while the unweighted channels represents the case where all channels were assigned equal weights. The 1st prior here represents the distribution of the k eff subscript 𝑘 eff k_{\rm eff} without the inclusion of experimental data while in the case of the 2nd prior, the information from the differential experimental data were used to exclude some random files based on a weight threshold (see Fig.  1 ).

	Unweighted channels			Weighted channels
Distributions	Mean $k_{\mathrm{eff}}$	\pbox20cmND uncertainty
(pcm)	ESS	$k_{\mathrm{eff}}$	\pbox20cmND uncertainty
(pcm)	ESS
1st prior	0.99841	1103$\pm$15	2700	0.99841	1103$\pm$15	2700
2nd Prior	0.99749	1065$\pm$17	2048	0.99707	1025$\pm$16	2038
Posterior (1st update)	0.99546	986$\pm$119	127	0.99558	1018$\pm$87	245
Posterior (2nd update)	0.99981	194$\pm$8	475	0.99979	195$\pm$8	492
Posterior (Combined)	0.99988	164$\pm$56	20	0.99941	203$\pm$55	73

Table 6. Table 6: The ratios of calculated to experimental values (C/E) for selected benchmarks for adjustments with differential experimental data only (1st update), the hmf57c1 benchmark (2nd update), and the combined (global) likelihood function (EXFOR + hmf57c1 benchmark). For the adjustments with hmf57c1, the ENDF/B-VII.0 was used as the reference (base) library for all isotopes, except for 208 Pb which was varied. In the case of ENDF/B-VII.0 (column 6), all isotopes were maintained as the ENDF/B-VII.0 nuclear data library. The benchmark experimental k eff subscript 𝑘 eff k_{\rm eff} values have been presented previously in Table  3 .

Benchmarks	\pbox20cm1st update
(unweighted
channels)	\pbox20cm1st update
(weighted
channels)	\pbox20cm2nd update
(hmf57c1)	\pbox20cmGlobal likelihood
(EXFOR + hmf57c1)	ENDF/B-VII.0
hmf57c1	1.01285	0.99629	1.00002	1.00042	0.98959
hmf57c2	1.01633	1.00439	1.00851	1.00201	0.99888
hmf57c3	1.04167	1.02459	1.03023	1.02918	1.01726
hmf57c4	1.00462	0.99293	0.99793	0.99565	0.98784
hmf57c5	1.04884	1.02934	1.03580	1.03437	1.02188
hmf57c6	1.02047	1.00322	1.00909	1.00769	0.99713
hmf27c1	1.01150	1.00363	1.00676	1.00558	1.00182
pmf35c1	1.00891	1.00015	1.00635	1.00290	0.99856
hmf64c1	1.02110	1.00237	1.00936	1.00263	0.99455
hmf64c2	1.02852	1.00532	1.01265	1.01058	0.99613

Equations

$$ESS=\frac{\big{(}\sum_{k=1}^{n}w_{k}\big{)}^{2}}{\sum_{k=1}^{n}w_{k}^{2}}$$

$$P^{post}({\sigma }_T,p|{\sigma }_E)\propto L({\sigma }_E|{\sigma }_T,p)\times P^{prior}({\sigma }_T,p)$$

$$\chi^{2}_{c(k)}=\frac{1}{N_{p}}\sum_{i=1}^{N_{p}}\bigg{(}\frac{\sigma^{i}_{T(k)}-\sigma^{i}_{E}}{\Delta\sigma^{i}_{E}}\bigg{)}^{2}$$

$${\chi }_{E(k)}^2=\frac{{\sum }_{c=1}^{N_c}{w}_c{\chi }_{c(k)}^2}{{\sum }_{c=1}^{N_c}{w}_c}$$

$$L({\sigma }_E|{\sigma }_T,p)\propto exp\left(-\frac{{\chi }_{E(k)}^2}{2}\right)$$

$$w_{E(k)}=e^{-\frac{1}{2}{\chi }_{E(k)}^2}$$

$$w_{E(k),N}=\frac{w_{E(k)}}{max(w_{E(k)})}$$

$$P^{post}(k_{\mathrm{eff},\mathrm{cal}}^B,{\sigma }_T|k_{\mathrm{eff},\exp }^B)\propto L(k_{\mathrm{eff},\exp }^B|k_{\mathrm{eff},\mathrm{cal}}^B,{\sigma }_T)\times P^{prior}(k_{\mathrm{eff},\mathrm{cal}}^B|{\sigma }_T)$$

$$L(k_{\mathrm{eff},\exp }^B|k_{\mathrm{eff},\mathrm{cal}}^B,{\sigma }_T)\propto exp\left(-\frac{{\chi }_{B(k,j)}^2}{2}\right)$$

$$\chi^{2}_{B(k,j)}=\Bigg{(}\frac{k^{B}_{\rm eff(k,j)}-k^{B}_{\rm eff,E(j)}}{\Delta k^{B}_{\rm eff,E(j)}}\Bigg{)}^{2}$$

$$\overline{k_{\mathrm{eff}(\mathrm{k},\mathrm{w})}^B}=\frac{\sum _{k=1}^n{w}_{B(k,j)}\cdot k_{\mathrm{eff}(\mathrm{k},\mathrm{j})}^B}{\sum _{k=1}^n{w}_{B(k,j)}}$$

$$Var(\overline{k_{\mathrm{eff}(\mathrm{k},\mathrm{w})}^B)}=\frac{\sum _{k=1}^n{w}_{B(k,j)}{k}_{\mathrm{eff}(\mathrm{k},\mathrm{j})}^B}{\sum _{k=1}^n{w}_{B(k,j)}}-{\overline{k_{\mathrm{eff}(\mathrm{k},\mathrm{w})}^B}}^2$$

$$L({\sigma }_E,k_{\mathrm{eff},\exp }^B|{\sigma }_T,p,k_{\mathrm{eff},\mathrm{cal}}^B)=L({\sigma }_E|{\sigma }_T,p)L(k_{\mathrm{eff},\exp }^B|k_{\mathrm{eff},\mathrm{cal}}^B,{\sigma }_T)$$

$$P_{comb}^{post}(k_{\mathrm{eff},\mathrm{cal}}^B|{\sigma }_E,k_{\mathrm{eff},\exp }^B)=P^{prior}(k_{\mathrm{eff},\mathrm{cal}}^B|{\sigma }_T)L({\sigma }_E|{\sigma }_T,p)L(k_{\mathrm{eff},\exp }^B|k_{\mathrm{eff},\mathrm{cal}}^B,{\sigma }_T)$$

$$w_{T,k}=w_{E(k)}\cdot w_{B(k,j)}$$

Full text

1 Introduction

Nuclear reactors like many other technical systems are complex in nature and therefore require the coupling of several mathematical models and sub-models in order to describe and analyse them. The models used for the analysis of nuclear reactor systems can broadly be classified into the following categories: (1) Basic nuclear physics models used for the computation of basic physical observables such as nuclear reaction cross sections and angular distributions, (2) models implemented in nuclear data processing codes such as NJOY [1] and PREPRO [2], (3) models for neutron transport including reactor kinetics, (4) models implemented in thermal-hydraulics and computational fluid dynamic codes, (5) reactor fuel mechanics models, and (6) reactor dosimetry models, among others. These models interact with each other since the output of some models (lower-level models) are normally used as input to other models (higher-level models). Also, feedbacks from higher-level models are normally given in order to improve lower-level models. Furthermore, different experimental data sets are utilized for the calibration and validation of these models at each level. This makes it particularly difficult and computationally expensive to integrate or combine the various activities over the entire calculation chain into a single process. With the increase in computational power and the improvements in nuclear reaction theory, a novel approach called ’Total Monte Carlo’ (TMC) was developed and presented in Ref. [3] with the goal of propagating uncertainties from basic physics parameters to applications. While this goal has been achieved to a large extent, the explicit inclusion of experimental data from both differential and integral benchmark data for combined data adjustments is still outstanding. This work seeks to combine these experimental data (differential and integral benchmark data) within a Bayesian framework, for nuclear data adjustments in the fast neutron energy region (below 20 MeV) using the Total Monte Carlo (TMC) method.

Within the TMC method, data adjustments to fit selected integral benchmark experimental data has been presented for example in Refs. [4, 5, 6, 7, 8] and used for uncertainty reductions in reactor applications [9, 10, 11, 12]. In Refs. [4, 5], a method that combines evaluations with a posteriori adjustments to integral measurements referred to as the ’Petten method’ was developed and utilized. A random search was performed to identify the best possible nuclear data file using a goodness of fit estimator. The ’Petten method’ was applied for the adjustment and evaluation of neutron induced reactions on 239Pu [13] and 63,65Cu [4] as well as for the improvement of H in $H_{2}O$ neutron thermal scattering data [14]. In Ref. [11], similar to Refs. [5, 6, 7, 8], only information from benchmark experiments were used; differential experimental data were not explicitly incorporated into model calculations.

Also, the inclusion of differential experimental data within the Total Monte Carlo (TMC) method has been presented previously in Refs. [15, 16, 17, 18, 19]. In Refs. [17, 18, 19], similar to this work, weights equal to the likelihood function were assigned to each random nuclear data file based on only differential experimental data. In Ref. [16], differential experimental data were incorporated into model calculations by computing a weighted $\chi^{2}$ for each reaction channel. In Refs. [15, 16, 17, 18, 19], only information from differential experimental data were utilized in the adjustments.

In this work, we present a method for combining experimental data from both differential and integral measurements within a Bayesian framework for data adjustments. The proposed method is applied for the adjustment of n+208Pb cross sections in the fast energy region.

2 Theory and Methods

2.1 Updating scheme

In Fig. 1, a diagram showing the proposed Bayesian updating scheme for multi-level uncertainty propagation within the TMC method is presented. The term uncertainty as used in this work is represented by an estimated one standard deviation of the distribution under consideration. The TMC method has been presented extensively in several references: [3, 4, 19, 20, 21] as well as applied to many applications: [10, 22, 23, 24, 25, 26, 27, 28, 29].

The first step in our updating scheme involves the identification of model input parameters (and the corresponding parameter widths) that play significant roles in defining or characterizing the reaction observables of interest. These parameters were determined through a combined use of expert judgement and sensitivity analysis. The next step is to determine the distribution for each model parameter of interest. Since our goal is to use non-informative (also known as ’objective’) priors as much as possible, we assumed here that, we have no prior knowledge of the parameters under consideration. Therefore the uniform distribution was used as the prior distribution for all selected model parameters. Also, the parameters were sampled in a non correlated way. A complete objective prior in our case as observed also in Refs.[18, 19], is however not possible since the evaluation in this work made use of a nominal file which was selected by comparing TALYS results with available differential experimental data. The nominal file is a complete evaluation with its corresponding input parameter set and represents the evaluator’s best effort before adjustment or assimilation.

The third step is to determine the nominal file with its parameter set for adjustment. A good starting point is to use reference input parameters from the RIPL database [30] for model calculations, if available. In this work however, the nominal file (with its parameter set) was taken from the TENDL-2015 evaluation [31]. It should be noted that, no significant changes have been made to the TENDL-2015 evaluation of 208Pb in the recent TENDL library release [32]. Next, all the model parameters were randomly varied within their assigned parameter widths using the TALYS code system (T6) [33] and a total of 2700 random nuclear data files were produced.

The fourth step involves carrying out a probabilistic data adjustment or assimilation using a Bayesian updating scheme (more details are presented in section 2.2). Two Bayesian updates are performed in this work: (1) using only differential experimental data (1st update), and (2) using integral experimental data (2nd update). The main objective of this work is to combine the individual likelihood functions from each update into a global (combined) likelihood function for each nuclear data (ND) file. From these global likelihood functions, the new ’best’ file which takes both differential and integral experimental data explicitly into account, is selected. The best file is then the file (with its parameter set) that maximises the combined (global) likelihood function. In the case of the 1st update, as shown in Fig. 1, we define an acceptance criterion such that any ND file with a weight less than an assigned threshold file weight ($w_{thres}$ = 2.06e-09) which corresponds to a critical $\chi^{2}$ value of 40, is discarded. The critical $\chi^{2}$ value correspond to the upper-tail critical value of the $\chi^{2}$ distribution with 28 degrees of freedom (i.e. varied model parameters; see Table 1) at 95% confidence level. The critical $\chi^{2}$ value was used here with the assumption that the $\chi^{2}$ distribution (from the 1st update) is in the form of a gamma distribution with the scale parameter equal to 2 and the shape parameter equal to N/2, where N is the number of degrees of freedom. The critical $\chi^{2}$ values can be found in the appendix of most statistical books, for example, in Appendix C, Table G of Ref. [34]. The relationship between the $\chi^{2}$ and the weight ($w_{k}$) is given in Eq. 6. If the $\chi^{2}$ value (or weight) computed for a particular random ND file is greater than the critical $\chi^{2}$ = 40 (or lower than the corresponding threshold weight ($w_{thres}$)), it gives an indication that the file under consideration did not fit the selected experimental data and is therefore rejected. In this way, files with very low weights were rejected which could lead to a savings in computational time. For example, in the case of this work, after setting a threshold weight (as presented in Fig. 1), a total of 2046 files were accepted; down from the initial 2700 files. This represents a 24% reduction in the number of random ND files used. It should be noted however that, this could introduce some bias into the calculation depending on the choice of the cut-off parameter. To address this, the use of Russian roulette for randomly rejecting ND files with insignificant weights was proposed in Refs. [17, 18]. This approach was however not used in this work.

Since the random samples were drawn from a rather large model parameter space (see section 2.2.1), a possibility could arise where a large number of the random ND files are assigned with very low or insignificant weights. Therefore an Effective Sample Size (ESS) (given by Eq. 1) has been computed for each update:

[TABLE]

The ESS was used in order to determine the number of random ND files with significant weights for both the prior and posterior distributions. The ESS gives an indication on the number of random ND files with significant impact on the distributions considered. Also, in order to determine if the prior and posterior $k_{\rm eff}$ distributions have converged, convergence plots (i.e. for the mean and the ND uncertainty as a function of random samples), are presented for each update.

2.2 Bayesian updating within TMC

2.2.1 Priors - model parameters and their distributions

In Table 1, the model parameters as well as the parameter widths (given as a fraction (%) of their absolute values), assigned to selected model parameters are presented. The parameter widths of $a$, $g_{\pi}$ and $g_{\nu}$ as presented in the Table are given in terms of the mass number $A$. Where $a$ is the real central diffuseness, $g_{\pi}$ and $g_{\nu}$ are the single-particle state densities as used in exciton model analyses [35]. The parameter widths were determined by means of comparisons of scattered random TALYS curves with differential experimental data from the EXFOR database [33]. In order to increase the parameter space, similar to Ref. [19], the model parameter widths (as presented in Table 1) were all multiplied by a factor of three. In addition, the parameters were sampled from uniform distributions. This was done in order to obtain non-informative priors as much as possible for each model parameter. However, as observed earlier in section 2.1, a complete state of non-informative prior can never be reached since experiments were used to fine tune the default model parameters used in the TALYS code as well as in determining the nominal file. All the parameters presented were then simultaneously varied within their parameter widths in a TMC fashion to produce a total of 2700 random 208Pb nuclear data files.

From Table 1, $r$ represents the real central radius, $V$ is the volume-central, $D$ and $SO$ respectively are the surface-central and spin-orbit potentials and $W$ is the imaginary depth of the optical model [36]. $\sigma^{2}$ is the spin cut-off parameter which represents the width of the angular momentum distribution of the level density [36]. The superscript $n$ denotes neutron induced reactions while the subscripts $v$ and $s$ respectively denote the real volume and real surface of the optical model. $\Gamma_{\gamma}$ is the average radiative capture width; $M^{2}$ is the average squared matrix element of the residual interaction as used in exciton model analyses [35]. The default optical model potentials (OMP) used in TALYS are the local and global parametrization of Koning and Delaroche (see Ref. [37]) and therefore were also used in this work.

2.2.2 1st Bayesian update: using differential experimental data

Selected differential experimental data from the EXFOR database [38] were used for the first Bayesian update. The selection of experimental data sets were based on a number of criteria: (1) We do not consider experiments that are inconsistent with other experimental sets and also deviates from the trend of other evaluations (if available) as well as our model calculations, (2) We assume 10% uncertainty for experimental data sets for which no uncertainties were reported but which however were reviewed and classified as good quality experiments in Ref. [39]. As noted in Ref. [19], these classifications are usually subjective in nature and could therefore introduce some bias into the evaluation process. However, (as observed also in Ref. [19]), blindly using all experiments in EXFOR without selection (or assigning weights) usually leads to very large $\chi^{2}$ values between model calculations and experiments. 3) If a particular author repeats a set of experiments for the same or similar energy range for the same cross section (and they are found to be inconsistent), we select the updated (or current) measurements, (4) In some cases, we also drop the first and/or last measured points (with respect to incident neutron energy) since these measurements can easily be dominated with background noise.

In Table 2, we present the experimental data that were selected for the 1st update showing the number of data points per reaction channel, the corresponding MT numbers (in ENDF terminology), the EXFOR ID, the first author and the year the measurements were carried out. It should be noted that only data points that fall between 5 and 20 MeV were used for the adjustments and therefore presented.

For the 1st Bayesian update we consider the following: given that, the probability distribution of the model parameters ($p$) before any data (cross sections in this case) were observed (also known here simply as the prior) is $P^{prior}(\sigma_{T},p)$, and $L(\sigma_{E}|\sigma_{T},p)$ is the probability of the experimental data occurring given that the model parameters are accurate (also known as the likelihood function); the posterior, which is the probability of the parameter values being accurate given differential experimental data ($P^{post}(\sigma_{T},p|\sigma_{E})$), can be given as:

[TABLE]

where $\sigma_{T}$ and $\sigma_{E}$ denote the TALYS calculated and experimental observables (cross sections in this case) respectively. If we assume that the off-diagonal elements of the experimental covariance matrix is zero, we can compute a reduced chi square ($\chi^{2}_{c(k)}$) per reaction channel and for the random ND file $k$, similar to Ref. [19] as follows:

[TABLE]

where $\sigma^{i}_{T(k)}$ is a vector of TALYS calculated observables as a function of incident neutron energy ($i$), found in the $k^{th}$ random ND file for the channel $c$ and $\sigma^{i}_{E}$ is a vector of experimental observables for the channel $c$ at incident neutron energy ($i$), $\Delta\sigma^{i}_{E}$ is the experimental uncertainty at an incident energy $i$ of channel $c$, and $N_{p}$ is the total number of experimental points per reaction channel considered. Where there were no matches in energy ($i$) between TALYS output and the experimental data set for the $c^{th}$ channel, TALYS data were interpolated to match corresponding experimental values. For the purpose of this work, the reduced $\chi^{2}$ and the $\chi^{2}$ are used interchangeably and refer to the same thing.

One major difference in the computation of the chi square presented in Eq. 3 and that used in Ref. [19] is that, while the average is taken per experimental data set in Ref. [19], the average is simply taken per channel in this work. By averaging per channel, we assign all experimental data points per channel with equal weights. Also, in order to give each channel equal weight, we then take an average over the considered channels. Additionally, in Ref. [19], the effects of model defects were taken into account by normalizing the $\chi^{2}_{E(k)}$ for each random ND file (as presented in Eq. 5) with the chi square obtained for the global TALYS calculation (usually TALYS run 0). A similar solution was presented also in Ref. [40] where instead of normalizing the $\chi^{2}_{E(k)}$ with the global TALYS calculation, the $\chi^{2}_{E(k)}$ was normalized with the minimum chi square ($\chi^{2}_{min}$). Since the underlying assumption for our Bayesian update is that the models used are ’good enough’, no modifications as in Refs. [19, 40], were made to the likelihood function. Since our models are not perfect, this assumption is entirely not true, however, we think that the quality of the models used are adequate for the purpose of this evaluation.

From Eq. 3, we can compute a weighted chi square ($\chi_{E(k)}^{2}$) given by:

[TABLE]

where $w_{c}$ is the channel weight and $N_{c}$ is the number of considered channels. Two types of channel weights were used in this work: (1) $w_{c}=\overline{\sigma_{c}}$ and (2) $w_{c}=1/N_{c}$ (also called unweighted channels since all channels have equal weights). $\overline{\sigma_{c}}$ is the average cross section for the channel $c$. By assigning channel weights equal to the average cross section as in (1), we give more weight to channels with relatively large cross sections such as the (n,tot), (n,inl) and the (n,2n) cross sections of 208Pb while channels such as the (n,$\gamma$) with smaller average cross sections are assigned with low channel weights. This is however not consistent for a ’general purpose’ library since the fit would favour some cross sections than others. In (2), all considered channels are assigned with equal weights and therefore all channels are of equal importance in the adjustment procedure. Case (2) is therefore referred to as the unweighted channel case in this work. From the $\chi_{E(k)}^{2}$ computed, the likelihood function can now be computed using an expression given within the Unified Monte Carlo approach (UMC-B) [41] and utilized also in the TMC + UMC-B method presented in Ref. [18]:

[TABLE]

The motivation for using Eq. 5 has been presented previously in Ref. [18]. From Eq. 5, each random ND file in the case of the 1st update was assigned a file weight based on only selected differential experimental data from EXFOR using Eq. 6:

[TABLE]

where $w_{E(k)}$ is the weight assigned to the $k^{th}$ ND file. In order to relate the file weight values to one, the weights were normalized with their maximum weight as follows:

[TABLE]

where $max(w_{E(k)})$ and $w_{E(k),N}$ are the maximum and the normalized weights from the 1st update respectively. The weights computed for the ’best’ files from this work have been compared with weights computed for the ENDF/B-VIII.0, JEFF-3.3, JENDL-4.0, TENDL-2017 and CENDL-3.1 nuclear data libraries using the same methodology and experimental data and presented in Table 4.

2.2.3 2nd Bayesian update: using integral experimental data

In Table 3, the benchmark experiments used for the 2nd update are presented showing the experimental benchmark $k_{\rm eff}$ and the experimental benchmark uncertainty. These benchmarks can be obtained from the International Handbook of Evaluated Criticality Safety Benchmark Experiments (ICSBEP) [42]. From the Table, HMF stands for Highly Enriched Uranium (HEU) Metallic Fast, PMF for Plutonium Metallic Fast, while LCT stands for Low Enriched Uranium (LEU) Compound Thermal benchmarks.

Given that $P^{prior}(k^{B}_{\rm eff,cal}|\sigma_{T})$ is the probability distribution of the $k_{\rm eff}$ given cross sections from the 1st update (i.e. before the inclusion of integral benchmark data) and $L(k^{B}_{\rm eff,exp}|k^{B}_{\rm eff,cal},\sigma_{T})$ is our likelihood function; the posterior distribution $(P^{post}(k^{B}_{\rm eff,cal},\sigma_{T}|k^{B}_{\rm eff,exp})$, can be expressed as:

[TABLE]

Where the likelihood function is given by:

[TABLE]

In the case of one benchmark (as utilized in this work), the chi square ($\chi^{2}_{B(k,j)}$) as a function of random ND ($k$) and benchmark $j$, can be expressed as:

[TABLE]

where $k^{B}_{\rm eff(k,j)}$ is a vector of calculated response parameters for the $k^{th}$ random nuclear data file and the $j^{th}$ benchmark, $k^{B}_{\rm eff,E(j)}$ is a vector of integral benchmark experimental observables with corresponding experimental uncertainty ($\Delta k^{B}_{\rm eff,E(j)}$) for the $j^{th}$ benchmark. Since only one benchmark was used in the adjustment, no correlations between benchmarks were considered. The $k^{B}_{\rm eff(k,j)}$ in Eq. 10 was computed by varying 208Pb nuclear data while holding constant, all other isotopes contained in the benchmark under consideration. In the computation of the $\chi^{2}$ as presented in Eq. 10, only benchmark experimental uncertainties were taken into account. However, as noted earlier in Ref. [11], because computer codes such as MCNP [43] are normally used for benchmark calculations and analyses, these benchmarks also contain calculation uncertainties. These uncertainties could come from (for example), our inability to model the benchmark geometry accurately with the code used, statistics (i.e. in the case where a Monte Carlo code is used) or from uncertainties due to nuclear data (since nuclear data are used as inputs to these codes). These calculation uncertainties were however not considered in this work.

Using the $\chi^{2}_{B(k,j)}$ presented in Eq. 10, weights similar to Eq. 6, were assigned to each random nuclear data based on experimental data from the hmf57c1 benchmark. The weights were also normalized with the maximum weight in order to keep the weight values between 0 and 1. As a proof of concept, we have only used one benchmark - the hmf57c1 benchmark for adjustment in this work; all other benchmarks (as presented in the Table 3) were used for testing the adjusted file’s performance. The motivation to use the hmf57c1 is because, aside being a fast neutron spectrum benchmark, it was observed to be sensitive to 208Pb nuclear data. Additionally, it has a relatively small benchmark uncertainty of 200 pcm and also comes with five different benchmark cases which can be used for verification purposes. A benchmark case as used in this work (as also used in Ref. [42]), is used to define a series of similar benchmarks where a limited number of parameters (e.g. geometrical parameters) have been varied.

After incorporating the weights in the 2nd Bayesian update, an updated (or weighted) mean and the corresponding weighted variance were computed for each random ND file using Eq. 11 and  12:

[TABLE]

[TABLE]

where $\overline{k^{B}_{\rm eff(k,w)}}$ and $Var(\overline{k^{B}_{\rm eff(k,w)}})$ are the weighted mean and weighted variance respectively, and $w_{B(k,j)}$ is the file weight assigned to the $k^{th}$ random ND file using experimental information from the benchmark $j$.

2.2.4 Global likelihood function: Combining weights from differential and integral experiments

According to Ref. [44], a statistically rigorous way to combine experimental data from different measurements is to combine the likelihood functions computed for the individual measurements. Similar to  [44], we combine here the likelihood functions computed from the 1st and 2nd updates into a global likelihood function for each random ND file. From probability theory, if two events are independent, the probability of both events occurring is given as the product of the probabilities of each event occurring. The likelihood function in this case can be considered as the probability of the experimental data occurring given that model parameters are accurate. Therefore, if we assume that the differential and integral experimental data are independent or uncorrelated, a global (combined) likelihood function (also known here as the combined weight) can be obtained as a product of the individual likelihood functions computed from differential experimental data ($L(\sigma_{E}|\sigma_{T},p)$), and integral benchmark data ($L(k^{B}_{\rm eff,exp}|k^{B}_{\rm eff,cal},\sigma_{T})$):

[TABLE]

where $L(\sigma_{E},k^{B}_{\rm eff,exp}|\sigma_{T},p,k^{B}_{\rm eff,cal})$ is the global (combined) likelihood function. From Bayes theorem, our new posterior distribution ($P^{post}_{comb}(k^{B}_{\rm eff,cal}|\sigma_{E},k^{B}_{\rm eff,exp}$)) which takes into consideration information from both differential and integral benchmark experiments, can be given as:

[TABLE]

From Eq. 13, a combined weight ($w_{T,k}$) which is equal to the global likelihood function, is computed for each $k$ ND file:

[TABLE]

where $w_{E(k)}$ and $w_{B(k,j)}$ are the weights for the $k^{th}$ random ND file computed using differential experimental data and integral benchmark experiment $j$ respectively. In order to compare $w_{E(k)}$ and $w_{B(k,j)}$, the weights computed for each update were both normalized using their maximum weights. In this way, the combined weight was constrained between 0 and 1. An algorithm was developed to automatically select the new ’best’ file i.e. the file with the parameter set that maximizes the global (combined) likelihood function. The performance of the ’best’ file has been compared against differential data from EXFOR and evaluations from the major nuclear data libraries using the reduced $\chi^{2}$ as the goodness of fit estimator and presented in Table 4 as well as with integral benchmarks (see Table 6).

2.3 Simulation and analysis

Since resonance structure extends well beyond 2 MeV for the (n,tot) and the (n,el) cross sections of 208Pb, the neutron energy range considered in the 1st update (using differential data only) was from 5 to 20 MeV. However, in order to create a complete ENDF file from 0 to 20 MeV as normally required for applications, the entire MF-2 (in ENDF terminology for the resonance parameters) was adopted from the JEFF-3.1 [45] library using the TARES code [46], while the elastic angular distributions (MF4-MT2) were adopted from ENDF/B-VII.1 [47] evaluation. The nuclear reaction code TALYS-1.6 [36] was used for the computation of all reaction observables while the TEFAL code [48] was used for converting TALYS output into the well known ENDF format.

The random ENDF ND files were then processed into the ACE format at 293K using NJOY99 version 336 [42]. The ACE formatted random files were used in the MCNPX code (version 2.5) for the computation of the $k_{\rm eff}$. For each benchmark, criticality calculations were carried out with 5000 neutron particles for 500 criticality cycles resulting in an average statistical uncertainty of 48 pcm. In order to speed up calculation time, the fast TMC as presented in Ref. [49] and utilized in Refs. [22, 25] for uncertainty propagation was also used in this work. Hence, the seed of the MCNPX code was changed for each random run by means of the DBCN card.

3 Results

3.1 1st Bayesian Update: EXFOR data

The distributions of cross sections as a function of incident neutron energy for the $(n,tot)$ and $(n,inl)$ cross sections are presented for 100 random 208Pb cross sections in Fig. 2. The spread observed in the cross sections can be attributed to the variation of model parameters using the TMC method. These cross sections, were used as the prior for the second Bayesian update.

In order to determine if the random nuclear data files for the prior $k_{\rm eff}$ distribution converged, the mean (of the prior $k_{\rm eff}$ distribution) with its corresponding ND uncertainty is plotted against the number of random files for a number of iterations and presented in Fig. 3. The $k_{\rm eff}$ distribution which is the 1st prior distribution obtained by varying 208Pb nuclear data in the hmf57c1 benchmark is also presented (right of Fig. 3). By 1st prior distribution, we refer to the $k_{\rm eff}$ distribution before the implementation of the threshold weight criterion presented in Fig. 1. It should be noted that, the $k_{\rm eff}$ distribution presented in Fig. 3 (as in the case for all other $k_{\rm eff}$ distributions in this work), is fitted to a normal distribution for the purpose of eye guidance only. From left of the figure (Fig. 3), even though some fluctuations can be observed in the convergence plots shown, the percentage change in the last five iterations (for both the mean and the ND uncertainty) was observed to be less than 1%. The error bars on the ND uncertainty represent the estimated uncertainty on the ND data uncertainty obtained for each iteration. More details on how to compute the uncertainty on the ND uncertainty was presented in Refs. [22, 24, 25]. It can be observed from the figure that, the uncertainty on the ND uncertainty estimated decreases with increasing sample size as expected.

From the $k_{\rm eff}$ distribution presented (right of Fig. 3), a mean, standard deviation and skewness of 0.99841, 1104 pcm and 0.55 respectively, were obtained. The positively (slightly) skewed $k_{\rm eff}$ distribution is in agreement with earlier results obtained for 208Pb nuclear data variation in fast criticality systems [50] and for the European Lead Cooled Training Reactor (ELECTRA) [22]. The skewness observed is due to the relatively high calculated $k_{\rm eff}$ values obtained for the hmf57c1 benchmark. According to Ref. [42], these high $k_{\rm eff}$ values give an indication that the current lead cross sections have some inaccuracies which could result in more neutron reflection than is warranted by the experimental results.

In this work, two types of computations of the posterior distributions have been performed: (1) using weighted channels and (2) unweighted channels. As mentioned earlier in section 3.1, in the case of weighted channels, each channel was assigned a weight equal to its average cross section over the energy range of interest while the unweighted channels represents the case where all channels were assigned equal weights. It should be noted however that, only the plots for the weighted case have been presented through out this work.

Similarly in Fig. 4, the convergence in the mean and the 208Pb ND uncertainty for the posterior $k_{\rm eff}$ distribution of the 1st update due to the variation of 208Pb nuclear data in the hmf57c1 benchmark is presented. Also presented is a scatter plot of the weights with its corresponding distribution, as well as the corresponding 2nd prior and posterior $k_{\rm eff}$ distributions. The 2nd prior distribution as presented in Fig. 4 represents the $k_{\rm eff}$ distribution after the implementation of the threshold weight as presented in Fig. 1(note: a total of 2046 random ND files were accepted). The weights were computed using selected differential experimental data from EXFOR presented ealier in Table 2. From the convergence plot, it can be seen that both the mean and the ND uncertainty converged after a number of iterations - the percentage change in the last five iterations was less than 1%. It can also be observed that the error bars (which represents the uncertainty on the ND uncertainty estimated) are large but reduces as the number of samples increases with the number of iterations. In all, a total of 2040 accepted random ND files (i.e. accepted files based on the weight threshold criterion presented in Fig. 1) were used which resulted in a mean and a ND uncertainty of 0.99707 and 1025$\pm$16 pcm respectively (i.e. in the case where the channels were weighted). The prior ND uncertainty reduced marginally to posterior ND uncertainty of 1018 pcm. This is however different in the case of unweighted channels, where a relatively larger reduction in the posterior ND uncertainty was achieved (ND uncertainty of 986 pcm was obtained). The reduction in ND uncertainty observed, gives an indication that, by including information from differential experimental data in the computation of integral observables such as the $k_{\rm eff}$, the spread in $k_{\rm eff}$ due to ND uncertainties could be reduced. From Fig. 4, it can be observed from the distribution of file weights that, a small number of random ND files were assigned with large weights, above 0.8 for example. This explains why a relatively small reduction in the posterior spread was obtained as seen in the Figure. The few ND files with significant weights obtained are however in agreement with earlier observations made in Refs. [17, 18]. As noted in Ref. [18], if we assume that the models are perfect, a possible explanation could be attributed to our sampling of model parameters from a rather wide parameter space which could lead to a combination of model parameters being drawn from a region of the parameter space where the likelihood is low. A possible solution is then to increase the sample size which could be computationally expensive. Another solution is to resample model parameters based on the ’best’ file but instead, with smaller parameter widths. This could lead to a more even distribution of file weights computed from EXFOR. But since the models are not perfect, another reason could be because of the presence of model defects. The effect of these model defects have not been considered in this work. There are however on-going efforts to include the effect of these defects in modern nuclear data evaluations (See Refs. [19, 51, 52, 53]). Another explanation could be due to the presence of resonance-like spikes in he (n,inl) cross section for example, which were not taken into account in our model calculations since the TALYS code only gives smooth curves for cross sections. Based on the file weights computed (i.e. in the case of weighted channels), an effective sample size (ESS) of 245 was obtained while an ESS of 127 was obtained for the unweighted channel’s case.

From the file weights computed for the 1st Bayesian update, a new ’best’ file (i.e. file number 0752) was selected. To determine the performance of the selected file, an average $\chi^{2}$ was computed between the selected file and selected experimental data from the following channels: ((n,tot), (n,non-el), (n,inl), (n,$\gamma$) and the (n,2n)), and compared with average $\chi^{2}$ values computed for different ND libraries (ENDF/B-VIII.0, JEFF-3.3, JENDL-4.0, CENDL-3.1 and TENDL-2017) using the same experimental data sets and methodology, and presented in Table 4. The nominal file as presented in the Table represents the prior file before adjustment. No data was available for the (n,non-el) cross section in the JEFF-3.3 and JENDL-4.0 libraries, therefore, the (n,non-el) cross sections for these libraries were obtained by subtracting the (n,el) from their total cross sections. In the case of the $\chi^{2}$ results presented for the ENDF/B-VIII.0, JEFF-3.3, JENDL-4.0, TENDL-2017, CENDL-3.1 libraries as well as our nominal file, the channels used were unweighted.

From Table 4, the 1st update (weighted channels) represents the adjustments with only differential data from EXFOR where each channel was assigned a weight equal to its average cross section while the (unweighted channels) represents adjustment where all channels were not weighted. In the case of the 2nd update (as presented in the Table), only the hmf57c1 was used for adjustment while in the case of ’This work (global likelihood)’, adjustments were carried out using selected experiments from EXFOR + the hmf57c1 benchmark. More results on the 2nd update and the global (combined) likelihood are presented in sections 3.2 and 3.3 respectively. It can be observed from Table 4 that the evaluation from this work (i.e. in the case of 1st update (unweighted) and adjustment with the global likelihood function) out performed the other libraries with respect to the average $\chi^{2}$. In the case where the 1st update (weighted channels) was compared with its unweighted, it was observed that, the weighted ’best’ file out performed the unweighted file for the (n,tot), (n,inl), (n,2n) channels but compared poorly in the case of the average $\chi^{2}$ and the (n,$\gamma$) channel. This is not surprising since channels with larger weights are given more importance in the random search for the best file than channels with smaller weights. Since the (n,tot), (n,inl), (n,2n) channels have relatively large average cross sections, they were assigned with larger channel weights while the (n,$\gamma$) channel with a relatively smaller average cross section was assigned a smaller channel weight. One possible reason for the relatively large average $\chi^{2}$ of 26.01 obtained for the 2nd update (as seen in Table 4) could be due to over-fitting to integral experimental data since the calculation uncertainties were not taken into account. Also, for a fast system such as the hmf57c1 benchmark, the (n,$\gamma$) channel for example would have minimal effect on the computation of the benchmark $k_{\rm eff}$, hence the large $\chi^{2}$ value obtained for this channel.

3.2 2nd Bayesian Update: integral benchmarks

The 2nd prior (as presented in the 1st update) was used as the prior distribution for the 2nd update. The prior distribution here then represents $k_{\rm eff}$ distributions using the ND files that were accepted (or passed the weight threshold test presented earlier in Fig. 1). In Fig. 5, the prior (2nd) and posterior $k_{\rm eff}$ distributions due to the variation of 208Pb nuclear data for the hmf57c1 benchmark with its corresponding weight distribution computed using only experimental information from the hmf57c1 benchmark is presented. Also presented is the convergence plot of the mean and ND uncertainty of the posterior $k_{\rm eff}$ distribution as well as a scatter plot of the weights.

Similar to Fig. 4, the fluctuations observed in the mean and the ND uncertainty were observed to be small (less than 1% (relative difference) was obtained between the last 5 iterations). Large error bars can be observed for the convergence of the ND uncertainty especially in the case of smaller sample sizes which reduces as the number of iterations (and sample size) increases. From Fig. 5, it can be observed that, the large prior uncertainty of 1062 pcm was reduced to a posterior uncertainty of 194 pcm representing an uncertainty reduction of 82%. Also, the mean of the posterior distribution of 0.99981 obtained is close to the experimental $k_{\rm eff}$ as desired. It should however be noted that, in this particular case, only the experimental benchmark uncertainty information was used in the computation of the file weights, hence the very narrow spread observed in the posterior distributions. As noted earlier, in Ref. [11] an additional uncertainty term (uncertainty in the calculation due to nuclear data) was included in the computation of file weights which resulted in much larger combined benchmark uncertainties and therefore larger posterior distributions. This approach was however not used in this work.

Additionally, it should be noted that, the $\chi^{2}$ values (and therefore the file weights) in the 2nd update were computed with only experimental data from the hmf57c1 benchmark and therefore, even though an updated $k_{\rm eff}$ = 1.00002 (see Table 6) was obtained for the new ’best’ file (i.e. for the 2nd update), the file performed poorly when compared back with differential experimental data from EXFOR (see $\chi^{2}$ values presented in Table 4). An average $\chi^{2}$ value of 26.01 was obtained for the new ’best’ file (2nd update), compared with 4.82 from the 1st update (unweighted channels), and 5.52, 6.11, 6.57, 5.08 and 6.64 for the ENDF/B-VIII.0, JEFF-3.3, JENDL-4.0, TENDL-2017 and CENDL-3.1 libraries respectively. The combined adjustment using the likelihoods from the two updates (i.e. the global likelihood (unweighted)), however, gave an average $\chi^{2}$ of 4.89. This gives an indication that by combining information from both differential and integral experiments, improvements on the nominal (prior) file can be achieved.

Table 5 presents a summary of results of the mean $k_{\rm eff}$ with corresponding ND uncertainties, and the Effective Sampling Sizes (ESS) for the 1st and 2nd prior distributions as well as for the posterior distributions of the 1st, 2nd and combined updates. The 1st prior represents the $k_{\rm eff}$ distribution of random nuclear data files before the implementation of the threshold weight (i.e. with a total of 2700 random ND files) while the 2nd prior represents the distribution after the implementation of the threshold weight (a total of 2046 random ND files were accepted) as presented in Fig. 1. From the Table, it can be observed that the ESS for the (unweighted case) are relatively smaller than the weighted case for the posterior distributions. This could be attributed to the large weights assigned to the channels with large average cross sections such as the (n,tot), (n,inl) and the (n,non-el) channels and which also, TALYS-1.6 was able to reproduce within experimental uncertainties. It was however observed that, TALYS-1.6 was not able to reproduce the (n,$\gamma$) cross sections of 208Pb (which also has a relatively smaller average cross section) for a large number of the random cross sections produced.

In the case of the posterior (combined) as presented in the Table, small ESS of 20 and 73 were obtained for the unweighted and weighted channels respectively. The small ESS obtained, could be attributed to the rather few files assigned with large weights which resulted in larger differences in file weights for some random ND files. Also, smaller ESS results in unconvergence of the distribution under consideration since it implies that very few files have significant impact on the posterior distribution. In order to achieve a more even weight distribution and hence higher ESS values, outlier ND files could have been discarded and the file weights renormalized with the maximum weight. This however, must be done with caution as extremely good ’outlier’ ND files with large weights may be discarded. Since the main objective of the combined adjustment is to identify the ND file that reproduces both differential and integral data, the outlier ND files were not discarded.

3.3 Using the global likelihood function: EXFOR + hmf57c1 benchmark

A final test of an adjustment is to compare the final adjusted file back with differential experimental data from EXFOR as well as with relevant benchmarks in order to determine its performance. Additionally, comparisons are made with other existing evaluations. In Figs. 6 and 7, the performance of the adjustments from this work are compared with available differential experimental data (between 5 to 20 MeV) from EXFOR for the (n,tot), (n,el), (n,inl), (n,2n), (non-el) and (n,$\gamma$) cross sections of 208Pb as well as with evaluations from the ENDF/B-VIII.0 and JEFF-3.3 nuclear data libraries. Since no experimental data were available in EXFOR for the (n,el) cross section, we have only compared our evaluations with the ENDF/B-VIII.0 and JEFF-3.3 nuclear data libraries in the case of the (n,el) cross section. The nominal (prior) file is the prior file around which the other random ND files were generated. As seen from the Figs. 6 and 7, the individual authors of the different experimental sets from EXFOR have not been listed, instead, they have been been lumped together and named as EXFOR. The weighted and unweighted as presented in the Figures, represents cases where channels were weighted with their average cross sections or where all channels carried equal weights respectively.

It can be seen from Fig. 6 (in the case of the (n,tot) cross section) that, the ’best’ file from the combined adjustment using the global likelihood function (EXFOR + hmf57c1) deviates slightly from experimental data between about 7 to 12 MeV; this explains the relatively large $\chi^{2}$ of 10.57 obtained for the (n,tot) cross section in Table 4. A deviation from experimental data is also observed in the (n,non-el) channel as can be seen in Fig. 7. Since the (n,tot) is a sum of the (n,el) with the (n,non-el) cross sections, the deviation observed in the (n,non-el) channel could be a contributory factor to the deviation observed in the (n,tot) channel. It should however be noted that, this interpretation must be made with caution since only one experimental measurement was available for the (n,non-el) channel in EXFOR at the time when this work was carried out. Our ’best’ file however outperforms the other evaluations for the (n,2n) and the (n,inl) cross sections as can be seen from the Fig. 6. The evaluations from JEFF-3.3 and JENDL-4.0 overlap each other for the (n,2n), (n,$\gamma$) and (n,tot) cross sections and therefore the evaluation from JENDL-4.0 has not be shown.

In order to determine how our adjusted files compared with integral experiments, the adjusted files (from this work) were inserted into the ENDF/B-VII.0 library and used for criticality calculations for selected benchmarks. The results obtained were compared with benchmark experimental data (see Table 3) as well as with calculational results obtained using only the ENDF/B-VII.0 ND library and presented in Table 6. ENDF/B-VII.0 was used because it was the library version that came with the version of the MCNPX code used for criticality calculations in this work. In Table 6, the ratios of calculated to experimental values (C/E) for the selected lead sensitive benchmarks are also presented. In the case of criticality calculations for each benchmark as presented in the Table, the the ENDF/B-VII.0 library was maintained as the reference (or base) library for all isotopes while only the nuclear data of 208Pb was varied. In the case of ENDF/B-VII.0 (column 6 of Table) however, the nuclear data for all isotopes were maintained as the ENDF/B-VII.0 nuclear data library. As mentioned earlier, only the hmf57c1 benchmark was used for adjustment in the 2nd update; all other benchmarks as presented in the Table were used for testing and validation purposes. This was done in order that the same benchmark used for adjustment was not also used for testing and validation. Consequently, three different ’best’ files were obtained: (1) adjustment with differential data only (1st update), (2) adjustment with the hmf57c1 benchmark only (2nd update), and (3) the combined adjustments using the combined (global) likelihood function.

From Table 6, it can be observed that the adjustment from the 1st update (weighted channels), performed relatively well when compared with its unweighted for most of the benchmarks. This could be because (as observed earlier in Table 4), relatively smaller $\chi^{2}$ values were obtained (i.e. in the case of the 1st update (weighted channels)) for the (n,tot), (n,inl) and (n,2n) cross sections, which are relatively important cross sections for fast systems. Fast systems such as the hmf57c1 benchmark have small capture reactions and hence, though the (n,$\gamma$) cross section ((in the case of weighted channels) performed badly with respect to the $\chi^{2}$ presented in Table 4, its impact on the calculated $k_{\rm eff}$ was minimal. This explains why even though the best file from the 2nd update (i.e. random file number 1835 with a $k_{\rm eff}$ = 1.00002 for the hmf57c1 benchmark) performed relatively better than the adjustment from the 1st update for most of the benchmarks, the file performed poorly when compared with differential experimental data (see Table 4). From the combined adjustments using the global likelihood function (EXFOR + hmf57c1), it can be observed that, the adjusted file performs relatively well against all the benchmarks as well as against the ENDF/B-VII.0 library. Furthermore, as seen previously from Table 4, the best file from the combined adjustment performed quite well against the other nuclear data libraries for the selected cross sections presented.

As mentioned earlier, in the 2nd update only one benchmark was used for the adjustment, however, a global adjustment with a single benchmark (as carried out in this work), could lead to a global fit to the experimental $k_{\rm eff}$ for the particular benchmark used for adjustment but could have discrepancies when compared with other benchmarks or with experimental data from EXFOR. Therefore, the use of multiple benchmarks (including their correlations) for adjustments has been proposed and is presented in a dedicated paper [54]. Further, a natural extension of this work is to utilized the global (combined) likelihood function computed for the reduction of ND uncertainty in applications as presented in Ref. [11]. Also, the inclusion of model defects in the adjustment procedure is proposed for future work.

4 Conclusion

A method was presented for combining differential and integral benchmark experimental data for nuclear data adjustments and multi-level uncertainty propagation within the TMC method. The method combines individual likelihood functions computed from two Bayesian updates into a global (combined) likelihood function. The proposed method was applied for the adjustment of neutron induced reactions on 208Pb in the fast energy region below 20 MeV. The results from the adjustments were compared with available experimental data from EXFOR, the major nuclear data libraries and against a selected number of criticality benchmarks and found to compare quite favourably. In order to take full advantage of relevant integral experiments available, the use of multiple benchmarks for adjustments (including their correlations) and the inclusion of benchmark calculation uncertainties as well as model defects, is proposed for future work.

ACKNOWLEDGMENTS

We would like to thank Petter Helgesson, Stephan Pomp and Georg Schnabel from Uppsala University, Sweden, and David R. Novog from McMaster University, Canada for insightful comments and discussions on this topic.