© P. Humbert and B. Méchitoua, published by EDP Sciences, 2017

1 Introduction

Real-time applications require fast and accurate calculation of the detected gamma-ray spectra produced by shielded sources. For this purpose, the SGaRD (Spectroscopy, Gamma rays, Rapid, Deterministic) code [1] that was used to calculate the leakage spectra of one-dimensional spherical assemblies has been updated in order to take into account the response function of various types of detectors and for the identification of gamma-shielded sources geometric characteristics.

SGaRD has two different transport solvers. The first one is a multigroup discrete ordinates S_N solver [2] for the integro-differential transport equation in one-dimensional spherical geometries using the spherical coordinates (r,μ). The second one is a “Method of Characteristics” or ray-tracing solver [2] for the integral transport equation along straight lines through the spherical geometry. The first solver noted S_N in the following is used to calculate the scattered component of the gamma leakage, and the second component is used to calculate the uncollided leakage of gamma lines.

In the first part, we recall the methods used in SGaRD to calculate the leakage spectra. We describe the precalculation of detector response functions (DRFs) using the Monte Carlo code MCNP5 [3], and we show some numerical results concerning the simulation of a natural uranium ball coated with plexiglass and measured using a NaI detector.

In the second part, we present the identification of nuclear radiation-source characteristics, which is a subject of interest for nonproliferation and nuclear safeguard applications. Gamma spectroscopy is used because of the sensitivity of these measurements to source isotopic composition and shielding materials properties. The determination of source characteristics using known signature measurements is an inverse transport problem. This subject has been studied by different authors (see for example [4–8]). In this paper, we present the determination of the unknown material interface positions from the measured uncollided gamma line spectrum obtained by processing high-precision gamma spectroscopy measurements. For this purpose, SGaRD code is used as a forward solver for iterative inverse transport calculations. The material thicknesses are computed using a nonlinear multidimensional iterative optimization algorithm that minimizes the error metric between calculated and measured signatures. The optimization is performed using the gradient-free Powell method [9,10]. For verification, numerical results are presented. A synthetic gamma lines spectrum is used as input to the inverse transport solver, and the obtained geometry is compared to the original one.

2 Pulse height spectra simulation

2.1 Gamma-ray spectra

The gamma-ray spectrum simulation is the result of successive steps. The first step is the computation of the primary gamma source emission rate and spectrum. The second step is the photon transport through the shielding materials up to the external surface of the detector in order to derive the leakage spectrum. The third step is the photon transport into the detector used to evaluate the distribution of the photon energy deposited inside the detector's sensitive volume.

When the source can be modeled as a set of concentric one-dimensional spherical shells, the second step can be handled very effectively using deterministic solvers [1,4]. The first and third steps can be performed using precalculated sources and DRF.

The SGaRD input/output flow diagram for gamma spectra simulation is presented in Figure 1.

Fig. 1 SGaRD input/output flow diagram for gamma spectra calculation.

2.1.1 Primary source

The primary photon source has several components: the gamma lines due to the radioactive decay of unstable isotopes, the gamma resulting from spontaneous fission, the bremsstrahlung radiation produced by charged particles slowing down and the neutron-induced gamma production.

The radioactive decay source line spectrum and emission rate are precalculated using DARWIN code [11]. The treatment with SGaRD of the other components of the primary source terms are under study. The contribution of these components may not be negligible. In particular, the continuum spectrum in depleted uranium is dominated by the bremsstrahlung radiations as shown in [4]. This component will be included in SGaRD for future comparisons with measured gamma spectra.

2.1.2 Leakage spectra

The leakage L is the number of particles leaving the external surface of the source per unit time and solid angle. (1)

The angular flux is solution of the transport equation. It has a discrete energy part due to the uncollided transport of gamma lines and a continuous part due to the scattering. The discrete component is calculated using a 1-D ray-tracing transport solver, presented in Section 3.1.2. The continuous component is transported using a 1-D multigroup discrete ordinates (S_N) solver.

The flux computation is performed using three transport calculations.

The first step is the uncollided transport of each source line with energy E using the ray-tracing algorithm with typically N = 64 directions. (2)

The second step is the total multigroup S_N transport taking into account the scattering term. The multigroup calculations are performed using the conventional discrete ordinates method in spherical curvilinear coordinates with Legendre polynomial expansion of the scattering source. The angular quadrature is an evenly spaced discretization of the direction cosine μ (equiweight quadrature), and the number of directions is typically N = 16. (3)

The multigroup source is obtained by condensing the gamma lines into groups (cf. Fig. 2). The source intensity in a given group is equal to the sum of the intensities of all the lines whose energies are included in the group energy boundaries.

The third step is the uncollided multigroup S_N transport of the multigroup source. (4)

The scattered component is given by subtracting the uncollided multigroup flux to the total multigroup flux. (5)

In Figure 3, the uncollided part has been added to the scattered background to produce the total multigroup leakage. The blue/red curve is the total leakage. The blue part of the curve is the contribution of the uncollided leakage, which has been added to the scattered background in red.

Finally, the total leakage spectrum is the superposition of the discrete uncollided and multigroup scattered components on a very fine multigroup grid (cf. Fig. 4).

Fig. 2 Gamma source lines (blue) are condensed into groups (red) for multigroup SN calculations.

Fig. 3 Multigroup SN leakage spectra for both scattered background (red) and uncollided component (blue).

Fig. 4 Discrete uncollided lines (blue) with multigroup scattered background (red).

2.1.3 Detector response function

The DRF is used to convert the leakage spectrum into a pulse height spectrum using equation (6) assuming that the pulse height and leakage spectra have the same energy discretization into n channels.

Notations

• Δt: Duration of the measurement.

• ΔΩ: Solid angle of the detector viewed from the center of the source.

• L_j: Leakage spectrum = number of gamma particles leaving the source with energy in channel j per unit time and solid angle.

• N_i: Pulse height spectrum = number of counts in channel i.

• R_ij: Detector response matrix = number of counts in channel i due to one gamma entering the detector with energy in channel j.

(6)

The response matrix is obtained by Monte Carlo simulations of the energy deposition within the detector using the same methodology as in [12]. A set of gamma ray incident energies E_j is chosen. For each of these energies, the detector response to a parallel beam normal to the detector is calculated using MCNP5. The DRF for a given energy channel is obtained by linear interpolation between the two nearest calculated responses.

The experimental spectra have a Gaussian distribution shape for the photons energy lines. We take this effect into account by modifying the MCNP5 simulation results. We use a fitting technique to approximate the resolution of the detector which is an experimentally measured data.

With MCNP5, we use an “FT8 Gaussian Energy Broadening” card in order to simulate the full width at half maximum (FWHM) around the peak.

The response matrix can then be used in conjunction with the SGaRD 1-D leakage using equation (6) to determine the number of counts in each detector channel. An example of NaI DRFs calculated with MCNP5 is presented in Figure 7.

Fig. 7 Scintillation NaI detector response functions. The different colors correspond to different incident energies of the gamma particles entering the detector.

2.2 Application to a scintillation NaI detector

As an example, we present the simulation of a pulse height spectra obtained with a scintillation sodium iodide detector and a spherical source made of a natural uranium ball surrounded by a plexiglass shell as shown in Figure 5.

The example geometry is detailed in Table 1. The distance between the center of the source and the detector is 10 cm. Even if this geometry is two-dimensional, SGaRD calculations stay one-dimensional because the detector is taken into account by using the precalculated DRF.

The leakage spectrum of gamma particles leaving the source is calculated using SGaRD. As shown in Figure 6, there is a good agreement with the corresponding MCNP5 calculation.

The pulse height spectrum is obtained by folding the SGaRD leakage spectrum with the DRF calculated with MCNP5 using the pulse height tally and Gaussian broadening to take into account the detector resolution.

As explained in [12], a nonlinear function is applied to estimate the values of the triplet (a, b, c) used to fit the FWHM in function of the photon energy.

Let where E is the incident gamma rays energy in MeV, and a, b and c are constants derived from FWHM measurements.

The triplet used is: (a, b, c) = (−0.002, 0.05, 2.86). These values are typical of a NaI detector and are close to those given in [12]. A set of response functions is presented in Figure 7 for different incident energy of the gamma particles entering the NaI detector. The incident energy range is between 40 keV and 1.5 MeV, and the energy intervals extend from 50 keV at low energy up to a few 100 keV above 1 MeV. The DRFs which are between two calculated ones are interpolated.

For verification purpose, the same calculation was performed using MCNP5 only. The bremsstrahlung gamma source is not simulated in the MCNP simulation in order to stay consistent with the SGaRD calculation.

Both spectra are presented in Figure 8. They show a very good agreement, although the SGaRD calculation is much faster and less prone to statistical artifacts. The computer time for SGaRD is of several seconds compared to 1 h for the MCNP simulation with 10⁹ particles studied.

Fig. 5 Detector and source geometry.

Fig. 6 Source leakage spectra calculated using MCNP5 and SGaRD.

Fig. 8 Pulse height spectra calculated using MCNP5 and SGaRD.

3 Material thickness identification

The identification of the source geometrical characteristics is performed by inverse transport using an optimization method which iterates on the ray-tracing simulations of the gamma lines leakage.

3.1 Uncollided leakage gamma current calculation

3.1.1 Gamma leakage

The observable is the uncollided gamma ray leakage line spectrum, produced by the radioactive decay gamma source.

The number L of gamma particles leaving a sphere of radius R per unit time and solid angle without collision is calculated using the uncollided angular flux ψ(R,μ). (7)

3.1.2 Ray tracing

The direct flux is solution of the transport equation without scattering. It is solved by SGaRD code using an accurate and fast ray-tracing method, also called Method of Characteristics (MoC [2]). When the source is a step function as in our application, this method gives the exact solution of the integral transport equation along discrete directions.

The shielded source is made of concentric spherical shells, each with a constant source intensity Q_i and total macroscopic cross section σ_i. The outgoing flux at the external surface boundary is discretized in N directions. The angular discretization is illustrated in Figure 9. A set of N_i characteristic rays is associated to each shell, from the most internal (i = 1) to the most external shell with a nonnull source (i = I).

In Figure 9, the lines are the boundaries of the angular discretization. The arrows show the directions of neutron travel. The dotted lines are the upper angular boundaries associated with each material shell. There is a different color for each material shell. The γ source is in the shaded shell. The direction of flight of the rays is characterized by the cosine of the angle θ between this direction and the spatial radial vector on the external surface with radius R. The maximum value of θ corresponding to a nonzero outgoing flux is θ_max.

The unscattered leakage is given by (8) where μ_ij characterizes the direction of flight associated with the jth ray crossing the ith shell.

μ_i is the cosine associated with the upper boundary of shell i and Δμ_i = μ_i−μ_i−1.

Each shell is evenly discretized, and the rays within a given shell have the same weight. The number of rays in a shell is proportional to its thickness with at least one ray by shell.

The angular flux ψ is calculated using the following transmission equation, taking into account a constant source approximation: (9)

The λ_i are the intersection lengths of the characteristics with the spherical shells. (10)

Fig. 9 1-D spherical ray tracing. The γ source is in the shaded shell.

3.1.3 Validation

The particle leakage from a homogeneous spherical source of radius R, opacity σ and intensity Q has an analytic expression [13]. (11)

The convergence of the ray-tracing algorithm with the number of directions is compared to the discrete ordinates S_N method in Table 2. The calculations use the following values: R = 5 cm, σ = 1.484673 cm⁻¹ and Q = 1948.531 cm⁻³/s.

The discrete ordinates S_N and ray-tracing calculations are both performed using SGaRD. The S_N solver is a one-dimensional spherical finite difference solver, using an equiweight angular quadrature and a 100-cell spatial discretization. Considering the ray tracing, there is no spatial discretization; one spherical shell is associated with each material region.

As shown in Table 2 and Figure 10, the ray-tracing method is faster and more accurate than the S_N discrete ordinates for unscattered transport.

Fig. 10 Convergence to the analytic solution of the ray tracing and SN leakage calculations with the number of directions.

3.2 Optimization − inverse transport

The goal of the inverse transport problem is to find the 1-D spherical assembly of source and shielding materials that produces the better agreement with the measured uncollided leakage gamma spectrum. Starting from a model composed of a given number of material atomic compositions and densities, the optimization algorithm automatically searches the shell thicknesses that give the best estimation of the measured line spectrum. The goal is then to minimize the error or distance between the measured and calculated spectrum. The quality of the estimation is quantified by the chosen error metric.

3.2.1 Error metric

We use a least square method, the error is defined as: (12)

• : Measured leakage in line i,

• : Calculated leakage in line i,

• M : Number of material shells,

• N : Number of spectral lines,

• x_m : Thickness of shell m,

• ω_i : Weight of line i.

Geometric or mass constraints can be incorporated in the error function. For instance, one can favor the external radius around R = R₀ using the following error function, where λ is the weight associated with the constraint. (13)

By default, the weight of the lines is all equal to one; however, it is possible to increase or decrease the importance of a line by using different weights.

3.2.2 Stochastically restarted Powell optimization method

The Powell method is a multidimensional deterministic optimization algorithm, which proceeds by a sequence of line minimizations along with well-chosen directions.

One advantage of the Powell method is that it needs to evaluate only the error function and does not necessitate the calculation of the gradient of the error with respect to the optimized variables for the choice of the successive minimization directions.

The Powell method discarding the direction of largest decrease, described in [9], is used in order to avoid the generation of linearly dependent set of directions.

The line minimizations are performed by successive bracketing of the minimum using Brent's method [10] which switches between inverse parabolic interpolation and golden section.

If the error function has local minima, the Powell method can miss the best solution. In order to move from a local minimum, the method can be stochastically restarted and another minimum is searched; the global minimum is estimated as the smallest local minimum after a given number of iterations.

We have implemented the following heuristic algorithm. The Powell optimization is embedded within a simple random optimization loop, which is periodically restarted (cf. Fig. 11). The external radius of each shell is randomly sampled, and the best configuration according to the error metric is retained and serves for a new Powell optimization. The best Powell's optimized geometry is retained.

In Figure 11, TEST1 means that the total number of iterations has reached the maximum (200 for instance). The algorithm ends when TEST1 is true. TEST2 means that the number of random iterations is equal to a maximum (10 for instance) or that the random configuration is better than the previous ones. The Powell's optimization is triggered only when TEST2 is true. The loop in Figure 11 is used to find a new random restart for the Powell's iterative algorithm, which is not detailed in the figure.

The method is based on the hypothesis that, when a solution exists and is unique, it will be found after a finite number of iterations. In applications, the algorithm has been found to be very robust; it always produces a solution which is close to the best one provided that the maximum number of iterations is high enough. This number is a parameter. In practical applications, it has been verified that two hundred total Powell iterations is a reasonably good value, but in theory, one has no certainty to handle the best minimum during the two hundred iteration search.

Fig. 11 Heuristic optimization algorithm: Powell optimization embedded within a simple random optimization loop.

3.3 Implementation with SGaRD code

The code has been adapted for material interface position determination using the direct gamma lines intensity, extracted from the measured gamma spectrum. For this purpose, the ray-tracing solver is embedded in the stochastically restarted Powell minimization iterative algorithm (cf. Fig. 12). The best solution after a prescribed maximum iteration number is the final optimized solution.

A typical computation time for the optimization process is of the order of seconds or less on a personal computer (DELL^® precision T5600).

Fig. 12 Inverse transport for geometry search using SGaRD ray tracing solver.

3.4 Numerical results

For illustration of the methodology, we present the numerical results obtained in the case of a spherical shielded source composed of an inner pellet of uranium surrounded by polyethylene, air and aluminum concentric shells. We suppose that the external radius, the materials isotopic composition and density and the direct line leakage spectrum are known. Here, the “measurements” are simulated with SGaRD uncollided direct transport calculations using the exact geometry. Five characteristic lines are chosen for the optimization (cf. Tab. 3).

The evolution of the error metric with the number of iterations is presented in Figure 13. Each jump corresponds to a stochastic restart of the Powell optimization algorithm. A first optimized configuration is found after 13 Powell iterations. The external radius of the uranium sphere is yet close to the solution, but the other interface positions are still approximate. In this case, an almost exact solution is found after 77 iterations and several stochastic restarts where the algorithm comes back to the same local minimum. Due to the stochastic restarts, if the random generator seed is changed, the best solution will be obtained after a different number of iterations. This solution is not guaranteed to be exactly the same as the previous one.

The results of the optimization process using Powell method are presented in Table 4.

The purpose of this numerical test is to verify the capacity of the algorithm to produce correct solutions. In practical applications, all the geometries with an error metric compatible with the experimental accuracy should be retained.

Fig. 13 Evolution of the error metric with the number of iterations.

4 Conclusion

SGaRD code used for fast calculation of uncollided and total leakage spectrum out of spherical shielded sources has been updated in order to simulate the detector pulse height spectra. For this purpose, the DRF calculated using MCNP5 is folded with the leakage spectrum. The methodology is illustrated on a test problem with a gamma source and a scintillation NaI detector. The results are in good agreement with a full Monte Carlo calculation. Future developments will include the extension of the source term with neutron-induced gamma rays and electron-bremsstrahlung radiation.

The second application is the fast material thicknesses identification performed using a ray-tracing transport solver embedded in a gradient-free stochastically restarted Powell iterative optimization loop.

The code has been verified to work well, using synthetic and measured data. As the method is based on the gamma ray measurements, it is affected by the optical thickness of the shielding materials when the gamma lines are badly or even not detected due to the attenuation.

This technique works when the material composition is known. It is worth noting that when the materials are not known a priori, information on the fissile and shielding materials can be inferred from the presence of specific gamma lines and the relative intensities of these lines.

Considering this application, future works will include the sensitivity analysis to measurement uncertainties, the use of other optimization parameters such as the source isotopic composition and the comparison with other optimization methods.

References

All Tables

All Figures

	Fig. 1 SGaRD input/output flow diagram for gamma spectra calculation.
In the text

	Fig. 2 Gamma source lines (blue) are condensed into groups (red) for multigroup SN calculations.
In the text

	Fig. 3 Multigroup SN leakage spectra for both scattered background (red) and uncollided component (blue).
In the text

	Fig. 4 Discrete uncollided lines (blue) with multigroup scattered background (red).
In the text

	Fig. 7 Scintillation NaI detector response functions. The different colors correspond to different incident energies of the gamma particles entering the detector.
In the text

	Fig. 5 Detector and source geometry.
In the text

	Fig. 6 Source leakage spectra calculated using MCNP5 and SGaRD.
In the text

	Fig. 8 Pulse height spectra calculated using MCNP5 and SGaRD.
In the text

	Fig. 9 1-D spherical ray tracing. The γ source is in the shaded shell.
In the text

	Fig. 10 Convergence to the analytic solution of the ray tracing and SN leakage calculations with the number of directions.
In the text

	Fig. 11 Heuristic optimization algorithm: Powell optimization embedded within a simple random optimization loop.
In the text

	Fig. 12 Inverse transport for geometry search using SGaRD ray tracing solver.
In the text

	Fig. 13 Evolution of the error metric with the number of iterations.
In the text