NARMER-1: a photon point-kernel code with build-up factors

This paper presents an overview of NARMER-1, the new generation of photon point-kernel code developed by the Reactor Studies and Applied Mathematics Unit (SERMA) at CEA Saclay Center. After a short introduction giving some history points and the current context of development of the code, the paper exposes the principles implemented in the calculation, the physical quantities computed and surveys the generic features: programming language, computer platforms, geometry package, sources description, etc. Moreover, specific and recent features are also detailed: exclusion sphere, tetrahedral meshes, parallel operations. Then some points about verification and validation are presented. Finally we present some tools that can help the user for operations like visualization and pre-treatment.


Introduction
Starting from the 60s, CEA launched the development of MERCURE, a gamma-ray shielding code using the pointkernel integration formalism. MERCURE-3 was published in 1967 [1], and MERCURE-4 in 1974 [2]. Later in the 80s, EDF and COGEMA (AREVA) became sponsors of MERCURE for its development and validation. MERCURE-5v3 released in 1995 [3] followed by MERCURE-6 [4]. MERCURE uses a Monte Carlo method for the point-kernel integration and build-up factors to take scattered photons into account. MERCURE-6 has enabled the build-up factor calculation for composed shields by using new algorithms developed by SERMA (Reactor Studies and Applied Mathematics Unit) [5]. Then, the development of NARMER-1 has been launched to integrate both gamma-ray attenuation and reflection features (implemented by MERCURE-6 and NARCISSE-3 [6] respectively). For this purpose, the code has been re-written in C++ allowing a modern programming based on object oriented design, with a great care to some extra-functional properties such as modularity, extensibility, maintainability and testability. Ongoing works are related to the improvement of reflection calculations, parallel computing, interoperability with CAD geometries, verification and validation, and integrability in the EDF, CEA software platforms. The development is also based on a components sharing strategy with other SERMA codes for example for the base libraries (geometry, composition, source term, etc.) NARMER-1 is part of a set of software developed by SERMA in the broad field of reactor studies. It is closely related to TRIPOLI-4 ® [7] as it uses the same geometry package allowing easier comparison to this reference Monte Carlo code for validation purpose. Depending on the period, up to 5 people have been contributing to MERCURE and now to NARMER-1, their activities covering development, integration & architecture, V&V, documentation, user support, distribution and licensing.

Straight-line attenuation
Considering an isotropic source represented by a Dirac in space and energy, the point-kernel of attenuation has the following classical expression: where μ ୧ (E) is the attenuation coefficient while t ୧ is the track length in the homogeneous medium i for the energy E.
The point-kernel represents the direct photon flow, without scattering. Diffusion contributions are taken into account by a build-up factor which depends on the computed physical quantity, the material and the track length in the medium. Indeed, the build-up factor is defined by the quotient of the response taking into account the diffusion contributions by the response without scattering.

Single-layer shield build-up factors
Single-layer shield build-up factors are calculated in a homogenous and infinite medium for a source defined as punctual, isotropic, mono-kinetic, placed in a spherical geometry. They have been calculated by running a discrete ordinates transport code, TWODANT [8], and following a method described by KITSOS, et al. [9]. Energy groups, angular mesh and space mesh have met the following features: 218 energy groups from 1 keV to 10 MeV, P7 for the Legendre expansion of the scattering cross sections, S96 for the Gauss-Legendre angular quadrature and at least 28 spatial meshes per mean free path. The build-up factors take the following phenomena into account: photoelectric effect, coherent and incoherent scattering, pair production, secondary sources of bremsstrahlung and fluorescence. Cross sections come from JEF2.2 evaluation files [10]. The build-up factors are available for each single chemical element from hydrogen (Z=1) to californium (Z=98), and are tabulated in 22 thicknesses (from 0.5 to 50 mean free path) and 195 energy groups (from 15keV to10MeV). Energy groups are identical to energy groups of total microscopic cross sections. For a homogeneous material composed of many chemical elements, the build-up factors are computed on the fly by the code which determines an equivalent atomic number of the mixture (ܼ ). For a given energy group g, the total macroscopic cross-section per electron of the mixture is calculated in equation (2) while equation (3) gives the total cross section per electron.
The equivalent atomic number of the mixture (ܼ , ) is obtained with respects to the corresponding computed cross-section (ߪ ). Thus, the build-up factors of the mixture correspond to the build-up for (ܼ ) obtained by a log-log interpolation of the ‫)ܼ(ܤ‬ function, by energy group, and as a function of the number of mean free path.

Multi-layer shield build-up factors
In this section layers are noted L i numbered increasingly from the source, Z i , X i (in unit of mfp) and B i are the corresponding atomic number, thickness and build-up, respectively.
Considering a double layer shield, NARMER-1 implements the SERMA formula (4) to compute the corresponding build-up factor.
For a multi-layer (more than two) shield, NARMER-1 uses an iterative process that replaces, from source to calculation point, two layers by an equivalent layer. First, the equivalent atomic number ܼ is determined by using a neural network trained on large training set. This equivalent atomic number ܼ is the closest to a solution to equation (8), where Z is the unknown quantity.
As shown by Table 1, the neural network is divided into five regions that are delimited with respect to physical considerations about sharp variations of ܼ induced by fluorescence [5]. Table 1. Description of the five neural networks.
Then, the equivalent thickness (ܺ ) of the equivalent material is calculated to reproduce, with the ܼ material, the double-layer build-up, by solving equations (9) where X is the unknown quantity.
The above requirements are summarized in equations (10).
Hence, considering a multi-layer shield composed of N layers, the iterative process starts with the first steps shown by Fig. 1. Fig. 1. Determination of a three layer shield build-up factor.

Source integration
For a source S distributed over energy and space, the response ܴ ఈ described by a response function ݂ ఈ ‫)ܧ(‬ to be calculated follows the equation (11).
NARMER-1 uses a Monte Carlo method to integrate points-kernel with build-up factor over the source domain. Considering the spatial and energy meshes defining the source, the code performs a stratified sampling as variance reduction technique. New sampling points are added in each cell proportionally to the variance of the cell, until the variance on the global result reaches the asked precision.

Computed physical quantities
The following response functions are available in NARMER- It is currently developed on Linux system but is maintained and available on both Linux and Windows.

Geometry package
NARMER-1 comes with the TRIPOLI-4 ® native geometry package, allowing for both a pure surface-based representation, and a combinatorial representation with predefined shapes and boolean operators (any combination of these two kinds of representations can be adopted).
As well as TRIPOLI-4 ® , NARMER-1 has been also made compatible with any geometry developed in the format of the ROOT software (Brun and Rademakers, 1997). This allows the user to directly use (without recompiling) an already built ROOT geometry. More generally, NARMER-1 may be linked to any thirdparty geometry module providing an Application Programming Interface for the current particle position and the next intersection in the direction of flight. Details concerning visualization tools for building and checking the geometry will be provided in Section 4.

Sources description
Sources are defined with an energetic spectrum and an emitting volume according to a mesh at a given position.

Energy Distribution
Line spectrum can be used as well as group spectrum based on a 195 bounds pre-defined mesh.

Spatial Distribution
Volume based sources can be described using the following type of meshes: Cartesian, Spherical, Cylindrical, Conic, Toroidal, General Conic, Constant Width Conic, and General Cylindrical.
Surface based sources can be described using the following type of meshes: Cartesian, Spherical Shell, Cylindrical Shell, Disk, Conic Shell, and Toroidal Shell.

Exclusion sphere
In order to perform computations for some detectors that can be very close to the emitting volume or inside the emitting volume, it is necessary to avoid a singularity due to the 1/r² factor that can lead to infinite variance. Therefore the code uses an exclusion sphere. The radius of the sphere can be set by the user.

Tetrahedral meshes
For specific uses, geometry and emitting volumes are modelled using computer-aided design software allowing easier definition of very complex 3D scenes. In this context the code has the possibility to define tetrahedral  As an illustration also allowing some considerations about verification and validation, the simple test case Fig.  2 shows a detector point belonging to a median plan defined by z=0 separating two identical emitting volumes having the form of a parallelepiped. The emitting volume in the region z>0 is described by a tetrahedral mesh and the one in the region z<0 is described by a cartesian mesh. NARMER-1 detailed results show that the geometrical intensity is identical for both emitting volumes, and also that the contribution to the global result is close to 50% for each source.

Double Energy Angle Differential Albedos
Lacunar media shielding studies can be dealt with albedo concept which is an alternative to the "exact" transport method calculations (SN, Monte Carlo). For multiple reflections on lacunar medium walls, one needs double energy-angle differential albedos because they are function of the incident energy. Reflection feature is a point of interest for EDF sponsor. R&D development and validation studies in this domain are still under progress at SERMA.

Parallel operations
Complex 3D scenes with a high number of gamma sources and a lot of geometrical volumes to cross can lead to substantial simulation times. Therefore, some developments and validation studies are currently under progress to use multi-threads algorithms.

Verification and Validation
Since the NARMER-1 code is attended to be used for radiation protection purposes, verification and validation are crucial issues. Verification is supported in particular by non-regression testing. Excellent results are observed about non-regression results with respect to MERCURE-6. Validation is based on comparison with calculations performed with the TRIPOLI-4 ® Monte Carlo code as well as with experimental benchmarks.
Here are presented four cases illustrating the validation aspect: two basic tests and two realistic cases in which the ambient dose equivalent rate H*(10) is calculated with NARMER-1 and confronted with TRIPOLI-4 ® results.
The definition of the geometry is common for both types of calculation.

Multi-layers configuration
In this case, a point-like J source is surrounded by six concentric spherical screens (from the centre to the outer blanket: Al -Cu -Al -Zr -Cu -Zr). Calculations were done for the following discrete values of J energy:   Table 2. Above 300 keV of J energy, the relative difference does not exceed 10%. However, below 300 keV, NARMER-1 results overestimate TRIPOLI-4 ® significantly. We note that 300 keV represents the cut energy in the definition of the neural networks domains used to determine the equivalent atomic number of the fictive layer for build-up calculations (see Table 1). Cu is replaced by Fe in the material list of the double-layer build-up formula. The back-scattering flux is overestimated in the third layer (Al) due the heavier material (Zr) in the fourth layer. This test case uses multiple times the iterative process (Fig. 1). A point-like J source is located at the pipe's centre with three different energies: 10 MeV, 1 MeV and 0.1 MeV.
The ambient dose equivalent rate H*(10) is calculated at three locations distant from the source of 5 cm (in the water), 8.05 cm (in the steel) and 9 cm (in the air). These locations are represented by dots in Fig. 5. The whole electromagnetic cascade has been taken into account for the 10 MeV while only J rays have been transported for the 1 MeV and 0.1 MeV sources in the TRIPOLI-4® calculations. Results of the H*(10) relative difference as defined by equation (12) are shown for different energy source as a function of the distance from the source in Fig. 6 whereas H*(10) numerical values of NARMER-1 calculations are shown in Table 3. It can be seen that the relative difference remains below about 10% in all cases which reflects the good agreement between the two codes.

Water pipe with protection
This case corresponds to a realistic water pipe which is made of a hollow steel cylinder (inner radius 30 cm, 2 cm thick) full of water which is covered by a 4 cm thick lead layer.
The source is distributed on the volume of the water and takes several discrete values (7 MeV, 3 MeV, 1 MeV and 0.6 MeV).
The ambient equivalent dose is calculated at different distances from the centre of the pipe: 30 cm (interface water/steel), 32 cm (interface steel/lead), 36 cm (interface lead/air), and in the air at 40 cm, 45 cm, 50 cm, 55 cm, 60 cm, 80 cm and 90 cm represented by dots in Fig. 7. The full electromagnetic cascade has been taken into account for the TRIPOLI-4® calculations for the 7 MeV and 3 MeV sources. Only J rays have been transported for the two others energy source. The difference between NARMER-1 and TRIPOLI-4 ® is between -7% and +18%. This difference does not vary linearly with the source energy because of different physical effects: Compton effect is dominant at low energy whereas at high energy pair creation can occur.

UO 2 -PuO 2 source
A cylinder of PuO 2 -UO 2 , in which the J source is equally distributed, is located in front of a series of screens as drawn on Fig. 9. The configuration of the series of 40 cm long screens is defined as following: x A 1 cm thick, 10 cm tall, iron screen is located at 28.5 cm from the center of the source, x A 5 cm thick, 17 cm tall, Kyowa-glass screen is located at 30 cm from the center of the source x A 0.9 cm thick, 17 cm tall, glass screen is located at 30 cm from the center of the source close to a 0.8 cm thick, 17 cm tall, leaded glass screen situated at 30.9 cm from the center of the source. The source energy varies from 15 keV to 3.67 MeV and its distribution is illustrated in Fig. 10. The atomic compositions are noted in the Table 5. The ambient dose equivalent rate H*(10) is calculated at three points, respectively at coordinates (50.,0.,-10.), (50.,0.,0.) and (50.,0.,10.), symbolized by dots in Fig. 9. The Table 6 presents the results of the calculations made with NARMER-1 and TRIPOLI-4® and their comparison.

Visualization and pre-treatment
Sharing the same geometry package as TRIPOLI-4 ® , NARMER-1 can be used with the pre-and-post treatment tool distributed with TRIPOLI-4 ® , known as "Visualizer T4G". The user also has possibility to check the definition of the geometry by using the last verification tools implemented in TRIPOLI-4 ® , including "Check Geom" and "Check Tracking" functionalities.

Summary and conclusion
NARMER-1 is a simplified code that allows to perform fast calculations to evaluate ambient dose equivalent H*(10) and KERMA for radiation protection studies.
A validation study has been done with comparison to TRIPOLI-4 ® considered as a reference code.
Basic configurations and sophisticated configurations show overall a good agreement between the two codes. In some cases, we observe an overestimation of NARMER-1 with respect to TRIPOLI-4 ® . The results are conservative for radiation protection studies.
NARMER-1 calculations are much faster than TRIPOLI-4 ® ones: NARMER-1 calculations last less than one minute and can reach several minutes for complicated configurations. These calculations times remain orders of magnitude below those met for TRIPOLI-4 ® .