A new unfolding code combining maximum entropy and maximum likelihood for neutron spectrum measurement

Abstract

We present a new spectrum unfolding code, the Maximum Entropy and Maximum Likelihood Unfolding Code (MEALU), based on the maximum likelihood method combined with the maximum entropy method, which can determine a neutron spectrum without requiring an initial guess spectrum. The Normal or Poisson distributions can be used for the statistical distribution. MEALU can treat full covariance data for a measured detector response and response function. The algorithm was verified through an analysis of mock-up data and its performance was checked by applying it to measured data. The results for measured data from the Joyo experimental fast reactor were also compared with those obtained by the conventional J-log method for neutron spectrum adjustment. It was found that MEALU has potential advantages over conventional methods with regard to preparation of a priori information and uncertainty estimation.

Keywords:

• neutron spectrum unfolding
• maximum entropy method
• maximum likelihood method
• underdetermined problem
• multi foil activation method

1. Introduction

The inverse problem of determining radiation source information from measured detector readings has been investigated by many researchers [1–5]. However, it remains unclear which methods achieve high accuracy for unfolding results, especially in underdetermined problems. In addition, most spectra unfolding codes require an a priori spectrum (i.e., a guess) to start the unfolding procedure for an unknown spectrum. The accuracy of the resulting spectrum strongly depends on the subjectively chosen guess spectrum.

A new spectrum unfolding code was developed, based on a combination of the maximum likelihood and maximum entropy methods. The code is called the Maximum Entropy and Maximum Likelihood Unfolding Code (MEALU). A critical advantage of the MEALU method is that it allows evaluation of the uncertainty and a determination of the neutron spectrum without an initial guess spectrum. Maximum likelihood combined with maximum entropy has been previously applied to the unfolding of neutron spectra, notably in the study of Itoh and Tsunoda [1]. However, the MEALU formulation is a different way of combing or using the maximum likelihood and maximum entropy methods from previous applications [1,5–6] in the way that it combines maximum likelihood and maximum entropy.

It is important to estimate the uncertainty for an unfolded neutron spectrum as well as the integral quantity. To improve the sensitivity and propagation of uncertainties of the MEALU solution from a previous study [7], we used analytical estimation of uncertainty propagation to consider the measurement uncertainty and the response function uncertainty. We checked the performance of the MEALU technique by analyzing mock-up data. The results from measured data were compared with those obtained by a conventional method for neutron spectrum unfolding.

2. Neutron spectrum measurement based on multi foil activation method

The M kinds of detector response c_i are related to the irradiated neutron spectrum φ(E) by a Fredholm integral equation of the first kind:

where R_i(E) is the response function of the i-th detector. The problem for neutron spectrum unfolding is to estimate φ(E) satisfying Equation (1) with the measured data c_i and the well evaluated data of R_i (E) within a reasonable uncertainty. The usual approach is to break the energy E into discrete intervals of N groups and rewrite Equation (1) in sum (or matrix) notation as

where R_ik is the average value of R_i (E) and φ_k is the flux in the k-th energy interval group. In general, since the number N of energy groups is larger than the number M of measured detector responses, the simultaneous linear Equation (2) results in an indefinite problem mathematically in the case of the foil activation method and multisphere neutron spectrometer (Bonner Ball).

In the underdetermined problem, the unfolding code and procedure have the following requirements:

•	Stability of solution (reducing the effect of measurement uncertainty and the response function uncertainty)
•	Uniqueness of solution (only one valid solution)
•	Rejection of unphysical solutions (guarantee of a positive solution)
•	Uncertainty estimation (uncertainty propagation from measurement, response function and numerical analysis).

To solve such problems, many researchers have investigated and developed numerous unfolding codes such as STAY’SL [2], FERRET [3], NEUPAC [4] and MAXED [5]. These codes, however, require the preparation of a group-wise initial guess spectrum with their covariance matrix as a priori input information, and sometimes involve very complicated calculations to obtain physically reasonable solutions when dealing with ill-conditioned matrices and/or assuming nonlinear probability functions such as lognormal distributions.

3. The MEALU algorithm

3.1. Definitions

Taking the expectation of both sides of Equation (2) and denoting the total number of neutrons incident on the detector by α, we have

where

The expectation of a variable q is denoted by 〈q〉.

3.2. Derivation of new method

For the vector p being a probability, the quantity

is known as Shannon's information entropy. Based on the maximum entropy method, we choose the energy distribution p which maximizes Equation (8).

On the other hand, the logarithmic likelihood relevant to the Poisson statistics of neutron detection [6] is given by the following:

In this study, we select p that maximizes the linear combination of the entropy and the likelihood. The combination of entropy and likelihood, Q, is given by

This study adopts other formulations to improve the reproducibility of unfolded neutron spectrum comparisons with conventional algorithms using maximum entropy, such as Itoh and Tsunoda's [1]. They used maximum entropy for the function formulation of p, and used the maximum likelihood method for the determination of the Lagrange multiplier in p. In our proposed method, the entropy and the likelihood were maximized together, that is the adoption of the maximum entropy method for regularization of the maximum likelihood equation.

To maximize Equation (10), Q is differentiated with respect to p_k , and the result is set equal to zero to give the following:

We define λ as

Finally, we obtain the set of (N+M) equations to be solved in terms of λ_i using Equation (7):

According to Equation (12), the number of λ _i matches the number of c_i for the underdetermined and overdetermined problems. Although this paper focuses on the underdetermined problem, the MEALU method can also be applied to the overdetermined problem without reformulation.

The problem now is to solve the nonlinear system (12) for the vector 〈c〉 included in λ. For this we apply Newton's method. The unfolded spectrum under the present method is not affected by an arbitrary choice of initial guess, but a starting point of λ is required in the iterative scheme. We employed all zeros for this starting value, so that the final result is independent of initial lambda.

The solution ⟨φ⟩ is calculated from the product of p and α:

On the other hand, the logarithmic likelihood relevant to the normal distribution statistics is given by the following:

where Σ_c and Σ_〈c〉 are the covariance matrices of c and 〈c〉, respectively. In this case, Equation (12) takes the following form:

By formula (17), since the number of the elements of lambda and the number of the elements of c_i are always in agreement, p becomes settled by formula (13) irrespective of the underdetermined or overdetermined problem. Moreover, it guarantees that p always becomes positive irrespective of the positivity/negativity of lambda.

Although this algorithm is not formulated for using an initial guess spectrum as a priori information, it is easy to consider a priori information for this algorithm. In this case, Equations (8), (13) and (14) are modified as follows:

3.3. Preparation of total number of neutrons

Determining the absolute value of the neutron spectrum is an important task in an unfolding problem.

Although the total neutron number α should be determined at the end of the experiment, it was not possible for us to know its true value. However, α can be measured when a detector for counting the total number of neutrons can be used.

In unfolding methods and codes that require an initial guess spectrum, the following two techniques are adopted:

1.	an absolute neutron spectrum is input as a priori information, equivalent to inputting α; or
2.	a relative neutron spectrum is input for the spectrum shape and α is determined by the least squares method as the maximum likelihood value.

In both cases, the uncertainty in α is evaluated from the covariance matrices of initial guess spectrum and measured response.

Since α appears in both λ and λ′, it is an important parameter for both scaling and determining the spectral shape in this algorithm. In the MEALU algorithm, α is input with the uncertainty as an input parameter. In other words, a priori information is required for MEALU. However, it is more convenient to prepare α than estimate the neutron spectral shape.

If there is no a priori information on α, a flat spectrum is adopted as the starting point of the nonlinear solution method in MEALU algorithm. For this technique, the following two methods are proposed for determining α:

1.	α is input with its uncertainty; or
2.	a flat spectrum is assumed and α is determined by the least squares method as the maximum likelihood value (if necessary, the MEALU algorithm including α is iterated to maximum H+lnL).

3.4. Uncertainty estimation

A previous study [7] applied the Monte Carlo method to estimate the uncertainty propagation. In this study, we adopt an analytical estimation of uncertainty propagation. For sensitivity analysis and propagation of uncertainties of the MEALU solutions, we have considered the effect of variations δc_i in the measurements, δR_ik in the response function and δα in the total number of neutrons incident on the detector. If the maximum entropy and maximum likelihood solution exist, it can be shown that we can find explicit solutions for the matrices δφ_k /δc_i , δφ_k /δR_ik and δφ_k /δα. We can then introduce the N × N total uncertainty matrix U, given by

provided the covariance matrices K, K ′ and K ^″ are available. We make the usual assumption that K is an M × M diagonal matrix with elements δc_i ² on the diagonal. For the N × N covariance matrix K ′, it is useful to consider the full covariance matrix in each detector response. K^″ is the uncertainty of α. The matrix U can be used to assign an uncertainty to any integral quantity H of the form H = Σ_ih_kφ_k . The uncertainty associated with H is given by

4. Performance test

We created a prototype MEALU program and checked its performance through analyses of mock-up data for verification and real measured data for validation for the application to real problems.

4.1. Mock up data

Inconsistencies between a measured detector response and a prepared response function make it difficult to evaluate the characteristics of an unfolding algorithm. To avoid this situation, a problem with five group neutron fluxes unfolded from three detector responses was prepared as shown in Table 1 . In this test, the total neutron flux α was assumed as a given parameter. It was assumed that the response function had no uncertainty. For the probability distribution for likelihood, a normal distribution was used. The true reaction rate data with no artificial noise was used as an input. The uncertainty value for the reaction rate input was set to 5%.

Table 1 Parameters for test problem

	Response Function
k	R 1,k	R 2,k	R 3,k	φk,exact
1	3	0	0	2
2	4	2	1	5
3	3	3	2	9
4	2	3	4	13
5	1	0	4	8
ci	87	76	107	37

4.2. Real measured data

The sample problem considered was dosimetry data analysis in the Joyo experimental fast reactor. Details of the Joyo dosimetry experiments are described elsewhere [8]. The unfolding calculations were carried out using five kinds of reaction rate data, where the input neutron cross section as a response function prepared by the NJOY code from JENDL Dosimetry File 99 [9] is given in 23 energy group numbers. The uncertainty value for the reaction rate input was set to the measurement uncertainty as a diagonal matrix, which consists mainly of calibration uncertainties of the γ-ray detector. The uncertainty of α was set to the square of 3.47% derived from the uncertainty of the reactor thermal output, which is equal to the total flux uncertainty. Initial α is calculated from the ²³⁵U fission rate dosimeter using a flat spectrum and final α is determined by iteration of the MEALU algorithm to maximize the probability of solution. In this problem, it is impossible to know the true results for an absolute neutron spectrum. The results of MEALU were compared with the conventional J-log method [4] to validate the algorithm, where the input initial guess spectrum was calculated by the MCNP code using the cross section based on JENDL-3.2 in 103 energy group numbers.

5. Results and discussion

5.1. Mock up data

Table 2 and Figure 1 show the unfolded spectrum with uncertainty estimation obtained from the MEALU prototype code. The results have no negative component, as predicted theoretically. The MEALU unfolding results are in good agreement with the true spectral shape and magnitude without an initial guess spectrum. The original spectrum was reproduced within 1σ at each energy bin. The detector response using the unfolded spectrum and results of χ ² tests are shown in Tables 3 and 4 . The re-folded detector responses using an unfolded spectrum also agree with the exact values.

Table 2 Result of uncertainty estimation by measurement data

k	φ unfold	Δφunfold (1σ)	Δφunfold (%)	φexact
1	2.5	2.2	85.8	2.0
2	4.9	4.5	91.8	5.0
3	8.1	3.6	45.0	9.0
4	13.9	1.3	9.1	13.0
5	7.5	1.3	17.8	8.0

Table 3 Unfolded response

i	c exact	c unfold
1	87.0	87.1
2	76.0	75.9
3	107.0	106.8

Table 4 Probability

χ^2 value	9.78 ×10^−4
Probability	100%

Figure 1 Unfolded spectrum

5.2. Real measured data

A comparison of the unfolded neutron spectrum with the uncertainty obtained by the MEALU and J-log methods from the same measured reaction rate data set is shown in Figure 2 , where 90% response range of each reaction is also presented. The 90% response range indicates the energy range containing 90% of the reaction rate energy integration. In this irradiation field, it is well known that the spectral shape is typically that of a fast reactor spectrum, and indeed the spectra are well reproduced by such a spectrum. Both results are in good quantitative agreement within the estimated uncertainties in the energy region of the 90% response range. In the energy region without sensitivity, only the J-log method shows a spectral shape due to the initial guess spectrum. Measurement of a hard spectrum such as that from a fast reactor fuel region requires a reaction rate with the 90% response range below 10⁴ eV.

Figure 2 Comparison of neutron spectra

Tables 5 – 7 show comparisons of the integral values and χ ² test. Both results are in good agreement quantitatively within the estimated uncertainties. The results of the χ² test show comparable unfolding between the MEALU and J-log methods.

Table 6 Comparison of reaction rate

	Neutron flux (n/cm^2/s)
Reaction	Measured	MEALU	J-log
^237Np(n,fission)	2.44 × 10^15	2.47 × 10^15	2.44 × 10^15
^235U(n,fission)	8.38 × 10^15	7.80 × 10^15	8.28 × 10^15
^45Sc(n,g)	1.46 × 10^14	1.46 × 10^14	1.49 × 10^14
^46Ti(n,p)	6.21 × 10^12	6.21 × 10^12	6.28 × 10^12
^58Ni(n,p)	7.07 × 10^13	7.07 × 10^13	7.15 × 10^13

Table 5 Comparison of integral value

	Neutron flux (n/cm^2/s)
Value	MEALU	J-log method
φtotal	(5.39 ± 0.21) × 10^15	(5.05 ± 0.31) × 10^15
φE>0.1MeV	(4.17 ± 0.42) × 10^15	(3.70 ± 0.31) × 10^15
φE>1.0MeV	(8.39 ± 1.77) × 10^14	(8.55 ± 0.91) × 10^14

Table 7 Comparison of probability

Value	MEALU	J-log method
χ ^2 value	0.70	0.16
Probability	95.2%	99.5%

6. Conclusions

A new spectrum unfolding code, MEALU, is currently under development based on a combination of the maximum likelihood and maximum entropy methods. We have demonstrated typical unfolding results of the MEALU algorithm using mock-up data and real measured data. The MEALU algorithm is found to accurately reproduce the neutron spectrum and agrees with the results of a conventional method within the estimated uncertainties. Compared with conventional methods, the new method has the following potential advantages:

•	It is not necessary to prepare an initial guess spectrum and covariance data as a priori information.
•	A positive solution is guaranteed.
•	Normal and Poisson distributions can be used for the statistical distribution.
•	The MEALU algorithm can treat full covariance data for a measured detector response and the response function for uncertainty estimation.

The MEALU algorithm can be applied to large problems with approximately 100 energy bins and other neutron fields with different spectral shapes. Although this algorithm was formulated for the multi activation foil method, it should be applicable to other detectors such as the Bonner counter system and the NE213 recoil proton detector system. It is also expected that MEALU could be used in other unfolding problems in nuclear engineering.

Appendix

A.1 Likelihood function based on Poisson distribution

It is common practise to use a Bonner ball for neutron measurements. Since the number of counts in each Bonner sphere follows the Poisson distribution, the probability P_i (c_i ) that we have counts c_i (i = 1, …,M) in the i-th Bonner sphere is given by

where λ_i is the expected value of c_i:

The expected value of a random variable x will be written as 〈x〉. In Equation (A.2), φ_j (j = 1, …,N) is the number of neutrons incident on the detector belonging to the j-th energy group, i.e.

and R_ij is the probability that a neutron of the j-th group produced a count in the i-th Bonner sphere. Then the joint probability distribution function of the c's, i.e. the likelihood function, is

by the mutual independence of c's.

The maximum likelihood method is a method to find φ with maximum L. At the point of maximum L for φ, the differential coefficient of the logarithmic likelihood function ℓnL also becomes zero. Therefore, we adopted the following function:

A.2 Likelihood function based on normal distribution

The probability density function of the normal distribution N(λ, σ) is

The logarithmic likelihood function ℓnL is