Interaction Mechanisms

The photoelectric absorption effect is an interaction mechanism in which the gamma ray provides all its energy to a bound electron, as depicted in Fig. 1(a). In this process, some of the energy overcomes the binding energy of a bound electron, and the free electron absorbs the remaining photon energy in the form of kinetic energy. Only a small portion of energy remains in the atom to conserve momentum [16,17]. Gamma rays can be scattered when they interact with a bound electron through Compton scattering, which consequently changes the photon direction and energy. Additionally, the bound electrons gain kinetic energy. After the interaction, the energy of the secondary photon due to scattering as a function of the scattering angle $\theta$ is computed using Eq. (1). $E_{\gamma \text{'}}=\frac{E_{\gamma }}{1+\frac{E_{\gamma }}{{\text{mc}}^2}\left(1-\text{cos}\theta \right)}$ (1)

Fig. 1Full-screen View

(a) Photoelectric effect; (b) Compton scattering; (c) pair production and annihilation of positrons and electrons.

The Compton scattering mechanism is illustrated schematically in Fig. 1(b). Pair production involves the total absorption of the primary gamma ray and the generation of an electron–positron pair with kinetic energies, as expressed in Eq. 2. $T_{{\text{e}}^{-}}+T_{{\text{e}}^{+}}=E_{\gamma }-1.022 \text{MeV},$ (2) where $T_{{\text{e}}^{-}}$ and $T_{{\text{e}}^{+}}$ are the kinetic energies of the electron and positron, respectively.

Based on Eq. (2), the threshold energy for pair production is above 1.022 MeV. After a pair is created, the positron annihilates to release two photons of 0.511 MeV energy each and traverses 180° opposite each other. This mechanism is illustrated in Fig. 1(c). Typically, the annihilated photons are followed individually while they are traversing 180° apart and are considered the classical model in this work. However, a strategy that can reduce computational time yet exhibits an extremely small error has been suggested previously; in this strategy, only one photon with an energy of 0.511 MeV is followed, and its weight is compensated by 2. This approach was considered a simple model [18]. The two strategies mentioned above were investigated to ascertain their accuracy.

Photon Path Length and Scoring

The photon is transported through a medium by sampling the path length to a collision site from the exponential distribution, which can be expressed as follows: $p\left(s\right)=\mu e^{-\mu s},0\le s\le \infty ,$ (3) where $e^{-\mu s}$ is the probability of a gamma ray traveling a distance s without interaction. Thus, the probability that a gamma ray will interact in the interval ds is μ ds. Therefore, the probability of a gamma ray interacting between distance s and s + ds by $\mu e^{\mu s}\text{ d}s$ can be expressed as $\mu \underset{0}{\overset{s}{{\displaystyle \int }}}{e}^{-\mu s\dot{\text{ d}s}}=\underset{0}{\overset{\xi }{{\displaystyle \int }}}\text{d}y\Rightarrow 1-e^{-\mu s}=\xi .$ (4)

Hence, the path length can be calculated as shown in Eq. (5). $s=-\frac{1}{\mu }\text{ln}\left(1-\xi \right)=-\frac{1}{\mu }\text{ln}\xi$ (5)

Here, ξ is a random number governed by the canonical distribution [p(ξ) = 1, 0 ≤ ξ ≤ 1]. Herein, $\mu$ is the total gamma photon attenuation coefficient that describes the probability of the photon undergoing an interaction; mathematically, it is the sum of the individual cross-sections of the interactions. For this code, all photon cross-sections were obtained using XCOM [19]. A graphical representation of the various gamma interaction cross-sections for the materials investigated in the study is shown in Fig. 2.

Fig. 2Full-screen View

(Color online) Gamma interaction cross-section variation with energy in (a) aluminum, (b) concrete, (c) lead, and (d) water.

The quantities of importance in the code, namely the dose rate, photon exposure buildup factor, and albedo, were calculated using a surface crossing estimator [20]. The estimator for a particle with energy ΔE_i in volume ΔV_i is expressed as $C\left(i,j\right)=C\left(i,j\right)+w\times l,$ (6) $C\left(i,j\right)=C\left(i,j\right)+w\times \frac{\text{Δ}{x}_i}{\left|\text{cos}\theta \right|},$ (7) where $C\left(i,j\right)$ is the counter, $w$ the weight of the particle, and $l$ the path length inside the volume ΔV_i. The projection along the thickness is expressed as $\frac{\Delta x_i}{\left|\text{cos}\theta \right|}$, as shown in Eq. 7. The absolute sign is considered to include all particles moving in all directions. When the absolute value approaches $\cos \theta =0,$ a value of 0.005 is assigned to prevent the estimator from becoming indeterminate. By specifying $\text{Δ}{x}_i$ as unity, the scalar flux for energy j outside a particular surface can be expressed as $\phi \left(E_j\right)=\phi \left(E_j\right)+w\times \frac{1}{\left|\text{cos}\theta \right|}.$ (8)

Translation and Rotation of Photon

After the tentative path length is calculated using Eq. 5, the photon traverses to the next interaction site via coordinate transformation and rotation. The new position of the gamma ray after an interaction ($x$, $y$, $z$) in the direction $\overrightarrow{u_0}=\left(u_0,v_0,w_0\right)=\left(\text{sin}\theta \text{ cos}\phi ,\text{ sin}\theta \text{ sin}\phi ,\text{ cos}\theta \right)$ can be computed using Eq. (9). $\overrightarrow{x}={\overrightarrow{x}}_0+{\overrightarrow{u}}_{0}s\Rightarrow \{\begin{array}{c}x=x_0+u_{0}s\\ y=y_0+v_{0}s\\ z=z_0+w_{0}s\end{array}$ (9)

Next, we introduce the rotation matrix $R\left(\theta ,\phi \right)$ (Eq. (10)) as follows: $R\left(\theta ,\phi \right)=\left(\begin{array}{lll}\text{cos}\theta \text{ cos}\phi \hfill & -\text{sin}\phi \hfill & \text{sin}\theta \text{ cos}\phi \hfill \\ \text{cos}\theta \text{ sin}\phi \hfill & \text{cos}\phi \hfill & \text{sin}\theta \text{ sin}\phi \hfill \\ -\text{sin}\theta \hfill & 0\hfill & \cos \theta \hfill \end{array}\right)$ (10)

The function of the rotation matrix $R\left(\theta ,\phi \right)$ is to transform the gamma-ray vector components from the local coordinates to the laboratory coordinate system (x, y, z) via the azimuthal angle $\phi$ about the positive $z$-axis, followed by the rotation of the polar angle $\theta$ about y′ [21]. After the transformation, the three components of the new direction cosine, namely $u,v$, and $w$ are calculated using Eq. (11). $\begin{array}{l}u=\text{sin}\theta \text{ cos}\phi =u_0\text{cos}\theta +\text{sin}\theta \left(w_0\text{ cos}\varphi \text{ cos}{\phi }_0-\text{sin}\varphi \text{ sin}{\phi }_0\right)\\ v=\text{sin}\theta \text{ sin}\phi =v_0\text{ cos}\theta +\text{sin}\theta \left(w_0\text{ cos}\varphi \text{ sin}{\phi }_0+\text{sin}\varphi \text{ cos}{\phi }_0\right)\\ w=\text{cos}\theta =w_0\text{ cos}\theta -\text{sin}\theta \text{ sin}{\theta }_0\text{ cos}\varphi \end{array}$ (11)

Using the aforementionedethodologies, the entire gamma-ray trajectory can be described. Each particle is monitored from birth until any cutoff condition terminates its life history or until it escapes from the medium. During particle transport, sampling is performed using known probability density functions for different physical processes. For double- and triple-layer systems, surface-to-surface tracking [22] accounts for boundary crossing and regional changes in an event in which a photon leaves its current region for a new one.

Variance Reduction (Biasing)

Variance is used to increase the efficiency and computational speed. The developed code applies the survival weight, exponential transform, and Russian roulette. The survival weight prevents the direct absorption of the particles by inducing forced scattering. This technique increases the number of photons that can penetrate important areas, thereby increasing the particle score in significant regions of the medium prior to premature absorption. The survival-weight technique is essential for shielding applications involving thick media. The introduced bias is mitigated by adjusting the particle weight after each collision by multiplying the photon survival weight by a factor, as shown in Eq. 12. $w=w_{\text{O}}\frac{{\text{σ}}_{\text{c}}}{\mu },$ (12) where the ratio $\frac{{\text{σ}}_{\text{c}}}{\mu }$ is the probability of a scattering event.

The photon is assigned a weight of unity from the source, i.e., $w_{\text{O}}=1$., which is adjusted as the simulation proceeds [14]. The exponential transform performs a function similar to the survival weight. However, this technique improves the particle score by increasing the path length of the photon when it approaches a region of interest (ROI), i.e., where scoring is desired. Conversely, the particle path length is reduced when the particle departs from the ROI. The exponential transformation is executed by replacing the actual total cross-section $\mu$ with a transformed total cross-section ${\mu }^{\prime }$ expressed as ${\mu }^{\prime }=\mu -cu$, which is dependent on the direction, or ${\mu }^{\prime }=\mu *c,$ which is the direction-independent transformed total cross-section [23]. In this code, the latter is employed because the directional form introduces fluctuations in the computation of the buildup factor. Subsequently, the path length is sampled using the transformed total cross-section as follows: $s=\frac{-\ln \left(\xi \right)}{{\mu }^{\prime }}\Rightarrow \frac{-\ln \left(\xi \right)}{\mu *c}.$ (13)

The weights of the particles are adjusted accordingly to remove the effect of bias as follows: $w=\frac{w_O}{c}{e}^{\left[-s\mu \left(1-c\right)\right]}.$ (14)

When c > 1, the path length is reduced and is typically used in situations where buildup zones exist at the medium entrance. If c = 1, then sampling is performed using the untransformed total cross-section; meanwhile, if 0 < c < 1, then the path length is stretched, which renders its applicable to thick media. This technique is useful for shielding problems in which the particle history is terminated before detection, which is characteristic of deep penetration problems.

The Russian roulette determines whether a particle’s life history must be terminated or continued when its weight is low. This technique is typically used to prolong the computational time. A particle is assigned a probability of survival, e.g., 1 out of 10. If a generated random number exceeds 0.1, then the history is terminated and the particle history continues [14].

ANN

As mentioned previously, obtaining an arbitrary parameter for the exponential transformation involves trial and error. Once data are available for parameters that require trial and error, such data can be trained using machine-learning techniques, such as ANNs, to predict such parameters. An ANN is based on the manner in which the human brain processes information using neurons. A typical ANN contains neurons, which are known as nodes, that are connected to form a web-like system. Additionally, an ANN comprises an input layer that receives input data and transmits it through a hidden layer, if such a layer exists. Finally, an output is produced after data processing is performed [24]. An ANN can be considered a perceptron if it does not contain any hidden layers, a multilayer perceptron if it comprises three layers, or a deep neural network if it comprises more than three hidden layers [25]. Deep learning is desirable owing to its high precision [26]. In this study, an ANN model, i.e., a deep neural network, is proposed to predict the exponential transformation parameter c. The proposed model is shown in Fig. 3.

Fig. 3Full-screen View

The proposed deep neural network architecture.

The input data are transmitted between the hidden layers using a transfer function, which is expressed as $H_j={\displaystyle \sum _{i=1}^n{w}_i^{\text{T}}{x}_i+b}$

Here, $H_j$ is the input to the next layer and $w_i^{\text{T}}$ is the weight of the ith $x_i$ input. The sum of the weights generates activation, which is the signal strength. After processing the output using the hidden layers, the output is generated using an activation function that can be a linear or nonlinear function, such as a sigmoidal hyperbolic function, expressed as $f\left(x\right)=1/\left[1+\text{exp}\left(-x_i\right)\right].$ (15)

Subsequently, the output returned in the output layer is $y_j=1/\left[1+\text{exp}\left(-H_j\right)\right]$ (16)

Because ANNs are machine-learning techniques, the network must be trained before deployment. In this context, the suitable learning schemes include the Newton, quasi-Newton, Conjugate gradient, gradient descent, and Levenberg–Marquardt schemes. The detailed algorithms of these schemes are not presented here as they are available in [27].

Dose Rate, Exposure Buildup factor, and Albedo Evaluation

After all histories have been simulated, the code computes quantities that characterize the shield's performance by evaluating all tallies corresponding to each. The dose rate provides a quantitative measure of the radiation risk at a point behind the shield and provides information regarding the hazards posed by radiation. The code evaluates the dose rate as follows: $\text{Does}={\displaystyle \sum _{g=1}^G\phi }\left(r,E\right)\cdot K\cdot {\overline{E}}_g\cdot \frac{{\mu }_{\text{en}}^{\text{air}}\left({\overline{E}}_g\right)}{\rho },$ (17) where $\phi \left(r,E\right)$ is the total photon flux at a distance $r$, K the dose conversion factor, ${\overline{E}}_g$ the average photon energy, and $\frac{{\mu }_{\text{en}}^{\text{air}}\left({\overline{E}}_g\right)}{\rho }$ the energy absorption coefficient of air at energy ${\overline{E}}_g$. The dose is the total dose at a point comprising both scattered and noncollided photons.

The exposure buildup factor is evaluated as the ratio of the total dose to the uncollided dose and is expressed as $EBF=\frac{{\displaystyle {\sum }_{g=1}^G{\phi }_g\cdot {\overline{E}}_g}\cdot \frac{{\mu }_{\text{en}}^{\text{μir}}\left({\overline{E}}_g\right)}{\rho }}{{\phi }_u\cdot E_0\frac{{\mu }_{\text{em}}^{\text{uir}}\left({\overline{E}}_g\right)}{\rho }}.$ (18)

Here, ${\phi }_u$ is the uncollided photon flux defined by Beer and Lambert, which is expressed as ${\phi }_u={\phi }_0\text{exp}\left(-\mu r\right)=\left(S/4\pi r^2\right)\text{exp}\left(-\mu r\right)$; and ${\phi }_0$ is the total photon history based on the initial photon flux impinging on the material surface.

The integration of the differential albedo generally computes the total albedo over energy and the solid angle for the directional domain as a function of both the total number and energy albedos as follows: $a_N\left(E_0,{\theta }_0\right)=\underset{0}{\overset{E_0}{{\displaystyle \int }}}\text{d}E\underset{0}{\overset{2\pi }{{\displaystyle \int }}}\text{d}\phi \underset{0}{\overset{\pi /2}{{\displaystyle \int }}}a\left(E_0,{\theta }_0;E,\theta ,\phi \right)\text{sin}\theta \text{d}\theta ,$ (19) $a_E\left(E_0,{\theta }_0\right)=\frac{1}{E_0}\underset{0}{\overset{E_0}{{\displaystyle \int }}}E\text{d}E\underset{0}{\overset{2\pi }{{\displaystyle \int }}}\text{d}\phi \underset{0}{\overset{\pi /2}{{\displaystyle \int }}}a\left(E_0,{\theta }_0;E,\theta ,\phi \right)\text{sin}\theta \text{d}\theta .$ (20)

Here, $a_N\left(E_0,{\theta }_0\right)$ and $a_E\left(E_0,{\theta }_0\right)$ are the number and energy of albedos; and $a\left(E_0,{\theta }_0;E,\theta ,\phi \right)$ is the differential albedo, which is integrated over the energy domain dE as well as over the polar and azimuthal domains ($d\theta$ and $d\phi ,$ respectively). For the Monte Carlo simulations, the total number and energy albedo are computed by replacing the integrals in Eqs. 21 and 22 with a finite sum and tallying for the particles in all directions and energies, which can be expressed as follows: ${\alpha }_N\left(E_0,{\theta }_0\right)=\frac{{\displaystyle {\sum }_{i=1}^{N_R}{w}_i}}{{\phi }_0},$ (21) ${\alpha }_E\left(E_0,{\theta }_0\right)=\frac{{\displaystyle {\sum }_{i=1}^{N_R}\frac{E_i}{E_0}}.w_i}{{\phi }_0},$ (22) where $w_i$ and $E_i$ are the weight and energy of the ith particle, respectively; and $E_0$ is the incident energy. A schematic representation of the reflection process in the material is shown in Fig. 4.

Fig. 4Full-screen View

Illustration showing gamma rays reflected from a surface. In most shielding applications, the interface between air and a shielding material is of interest, with the shield material being the reflective surface (adapted from [14]).

Comparison of EJUSTCO simulation with literature

The aim of this study is to develop a code to investigate gamma-ray transport through different shielding materials and calculate the quantities used to evaluate the capability of gamma-ray shields. In the initial study, monoenergetic photons with an energy of 0.662 MeV were simulated using 20001 photons/cm²/s as the number of histories in a lead medium with a thickness of 1.7 cm (two mean free paths). The photon distribution, deposited energy, and final positions of the photons after each cycle are shown in Fig. 6. This plot allows the particle behavior within a material to be visualized and understood easily.

Fig. 6Full-screen View

(Color online) Distribution of final position of each photon after each history in lead media using EJUSTCO.

The fates of the photons after the simulations as fractions of the normalized photon history are presented in Table 1. In addition, the dose rate and buildup factor at material thicknesses of 1, 1.5, and 2 mfp were computed (see Table 2). The results shown in both tables were compared with those in the literature, and the percentage errors were computed.

Fig. 7 shows the energy spectrum of gamma rays transmitted onto the surface of a 2-mfp-thick water medium, which agreed well with that obtained from the wood MONTERAY Mark (I) code.

Fig. 7Full-screen View

(Color online) Energy spectrum of gamma-ray (6 MeV) transmitted to 2-mfp-thick water.

In another simulation using the same number of particle histories, the effect of pair production inclusion on the fate of photons, buildup factor, and dose rate was investigated using lead with thicknesses of 1 and 3 mfp and a monoenergetic photon source of 10 MeV. The results without pair production are presented in Tables 3 and 4, whereas those with pair production are presented in Tables 5 and 6. In both cases, the results were compared with those reported in the literature. The results of this study show that at energies above 1.022 MeV, the effect of pair production is not negligible as it increases the buildup factor and hence the number of absorption processes in the media. The EJUSTCO simulation results agreed well to the results from literature for both cases. The energy spectrum for the case of 3-mfp-thick lead is presented in Fig. 8 for comparison.

Fig. 8Full-screen View

(Color online) Energy spectrum of gamma-ray (10 MeV) transport in 3-mfp-thick lead.

According to the literature, annihilated photons can be treated via two approaches, as explained in Section 2.1. Thus, to understand the effect of each model on the simulation results and computational time, both models were tested. Their effects are listed in Tables 7 and 8. As shown in Table 7, the simple model can be safely employed because its figure-of-merit (FOM) is higher than that of the classical model, thus confirming the assertions presented in the literature. Meanwhile, Table 8 shows that employing a simple model significantly reduces the computational time by more than one-half of that required by the classical model. Therefore, all other analyses were performed using a simple model.

Comparison of exposure buildup factor results in single materials with standard codes and data

To fully ascertain the usefulness and robustness of the code, it must be validated by comparing the results yielded by it with standard codes and data for different materials with varying thicknesses and photon energies. Thus, exposure buildup factors computed using the EJUSTCO code were validated for water (Table 9), lead (Table 10), iron (Table 11), aluminum (Table 12), and concrete (Table 13) based on results obtained from the EGS4 [29] and MCNP [30] standard codes and the ANS-6.4.3 standard [31] for material thicknesses up to 8 mfp and energies of 1, 2, 3, and 4 MeV. In all cases, a particle history of 1 × 10⁶ was used in the validation, which resulted in an average maximum standard error of approximately 5%.

Exponential transform parameter prediction using ANN

As stated in Section 2, determining an arbitrary parameter for use in the simulation when the exponential transformation technique is active can be time consuming and may require several trials. Therefore, some knowledge regarding the expected value is required. Hence, an ANN model was trained to predict these parameters. The model was trained using the Levenberg–Marquardt algorithm, which achieved an MSE of 0.00076752, as shown in Fig. 9, 10 shows the regression plots of the training and validation data. In the regression, the accuracy of data fitting for training was indicated by an R-value of 0.99527.

Fig. 9Full-screen View

Mean square error analysis of ANN training.

Fig. 10Full-screen View

Regression plots showing fitting performance of the ANN model.

The error and fitting accuracy confirmed that the model can be used to predict the parameter c. The values predicted using the ANN model were compared with those obtained using the EJUSTCO code (see Table 14). The results confirmed the viability of using ANNs to predict the exponential transform parameter. However, because of the few data points employed, predictions beyond the range of the dataset would be inaccurate, which is a significant disadvantage in employing ANNs in regression analysis. Hence, increasing the dataset would extend its applicability beyond current data ranges. To demonstrate the representativeness of the training dataset, Fig. 11 shows the distribution of the data in terms of the input parameters. Fig. 11(a) shows that the energy range extends from low to high energies in a range between 0.5 and 10 MeV. Similarly, Fig. 11(b) shows the distribution of the densities used in the training. Based on the distribution, the density spans from less dense to highly dense materials, thus capturing materials of diverse densities. Further analyses were conducted to understand the effect of incorporating the ANN model on the computational time required by EJUSTCO. The computational times required when the ANN model was used and not used were 88.776 and 89.673 s, respectively, based on a photon history of 20000. Clearly, incorporating the ANN model did not significantly affect the EJUSTCO code.

Fig. 11Full-screen View

Distribution of training sample based on input parameters.

In Table 15, the total times required to predict parameter C using trial and error and the ANN model are listed. Using trial and error, an average of four trials was required, which corresponded to a total time of 1287.44 s. However, no trials were necessitated for the ANN and hence a short time of 273.47 s was required. The results clearly indicate that the ANN is advantageous because it does not require simulation trials to determine the optimum value. Additionally, the effect of incorporating the ANN inclusion on the efficiency of the simulation must be ascertained. Thus, the FOMs of different simulation trials involving aluminum, lead, and carbon at 3 MeV were computed (see Table 16). For the simulation involving aluminum and lead, the FOMs were extremely high in the simulation employing the ANN compared with that without the ANN. However, the FOM for the carbon simulation with the ANN was lower than that without the ANN. This shows that incorporating the ANN can significantly affect the simulation efficiency, which is primarily governed by the prediction accuracy, as listed in Table 14.

Dose rate analysis in single- and triple-layer material systems

This section presents the dose rate calculated using the EJUSTCO code for single- and triple-layer systems. The dose rates in aluminum, lead, water, concrete, and iron are shown in Fig. 12 for the 1 and 3 MeV energy sources. As shown in the plot, the dose reduction reflects an exponential decay, which characterizes gamma-ray attenuation. Furthermore, the attenuation reduction for both energies shows that as the energy increases, the ability of the materials to shield gamma rays decreases owing to the decrease in the attenuation coefficient of the material.

Fig. 12Full-screen View

Dose rate for iron, lead, and water at 1 MeV (a) and 3 MeV (b) with 8 mfp thickness, as predicted by EJUSTCO.

Similarly, the dose rates within layered materials are shown in Fig. 13 and Fig. 14, which have been validated via comparison with results in the literature [14,32]. Based on the figures, the behavior of photon attenuation can be visualized, i.e., the attenuation process involves a combination of individual materials within the system. To demonstrate the attenuation performance of different single materials and their combinations in triple-layer systems, the efficiencies of these systems were computed, as shown in Fig. 15.

Fig. 13Full-screen View

EJUSTCO prediction of relative dose rate variation for different triple-layer shield systems: (a) water/lead/iron and (b) lead/water/iron for 1 MeV gamma photon.

Fig. 14Full-screen View

EJUSTCO prediction of relative dose rate variation for different triple-layer shield systems: (a) water/lead/iron and (b) lead/water/iron for 10 MeV gamma photon.

Fig. 15Full-screen View

Efficiencies of single and triple-layer materials with 5 mfp thickness in attenuating gamma radiation, as predicted by EJUSTCO.

Comparison of exposure buildup factor in double layer system with standard codes and data

The transmission buildup factors for double-layer systems were investigated. Specifically, the buildup factor values for an iron–water system (medium- and high-atomic-number materials) were investigated for two cases. In the first case (as presented in Table 17), the iron thickness fixed while the thickness of water was varied at energies of 1 and 3 MeV. Meanwhile, the converse was applied for the second case, as shown in Table 18. The results from both cases were compared with the results yielded by the EGS4 and MCBLD codes in the literature. Finally, the transmission buildup factors for double layers of lead, aluminum, and iron were validated at energies of 0.662 and 1.25 MeV via comparison with experimental results. In all cases, the results showed good agreement, as presented in Table. 19. The exposure buildup factor of the triple-layer materials was computed and compared with the results yielded by the EGS and PENELOPE codes and an empirical formula derived by Lin and Jiang [29] (see Table 20).

Total number and energy albedo computation using EJUSTCO for perpendicular incidence

The concept of albedo is provided in the Inoduction section and in Section 2. This section presents the numerical results obtained via EJUSTCOimulations. Figures 16 and 17 show the total number albedo values for aluminum and water at varying energies. The results were validated using the Monte Carlo codes MCNP, FOTELP, and PENELOPE. As shown in the figures, all the codes agreed well with one another, thus indicating the ability of the developed code to accurately compute albedo values. The general trend shows that the albedo increases with energy and is lower for heavier materials because of the increased absorption rate in such materials. These observations are consistent with those reported previously.

Fig. 16Full-screen View

Comparison of total number albedo as a function of energy computed using EJUSTCO with general purpose codes MCNP, FOTELP, and PENELOPE for aluminum.

Fig. 17Full-screen View

Comparison of total number albedo as a function of energy computed using EJUSTCO with general purpose codes MCNP, FOTELP, and PENELOPE for water.

Similarly, Fig. 18 shows the energy albedo for both water and aluminum, where the same trend is shown. However, a slight error exists between our code and that of MCNP, which is attributable to the nonconsideration of coherent scattering and secondary photon emissions [33] in our code. Nonetheless, the results show good agreement.

Fig. 18Full-screen View

Comparison of total energy albedo values computed using EJUSTCO and general-purpose code MCNP.