https://doi.org/10.1051/epjn/2017006
Regular Article
Gamma ray transport simulations using SGaRD code
CEA DAM IledeFrance, BruyèresleChâtel,
91297
Arpajon cedex, France
^{⁎} email: philippe.humbert@cea.fr
Received:
18
November
2016
Received in final form:
7
February
2017
Accepted:
15
February
2017
Published online: 21 March 2017
SGaRD (Spectroscopy, Gamma rays, Rapid, Deterministic) code is used for the fast calculation of the gammaray spectrum, produced by a spherical shielded source and measured by a detector. The photon source lines originate from the radioactive decay of the unstable isotopes. The leakage spectrum is separated in two parts: the uncollided component is transported by ray tracing, and the scattered component is calculated using a multigroup discrete ordinates method. The pulse height spectrum is then simulated by folding the leakage spectrum with the detector response function, which is precalculated for each considered detector type. An application to the simulation of the gamma spectrum produced by a natural uranium ball coated with plexiglass and measured using a NaI detector is presented. The SGaRD code is also used to infer the dimensions of a onedimensional model of a shielded gamma ray source. The method is based on the simulation of the uncollided leakage current of discrete gamma lines that are produced by nuclear decay. The material thicknesses are computed with SGaRD using a fast raytracing algorithm embedded in a nonlinear multidimensional iterative optimization procedure that minimizes the error metric between calculated and measured signatures.
© P. Humbert and B. Méchitoua, published by EDP Sciences, 2017
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Realtime applications require fast and accurate calculation of the detected gammaray 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 onedimensional spherical assemblies has been updated in order to take into account the response function of various types of detectors and for the identification of gammashielded sources geometric characteristics.
SGaRD has two different transport solvers. The first one is a multigroup discrete ordinates S_{N} solver [2] for the integrodifferential transport equation in onedimensional spherical geometries using the spherical coordinates (r,μ). The second one is a “Method of Characteristics” or raytracing 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 radiationsource 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 highprecision 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 gradientfree 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 Gammaray spectra
The gammaray 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 onedimensional 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 neutroninduced 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 1D raytracing transport solver, presented in Section 3.1.2. The continuous component is transported using a 1D 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 raytracing 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 S_{N} calculations. 
Fig. 3 Multigroup S_{N} 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.
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 1D 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 twodimensional, SGaRD calculations stay onedimensional 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^{9} particles studied.
Fig. 5 Detector and source geometry. 
Example specification.
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 raytracing 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 raytracing 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 1D 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 raytracing 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^{−1} and Q = 1948.531 cm^{−3}/s.
The discrete ordinates S_{N} and raytracing calculations are both performed using SGaRD. The S_{N} solver is a onedimensional spherical finite difference solver, using an equiweight angular quadrature and a 100cell 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 raytracing method is faster and more accurate than the S_{N} discrete ordinates for unscattered transport.
Convergence of the raytracing and S_{N} solvers with the number of directions.
Fig. 10 Convergence to the analytic solution of the ray tracing and S_{N} leakage calculations with the number of directions. 
3.2 Optimization − inverse transport
The goal of the inverse transport problem is to find the 1D 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_{0} 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 wellchosen 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 raytracing 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.
γ line characteristics and simulated leakage measurements.
Fig. 13 Evolution of the error metric with the number of iterations. 
Results of the optimization, external radius.
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 neutroninduced gamma rays and electronbremsstrahlung radiation.
The second application is the fast material thicknesses identification performed using a raytracing transport solver embedded in a gradientfree 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
 Ph. Humbert, B. Méchitoua, Fast gamma ray leakage spectra simulation, in Proc. M&C 2009, Saratoga Springs, New York, May 3–7, 2009 (American Nuclear Society, 2009) (CDROM) (In the text)
 A. Hébert, Applied reactor physics (Presses Internationales Polytechnique, 2009) (In the text)
 X5 Monte Carlo Team 2003, “MCNP A General Monte Carlo NParticle Transport Code Version 5,” LAUR031987 (2003) (In the text)
 D.J. Mitchell, J. Mattingly, Rapid computation of gammaray spectra for onedimensional source models, Trans. Am. Nucl. Soc. 98, 565 (2008) (In the text)
 J.A. Favorite, Using the Schwinger variational functional for solution of inverse transport problems, Nucl. Sci. Eng. 146, 51 (2004) [CrossRef]
 K.C. Bledsoe, J.A. Favorite, T. Aldemir, Using the LevenbergMarquardt method for solutions of inverse transport problems in one and twodimensional geometries, Nucl. Tech. 176, 106 (2011) [CrossRef]
 J.C. Armstrong, J.A. Favorite, Identification of unknown interface locations in source/shield system using the mesh adaptive direct search method, Trans. Am. Nucl. Soc. 106, 375 (2012)
 J. Mattingly, D.J. Mitchell, A framework for the solution of inverse radiation transport problems, IEEE Trans. Nucl. Sci. 57, 3734 (2010) (In the text)
 W.H. Press et al., Numerical recipes (University Press, Cambridge, 1994), Chapter 10 (In the text)
 R. Brent, Algorithms for minimization without derivatives (Prentice Hall, Englewood Cliffs, N.J., 1973), Chapter 7 (In the text)
 A. Tsilanizara et al., DARWIN: an evolution code system for a large range of applications, in Proc. Of the 9th Int. Conf. on Radiation Shielding, Tsukuba, Japan, Oct. 17–22 (1999), p. 845. (In the text)
 C.M. Salgado et al., Validation of a NaI(Tl) detector's model developed with MCNPX code, Prog. Nucl. Energy 59, 19 (2012) [CrossRef] (In the text)
 K.M. Case, F. de Hoffmann, G. Placzek, Introduction to the theory of neutron diffusion, Los Alamos Scientific Laboratory report (USAEC, 1953), Vol. 1 (In the text)
Cite this article as: Philippe Humbert, Boukhmès Méchitoua, Gamma ray transport simulations using SGaRD code, EPJ Nuclear Sci. Technol. 3, 9 (2017)
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 S_{N} calculations. 

In the text 
Fig. 3 Multigroup S_{N} 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 1D 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 S_{N} 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 