Deconvolution methods used for the development of a neutron spectrometer

To fulfill the requirements of the industry, the feasibility of a transportable neutron spectrometer is under study by our collaboration. Preliminary studies have led to a solution based on a unique-multi-detectors-Bonner sphere, which has the advantage to simultaneously perform many neutrons-measurements through thermal neutron detectors placed at different depths inside a polyethylene sphere. The incident neutron energy is then reconstructed using unfolding methods. The optimization of the layout of the detectors in the sphere and of the diameter of the sphere was performed thanks to GEANT4 simulations, coupled to unfolding methods such as least-squares, maximum of entropy or maximum of likelihood. Different unfolding methods have been tested. For the time being, the best unfolding results were obtained by the computer codes M.A.X.E.D. and G.R.A.V.E.L., formalized by the Nuclear Energy Agency (N.E.A.). A “personal” method based on the maximum of entropy coupled to the maximum of likelihood gives good results as well, but a few convergence parameters have still to be optimized. In the present paper, a solution of a multi-detectors-Bonner sphere is presented as well as results of unfolded neutron energy spectra. The results obtained show a good agreement with unmoderated and moderated Am/Be and Cf spectra.


I. INTRODUCTION
CCURATE measurements of neutron dose rates is of great importance for radiation protection applications. To fulfill these objectives, Bonner spheres [1] were developed more than 40 years ago and are still used, due to continuous optimizations of the detectors efficiencies and new patterns of use, from one detector per Bonner sphere at the beginning to multiple detectors at different depths in only one sphere [1,2,3].
This report is about a project in progress, to develop a neutron spectrometer with the following general specifications: it has to be transportable, to return a 0 to 20 MeV neutrons spectrum and be able to measure an ambient dose equivalent (H*(10)) from 1 µSv/h to ~10 mSv/h in less than 5 to 10 minutes.
To fulfill theses specifications, a multi-detectors Bonner sphere setup was chosen. It is composed of a sphere made of a solid polyethylene moderator, with multiple thermal-neutron detectors inside. For this prototype, the 6 Li(n,) 3 H reaction was chosen for the detectors to optimize the neutron gamma discrimination thanks to the high Q of reaction, equal to 4.78 MeV. An important work of optimization has been done with GEANT4 [4,5,6] Monte-Carlo simulations coupled to unfolding methods. This optimization concerns not only the size of the sphere, but also its thermal-neutron-detectors layout, number and geometry. Numerous simulations were performed to get information on the sensitivity of the diameter of the sphere and of the detectors depth and number in order to get an adequate layout compatible with the previously mentioned specifications. Then, different unfolding methods were tested, leading to results, with a good agreement, for unmoderated and moderated Am/Be and 252 Cf neutron energy spectra.
In the following, the geometry of the prototype and the reference spectra will be presented. Then, the unfolding methods used in this project will be exposed. Finally, the results of GEANT4 optimization and of unfolding will be exposed.

II. GEOMETRY OPTIMIZATION
Monte-Carlo simulations with the GEANT4 environment have been processed to optimize the geometry of the multidetectors-Bonner sphere. The version 4.10.5 of GEANT4 was selected, and the QGSP_BIC_HP physics list was loaded to allow a high precision for neutrons energies below 20 MeV [7].

A. Geometry in GEANT4
GEANT4 simulations were performed with Bonner sphere diameters from 20 to 30 centimeters, in order to determine the most appropriate diameter to the specifications. In these spheres made of polyethylene, one thermal-neutrons detector was placed each centimeter along the 3 axis (x, y and z). The detectors were composed of lithium fluoride, LiF, with an isotopic ratio of 14.3% for 6 Li and 85.7% for 7 Li, to get 8.81 atoms/cm 3 of 6 Li for a LiF density of 2.64 g/cm 3 . The dimensions of the detectors were 1.0 cm   cm   cm. Fig. 1 presents an example of the 25 cm diameter sphere, which makes a total of 73 detectors. Fig. 1. Unique-multi-detectors Bonner sphere modelized with GEANT4, with a parallel beam in the direction : (-1 0 0). Here, a 25 cm diameter polyethylene sphere is represented in red, with one LiF detector placed each centimeter along x, y and z axis.

B. Neutron source
The first step of the study was to get a maximum of information on the Bonner sphere with GEANT4. Two types of sources were necessary to reach this goal : ̶ Monoenergetic sources vectors named Xj, where j is the number of the energy bin considered, j∈{0,…,m-1} with m the total number of energy bins, which corresponds to the number of components of the X source vector. These sources were used to obtain the response function of the setup. ̶ Reference sources vectors, corresponding to unmoderated and moderated Am/Be and 252 Cf energy spectra, named XAmBe, XCf, etc. These spectra were used to validate the unfolding. For the moment, only parallel sources were sent to the sphere, with the direction : (-1 0 0). Other directions, parallel or isotropic, will be studied in the future.

1) Monoenergetic sources
Monoenergetic sources of neutrons were sent to the Bonner sphere to build a response matrix of the setup, named R. For monoenergetic-parallel sources, the detector responses will be stored in a rectangular response matrix named Rparallel, with n lines and m columns : the Rij scalar represents the response of the detector or group of equidistant from the center detectors i, i∈{0,…,n-1} to a mono-energetic-parallel source of neutron Xj with the energy bin j, j∈{0,…,m-1}. In our case, 56 monoenergetic parallel beams were simulated with the following energy values: Where j∈{0,…,55} to cover the energy range from 10 -9 MeV to 10 2 MeV. For a given mono-energetic source j, an associated logarithmic bin of energy is then defined, whose bin edges are: So that Ej is the geometrical mean of these two bin edges. A beam surface fluence (Φ) of 1.41×10 5 cm -2 , which corresponds to 10 8 neutrons sent into a 30 cm diameter beam, was used.

2) Reference spectra
Four reference-input-parallel spectra were tested, with the same surface fluence : -Unmoderated and moderated Am/Be (XAmBe and XAmBe_mod). -Unmoderated and moderated 252 Cf (XCf and XCf_mod). The Am/Be source was defined according to the ISO 8529 standard spectrum [8]. The 252 Cf spectrum was defined following a normalized Maxwellian energy distribution [9]: With E the energy in MeV and EM the temperature parameter equal to 1.42 MeV.
The moderated spectra were obtained by GEANT4 simulations of the previous spectra moderated by 6 cm of polyethylene. Fig. 2 shows an example of the reference spectrum, XAmBe_mod, corresponding to the moderated Am/Be source.

III. UNFOLDING METHODS
Whereas the response-matrix building, the spectra initialization and the prototype-geometry optimization were done with the GEANT4 environment, the aim of the unfolding is to retrieve the neutron source energy spectra X by solving the following equation: With X, the unknown source vector (in cm -2 ) to be reconstructed, M the measured vector (no unit) associated to X, and R the matrix corresponding to the response function of the setup (in cm 2 ). In our case, the X vector has 56 components which correspond to 56 groups of energies from 10 -9 MeV to 10 2 MeV. In general, X is normalized to 1, image of the total probability, and then multiplied by the calculated surface fluence Φ. If the unfolding perfectly works, X has to be equal to the energy spectra of the incident source.

A. General principles
The method is composed of two parts : at first, getting M, the measure vector M associated to one of the four reference spectra mentioned above, thanks to GEANT4 simulations. Then, getting X from M, thanks to an iterative unfolding method.
In the first step, the measure vector M has n components corresponding to the n detectors or groups of equidistant from the center detectors. R=Rparallel is the response matrix obtained thanks to the parallel monoenergetic beams, whose dimensions are n×m, with m=56 (number of groups of energies). R is constant, whereas the measure vector M depends on the input energy spectrum X.
In the second step, the aim is to get X from the M vector with an iterative deconvolution method. Indeed, as R is not square in general, it is not possible to simply invert it to get the unknown spectrum X, and even if R is invertible, the slightest disturbance in the system implies unstable solutions, even negative ones for X.

B. Unfolding methods
Three unfolding methods have been experimented in this project: -Two methods, M.A.X.E.D. [10] and G.R.A.V.E.L. [11], formalized by the N.E.A. within the U.M.G. package [12]. M.A.X.E.D. is a computer code based on the maximum of entropy. G.R.A.V.E.L. is based on a nonlinear least-squares method. -One "personal method" inspired by [13,14] and based on the maximum of entropy and likelihood.
Concerning the "personal method", the general principles are based on the maximization of the likelihood (L) and of the entropy (H) simultaneously.
The natural logarithm of the likelihood is described in equation (6), considering a Poisson distribution: With: -the total number of detectors. -the measure of the detector . -〈 〉 the mathematical expectation of .
Then, according to Shannon, the entropy follows the equation (7): the total number of energy groups. -the probability for a neutron to belong to the energy group .
If considering a scalar equation, the differentiation of the function = ln + and the use of the Newton's method imply the following iterative equation: With the iteration number. In this case, there are no convergence issues in general, unless the initial value ( 0 ) is far from the solution.
Though, if considering a matrix equation, the Newton's method leads to the following equation, where x becomes a vector (X), becomes a vector ( ) and ′ becomes a Jacobian matrix (J): This Jacobian matrix, even invertible, can rapidly become non-invertible with the slightest perturbation of one of its components. This issue can be solved by the Cholesky method [13,15]. Moreover, to help the Newton's method to converge, an adjustment of the size of = ( + − ) can be necessary, by multiplying by a factor after each iteration, thanks to Armijo's method [13,16]. Despite Cholesky and Armijo's methods, optimizations still remain to be applied to optimize the "personal" method, the results obtained are in good agreement with the source spectra.
The results obtained with these three methods will be discussed in the following part.

IV. RESULTS
The result which have led to the geometry optimization will be presented. Then, the results of unfolding will be exposed, in a first step with high statistics (i.e. a beam of 10 8 neutrons) and then with lower statistics (i.e. a beam of 10 6 neutrons) to validate the performances of the setup.

A. Geometry optimization
The first step of the optimization of the system consists in optimizing both the size of the polyethylene sphere and the number and layout of the detectors. The diameter of the polyethylene sphere has to be a compromise between the capability to moderate the neutrons and the weight of the spectrometer. The mean free path of a 20 MeV neutron in a moderator like polyethylene is about 25 centimeters, by using the Fermi age method [17]. This length should approximately correspond to the most appropriate diameter of the polyethylene sphere. To confirm this rough estimate, GEANT4 simulations were performed, with 56 groups of energies, for 3 different Bonner sphere diameters : 20, 25 and 30 centimeters, with one thermal neutrons detector each centimeter along x, y and z, in accordance with Fig. 1.
Then, the G.R.A.V.E.L. unfolding code was used on the 3 response matrices related to these 3 different diameters. Fig. 3 presents an example of the results obtained, for the three diameters, for a moderated Am/Be source. For a 20 centimeters diameter, the low energies components seem to be correctly returned. But the high energies slightly underestimated, probably due to the small sphere diameter, not sufficient to correctly moderate the highest energy neutrons. For a 25 or a 30 centimeters diameter, high energies are better reconstructed, despite the fact that the energy of the moderated component was underestimated. In any case, for the calculus of the ambient dose equivalent H*(10), to have an accurate estimation of the high energy values is very important because the ambient dose equivalent per unit neutron fluence coefficients ℎ Φ * are much higher above 0.1 MeV, according to the I.C.R.U. 66 [18]. Consequently, to minimize the weight of the spectrometer, a 25 centimeters diameter sphere has been retained in the following.
The second step is the optimization of the detectors layout within the sphere. The analysis of the response matrix Rparallel obtained for the 56 monoenergetic beams with detectors every centimeters is represented in Fig. 4 for a selection of energies. This figure emphasized that close to the surface (~10 to 12 cm to the center), the variation of the curves is more pronounced than next to the center. Consequently, a more important detector granularity is required close to the surface than in the depth of the sphere. An optimization, by removing some detectors, particularly next to the center, leads to an optimized LiF detectors layout based on 6 depths : 1, 6, 8, 10, 11 and 11.5 cm from the center. With this new layout, thermal-neutrons detectors are concentrated next to the surface.  This layout is the result of a first optimization and will may be subject to more refinements in the future [19].

B. Unfolding results
The aims of these studies are to confirm the optimized geometry and to study the most appropriate unfolding method for different energy spectra (moderated or not) and neutron flux.

1) With high statistics
The first unfolding tests were done with a surface fluence of 1.41×10 5 neutrons per square centimeter, which approximately corresponds to an ambient dose equivalent of about 50 µSv for Am/Be or 252 Cf sources.
For each unfolding method, the first iteration has to be initialized with a guess spectrum. Three kinds of guess spectra were tested : a flat one, an Am/Be or 252 Cf spectrum and a moderate Am/Be or 252 Cf spectrum. Fig. 6 presents the results obtained with a flat guess spectrum. Fig. 6 (a) shows the unfolding results for an Am/Be source and Fig. 6 (b) for a moderated Am/Be source. The blue dotted line represents the source energy spectrum. An offset in the high energies with the 3 tested methods can be observed for both sources. This offset is much more significant for a moderated source. This offset will lead to an overestimation of the ambient dose equivalent (H*(10)). Thus, a flat guess spectrum does not seem to be adapted to initialize the unfolding.   Fig. 7 (a), where the unfolded spectra and the source spectrum are quite equivalent. With the personal method, the mean value is correct, but not the shape of the spectrum, showing that further optimizations in this algorithm are to be performed as expected. In Fig. 7 (b), the results totally diverge for M.A.X.E.D. and G.R.A.V.E.L., because the firstiteration-spectrum components are equal to zero for the thermal and epithermal energies. In the personal method, first-iteration components are never equal to zero, they are at least equal to an infinitesimal value, that gives more degrees of freedom to the algorithm and leads to a more appropriate result.
Thus, if a thermal component is present in the source spectrum, which will be usually the case in the industrial environment, a thermal component has probably to be introduced in the guess spectrum to reconstruct thermalized sources.   As a partial conclusion, when the awaited spectrum is totally unknown, it is better to have a moderated guess spectrum, which covers a wide range of energies and confers to the algorithm more degrees of freedom. Indeed, good results were obtained, regardless of the source spectra. Similar results were obtained for 252 Cf and moderated 252 Cf sources.

2) With lower statistics
The previous results with high statistics correspond to the following ambient dose equivalent : H*(10) ~ 50 µSv, which can be converted into the following ambient dose equivalent rate : Ḣ*(10)) ~ 300 µSv/h measured in 10 minutes. The aim in this subsection is to test the unfolding methods with one hundred times less statistics (H*(10) ~ 0.5 µSv), in order to validate our setup at lower dose rates (Ḣ*(10) ~ 3 µSv/h in 10 minutes), to approach the specifications discussed in the introduction. Fig. 9 presents an example of the results obtained for a moderated Am/Be source and an unfolding initialization with a moderated-guess 252 Cf spectrum. With a dose divided by one hundred, the results still show a good agreement. Similar results were obtained for moderated or not Am/Be and 252 Cf sources.

3) Conversion into ambient dose equivalent (H*(10))
One of the main objectives of the ratemeter prototype in development is to be able to measure an ambient dose equivalent rate. To check the accuracy of the results obtained with the different unfolding methods, the neutron energy spectra will now be converted into dose equivalent. This was achieved by applying to the previous spectra the ambient dose equivalent per unit neutron fluence coefficients (ℎ Φ * ) from the I.C.R.U. 66 [18]. Fig. 10 presents an example of the results obtained when converting the spectra of the Fig. 8 (b) and Fig. 9 into ambient dose equivalent H*(10). The results obtained for a moderated Am/Be source unfolded with a moderated-guess 252 Cf spectrum, for 2 different doses are presented : in (a), with high statistics, i.e. 1.41×10 5 cm -2 or a dose rate (Ḣ*(10)) of approximately 300 µSv/h measured in 10 minutes, and in (b), with one hundred times lower statistics, i.e. 1.41×10 3 cm -2 or a dose rate (Ḣ*(10)) of approximately 3 µSv/h measured in 10 minutes. The doses obtained after deconvolution are in good agreement with the doses calculated for the source spectra. Errors, below 15% from the reference doses are obtained even in the worst case, i.e. for low statistics and an input spectrum which is not the one of the source. This accuracy corresponds to the one obtained usually with Bonner spheres [19], validating the detector setup and the optimized detectors layout.

V. CONCLUSIONS
Monte-Carlo simulation with GEANT4 associated to unfolding methods of U.M.G. 3.3 and to one "personal" method have been processed, in order to optimize the geometry of a neutron dose ratemeter based on a multi-detectors-Bonner sphere. Results obtained are in good agreement with initial sources energy spectra, thanks to the optimization of the sphere with a diameter of 25 centimeters and the detectors layout composed of 6 groups of detectors.
Concerning the unfolding initialization, guess spectra have to be pre-defined with the following components based on real spectra: -A rapid component.
-An epithermal component, at least equal to non-zero values. -A thermal component. The adequacy of the guess spectra could be determined by calculating the 2 per degree of freedom after the first iteration.
In this work, 4 simulated energy spectra of moderated and unmoderated Am/Be and 252 Cf sources have been tested, with high and low statistics. Results have shown that even with low statistics, the errors are less than 15% on the dose calculation.
This study was performed with parallel-monoenergetic beams in only one direction (-1 0 0). Thus, two other types of sources will have to be tested: -Parallel sources with a 45 degrees-shift-angle from the initial beam (-1 -1 0). -Isotropic sources.