Open-Source Optimization of Hybrid Monte Carlo Methods for Fast Response Modeling of NaI (Tl) and HPGe Gamma Detectors

Abstract

Modeling the response of gamma detectors has long been a challenge within the nuclear community. Significant research has been conducted to digitally replicate instruments that can cost over USD 100,000 and are difficult to operate outside of a laboratory setting. The large cost and availability prevent some from making use of such equipment. Subsequently, there have been multiple attempts to create cost-effective codes that replicate the response of sodium-iodide and high-purity germanium detectors for data derivation related to gamma-ray interaction with matter. While robust programs do exist, they are often subject to export controls and/or they are not intuitive to use. Through the use of hybrid Monte Carlo methods, MATLAB can be used to produce a fast first-order response of various gamma-ray detectors. The combination of a graphical user interface with a numerical-based script allows for open-source and intuitive code. When benchmarked with experimental data from Co-60, Cs-137, and Na-22, the code can numerically calculate a response comparable to experimental and industry-standard response codes. Evidence supports both savings in computational requirements and the inclusion of an intuitive user experience that does not heavily compromise data when compared to other standard codes, such as MCNP and GADRAS, or experimental results. When the application is installed on a Dell Intel i7 computer with 16 cores, the average time to simulate the benchmarked isotopes is 0.26 s. Installation on an HP Intel i7 four-core machine runs the same isotopes in 1.63 s. The results indicate that simple gamma detectors can be modeled in an open-source format. The anticipation for the MATLAB application is to be a tool that can be easily accessible and provide datasets for use in an academic setting requiring gamma-ray detectors. Ultimately, this article provides evidence that hybrid Monte Carlo codes in an open-source format can benefit the nuclear community in both computational time and up-front cost for access.

1. Introduction

Since the 1970s, non-proliferation efforts have been the cornerstone of the nuclear community. A significant portion of these efforts has been dedicated to the detection of gamma-ray radiation and the isotope identification of radioactive materials. Since then, significant progress has been made on detector technology. Advanced materials and algorithms have improved the energy resolution of most detectors to the point where the identification of isotopes can be made with a high level of accuracy. However, for most of the nuclear community, detection remains a costly endeavor. As a result, significant research has been conducted in an attempt to successfully replicate gamma detector response for use in academia and the production of large datasets.

The United States Department of Energy and Oakridge National Laboratory maintains a repository of transport codes that have long been regarded as the industry standard for modeling the transport of gammas and neutrons. Among the long list of specialized codes, two stand out in particular for modeling detector response. The well-known Monte Carlo Neutron Particle Transport (MCNP 6.2) can be used to track gammas and neutrons throughout a user-defined universe of various materials and geometries. The second, Gamma Detector Response and Analysis Software (GADRAS 19.4.0), aims to replicate both gamma and neutron detectors based on user-defined scintillating material and instrument parameters. While both codes provide high-fidelity results when compared to experimental data, three significant issues arise when attempting to use these programs to build models.

A combination of export controls and the undermanning of the Oakridge website creates a lead time of several months for requesting software. Additionally, the export controls prevent some users from applying at all due to their country of nationality. For many settings, waiting months to receive a copy of the codes creates barriers to fully using the software. Despite the usefulness and quality of the software, the lack of timely distribution complicates widespread use.

Additionally, user-friendliness is a characteristic that is lacking with both the MCNP and GADRAS. The MCNP is a terminal-based program that takes a considerable amount of time to learn. For a user who is not well versed with the Windows command terminal, the MCNP can be difficult to utilize to its full potential. Moreover, the input files require a highly specialized language that takes additional time to learn. Although, there is online documentation and user manuals, modeling a detector’s response may offer additional challenges to users over using an actual detector. Conversely, GADRAS makes use of a graphical user interface but lacks the robust documentation that exists for the MCNP. To appropriately model a detector requires trial and error, which may be a deterrent for some users.

The final issue that is present with the codes, particularly with the MCNP, is the computational requirements for tracking particles through the traditional Monte Carlo methods. As a result, simple simulations can take hours or days to produce useful results. In a worst-case scenario, the computing requirement can be beyond that of which common computers can handle. The solution requires super-computers that are not practical for most applications.

While these three issues are not all-inclusive to every user, they present a challenge that reduces the availability of modeling gamma detector response. As a result, they are better off making use of the costly detectors. Recently, research has been conducted with hybrid Monte Carlo methods. There is strong evidence supporting reduced computational requirements to achieve similar results when compared to traditional transport codes. Although the hybrid methods most likely cannot replace traditional Monte Carlo methods, making use of these methods in detector modeling may eliminate computational limitations. Moreover, an open-source code that creates a clean user interface with intuitive inputs eliminates the export control problem and the need to learn complex coding for simple applications.

Traditional codes will remain the industry standard for the foreseeable future. However, an open-source code that implements a hybrid Monte Carlo approach along with an intuitive interface will minimize the problems that commonly arise when modeling gamma detectors. The following is a description of a MATLAB 2023b application integrating these solutions into a first-order approximation of thallium-dopped sodium-iodide (NaI (Tl)) and high-purity germanium (HPGe) detectors.

2. Background

Previous work has been conducted in this area with an attempt to model Compton scattering with the Python language; while useful, this model fails to address all of the interactions that occur within detectors [1]. An additional project addresses the use of parametric equations to model NaI response. However, this model makes use of the MCNP for the production of the models [2].

A simulation that removes the dependence on the MCNP and represents a majority of radiation–matter interactions is ideal. While work has been conducted on neutron detectors, the team focused on the highly specific parameters of BF₃ detectors rather than other detector types [3]. A detector suite that offers simulations for three systems and two common radiation–particle interactions will offer a more comprehensive contribution to radiation modeling.

Radiation effects with matter can be treated deterministically in some specific applications. A common example is the dose rate and exposure. Deterministic methods are often conservative because they return the same results for a given set of inputs. Due to the use in safety protocols, deterministic results are sought after. However, radiation effects are highly probabilistic and most applications require a different approach. Traditionally, this is conducted via direct Monte Carlo modeling, which applies randomness to a system by statistical sampling over time [4]. Direct Monte Carlo methods require significant computational power and are a generalized approach to randomness, but the hybrid Monte Carlo method was created to apply random sampling with reduced computational requirements for complex modeling problems. While the decrease in computational requirement varies from the problem to the computer used, estimates are that the hybrid Monte Carlo methods are around 10 times faster [5].

An attempt to utilize hybrid Monte Carlo methods for calculating the counting efficiency of a NaI detector was conducted in 2007. The team integrated analytical relationships with traditional Monte Carlo methods to simulate point-source and planer radiation fields in a 3″ × 3″ NaI detector. Their findings support that this approach produces efficiencies similar to the results produced by earlier accepted models for three discrete energies of Cs-137, Co-60, and Na-22. The overall result is that the time to compute was 0.2–0.8 s on a computer running Windows XP x86. The results of their findings highlight the computational advantage of integrating analytical relationships with traditional Monte Carlo methods [6].

2.1. Overview of the Hybrid Monte Carlo Method

The hybrid Monte Carlo method is a means of methodically sampling random distributions. For instance, Compton scattering can be efficiently modeled because the input energies are known and the approximate gamma interactions adhere to mathematical relations based on the density of the material. While scattering angles are still computationally random, the distribution of each scatter is statistically determined based on incoming energy. The pre-allocated distribution of scattering angles allows for a reduced computational requirement because it is calculated independently of each time-dependent event [7]. The process follows a 5-step process, as illustrated in Figure 1 [8].

The distribution is generated in the initialization block. The empirical correlations of the gamma interactions are applied in the relationship block. The algorithm then conditionally accepts or rejects the outcome of the inputs. Finally, randomness is applied to the inputs based on the initialized distribution. Iteration is applied to represent each event in a sample.

2.2. Gamma Detection Simulation Procedure

The NaI and the HPGe detectors feature systems capable of sensing gamma and high-energy X-ray radiation. While several different detector materials are on the market, we focused on two of the most readily available. The simulation procedure is shown in the following block diagram in Figure 2. It consists of 5 “If statements” that compare the input gammas’ energies to MATLAB-generated random numbers. A version of the pseudo-code is provided in Appendix A.

The gamma photon detector interactions are governed by the following three gamma interactions with matter: the photoelectric effect, Compton scattering, and pair production. These physical interactions determine the response of the detector and the output that is displayed in the graphical user interface. The fraction of energy distribution is based on the probability that all photons interacting with the matter inside the detector undergo one of the aforementioned phenomena. Equation (1) is the total probability of interaction for any single source [9]. The derivation for Equation (2) through (4) was produced by digitizing a plot of the energy dependence of the various gamma-ray interaction processes in sodium-iodide, then using MATLAB curve fitting to produce correlations [10].

$${\mu }_T={\mu }_{pe}+{\mu }_{cs}+{\mu }_{pp}$$

(1)

$${\mu }_{pe }=(0.000243+0.72010)[{\displaystyle \frac{{\mathrm{cm}}^2}{\rho }}]\ast e^{\left(-7.36405\ast \mathrm{E}\left[\mathrm{M}\mathrm{e}\mathrm{V}\right]\right)}$$

(2)

$${\mu }_{cs }=(0.0319+0.09681)[{\displaystyle \frac{{\mathrm{cm}}^2}{\rho }}]\ast e^{\left(-1.2035\ast \mathrm{E}\left[\mathrm{M}\mathrm{e}\mathrm{V}\right]\right)}$$

(3)

$${\mu }_{pp }=(0.000595)[{\displaystyle \frac{{\mathrm{cm}}^2}{\rho }}]\ast e^{(0.597\ast \mathrm{E}[\mathrm{M}\mathrm{e}\mathrm{V}\left]\right)}$$

(4)

where µ_T is the total probability of gamma interaction with matter, µ_pe is the probability of the photoelectric effect, µ_CS is for Compton scattering, and µ_pp is for pair production. E is the gamma energy in MeV. The values were estimated by digitization and curve fitting using energy-dependent relations for sodium-iodide [10,11]. While a good first approximation for the model, the values for the photoelectric effect and Compton scattering required modification to better match the benchmarked data. A comparison of the NaI-averaged values shown in Equations (2)–(4) and their modification from the above sources are shown below in Figure 3.

Compton scattering remains the most difficult to model because it is the most dominant effect for the energies used in the model and the additional requirements to calculate the scattering angle. This is partially due to the random nature of Compton scattering. While it is relatively easy to determine the energy following a scatter, it is much more difficult to predict the angle for scattering. Not only are the scattering angles dependent upon the energy of the incoming gamma, but there is a bias for forward scattering with higher energies and backscattering at lower energies. The energy of the gamma following scattering can be found from Equation (5) and the bias is based on the Klein–Nishina formula in Equation (6) [12], as follows:

$${\lambda }^{\prime }-\lambda ={\displaystyle \frac{h}{m_{0}c}}(1-cos(\theta ))$$

(5)

$${\displaystyle \frac{\delta \sigma }{\delta \Omega }}={\displaystyle \frac{1}{2}}{r_e}^2({\displaystyle \frac{\lambda }{{\lambda }^{\prime }}})({\displaystyle \frac{\lambda }{{\lambda }^{\prime }}})+({\displaystyle \frac{{\lambda }^{\prime }}{\lambda }})-si{n}^2(\theta )$$

(6)

where λ is the initial gamma wavelength, λ′ is the scattered gamma wavelength, h is the Plank’s constant, m₀c² is the rest mass of an electron, ϴ is the scattering angle, r_e is the classical electron radius, and ${\displaystyle \frac{\delta \sigma }{\delta \Omega }}$ is the differential cross-section of scattered gammas. The sample angles for scattering are generated during the initialization of the hybrid Monte Carlo method and based on the Klein–Nishina formula for an initial energy of 0.5 MeV. Figure 4 shows this distribution and the probability of occurrence.

Pair production is a binary response in the MATLAB script where the energy of the current gamma is compared to 1.02 MeV. If the current energy is greater than 1.02 MeV, a positron–electron production and subsequent positron annihilation occur if a variable set to a random number between 0 and 1 is greater than another random number of the same characteristic. The relationship is shown in Equations (7) and (8). Both values are produced through the Mersenne Twister method as MATLAB’s standard random number generator [13]. A single-escape or double-escape event is then simulated by another variable defined by MATLAB’s standard random number generator.

$$\gamma \to e^{+}+e^{-}$$

(7)

$$e^{+}+e^{-}\to 2\gamma$$

(8)

The photoelectric effect is simulated by creating a photo peak based on a user-inputted resolution at the initial gamma energies that each radioactive source emits, e.g., Co-60 at 1.17 and 1.33 MeV. The resolution is calculated by assuming that the resolution for a NaI at Cs-137 is 8% and HPGe is 0.8%. Using the NaI, the resolution was established by using the FWHM of the three isotopes. Then, a linear fit was applied to produce Equation (9). The resolution requires a user input for a detector of a higher resolution than NaI. The correction factor ranges from 0.1 to 1 and corresponds to a user input of 100% and 0%, respectively.

$$Resolution=Correction Factor\ast (-0.0298\ast Energy+0.08972)$$

(9)

While the resolution can vary for each source, using Cs-137, where the peak gamma is near the center of the ranges displayed for these models, the assumption is that the resolution is subjectively sufficient for all gammas between 0.5 MeV and 2 MeV. It is important to note that the user-specified input when at 0% is 0.01% of the correction factor. This prevents a null resolution value and is a result of MATLAB requiring the input to be between 0 and 100. When a user selects 0%, the model simulates a perfect 3″ × 3″ NaI detector and 100% results in a perfect 3″ × 3″ HPGe detector. Figure 5 shows the linear estimated resolution for a 0% user input compared to the benchmarked data used to derive Equation (9). Note that the line denoted detector resolution is an estimate of the NaI detector corrected for the inverse square root relationship.

An analysis of the detector energy resolution and absolute efficiency was conducted to ensure that the estimate presented in Figure 5 is a reasonable assumption. The analysis was conducted with 10-microcurie sources verified on 1 July 2003. By measuring the counted activity and comparing the value with the source-corrected activity for the date of measurement, the absolute efficiency can be determined. Equation (10) presents the relationship for the absolute efficiency and Equation (11) presents the relationship for the energy resolution as a function of the gamma energy, as follows:

$$\epsilon ={\displaystyle \frac{A_0}{A_f\ast BR}}$$

(10)

$$Energy Resolution={\displaystyle \frac{FWHM}{Gamma Energy \left(keV\right)}}$$

(11)

where A₀ is the initial activity, A_f is the final activity, BR is the branching ratio, and FWHM is the full-width half-maximum of the gamma photo peak.

Table 1. Experimental materials for the absolute efficiency and energy resolution.

Description	Make	Model
10-microcurie button sources	-	Validated 01 JUL 2003
1.5 × 2.5″ NaI Detector	Ortec	905-4
2.25 × 2.25″ HPGe Detector	Ortec	GEM25-70-ICS
Lead Well	Canberra	-
Lead Well	-	Lead Bricks
Maestro	-	Windows XP

presents the materials for the experimental setup. Each source was placed 3″ from the end of the detector face using a generic 3D-printed ABS stand. Additionally, each setup took place inside of at least 2″ of lead to minimize background radiation. Appendix B provides the tabulated results for the isotopes Co-60, Ba-133, and Cs-137. Figure 6 displays the result of the absolute efficiency, Figure 7 shows the energy resolution for both detectors, and Figure 8 provides a rescaled plot of the HPGe energy resolution.

3. Benchmarking

The benchmarking of the algorithm was conducted with NaI and HPGe spectra of approximately 10 ± 0.234 microcuries each of Co-60, Cs-137, and Na-22 at 100 s, 300 s, and 600 s. All sources were sealed cylindrical buttons with a radius of 0.5 inches. Additionally, the sources were validated on 1 July 2003. To ensure the spectra collected for benchmarking occupied approximately two-thirds of the Maestro viewing window, a 795 bias Voltage/2.4 Gain was selected for NaI and 2000 Volts/2.4 Gain for the HPGe. These settings are subjective and not used in the model, but allow for the user to adequately observe all photo peaks and Compton events in benchmarking the model. Furthermore, the spectra collected for benchmarking the simulation required a collimated beam of gammas to ensure even energy deposition and the loss of gammas between the source and scintillating material. An uncollimated beam of gammas allowed the gammas to spread isotropically. The collimator is a 2-inch × 4-inch × 4-inch block of lead with a 0.25-inch diameter hole bored through the entirety and a 1-inch diameter hole centered on the smaller part and bored 1 inch deep to hold the source, as shown in Figure 9. A similar approach was taken for the HPGe detector; however, due to the orientation of the detector, a longer neck was required on the 3D-printed source/collimator housing, as shown in Figure 10. It is critical to note that the experimental setup for both detector types took place within a lead well to minimize the detection of background radiation in the benchmark spectra. Figure 11 displays the dimensions of the 3D-printed ABS collimator holder, shown in yellow. The 2″ × 4″ × 4″ lead collimator fits within the holder. The source and the detectors are placed on opposite ends of the holder.

Figure 12 compares the difference between the simulation and a 600-s benchmark spectrum for a 20-year-old, 10-microcurie sample of Co-60. Figure 13 provides a relative error between the experimental results and the Purdue model. The upper panel displays the overall relative error for all channels, the middle panel shows the channels relevant to the Compton continuum, and the lower panel displays the channels relevant to the photo peaks.

At the surface, the model overlaps with the experimental data when considering the main architecture of the spectra. The Compton continuum and the photo peaks are nearly aligned when viewing all channels. There is some variation in the counts between the model and experimental data, indicating that the total number of photons generated in the code is not fully accurate compared to the activity of the actual Co-60 sample. This is to be expected, as the sample activity follows an exponential decay trend over time. However, the exact sample activity cannot be predicted with absolute certainty due to the random behavior of decay. The closest approach to knowing the exact activity at any given time is to measure it for as long as reasonably achievable to minimize the uncertainty in the measurement. Currently, the decay constant variables set in the code are two significant figures. The lack of high precision in the variables produces variation between the generated and sample activities, especially given that the samples are over 20 years old. In future editions of the code, this can be corrected by modifying the time constants for the decay by including the highest precision measurement available. Although, there will never fully be a perfect match.

Considering the relative error between the data points, there is greater error in the low and high channels. In the channels under 150, this error can be accounted for by the lack of background radiation and X-rays in the code. This is again apparent in Figure 13, where these features are present in the experimental data but not the model. The higher channels experience a much larger relative error, with values of up to 400%. This occurrence can be explained by the low counts at a higher energy in the experimental data, and the noise generated in the code. In this particular case, the large value comes from the relative error between 4 counts from the background and 20 artificially generated counts. A different mode of presenting the error, such as the root mean square, may suppress the extremes. However, using the simple relative error shows the important features such as the Compton continuum and photo peaks appropriately when comparing the data spectra, as shown in Figure 12. In the lower panels of Figure 13, there is a relative error average of about 45%. For two signals with random noise, a 45% relative error is fair. A comparison between the two spectra of the experimental data produces relative errors on the same order of magnitude.

A similar comparison can be made for the HPGe model setting and experimental data. Figure 14 compares the difference between the simulation and a 100-s benchmark spectrum for a 20-year-old, 10-microcurie sample of Co-60. Figure 15 shows the relative error between the model and experimental data.

Again, considering the relative error between the model and experimental data, the average is approximately 80%. However, different features occur when comparing Figure 13 and Figure 15. The data overall have less variance, except for the photo peaks. Given that the channels are not exactly aligned and given the high resolution of the HPGe detector, any misalignment between the channels of each dataset will lead to large errors. This is apparent in the last panel of Figure 15, where the two datasets do not perfectly line up. Counter to the photo peaks, the significantly larger number of channels in an HPGe detector reduces the variance in the data for the Compton continuum.

It is critical to note that there is noise present in both the experimental and model data, which leads to what would normally be considered a high relative error. As a result, a comparison of how the two signals align with their spectra architecture and relative error must be considered together. By doing so, the evidence supports that the models in both the HPGe and NaI modes align reasonably well with their respective experimental detectors.

4. Graphical User Interface

The graphical user interface (GUI) offers a clean and organized means of adjusting the algorithm’s input variables. The user can select up to eight pre-loaded 1-microcurie gamma-emitting isotopes or create their own between 1 and 2 microcuries with an option of two initial gamma energies.

Table 2. List of isotopes pre-loaded in Purdue code.

Isotope Selection
Co-60	Th-228
Co-57	Am-241
Na-22	Au-198
Ba-133	Cs-137

provides a list of isotopes loaded into the code. Then, they can select the number of MCA channels available. Also included is the ability to select the isotope age in years or days depending on the initial isotope selected. The count time is adjustable from 1 to 3000 s. Furthermore, the user can select the detector’s resolution from 0 to 100 percent. This allows for the response of several detectors or to display the three-gamma interaction without resolution. Perhaps the most useful feature is the ability to download the spectrum of the output in a comma-separated value file. Allowing the user to manipulate the values allows for a similar experience and universal utility that occurs with an actual detector. User input variables increase the utility of the model for various applications seen fit by the user.

Earlier versions of the algorithm relied on single-stream variable computation. This significantly slowed the computation of several combinations of user input variables. The solution was to add parallelization to the script. An example of this was running a 20-year-old Cs-137 for a count time of 600 s. Without parallelization, the algorithm required 2 min and 27 s to return a response. The addition of multiple cores reduced that time to just under 3 s. However, a consequence of running multiple cores is the requirement of a longer startup time. For newer systems, it can take 10 s to start multiple cores, and up to 60 s for older systems running MATLAB. This is performed in the background of the GUI when the power button is pressed, but this “dead time” can be analogous to the application of voltage to a real detector. Figure 16 shows the options available to the end user.

Currently, the model does not allow for the manipulation of the geometric features of the detector, such as distance from the source or the angle relative to the crystal. The only user input as of now related to the geometry is the size of the scintillating crystal. The options allow for a 1″ × 1″ crystal, 2″ × 2″ crystal, 3″ × 3″ crystal, or an “infinite” crystal. Note that the infinite crystal in the algorithm is set at 100″ × 100″. This input is based on volumetric differences in the material, where the 3″ × 3″ crystal is the benchmark standard.

5. Results

The goal of the model was to provide flexibility to the user to simulate a wide range of settings that are present on actual detectors. The resolution feature can be seen between different detector brands and different scintillating materials.

An example of the versatility can be expressed for a 20-year-old, 10-microcurie sample of Cs-137 counted for 600 s for different resolutions. Critical to note, given the two extreme ends of the resolution slider representing the “perfect” version of either NaI or HPGe, is that any value in-between can be selected to match that of an actual detector. In testing, it was found that the NaI detector closely resembled the code at a 20% resolution. Subsequently, Figure 17 shows the result at a 100% resolution (perfect HPGe) and Figure 18 shows the result at a 20% resolution (testing NaI).

Another example expressing the flexibility of the model is the ability to modify the user-defined sample age. Take Ba-133, which has a half-life of 10.51 years. By modifying the user-defined sample age from 5 years to 50 years (approximately five half-lives), the code represents the decline in the simulated sample activity proportional to how two actual samples of 5-year-old and 50-year-old Ba-133 would respond. The responses of the code for the user input are displayed in Figure 19 and Figure 20, respectively.

The model is also capable of producing a gamma spectrum with a combination of up to three isotopes and an option to include background radiation. The background spectra are difficult to produce because of the large variety of naturally occurring isotopes for given locations, building types, or even altitudes from ground level. However, a basic spectrum representing the background counts from U-235 and Ra-226 decay chains, as well as naturally occurring K-40, Cu-63, and Cu-65, is present. The model assumes a 0.25 µCi activity of all the isotopes with a half-life of 1000 years. While none of these values are truly accurate, the background spectrum produced is consistent, with many background spectra [14]. Figure 21 shows a combination spectrum of 1 microcurie each of Ba-133, Co-60, and Na-22 at an age of 1 year and a count time of 600 s. Conversely, Figure 22 displays the same combination at 100 years, which is well beyond the five half-lives of all the isotopes. The user interface allows for the inclusion of the background counts with a selection of a checkbox. For most applications, the background is much lower than the counts and does not affect the results of the spectrum, but it is evident that it plays a role when all of the isotopes have passed five half-lives and have decayed to nearly null activity.

The advantage of making use of hybrid Monte Carlo methods is the reduced computational times required to conduct calculations. The model was run on two systems and timed to highlight the low computational times. The first system was a Dell OptiPlex with an i7 14700 vPro processor and 16 computational cores. The other was an HP EliteBook running an i7-8550U with four computational cores. The results of the different combinations of the model running for Co-60 are displayed in

Table 3. Computation time in seconds for variations in Co-60 between two processors.

	Dell i7 14700 vPro	HP i7-8550U
1 µCi/1 Year/100 s/256 Channels Co-60	0.47466	1.64626
1 µCi/1 Year/100 s/2048 Channels Co-60	0.485047	2.095934
1 µCi/1 Year/1 s/2048 Channels Co-60	0.244663	0.827295
1 µCi/1 Year/3000 s/2048 Channels Co-60	5.304477	15.77539
1 µCi/0.1 Years/100 s/2048 Channels Co-60	0.353323	1.400918
1 µCi/100 Years/100 s/2048 Channels Co-60	0.134576	0.767599
0.1 µCi/1 Year/100 s/2048 Channels Co-60	0.102322	1.072231
10 µCi/1 Year/100 s/2048 Channels Co-60	1.042232	6.820869

. Different combinations of isotopes are shown in

Table 4. Computation time in seconds for variations in isotopes between two processors.

	Dell i7 14700 vPro	HP i7-8550U
1 µCi/1 Year/100 s/2048 Channels Co-60	0.288089	1.437798
1 µCi/1 Year/100 s/2048 Channels Cs-137	0.34729	1.593292
1 µCi/1 Year/100 s/2048 Channels Na-22	0.131442	1.009565
1 µCi/1 Day/100 s/2048 Channels Au-198	0.127421	0.944222
1 µCi/1 Year/100 s/2048 Channels Th-228	0.166742	1.339904
1 µCi/1 Year/100 s/2048 Channels Co-60, Cs-137	0.347928	2.203138
1 µCi/1 Year/100 s/2048 Channels Co-60, Cs-137, Na-22	0.405756	2.93693

.

6. Comparison with GADRAS

One of the Oakridge National Lab Codes that attempts to solve the detector response is Gamma Detector Response and Analysis (GADRAS) [15]. Similar to the MCNP, GADRAS requires a vetting process that often takes months and has strict export control parameters. Additionally, the popularity of the software is much less than for the MCNP and has sparse third-party documentation. As a result, the user manual is the only means to understand the software.

Objectively, GADRAS is much more robust than the Purdue model. The user is allowed to create unique detectors with specific parameters, the background radiation of multiple locations globally can be considered, and there is a large library of sources to include neutrons. However, GADRAS is an analytical solution that produces a solution that does not vary between runs and does not replicate detector noise. The Purdue model has an energy resolution option, which can theoretically be set to match any scintillating material provided in the detector selection and parameter page. At this point in development, the user can only adjust the size of the scintillating crystal in the Purdue code. Another limiting difference in the Purdue model is the small number of isotopes. Currently, there is an option to combine up to 3 of the isotopes combined at any given run, which effectively provides 336 combinations of isotopes.

When comparing the two models for a given isotope, there are a few critical differences. GADRAS allows for the energy calibration of the detector, while the Purdue code does not. As a result, the channels in which the energy peaks are displayed are distorted. The Purdue model is slightly expanded when compared to GADRAS. This effect is a direct effect of the bias voltage used when the samples were collected for benchmarking and built into the code. GADRAS cannot vary the counting time and defaults to 10 s of live time. To make a useful comparison between GADRAS and the Purdue code using a 3″ × 3″ NaI crystal, the Purdue code was set to the following: 1 µCi Co-60, a 1-year-old sample, and a counting time of 10 s. The raw results are shown in Figure 23, and the results to highlight the differences in the spectrum compression are displayed in Figure 24 by aligning the 1.17 MeV peaks. This was performed by finding the channel with the peak in the GADRAS code, then artificially setting the same peak of the Purdue code at that channel and adjusting the rest accordingly. Due to this, the Purdue code has negative channel numbers, but the absolute value of the channels remains 512.

To further highlight the similarities between the outputs, the standard HPGe detector was set up in GADRAS for the same Co-60 sample parameters listed previously. However, the energy resolution in the Purdue code was set to 100%. Figure 25 shows the raw output of both files and Figure 26 displays an adjusted channels aligned at the 1.17 MeV peaks. The same apparent data compression is apparent with the energy resolution, showing that channel number is not a function of the resolution in this code.

An additional comparison that can be made between the codes is the performance with multiple isotopes. It should be noted that GADRAS does not combine the spectrum of each isotope on the displayed plot. Rather, this was conducted in a spreadsheet by summing the output of each isotope from the CSV file. The benefit of the Purdue code is that this process is performed automatically in the GUI. To highlight the similarities, 1-year-old samples of Co-60, Cs-137, and Na-22 at 1 microcurie were counted for 10 s with the standard GADRAS NaI detector. The results are shown below in Figure 27. Similar to before, Figure 28 shows an alignment of one select peak for better reference. The Cs-137 662 keV peaks are aligned.

7. Conclusions

The function of the Purdue code allows users to gain an understanding of radiation and the underlying physics without the requirement of having access to expensive equipment or the MCNP/GADRAS. The use of hybrid Monte Carlo methods decreases the computation requirements. A code that could take minutes to hours to run in traditional Monte Carlo codes can take a matter of seconds or less in the Purdue code. This is dependent upon the type of processor used, but the correlation between specific tasks on different machines produces a time change that is proportional between the two processors.

Future work for the Purdue code includes the addition of a robust isotope library that can match that of GADRAS. Currently, a user can apply up to three isotope combinations of the pre-loaded eight-isotope library, but this can be expanded to allow for more complex spectra to be created. Additionally, the inclusion of calibration or bias voltage-like features would help create spectra that can be expanded to match the calibration of any detector. Furthermore, the inclusion of neutrons via gas-filled detectors would make the Purdue code a more competitive model. However, there is no claim that the Purdue code can replace widely used transport codes. It may offer an alternative that can be used in an educational setting by more people, with comparable results to experimental data and data produced by industry-leading gamma detector response codes.