MAXIMUM LIKELIHOOD ESTIMATION USING EXPECTATION MAXIMIZATION APPLIED TO AMBIENT DOSE EQUIVALENT MEASUREMENTS

MAXIMUM LIKELIHOOD ESTIMATION USING EXPECTATION MAXIMIZATION APPLIED TO AMBIENT DOSE EQUIVALENT... Abstract The use of spectrometric detectors for ambient dose equivalent measurements requires the development of calculation procedures to obtain the incident photon fluence spectra from the measured pulse-height spectra. In this work, the implementation of a methodology based on the maximum likelihood estimation using expectation maximization is described. The method has been validated using simulated and reference photon beams. The results are in excellent agreement with the reference values in the framework of the involved uncertainties. The presented method allows the consideration of the detector resolution, improving thus the energy resolution of the obtained incident photon spectra. INTRODUCTION The appropriate quantity to estimate the effective dose in area monitoring is the ambient dose equivalent, H*(10). This operational quantity was introduced by the International Commission on Radiation Units and Measurements(1) and its use is mandatory in the European Union according to EU Directive 2013/59 EURATOM(2). There are thousands of radiological early warning network stations in Europe in which typically Geiger Müller detectors are being used for ambient dose equivalent rate measurements. Moreover, some detector systems designed for other quantities, such as exposure or air kerma are also in use. Due basically to their inability to provide information about the incident photon spectra and to the pronounced energy and angular response to gamma radiation of these detectors in terms of H*(10), new compact and handy spectrometric systems based on scintillation detectors like LaBr3, CeBr3 and SrI2 or on room-temperature semiconductor detectors like CdZnTe are installed instead of or in addition to the traditional dose rate meters. These spectrometers will considerably improve the information provided by early warning networks in real-time, in particular concerning radionuclide identification(3). The use of spectrometric detectors requires the development of calculation procedures to obtain the H*(10) rates from the measured photon pulse-height spectra. A comprehensive description of the so far applied methods, including their advantages and disadvantages, has been provided by Dombrowski(4). A different approach, based on the maximum likelihood estimation using expectation maximization (ML-EM), has been recently implemented by the author in the framework of the project: ‘Metrology for Radiological Early Warning Networks in Europe’ (MetroERM project(5)), as part of the tasks related to the investigations of detector features and of spectra evaluation and deconvolution methods for new and improved measurement systems based on novel spectrometric detectors. The ML-EM is a statistical method initially developed for the use in the reconstruction of images in positron emission tomography(6, 7) and successfully tested for spectra deconvolution in gamma-ray spectrometry(8). The implementation of the ML-EM method for H*(10) rate calculations from measured gamma-ray spectra is explained in the present work. The associated mathematical formalism and the general features of the specific computer code DET-H10R, developed for spectra deconvolution and H*(10) rate calculations using the ML-EM algorithm, are presented. The incorporation of the detector resolution during the unfolding of measured photon spectra using the code DET-H10R is described and the results of theoretical and experimental verifications of the implemented method are also presented. MATERIALS AND METHODS General concepts After background subtraction, the measured pulse-height spectrum O(E) can be represented by the integral: O(E)=∫0∞R(E,E0)⋅ΦE(E0)⋅dE0 (1) where R(E, E0) is the detector response function and ΦE(E0) is the distribution of the incident photon fluence in energy. The discrete form of this equation is used in practical applications: Oi=Ri,j⋅Φj{i=1,…,n;j=1,…,m} (2) where Oi is the number of counts in channel i of the measured pulse-height spectrum, Ri,j is the number of counts expected in channel i of the measured pulse-height spectrum per unit of photon fluence with energy within the energy bin j, and Φj is the unknown fluence of photons with energy in bin j. The number of channels in the measured spectrum is denoted as ‘n’, the number of energy bins of the photon fluence spectrum as ‘m’. The summation convention has been applied, i.e. the repeated indices in the right term are implicitly summed over. The matrix Rˆ(n,m) in equation (2) represents the detector response operator and is also known as ‘detector response matrix’. The ‘noise’ contribution due, for instance, to the counting statistics, the lack of accuracy in the determination of the detector response operator and its finite resolution, has been omitted in equation (2) and is therefore neglected in the analysis. The unknown fluence vector Φ→ in equation (2) generates, and is to be calculated from, the observed pulse-height spectrum O→ by using a known detector response matrix. The ‘noise’ referred in the previous paragraph contributes to the uncertainty of these calculations. Once the photon fluence spectrum is determined, the air kerma, Kair, can be calculated using the following equation: Kair=∑j=1mΦj⋅Ej⋅(μtrρ)air,j (3) where Φj is the fluence of photons with energy in bin j, Ej is the mean energy of the energy bin j and (μtrρ )air,j is the mass energy transfer coefficient in air for photons with energy Ej. The radiation protection quantities can be also calculated from the photon fluence spectrum. In particular, the ambient dose equivalent can be obtained as follows: H*(10)=∑j=1mΦj⋅Ej⋅(μtrρ)air,j⋅hk*(10,Ej) (4) where hk*(10,Ej) is the air kerma to ambient dose equivalent conversion coefficient corresponding to the energy Ej. The rest of the parameters are the same as in equation (3). One of the requirements of the ML-EM method is that the response matrix Rˆ(n,m) in equation (2) is needed for a given detector. It could be obtained by Monte Carlo simulation after a detailed characterization of the detection system, i.e. the geometry and material composition of the elements that constitute the detection probe are to be well known. Then, the response matrix is filled with m detector specific pulse-height spectra of n channels, where each column j contains the pulse-height spectrum produced by incident photons with an energy within the energy bin j. In what follows, some simple transformations of equation (2) are introduced in order to apply the ML-EM method(6, 7) to this particular case, in which the incident photon fluence spectrum and the corresponding ambient dose equivalent are calculated from the measured pulse-height spectrum. In order to simplify the formulations, the summation convention is applied for repeating indices ‘i’ or ‘s’ in matrix products. Basic transformations required An important requirement for the ML-EM method is that the sum of the elements in each column of the response operator should be lower than or equal to one. In order to achieve this, a total efficiency operator RT^(m,m){RTj,l=0∀(j≠l);RTj,j=RTj} is derived from the detector response matrix Rˆ(n,m), and the vector of incident photons which are detected, I→, is defined according to the following equation: Ij=RTj⋅Φj{j=1,…,m} (5) where RTj is the total number of counts expected in the measured spectrum per unit photon fluence in energy bin j and Φj is the fluence of photons with energy in bin j. The RTj values are calculated from the detector response matrix by summing the elements of each column: RTj=∑i=1nRi,j{j=1,…,m} (6) Combining equations (2) and (5), the following relationship is obtained: Oi=Ci,s⋅Is{i=1,…,n;s=1,…,m} (7) where Ci,s represents the fraction of detected photons, incident with energy s, that are detected in channel i: Ci,j=Ri,jRTj{i=1,…,n;j=1,…,m} (8) As can be noted from equations (6) and (8), the total number of incident photons that are detected is equal to the number of counts in the measured spectrum, i.e.: ∑i=1nCi,j=1 (9) This is true when coincidence summing effects can be neglected, which is the case when in-situ measurements are performed. This is also true for low ambient dose equivalent rates for which the dead time corrections of the detection systems are not significant. The ML-EM method will be applied to solve equation (7), i.e. to calculate the vector I→ from the measured pulse-height spectrum O→ by using the normalized detector response matrix, Cˆ(n,m), whose elements are defined by equation (8). The calculation of the incident fluence spectrum is then straightforward from equation (5). The ML-EM algorithm For a given detection system and for each possible vector I→, the measured spectrum O→ has a probability or likelihood: L(I→)=P(O→|I→) (10) The idea behind the implemented method is to find the mathematical expectation of the vector Í→ which maximizes the probability of measuring the actual pulse-height spectrum O→ over all possible incident photon spectra, i.e. max(L(I→))=L(Í→). It is obvious that for a given detector response matrix Rˆ(n,m) this solution also provides the maximization of the likelihood L(Φ→)=P(O→|Φ→) because: Φ′j=ÍjRTj{∀j∈mˆ} (11) The likelihood function deduced by Shepp and Vardi(6) has been chosen for L(I→), based on the assumption that the components of the involved photon spectra, representing a certain number of counts in a given time interval, can be considered independent Poisson variables. The iterative scheme based on the expectation maximization (EM) algorithm(6, 7) is thus applied to find the maximum value of this likelihood function in equation (10), i.e. to find the vector Í→ for which the probability of having measured the spectrum O→ reaches the maximum value. Mathematical demonstrations of the method are discussed in detail in the pioneering work of Shepp and Vardi(6) and are beyond the scope of this work, in which the efforts have been aimed at the achievement of a compliance with the necessary requirements for its application. The equality according to equation (9) by using the normalized response operator defined by equations (7) and (8) makes the following iteration procedure simpler and more efficient. Starting with an initial guess of Í→: Íjr=0≥0{j=1,…,m} for which: ∑j=1mÍjr=0=∑i=1nOi (12) a new estimate of the spectrum Í→ is calculated at any step r as follows: Íjr=Íjr−1⋅[OiCi,s⋅Ísr−1⋅Ci,j]{i=1,…,n;s,j=1,…,m;n≥m}, (13) where, the likelihood function increases: P(O→|Ír→)≥P(O→|Ír−1→){r≥1} with equality if and only if L(Ír−1→)=max(L(I→)). As can be deduced from equation (12), the measured spectrum O→ could be taken as an initial estimate of the incident spectrum when n = m. The incident photon fluence spectrum is obtained at each iteration step using equation (11) and the corresponding value of the ambient dose equivalent is calculated according to equation (4), i.e. as follows: H*(10)r=∑j=1mΦjr⋅Ej⋅(μtrρ)air,j⋅hk*(10,Ej) (14) where H*(10)r and Φjr are respectively, the ambient dose equivalent and the fluence of photons with energy in bin j, obtained at iteration step r. The rest of the parameters have the same meaning as in equations (3) and (4). The iterative process stops when Δ=|H*(10)r−H*(10)r−1H*(10)r−1|≤δ (15) where δ is the accepted relative deviation for the two consecutive estimated values of H*(10). The ambient dose equivalent rate can be calculated using the live collection time of the spectrometric system. Implementation of the ML-EM method Although the ML-EM method is more computer-power-demanding than other simpler procedures(4), the involved algorithm, i.e. the iterative calculation given by equation (13), is computationally of low complexity and could be implemented on a low-power, industrial CPU, of a type which can be integrated in a hand-held instrument or measurement probe. It should be also considered that field instruments, for instance in early warning stations, could easily send the measured pulse-height spectra (few tens of kB for 1024 channels, in the most used formats) in real-time to an external central processing system. Thus, these spectra could be stored in a central memory unit, from which they could be read by the main computer with a predefined frequency for dose calculations, radionuclide identification, alarms triggering, etc. This could be the better option, even for apparently simpler or less-demanding H*(10) calculation procedures, taking into account that the CPU integrated in the measuring devises should be used rather to produce reliable pulse-height spectra, corrected for instrument instabilities, in particular for energy shifts in the measured peaks due to thermal instability of the electronic components. General features of the code DET-H10R The presented ML-EM method has been implemented in the code DET-H10R, written in free PASCAL for the calculation of air kerma rates and H*(10) rates from measured photon pulse-height spectra in PC platforms. The calculated spectra of the incident photon fluence rate, air kerma rate and H*(10) rate are also provided by the code and can be saved in text files. The values of the mass energy transfer coefficient in air have been obtained from the publication of Hubbell and Seltzer(9) and the values of the air kerma to ambient dose equivalent conversion coefficient hk*(10,E) have been taken from the ICRU Report 57(10). To facilitate the calculations, the discrete values of the provided coefficients were fitted with empirical analytical functions. The maximum relative deviation between theoretical and calculated values is below 0.3% with 1 − P ≤ 10−5, where, P is the Pearson r-moment. The supported file formats for the measured spectra are the ISO N42.42 and the binary CAM CNF format from CANBERRA. A code-specific text format has been also defined for users to upload those measured pulse-height spectra with formats not supported by the code. Considering that the resolution, i.e. the number of keV per channel, of measured pulse-height spectra should be equal to or better than the resolution of the constructed detector response matrix, a straightforward procedure has been incorporated into the code DET-H10R to fit the energy bin width of measured spectra with the energy bin width of the corresponding detector response matrix. Starting from the lower energy, this procedure gradually groups the energy channels of measured spectra according to the resolution of the response matrix. For instance, if a measured spectrum has energy channels of 3 keV width and the energy bin width of the response matrix is 10 keV, the first energy bin of the adjusted spectrum will contain the sum of the counts from the first three channels plus one-third of the counts from the fourth channel of the measured spectrum. The second energy bin of the adjusted spectrum will contain two-thirds of the counts from the fourth channel plus the counts from the fifth and the sixth channels, plus two-thirds of the counts from the seventh channel, and so on. A graphical user interface, including a user help, facilitates the selection of the detector response matrix, the uploading of the measured pulse-height spectra and the definition of the relevant calculation parameters, like the ‘iteration resolution’ and the ‘cut-off energy’. The iteration resolution is the accepted relative deviation between two successively calculated values of the H*(10) rate, as defined by equation (15). The cut-off energy is the lower energy considered during the calculations of air kerma and H*(10) rates. This option could be useful to discard the contribution of low-energy incident photons in those cases where the energy bin width (resolution) of the detector response matrix could be insufficient to properly consider the interaction processes of these photons or to adequately reflect the pronounced energy dependence of the involved photon interaction coefficients(11) and the air kerma to ambient dose equivalent conversion coefficient, in the low-energy region. In these cases, the response matrix should be constructed with the appropriate resolution, after a more detailed characterization of the detector geometry and materials. The code interface also incorporates tools to analyse the measured pulse-height spectra and the calculated H*(10) rate spectra (see the example in Figures 1 and 2) facilitating the verification of the correct energy calibration of spectrometers, the detection of possible energy shifts due to thermal instability of the electronic components and the radionuclide identification. Figure 1. View largeDownload slide Measured gross pulse-height spectrum using a LaBr3 (1″ × 1″) detector in a low dose rate 60Co collimated beam. Figure 1. View largeDownload slide Measured gross pulse-height spectrum using a LaBr3 (1″ × 1″) detector in a low dose rate 60Co collimated beam. Figure 2. View largeDownload slide Calculated H*(10) rate spectrum from the 60Co pulse-height spectrum shown in Figure 1. The detector response matrix with 10 keV energy bin width has been obtained by Monte Carlo simulation. The small fraction of scattered photons produced by the collimator of the 60Co source can be seen as well as a small broad peak around 1470 kV, corresponding to the external 40K background and the internal detector contamination with 138La(12). Figure 2. View largeDownload slide Calculated H*(10) rate spectrum from the 60Co pulse-height spectrum shown in Figure 1. The detector response matrix with 10 keV energy bin width has been obtained by Monte Carlo simulation. The small fraction of scattered photons produced by the collimator of the 60Co source can be seen as well as a small broad peak around 1470 kV, corresponding to the external 40K background and the internal detector contamination with 138La(12). Incorporation of the detector resolution One advantage of the ML-EM method is that a triangular matrix is not required in equation (2). Therefore, the detector resolution could be taken into account during the unfolding of measured pulse-height spectra. A convolution algorithm included in DET-H10R facilitates the incorporation of the detector energy resolution into the detector response matrix previously constructed by Monte Carlo simulation. The resulting response operator allows the calculation of more realistic incident photon fluence spectra, with improved resolution. A brief description of the implemented procedure is presented below. During the convolution procedure the detector response function in equation (1) is modified accordingly to the following equation: R′(E,E0)=∫0∞G(E−E′)⋅R(E′,E0)⋅d(E′) (16) where, R′(E,E0) is the response function for incident energy E0 considering the resolution of the detector system, G(E−E′) represents the probability per unit pulse height that a deposited energy E′ will produce counts on the energy channel E, and R(E′,E0) is the ideal or mathematical detector response function for incident energy E0, obtained by Monte Carlo simulation. The resolution operator G(E−E′) is approximated by a Gaussian distribution: G(E−E′)=12πσ2e−((E−E′)2/2σ2){σ≡σ(E′)} (17) where the standard deviation σ(E′) is calculated from the peak full width at half maximum (FWHM) corresponding to energy E′, according to the following relationship: σ(E′)=FWHM(E′)22ln2 (18) As in the case of equation (1), the discrete form of equation (16) is used: R′i,j=Gi,s⋅Rs,j{i,s=1,…,n;j=1,…,m} (19) Each column of the spreading operator Gˆ(n,n) is a Gaussian probability density function in form of histogram, with the most probable values in the diagonal elements. These columns are constructed by generating energy values, E, using the Box–Müller algorithm(13): E=E′+σ(E′)⋅sin(2π⋅γ1)⋅−2ln⋅γ2 (20) where E′ is an energy value sampled uniformly in the energy bin corresponding to the diagonal element in a given column. γ1 and γ2 are the independent random variables uniformly distributed in the interval (0,1). One event is counted in the corresponding energy bin per each energy value E generated. To obtain the probability density function, the counts in each energy bin are then divided by the total energy values generated for a given column. It can be noted that this normalization of the spreading operator implies that the counting statistics in the new response matrix is not altered. The number of values with an energy, E, to be generated in each column can be selected by the user according to the required statistical uncertainty. A graphical interface has been included to facilitate the input of the detector specific resolution function: FWHM≡FWHM(E′) determined from the energy resolution calibration. The following two general equations are provided to facilitate the definition of a wide variety of possible FWHM functions(14) by selecting the adequate fitting parameters. FWHM(E′)=a⋅b+c⋅E′+d⋅E′2 (21) FWHM(E′)=a+b⋅E′+c⋅E′2+d⋅E′ (22) The convolution process indicated by equation (19) readjusts the distribution of the deposited energy in each pulse-height spectrum that constitutes the detector response matrix, producing the broadening of this pulse-height spectrum. RESULTS AND DISCUSSION Calculations with a simulated spectrum The ML-EM method has been applied to a relatively complex photon pulse-height spectrum generated by Monte Carlo simulation, using the model of a LaBr3 (5% Ce) (1.0″ × 1.0″) detector. A parallel photon beam, perpendicular to the main detector axis and constituted by 10 gamma-energy lines was simulated. The beam covered the entire probe volume which was exposed to a total ambient dose equivalent rate of 500 nSv/h. The photon fluence rate values were selected so that the ambient dose equivalent rate corresponding to each energy line was 50 nSv/h (Table 1). The number of simulated photons results from the adopted acquisition time of the spectrometer and the area of the beam cross section. The simulated pulse-height spectrum (Figure 3) does not consider the statistics associated with the process between the deposition of the photon’s energy and the counting of the corresponding event in the spectrum. Therefore, in order to obtain a more realistic spectrum, the simulated spectrum was convolved with the resolution operator constructed from the specific detector resolution function shown in Figure 4. The resulting spectrum, to which the ML-EM method was applied, is shown in Figure 5. Table 1. Simulated gamma lines, fluence rates and calculated H*(10) rates. Energy Fluence rate Ḣ*(10)ML−EM keV cm−2 s−1 nSv/h 59,54 (241Am) 27.23 50.5 ± 1.6 134 (144Ce) 17.43 49.5 ± 1.1 364 (131I) 6.54 49.7 ± 1.1 609 (214Po) 3.96 49.9 ± 1.1 662 (137Cs) 3.68 48.9 ± 1.2 796 (134Cs) 3.16 49.2 ± 1.2 1173 (60Co) 2.34 49.2 ± 1.6 1332 (60Co) 2.13 49.0 ± 1.6 1460 (40K) 2.00 51.1 ± 2.1 1764 (214Po) 1.75 49.5 ± 2.0 Total — 496,5 ± 4,8 Energy Fluence rate Ḣ*(10)ML−EM keV cm−2 s−1 nSv/h 59,54 (241Am) 27.23 50.5 ± 1.6 134 (144Ce) 17.43 49.5 ± 1.1 364 (131I) 6.54 49.7 ± 1.1 609 (214Po) 3.96 49.9 ± 1.1 662 (137Cs) 3.68 48.9 ± 1.2 796 (134Cs) 3.16 49.2 ± 1.2 1173 (60Co) 2.34 49.2 ± 1.6 1332 (60Co) 2.13 49.0 ± 1.6 1460 (40K) 2.00 51.1 ± 2.1 1764 (214Po) 1.75 49.5 ± 2.0 Total — 496,5 ± 4,8 Table 1. Simulated gamma lines, fluence rates and calculated H*(10) rates. Energy Fluence rate Ḣ*(10)ML−EM keV cm−2 s−1 nSv/h 59,54 (241Am) 27.23 50.5 ± 1.6 134 (144Ce) 17.43 49.5 ± 1.1 364 (131I) 6.54 49.7 ± 1.1 609 (214Po) 3.96 49.9 ± 1.1 662 (137Cs) 3.68 48.9 ± 1.2 796 (134Cs) 3.16 49.2 ± 1.2 1173 (60Co) 2.34 49.2 ± 1.6 1332 (60Co) 2.13 49.0 ± 1.6 1460 (40K) 2.00 51.1 ± 2.1 1764 (214Po) 1.75 49.5 ± 2.0 Total — 496,5 ± 4,8 Energy Fluence rate Ḣ*(10)ML−EM keV cm−2 s−1 nSv/h 59,54 (241Am) 27.23 50.5 ± 1.6 134 (144Ce) 17.43 49.5 ± 1.1 364 (131I) 6.54 49.7 ± 1.1 609 (214Po) 3.96 49.9 ± 1.1 662 (137Cs) 3.68 48.9 ± 1.2 796 (134Cs) 3.16 49.2 ± 1.2 1173 (60Co) 2.34 49.2 ± 1.6 1332 (60Co) 2.13 49.0 ± 1.6 1460 (40K) 2.00 51.1 ± 2.1 1764 (214Po) 1.75 49.5 ± 2.0 Total — 496,5 ± 4,8 Figure 3. View largeDownload slide Simulated spectrum in a LaBr3 (1.0″ × 1.0″) detector: 1024 channels, acquisition time = 600 s, total H*(10) rate = 500 nSv/h. Figure 3. View largeDownload slide Simulated spectrum in a LaBr3 (1.0″ × 1.0″) detector: 1024 channels, acquisition time = 600 s, total H*(10) rate = 500 nSv/h. Figure 4. View largeDownload slide Resolution function applied in the simulations of the test spectrum with the LaBr3 (1.0″ × 1.0″) detector. Figure 4. View largeDownload slide Resolution function applied in the simulations of the test spectrum with the LaBr3 (1.0″ × 1.0″) detector. Figure 5. View largeDownload slide Test spectrum obtained by convolving the simulated spectrum (Figure 3) with the resolution operator corresponding to experimental data of the studied LaBr3 (1.0″ × 1.0″) detector. Number of channels = 1024, acquisition time = 600 s, total H*(10) rate = 500 nSv/h. Figure 5. View largeDownload slide Test spectrum obtained by convolving the simulated spectrum (Figure 3) with the resolution operator corresponding to experimental data of the studied LaBr3 (1.0″ × 1.0″) detector. Number of channels = 1024, acquisition time = 600 s, total H*(10) rate = 500 nSv/h. The detector response matrix was obtained by Monte Carlo simulation using an energy bin width of 10 keV and assuming the same geometry, i.e. perpendicular beam incidence to the longitudinal detector axis. According to the work of Camp and Vargas(15), a bin width of 10 keV constitutes a good trade-off between energy resolution and counting statistics during ambient dose equivalent measurements with spectrometric systems based on scintillation detectors. The energy span considered was 3000 keV, therefore 300 spectra were generated. For each spectrum 109 photon histories were simulated in order to achieve a statistical standard uncertainty below 1% in each energy bin. The cut-off energies were 1 keV for photons and 10 keV for electrons. No additional variance reduction techniques have been applied in order to minimize the possibility of biased results, therefore 64 processors (cores) have been used (3 GHz, 2 GB ram/core). The generation of the test spectrum and the detector response matrix were carried out with the code MCNPX(16), using the photon library mcplib04 and the electron library el03. The incident ambient dose equivalent rate spectrum was calculated by unfolding the test spectrum (Figure 5) with the ML-EM algorithm. As expected, the resolution of the calculated incident spectrum using the ideal detector response matrix, i.e. without including the detector energy resolution, was similar to the resolution of the pulse-height spectrum, as can be seen in Figure 6. In this case, the peaks related to the 609 and 662 keV gamma-energy lines are not well resolved and consequently the ambient dose equivalent rate corresponding to each of these energies cannot be determined accurately. Figure 6. View largeDownload slide H*(10) rate spectrum calculated by unfolding the test spectrum (Figure 5) with the ideal detector response matrix, i.e. without considering the energy resolution of the detector. Figure 6. View largeDownload slide H*(10) rate spectrum calculated by unfolding the test spectrum (Figure 5) with the ideal detector response matrix, i.e. without considering the energy resolution of the detector. Nevertheless, the resolution of the recalculated incident spectrum was considerably improved when the energy resolution of the detector (Figure 4) was incorporated into its response matrix using the previously described convolution algorithm. The obtained spectrum is shown in Figure 7. In this case, the ambient dose equivalent rate corresponding to each incident photon-energy line can be easily obtained as the area of the related peak in the spectrum. It should be indicated that the energy resolution in the obtained spectra could be further improved by using detector response matrixes with lower energy bin widths, i.e. with better resolution, in detriment, however, of the counting statistics. Figure 7. View largeDownload slide H*(10) rate spectrum calculated by unfolding the test spectrum (Figure 5) with the spread detector response matrix, i.e. taking into account the energy resolution of the detector. Figure 7. View largeDownload slide H*(10) rate spectrum calculated by unfolding the test spectrum (Figure 5) with the spread detector response matrix, i.e. taking into account the energy resolution of the detector. The calculated ambient dose equivalent rates ( Ḣ*(10)ML−EM) are shown in Table 1. The indicated uncertainties have been expressed as expanded uncertainties with a coverage factor k = 2 and have been estimated considering the contributions of the counting statistics in the simulated 1024 channel spectrum and the finite resolution (10 keV per channel) of the detector response matrix. The differences of the calculated ambient dose equivalent rates and the expected values, i.e. 50 nSv/h for each photon energy and 500 nSv/h for the total spectrum, are consistent with the involved uncertainties. Measurements in reference photon beams As part of the validation studies, the implemented method has been applied to experimental spectra obtained at the Metrology Laboratory for Ionizing Radiation (LMRI) of CIEMAT. Two probes with spectrometric detectors (LaBr3 (1.5″ × 1.5″) and CeBr3 (1.5″ × 1.5″)) were characterized and their response matrixes were calculated by Monte Carlo simulation assuming a perpendicular beam incidence to the longitudinal detector axis and applying the same procedures and parameters as in the previous calculations. As a result, for each detector system a 300 × 300 response matrix with a statistical standard deviation better than 1% in each energy bin was obtained. Both detectors were exposed to known values of H*(10) rates in 137Cs and 60Co reference calibration beams at LMRI. The longitudinal axes of the detector probes were placed perpendicular to the axes of the photon beams. At each calibration point, five gross spectra and five background spectra were collected during 120 s each. The ML-EM method was applied to each measured spectrum to calculate the corresponding H*(10) rates. The calculations were performed in fractions of a second using a personal computer with an i7-4790 CPU @ 3.60 GHz processor for an iteration resolution better than 10−3%. The net measured H*(10) rate was obtained by subtracting the mean value of the background H*(10) rates from the mean value of the gross H*(10) rates. The results are shown in Table 2. The laboratory reference values Ḣ*(10)REF have a relative standard uncertainty of 2.2% and the relative statistical standard uncertainty of the measured H*(10) rate values Ḣ*(10)ML−EM is lower than 1%. The results provided by the implemented ML-EM algorithm are in excellent agreement with the reference values, considering the involved uncertainties. Table 2. Reference and calculated H*(10) rates at the LMRI radiation beams. Detector Source Ḣ*(10)REF Ḣ*(10)ML−EM Relative difference % μSv h−1 μSv h−1 LaBr3 137Cs 1.00 1.01 1.0 2.00 1.99 −0.5 5.00 4.99 −0.2 10.0 9.99 −0.1 60Co 2.00 2.03 1.5 CeBr3 137Cs 0.842 0.847 0.6 2.00 2.01 0.5 5.00 4.92 −1.6 10.0 9.93 −0.7 60Co 5.00 4.99 −0.2 Detector Source Ḣ*(10)REF Ḣ*(10)ML−EM Relative difference % μSv h−1 μSv h−1 LaBr3 137Cs 1.00 1.01 1.0 2.00 1.99 −0.5 5.00 4.99 −0.2 10.0 9.99 −0.1 60Co 2.00 2.03 1.5 CeBr3 137Cs 0.842 0.847 0.6 2.00 2.01 0.5 5.00 4.92 −1.6 10.0 9.93 −0.7 60Co 5.00 4.99 −0.2 Table 2. Reference and calculated H*(10) rates at the LMRI radiation beams. Detector Source Ḣ*(10)REF Ḣ*(10)ML−EM Relative difference % μSv h−1 μSv h−1 LaBr3 137Cs 1.00 1.01 1.0 2.00 1.99 −0.5 5.00 4.99 −0.2 10.0 9.99 −0.1 60Co 2.00 2.03 1.5 CeBr3 137Cs 0.842 0.847 0.6 2.00 2.01 0.5 5.00 4.92 −1.6 10.0 9.93 −0.7 60Co 5.00 4.99 −0.2 Detector Source Ḣ*(10)REF Ḣ*(10)ML−EM Relative difference % μSv h−1 μSv h−1 LaBr3 137Cs 1.00 1.01 1.0 2.00 1.99 −0.5 5.00 4.99 −0.2 10.0 9.99 −0.1 60Co 2.00 2.03 1.5 CeBr3 137Cs 0.842 0.847 0.6 2.00 2.01 0.5 5.00 4.92 −1.6 10.0 9.93 −0.7 60Co 5.00 4.99 −0.2 In order to minimize the influence of the detector angular response in the evaluation of the implemented method the irradiations were performed for a similar geometry to the one used during the calculations of the response matrixes, i.e. for perpendicular beam incidence. According to Camp and Vargas(15) a nearly isotropic response of this kind of instruments, i.e. cylindrical scintillation detectors with the same diameter and length, could be assumed in environmental monitoring. In general, depending also on the probe materials and their thicknesses, the angular dependence of the response is not negligible at lower energies. For a given detector and energy interval, the evaluation of the uncertainty component due to the angular dependence of the response could be carried out by Monte Carlo simulations. An average response matrix, weighted with the relevant angles of incidence, could then be used for H*(10) rate calculations. CONCLUSIONS In this paper the ML-EM method has been implemented for ambient dose equivalent measurements with spectrometric detectors. The results obtained during the verifications carried out with a simulated pulse-height spectrum and in reference photon beams are in excellent agreement with the expected values, considering the involved uncertainties. The incident gamma-ray fluence spectrum can be obtained in fractions of a second by unfolding the measured spectrum and the ambient dose equivalent can be accurately calculated using the relevant energy dependent coefficients, i.e. the mass energy transfer coefficient in air and the air kerma to ambient dose equivalent coefficient. A flat-energy response with respect to the quantity H*(10) is thus expected when using well-characterized spectrometric detectors. The application of the ML-EM method, like that of any other H*(10) calculation procedure using measured pulse-height spectra, requires the knowledge of the detector response matrix, as described above. Unfortunately, detector response matrixes are not provided by the manufacturers of spectrometric systems, and therefore they must be created by the users. The detector response matrix may be determined very accurately by Monte Carlo simulations, considering that a detailed characterization of the detection probe geometry and materials is possible. Although these simulations could be time demanding, they need to be carried out only once for a particular detection system. Well-characterized spectrometric detectors, with an isotropic response, can be used for ambient dose equivalent measurements without any knowledge of the incident photon field. Furthermore, the possibility of incorporating the detector resolution into the response matrix allows the improvement of the energy resolution in the calculated photon spectra, facilitating the identification of the involved radionuclides and the evaluation of their contributions to the measured ambient dose equivalent rate. Due to the diversity of currently available methods for ambient dose rate calculations from measured pulse-height spectra, the analysis of the advantages and disadvantages of the implemented ML-EM method in comparison with each of these existing procedures is beyond the scope of this work. The readers have been provided with the information required to evaluate the adequacy of the described approach for their particular applications or in their particular cases. Moreover, several comparison exercises that have been carried out in the framework of the MetroERM project among the implemented H*(10) calculation methods will be subject of further publications. ACKNOWLEDGEMENTS The author would like to thank CIEMAT for providing the working conditions required to work out this article. The author is also grateful to Dr José Antonio Suárez Navarro for his comments and suggestions. FUNDING This work was supported by EMRP within the JRP-Env57 (MetroERM) project. The EMRP is jointly funded by the EMRP participating countries within EURAMET and the European Union. REFERENCES 1 International Commission on Radiation Units and Measurements . Determination of dose equivalents resulting from external radiation sources. ICRU Report 39. ICRU Publications ( 1985 ). 2 European Council Directive 2013/59/EURATOM on basic safety standards for protection against the dangers arising from exposure to ionising radiation and repealing. Directives 89/618/Euratom, 90/641/Euratom, 96/29/Euratom, 97/43/Euratom and 2003/122/Euratom. OJ of the EU. L13; 57: 1–73 (2014). 3 Neumaier , S. , Dombrowski , H. and Kessler , P. Metrology for radiological early warning networks in Europe (‘METROERM’)—a join European metrology research project . Health Phys. 111 ( 2 ), 100 – 105 ( 2016 ). Google Scholar CrossRef Search ADS PubMed 4 Dombrowski , H. Area dose rate values derived from NaI or LaBr3 spectra . Radiat. Prot. Dosim. 160 ( 4 ), 269 – 276 ( 2014 ). Google Scholar CrossRef Search ADS 5 EUROPEAN ASSOCIATION OF NATIONAL METROLOGY INSTITUTES . Final Publishable Join Research Project Summary Report for ENV57 MetroERM—Metrology for Radiological Early Warning Networks in Europe. Available on https://www.euramet.org/research-innovation/search-research-projects. Issued in September ( 2017 ). 6 Shepp , L. A. and Vardi , Y. Maximum likelihood reconstruction for emission tomography . IEEE Trans. Med. Imaging MI-1 ( 2 ), 113 – 122 ( 1982 ). Google Scholar CrossRef Search ADS 7 Shepp , L. A. , Vardi , Y. and Kaufman , L. A statistical model for positron emission tomography . J. Am. Stat. Assoc. 80 ( 389 ), 8 – 20 ( 1985 ). Google Scholar CrossRef Search ADS 8 Meng , L. J. and Ramsden , D. An inter-comparison of three spectral-deconvolution algorithms for gamma-ray spectroscopy . IEEE Trans. Nucl. Sci. 47 ( 4 ), 1329 – 1336 ( 2000 ). Google Scholar CrossRef Search ADS 9 Hubbell , J. H. and Seltzer , S. M. Tables of X-Ray Mass Attenuation Coefficients and Mass Energy-Absorption Coefficients (version 1.4) ( Gaithersburg, MD : National Institute of Standards and Technology ) ( 2004 ) [Online] Available: http://physics.nist.gov/xaamdi. 10 International Commission on Radiation Units and Measurements . Conversion coefficients for use in radiological protection against external radiation. Bethesda, MN: ICRU Report 57 ( 1998 ). 11 International Commission on Radiation Units and Measurements . Fundamental quantities and units for ionizing radiation. ICRU Report 85. ICRU Publications ( 2011 ). 12 Camp , A. , Vargas , A. and Fernández-Varea , J. M. Determination of LaBr3(Ce) internal background using a HPGe detector and Monte Carlo simulations . Appl. Radiat. Isot. 109 , 512 – 517 ( 2016 ). Google Scholar CrossRef Search ADS PubMed 13 Box , G. E. P. and Müller , M. E. A Note on the generation of random normal deviates . Ann. Math. Stat. 29 ( 2 ), 610 – 611 ( 1958 ). Google Scholar CrossRef Search ADS 14 Casanovas , R. , Morant , J. J. and Salvadó , M. Energy and resolution calibration of NaI(Tl) and LaBr3(Ce) scintillators and validation of an EGS5 Monte Carlo user code for efficiency calculations . Nucl. Instrum. Methods Phys. Res. A 675 , 78 – 83 ( 2012 ). Google Scholar CrossRef Search ADS 15 Camp , A. and Vargas , A. Ambient dose estimation H*(10) from LaBr3(Ce) spectra . Radiat. Prot. Dosim. 160 ( 4 ), 264 – 268 ( 2014 ). Google Scholar CrossRef Search ADS 16 Hendricks , J. S. , McKinney , G. W. , Fensin , M. L. , James , M. R. , Johns , R. C. , Durkee , J. W. , Finch , J. P. , Pelowitz , D. B. , Waters , L. S. and Johnson , M. W. ‘MCNPX 2.6.0 Extensions’, Los Alamos National Laboratory. LA-UR-08-2216 ( 2008 ). © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Radiation Protection Dosimetry Oxford University Press

MAXIMUM LIKELIHOOD ESTIMATION USING EXPECTATION MAXIMIZATION APPLIED TO AMBIENT DOSE EQUIVALENT MEASUREMENTS

Loading next page...
 
/lp/ou_press/maximum-likelihood-estimation-using-expectation-maximization-applied-rJVZUujsQ0
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com
ISSN
0144-8420
eISSN
1742-3406
D.O.I.
10.1093/rpd/ncy064
Publisher site
See Article on Publisher Site

Abstract

Abstract The use of spectrometric detectors for ambient dose equivalent measurements requires the development of calculation procedures to obtain the incident photon fluence spectra from the measured pulse-height spectra. In this work, the implementation of a methodology based on the maximum likelihood estimation using expectation maximization is described. The method has been validated using simulated and reference photon beams. The results are in excellent agreement with the reference values in the framework of the involved uncertainties. The presented method allows the consideration of the detector resolution, improving thus the energy resolution of the obtained incident photon spectra. INTRODUCTION The appropriate quantity to estimate the effective dose in area monitoring is the ambient dose equivalent, H*(10). This operational quantity was introduced by the International Commission on Radiation Units and Measurements(1) and its use is mandatory in the European Union according to EU Directive 2013/59 EURATOM(2). There are thousands of radiological early warning network stations in Europe in which typically Geiger Müller detectors are being used for ambient dose equivalent rate measurements. Moreover, some detector systems designed for other quantities, such as exposure or air kerma are also in use. Due basically to their inability to provide information about the incident photon spectra and to the pronounced energy and angular response to gamma radiation of these detectors in terms of H*(10), new compact and handy spectrometric systems based on scintillation detectors like LaBr3, CeBr3 and SrI2 or on room-temperature semiconductor detectors like CdZnTe are installed instead of or in addition to the traditional dose rate meters. These spectrometers will considerably improve the information provided by early warning networks in real-time, in particular concerning radionuclide identification(3). The use of spectrometric detectors requires the development of calculation procedures to obtain the H*(10) rates from the measured photon pulse-height spectra. A comprehensive description of the so far applied methods, including their advantages and disadvantages, has been provided by Dombrowski(4). A different approach, based on the maximum likelihood estimation using expectation maximization (ML-EM), has been recently implemented by the author in the framework of the project: ‘Metrology for Radiological Early Warning Networks in Europe’ (MetroERM project(5)), as part of the tasks related to the investigations of detector features and of spectra evaluation and deconvolution methods for new and improved measurement systems based on novel spectrometric detectors. The ML-EM is a statistical method initially developed for the use in the reconstruction of images in positron emission tomography(6, 7) and successfully tested for spectra deconvolution in gamma-ray spectrometry(8). The implementation of the ML-EM method for H*(10) rate calculations from measured gamma-ray spectra is explained in the present work. The associated mathematical formalism and the general features of the specific computer code DET-H10R, developed for spectra deconvolution and H*(10) rate calculations using the ML-EM algorithm, are presented. The incorporation of the detector resolution during the unfolding of measured photon spectra using the code DET-H10R is described and the results of theoretical and experimental verifications of the implemented method are also presented. MATERIALS AND METHODS General concepts After background subtraction, the measured pulse-height spectrum O(E) can be represented by the integral: O(E)=∫0∞R(E,E0)⋅ΦE(E0)⋅dE0 (1) where R(E, E0) is the detector response function and ΦE(E0) is the distribution of the incident photon fluence in energy. The discrete form of this equation is used in practical applications: Oi=Ri,j⋅Φj{i=1,…,n;j=1,…,m} (2) where Oi is the number of counts in channel i of the measured pulse-height spectrum, Ri,j is the number of counts expected in channel i of the measured pulse-height spectrum per unit of photon fluence with energy within the energy bin j, and Φj is the unknown fluence of photons with energy in bin j. The number of channels in the measured spectrum is denoted as ‘n’, the number of energy bins of the photon fluence spectrum as ‘m’. The summation convention has been applied, i.e. the repeated indices in the right term are implicitly summed over. The matrix Rˆ(n,m) in equation (2) represents the detector response operator and is also known as ‘detector response matrix’. The ‘noise’ contribution due, for instance, to the counting statistics, the lack of accuracy in the determination of the detector response operator and its finite resolution, has been omitted in equation (2) and is therefore neglected in the analysis. The unknown fluence vector Φ→ in equation (2) generates, and is to be calculated from, the observed pulse-height spectrum O→ by using a known detector response matrix. The ‘noise’ referred in the previous paragraph contributes to the uncertainty of these calculations. Once the photon fluence spectrum is determined, the air kerma, Kair, can be calculated using the following equation: Kair=∑j=1mΦj⋅Ej⋅(μtrρ)air,j (3) where Φj is the fluence of photons with energy in bin j, Ej is the mean energy of the energy bin j and (μtrρ )air,j is the mass energy transfer coefficient in air for photons with energy Ej. The radiation protection quantities can be also calculated from the photon fluence spectrum. In particular, the ambient dose equivalent can be obtained as follows: H*(10)=∑j=1mΦj⋅Ej⋅(μtrρ)air,j⋅hk*(10,Ej) (4) where hk*(10,Ej) is the air kerma to ambient dose equivalent conversion coefficient corresponding to the energy Ej. The rest of the parameters are the same as in equation (3). One of the requirements of the ML-EM method is that the response matrix Rˆ(n,m) in equation (2) is needed for a given detector. It could be obtained by Monte Carlo simulation after a detailed characterization of the detection system, i.e. the geometry and material composition of the elements that constitute the detection probe are to be well known. Then, the response matrix is filled with m detector specific pulse-height spectra of n channels, where each column j contains the pulse-height spectrum produced by incident photons with an energy within the energy bin j. In what follows, some simple transformations of equation (2) are introduced in order to apply the ML-EM method(6, 7) to this particular case, in which the incident photon fluence spectrum and the corresponding ambient dose equivalent are calculated from the measured pulse-height spectrum. In order to simplify the formulations, the summation convention is applied for repeating indices ‘i’ or ‘s’ in matrix products. Basic transformations required An important requirement for the ML-EM method is that the sum of the elements in each column of the response operator should be lower than or equal to one. In order to achieve this, a total efficiency operator RT^(m,m){RTj,l=0∀(j≠l);RTj,j=RTj} is derived from the detector response matrix Rˆ(n,m), and the vector of incident photons which are detected, I→, is defined according to the following equation: Ij=RTj⋅Φj{j=1,…,m} (5) where RTj is the total number of counts expected in the measured spectrum per unit photon fluence in energy bin j and Φj is the fluence of photons with energy in bin j. The RTj values are calculated from the detector response matrix by summing the elements of each column: RTj=∑i=1nRi,j{j=1,…,m} (6) Combining equations (2) and (5), the following relationship is obtained: Oi=Ci,s⋅Is{i=1,…,n;s=1,…,m} (7) where Ci,s represents the fraction of detected photons, incident with energy s, that are detected in channel i: Ci,j=Ri,jRTj{i=1,…,n;j=1,…,m} (8) As can be noted from equations (6) and (8), the total number of incident photons that are detected is equal to the number of counts in the measured spectrum, i.e.: ∑i=1nCi,j=1 (9) This is true when coincidence summing effects can be neglected, which is the case when in-situ measurements are performed. This is also true for low ambient dose equivalent rates for which the dead time corrections of the detection systems are not significant. The ML-EM method will be applied to solve equation (7), i.e. to calculate the vector I→ from the measured pulse-height spectrum O→ by using the normalized detector response matrix, Cˆ(n,m), whose elements are defined by equation (8). The calculation of the incident fluence spectrum is then straightforward from equation (5). The ML-EM algorithm For a given detection system and for each possible vector I→, the measured spectrum O→ has a probability or likelihood: L(I→)=P(O→|I→) (10) The idea behind the implemented method is to find the mathematical expectation of the vector Í→ which maximizes the probability of measuring the actual pulse-height spectrum O→ over all possible incident photon spectra, i.e. max(L(I→))=L(Í→). It is obvious that for a given detector response matrix Rˆ(n,m) this solution also provides the maximization of the likelihood L(Φ→)=P(O→|Φ→) because: Φ′j=ÍjRTj{∀j∈mˆ} (11) The likelihood function deduced by Shepp and Vardi(6) has been chosen for L(I→), based on the assumption that the components of the involved photon spectra, representing a certain number of counts in a given time interval, can be considered independent Poisson variables. The iterative scheme based on the expectation maximization (EM) algorithm(6, 7) is thus applied to find the maximum value of this likelihood function in equation (10), i.e. to find the vector Í→ for which the probability of having measured the spectrum O→ reaches the maximum value. Mathematical demonstrations of the method are discussed in detail in the pioneering work of Shepp and Vardi(6) and are beyond the scope of this work, in which the efforts have been aimed at the achievement of a compliance with the necessary requirements for its application. The equality according to equation (9) by using the normalized response operator defined by equations (7) and (8) makes the following iteration procedure simpler and more efficient. Starting with an initial guess of Í→: Íjr=0≥0{j=1,…,m} for which: ∑j=1mÍjr=0=∑i=1nOi (12) a new estimate of the spectrum Í→ is calculated at any step r as follows: Íjr=Íjr−1⋅[OiCi,s⋅Ísr−1⋅Ci,j]{i=1,…,n;s,j=1,…,m;n≥m}, (13) where, the likelihood function increases: P(O→|Ír→)≥P(O→|Ír−1→){r≥1} with equality if and only if L(Ír−1→)=max(L(I→)). As can be deduced from equation (12), the measured spectrum O→ could be taken as an initial estimate of the incident spectrum when n = m. The incident photon fluence spectrum is obtained at each iteration step using equation (11) and the corresponding value of the ambient dose equivalent is calculated according to equation (4), i.e. as follows: H*(10)r=∑j=1mΦjr⋅Ej⋅(μtrρ)air,j⋅hk*(10,Ej) (14) where H*(10)r and Φjr are respectively, the ambient dose equivalent and the fluence of photons with energy in bin j, obtained at iteration step r. The rest of the parameters have the same meaning as in equations (3) and (4). The iterative process stops when Δ=|H*(10)r−H*(10)r−1H*(10)r−1|≤δ (15) where δ is the accepted relative deviation for the two consecutive estimated values of H*(10). The ambient dose equivalent rate can be calculated using the live collection time of the spectrometric system. Implementation of the ML-EM method Although the ML-EM method is more computer-power-demanding than other simpler procedures(4), the involved algorithm, i.e. the iterative calculation given by equation (13), is computationally of low complexity and could be implemented on a low-power, industrial CPU, of a type which can be integrated in a hand-held instrument or measurement probe. It should be also considered that field instruments, for instance in early warning stations, could easily send the measured pulse-height spectra (few tens of kB for 1024 channels, in the most used formats) in real-time to an external central processing system. Thus, these spectra could be stored in a central memory unit, from which they could be read by the main computer with a predefined frequency for dose calculations, radionuclide identification, alarms triggering, etc. This could be the better option, even for apparently simpler or less-demanding H*(10) calculation procedures, taking into account that the CPU integrated in the measuring devises should be used rather to produce reliable pulse-height spectra, corrected for instrument instabilities, in particular for energy shifts in the measured peaks due to thermal instability of the electronic components. General features of the code DET-H10R The presented ML-EM method has been implemented in the code DET-H10R, written in free PASCAL for the calculation of air kerma rates and H*(10) rates from measured photon pulse-height spectra in PC platforms. The calculated spectra of the incident photon fluence rate, air kerma rate and H*(10) rate are also provided by the code and can be saved in text files. The values of the mass energy transfer coefficient in air have been obtained from the publication of Hubbell and Seltzer(9) and the values of the air kerma to ambient dose equivalent conversion coefficient hk*(10,E) have been taken from the ICRU Report 57(10). To facilitate the calculations, the discrete values of the provided coefficients were fitted with empirical analytical functions. The maximum relative deviation between theoretical and calculated values is below 0.3% with 1 − P ≤ 10−5, where, P is the Pearson r-moment. The supported file formats for the measured spectra are the ISO N42.42 and the binary CAM CNF format from CANBERRA. A code-specific text format has been also defined for users to upload those measured pulse-height spectra with formats not supported by the code. Considering that the resolution, i.e. the number of keV per channel, of measured pulse-height spectra should be equal to or better than the resolution of the constructed detector response matrix, a straightforward procedure has been incorporated into the code DET-H10R to fit the energy bin width of measured spectra with the energy bin width of the corresponding detector response matrix. Starting from the lower energy, this procedure gradually groups the energy channels of measured spectra according to the resolution of the response matrix. For instance, if a measured spectrum has energy channels of 3 keV width and the energy bin width of the response matrix is 10 keV, the first energy bin of the adjusted spectrum will contain the sum of the counts from the first three channels plus one-third of the counts from the fourth channel of the measured spectrum. The second energy bin of the adjusted spectrum will contain two-thirds of the counts from the fourth channel plus the counts from the fifth and the sixth channels, plus two-thirds of the counts from the seventh channel, and so on. A graphical user interface, including a user help, facilitates the selection of the detector response matrix, the uploading of the measured pulse-height spectra and the definition of the relevant calculation parameters, like the ‘iteration resolution’ and the ‘cut-off energy’. The iteration resolution is the accepted relative deviation between two successively calculated values of the H*(10) rate, as defined by equation (15). The cut-off energy is the lower energy considered during the calculations of air kerma and H*(10) rates. This option could be useful to discard the contribution of low-energy incident photons in those cases where the energy bin width (resolution) of the detector response matrix could be insufficient to properly consider the interaction processes of these photons or to adequately reflect the pronounced energy dependence of the involved photon interaction coefficients(11) and the air kerma to ambient dose equivalent conversion coefficient, in the low-energy region. In these cases, the response matrix should be constructed with the appropriate resolution, after a more detailed characterization of the detector geometry and materials. The code interface also incorporates tools to analyse the measured pulse-height spectra and the calculated H*(10) rate spectra (see the example in Figures 1 and 2) facilitating the verification of the correct energy calibration of spectrometers, the detection of possible energy shifts due to thermal instability of the electronic components and the radionuclide identification. Figure 1. View largeDownload slide Measured gross pulse-height spectrum using a LaBr3 (1″ × 1″) detector in a low dose rate 60Co collimated beam. Figure 1. View largeDownload slide Measured gross pulse-height spectrum using a LaBr3 (1″ × 1″) detector in a low dose rate 60Co collimated beam. Figure 2. View largeDownload slide Calculated H*(10) rate spectrum from the 60Co pulse-height spectrum shown in Figure 1. The detector response matrix with 10 keV energy bin width has been obtained by Monte Carlo simulation. The small fraction of scattered photons produced by the collimator of the 60Co source can be seen as well as a small broad peak around 1470 kV, corresponding to the external 40K background and the internal detector contamination with 138La(12). Figure 2. View largeDownload slide Calculated H*(10) rate spectrum from the 60Co pulse-height spectrum shown in Figure 1. The detector response matrix with 10 keV energy bin width has been obtained by Monte Carlo simulation. The small fraction of scattered photons produced by the collimator of the 60Co source can be seen as well as a small broad peak around 1470 kV, corresponding to the external 40K background and the internal detector contamination with 138La(12). Incorporation of the detector resolution One advantage of the ML-EM method is that a triangular matrix is not required in equation (2). Therefore, the detector resolution could be taken into account during the unfolding of measured pulse-height spectra. A convolution algorithm included in DET-H10R facilitates the incorporation of the detector energy resolution into the detector response matrix previously constructed by Monte Carlo simulation. The resulting response operator allows the calculation of more realistic incident photon fluence spectra, with improved resolution. A brief description of the implemented procedure is presented below. During the convolution procedure the detector response function in equation (1) is modified accordingly to the following equation: R′(E,E0)=∫0∞G(E−E′)⋅R(E′,E0)⋅d(E′) (16) where, R′(E,E0) is the response function for incident energy E0 considering the resolution of the detector system, G(E−E′) represents the probability per unit pulse height that a deposited energy E′ will produce counts on the energy channel E, and R(E′,E0) is the ideal or mathematical detector response function for incident energy E0, obtained by Monte Carlo simulation. The resolution operator G(E−E′) is approximated by a Gaussian distribution: G(E−E′)=12πσ2e−((E−E′)2/2σ2){σ≡σ(E′)} (17) where the standard deviation σ(E′) is calculated from the peak full width at half maximum (FWHM) corresponding to energy E′, according to the following relationship: σ(E′)=FWHM(E′)22ln2 (18) As in the case of equation (1), the discrete form of equation (16) is used: R′i,j=Gi,s⋅Rs,j{i,s=1,…,n;j=1,…,m} (19) Each column of the spreading operator Gˆ(n,n) is a Gaussian probability density function in form of histogram, with the most probable values in the diagonal elements. These columns are constructed by generating energy values, E, using the Box–Müller algorithm(13): E=E′+σ(E′)⋅sin(2π⋅γ1)⋅−2ln⋅γ2 (20) where E′ is an energy value sampled uniformly in the energy bin corresponding to the diagonal element in a given column. γ1 and γ2 are the independent random variables uniformly distributed in the interval (0,1). One event is counted in the corresponding energy bin per each energy value E generated. To obtain the probability density function, the counts in each energy bin are then divided by the total energy values generated for a given column. It can be noted that this normalization of the spreading operator implies that the counting statistics in the new response matrix is not altered. The number of values with an energy, E, to be generated in each column can be selected by the user according to the required statistical uncertainty. A graphical interface has been included to facilitate the input of the detector specific resolution function: FWHM≡FWHM(E′) determined from the energy resolution calibration. The following two general equations are provided to facilitate the definition of a wide variety of possible FWHM functions(14) by selecting the adequate fitting parameters. FWHM(E′)=a⋅b+c⋅E′+d⋅E′2 (21) FWHM(E′)=a+b⋅E′+c⋅E′2+d⋅E′ (22) The convolution process indicated by equation (19) readjusts the distribution of the deposited energy in each pulse-height spectrum that constitutes the detector response matrix, producing the broadening of this pulse-height spectrum. RESULTS AND DISCUSSION Calculations with a simulated spectrum The ML-EM method has been applied to a relatively complex photon pulse-height spectrum generated by Monte Carlo simulation, using the model of a LaBr3 (5% Ce) (1.0″ × 1.0″) detector. A parallel photon beam, perpendicular to the main detector axis and constituted by 10 gamma-energy lines was simulated. The beam covered the entire probe volume which was exposed to a total ambient dose equivalent rate of 500 nSv/h. The photon fluence rate values were selected so that the ambient dose equivalent rate corresponding to each energy line was 50 nSv/h (Table 1). The number of simulated photons results from the adopted acquisition time of the spectrometer and the area of the beam cross section. The simulated pulse-height spectrum (Figure 3) does not consider the statistics associated with the process between the deposition of the photon’s energy and the counting of the corresponding event in the spectrum. Therefore, in order to obtain a more realistic spectrum, the simulated spectrum was convolved with the resolution operator constructed from the specific detector resolution function shown in Figure 4. The resulting spectrum, to which the ML-EM method was applied, is shown in Figure 5. Table 1. Simulated gamma lines, fluence rates and calculated H*(10) rates. Energy Fluence rate Ḣ*(10)ML−EM keV cm−2 s−1 nSv/h 59,54 (241Am) 27.23 50.5 ± 1.6 134 (144Ce) 17.43 49.5 ± 1.1 364 (131I) 6.54 49.7 ± 1.1 609 (214Po) 3.96 49.9 ± 1.1 662 (137Cs) 3.68 48.9 ± 1.2 796 (134Cs) 3.16 49.2 ± 1.2 1173 (60Co) 2.34 49.2 ± 1.6 1332 (60Co) 2.13 49.0 ± 1.6 1460 (40K) 2.00 51.1 ± 2.1 1764 (214Po) 1.75 49.5 ± 2.0 Total — 496,5 ± 4,8 Energy Fluence rate Ḣ*(10)ML−EM keV cm−2 s−1 nSv/h 59,54 (241Am) 27.23 50.5 ± 1.6 134 (144Ce) 17.43 49.5 ± 1.1 364 (131I) 6.54 49.7 ± 1.1 609 (214Po) 3.96 49.9 ± 1.1 662 (137Cs) 3.68 48.9 ± 1.2 796 (134Cs) 3.16 49.2 ± 1.2 1173 (60Co) 2.34 49.2 ± 1.6 1332 (60Co) 2.13 49.0 ± 1.6 1460 (40K) 2.00 51.1 ± 2.1 1764 (214Po) 1.75 49.5 ± 2.0 Total — 496,5 ± 4,8 Table 1. Simulated gamma lines, fluence rates and calculated H*(10) rates. Energy Fluence rate Ḣ*(10)ML−EM keV cm−2 s−1 nSv/h 59,54 (241Am) 27.23 50.5 ± 1.6 134 (144Ce) 17.43 49.5 ± 1.1 364 (131I) 6.54 49.7 ± 1.1 609 (214Po) 3.96 49.9 ± 1.1 662 (137Cs) 3.68 48.9 ± 1.2 796 (134Cs) 3.16 49.2 ± 1.2 1173 (60Co) 2.34 49.2 ± 1.6 1332 (60Co) 2.13 49.0 ± 1.6 1460 (40K) 2.00 51.1 ± 2.1 1764 (214Po) 1.75 49.5 ± 2.0 Total — 496,5 ± 4,8 Energy Fluence rate Ḣ*(10)ML−EM keV cm−2 s−1 nSv/h 59,54 (241Am) 27.23 50.5 ± 1.6 134 (144Ce) 17.43 49.5 ± 1.1 364 (131I) 6.54 49.7 ± 1.1 609 (214Po) 3.96 49.9 ± 1.1 662 (137Cs) 3.68 48.9 ± 1.2 796 (134Cs) 3.16 49.2 ± 1.2 1173 (60Co) 2.34 49.2 ± 1.6 1332 (60Co) 2.13 49.0 ± 1.6 1460 (40K) 2.00 51.1 ± 2.1 1764 (214Po) 1.75 49.5 ± 2.0 Total — 496,5 ± 4,8 Figure 3. View largeDownload slide Simulated spectrum in a LaBr3 (1.0″ × 1.0″) detector: 1024 channels, acquisition time = 600 s, total H*(10) rate = 500 nSv/h. Figure 3. View largeDownload slide Simulated spectrum in a LaBr3 (1.0″ × 1.0″) detector: 1024 channels, acquisition time = 600 s, total H*(10) rate = 500 nSv/h. Figure 4. View largeDownload slide Resolution function applied in the simulations of the test spectrum with the LaBr3 (1.0″ × 1.0″) detector. Figure 4. View largeDownload slide Resolution function applied in the simulations of the test spectrum with the LaBr3 (1.0″ × 1.0″) detector. Figure 5. View largeDownload slide Test spectrum obtained by convolving the simulated spectrum (Figure 3) with the resolution operator corresponding to experimental data of the studied LaBr3 (1.0″ × 1.0″) detector. Number of channels = 1024, acquisition time = 600 s, total H*(10) rate = 500 nSv/h. Figure 5. View largeDownload slide Test spectrum obtained by convolving the simulated spectrum (Figure 3) with the resolution operator corresponding to experimental data of the studied LaBr3 (1.0″ × 1.0″) detector. Number of channels = 1024, acquisition time = 600 s, total H*(10) rate = 500 nSv/h. The detector response matrix was obtained by Monte Carlo simulation using an energy bin width of 10 keV and assuming the same geometry, i.e. perpendicular beam incidence to the longitudinal detector axis. According to the work of Camp and Vargas(15), a bin width of 10 keV constitutes a good trade-off between energy resolution and counting statistics during ambient dose equivalent measurements with spectrometric systems based on scintillation detectors. The energy span considered was 3000 keV, therefore 300 spectra were generated. For each spectrum 109 photon histories were simulated in order to achieve a statistical standard uncertainty below 1% in each energy bin. The cut-off energies were 1 keV for photons and 10 keV for electrons. No additional variance reduction techniques have been applied in order to minimize the possibility of biased results, therefore 64 processors (cores) have been used (3 GHz, 2 GB ram/core). The generation of the test spectrum and the detector response matrix were carried out with the code MCNPX(16), using the photon library mcplib04 and the electron library el03. The incident ambient dose equivalent rate spectrum was calculated by unfolding the test spectrum (Figure 5) with the ML-EM algorithm. As expected, the resolution of the calculated incident spectrum using the ideal detector response matrix, i.e. without including the detector energy resolution, was similar to the resolution of the pulse-height spectrum, as can be seen in Figure 6. In this case, the peaks related to the 609 and 662 keV gamma-energy lines are not well resolved and consequently the ambient dose equivalent rate corresponding to each of these energies cannot be determined accurately. Figure 6. View largeDownload slide H*(10) rate spectrum calculated by unfolding the test spectrum (Figure 5) with the ideal detector response matrix, i.e. without considering the energy resolution of the detector. Figure 6. View largeDownload slide H*(10) rate spectrum calculated by unfolding the test spectrum (Figure 5) with the ideal detector response matrix, i.e. without considering the energy resolution of the detector. Nevertheless, the resolution of the recalculated incident spectrum was considerably improved when the energy resolution of the detector (Figure 4) was incorporated into its response matrix using the previously described convolution algorithm. The obtained spectrum is shown in Figure 7. In this case, the ambient dose equivalent rate corresponding to each incident photon-energy line can be easily obtained as the area of the related peak in the spectrum. It should be indicated that the energy resolution in the obtained spectra could be further improved by using detector response matrixes with lower energy bin widths, i.e. with better resolution, in detriment, however, of the counting statistics. Figure 7. View largeDownload slide H*(10) rate spectrum calculated by unfolding the test spectrum (Figure 5) with the spread detector response matrix, i.e. taking into account the energy resolution of the detector. Figure 7. View largeDownload slide H*(10) rate spectrum calculated by unfolding the test spectrum (Figure 5) with the spread detector response matrix, i.e. taking into account the energy resolution of the detector. The calculated ambient dose equivalent rates ( Ḣ*(10)ML−EM) are shown in Table 1. The indicated uncertainties have been expressed as expanded uncertainties with a coverage factor k = 2 and have been estimated considering the contributions of the counting statistics in the simulated 1024 channel spectrum and the finite resolution (10 keV per channel) of the detector response matrix. The differences of the calculated ambient dose equivalent rates and the expected values, i.e. 50 nSv/h for each photon energy and 500 nSv/h for the total spectrum, are consistent with the involved uncertainties. Measurements in reference photon beams As part of the validation studies, the implemented method has been applied to experimental spectra obtained at the Metrology Laboratory for Ionizing Radiation (LMRI) of CIEMAT. Two probes with spectrometric detectors (LaBr3 (1.5″ × 1.5″) and CeBr3 (1.5″ × 1.5″)) were characterized and their response matrixes were calculated by Monte Carlo simulation assuming a perpendicular beam incidence to the longitudinal detector axis and applying the same procedures and parameters as in the previous calculations. As a result, for each detector system a 300 × 300 response matrix with a statistical standard deviation better than 1% in each energy bin was obtained. Both detectors were exposed to known values of H*(10) rates in 137Cs and 60Co reference calibration beams at LMRI. The longitudinal axes of the detector probes were placed perpendicular to the axes of the photon beams. At each calibration point, five gross spectra and five background spectra were collected during 120 s each. The ML-EM method was applied to each measured spectrum to calculate the corresponding H*(10) rates. The calculations were performed in fractions of a second using a personal computer with an i7-4790 CPU @ 3.60 GHz processor for an iteration resolution better than 10−3%. The net measured H*(10) rate was obtained by subtracting the mean value of the background H*(10) rates from the mean value of the gross H*(10) rates. The results are shown in Table 2. The laboratory reference values Ḣ*(10)REF have a relative standard uncertainty of 2.2% and the relative statistical standard uncertainty of the measured H*(10) rate values Ḣ*(10)ML−EM is lower than 1%. The results provided by the implemented ML-EM algorithm are in excellent agreement with the reference values, considering the involved uncertainties. Table 2. Reference and calculated H*(10) rates at the LMRI radiation beams. Detector Source Ḣ*(10)REF Ḣ*(10)ML−EM Relative difference % μSv h−1 μSv h−1 LaBr3 137Cs 1.00 1.01 1.0 2.00 1.99 −0.5 5.00 4.99 −0.2 10.0 9.99 −0.1 60Co 2.00 2.03 1.5 CeBr3 137Cs 0.842 0.847 0.6 2.00 2.01 0.5 5.00 4.92 −1.6 10.0 9.93 −0.7 60Co 5.00 4.99 −0.2 Detector Source Ḣ*(10)REF Ḣ*(10)ML−EM Relative difference % μSv h−1 μSv h−1 LaBr3 137Cs 1.00 1.01 1.0 2.00 1.99 −0.5 5.00 4.99 −0.2 10.0 9.99 −0.1 60Co 2.00 2.03 1.5 CeBr3 137Cs 0.842 0.847 0.6 2.00 2.01 0.5 5.00 4.92 −1.6 10.0 9.93 −0.7 60Co 5.00 4.99 −0.2 Table 2. Reference and calculated H*(10) rates at the LMRI radiation beams. Detector Source Ḣ*(10)REF Ḣ*(10)ML−EM Relative difference % μSv h−1 μSv h−1 LaBr3 137Cs 1.00 1.01 1.0 2.00 1.99 −0.5 5.00 4.99 −0.2 10.0 9.99 −0.1 60Co 2.00 2.03 1.5 CeBr3 137Cs 0.842 0.847 0.6 2.00 2.01 0.5 5.00 4.92 −1.6 10.0 9.93 −0.7 60Co 5.00 4.99 −0.2 Detector Source Ḣ*(10)REF Ḣ*(10)ML−EM Relative difference % μSv h−1 μSv h−1 LaBr3 137Cs 1.00 1.01 1.0 2.00 1.99 −0.5 5.00 4.99 −0.2 10.0 9.99 −0.1 60Co 2.00 2.03 1.5 CeBr3 137Cs 0.842 0.847 0.6 2.00 2.01 0.5 5.00 4.92 −1.6 10.0 9.93 −0.7 60Co 5.00 4.99 −0.2 In order to minimize the influence of the detector angular response in the evaluation of the implemented method the irradiations were performed for a similar geometry to the one used during the calculations of the response matrixes, i.e. for perpendicular beam incidence. According to Camp and Vargas(15) a nearly isotropic response of this kind of instruments, i.e. cylindrical scintillation detectors with the same diameter and length, could be assumed in environmental monitoring. In general, depending also on the probe materials and their thicknesses, the angular dependence of the response is not negligible at lower energies. For a given detector and energy interval, the evaluation of the uncertainty component due to the angular dependence of the response could be carried out by Monte Carlo simulations. An average response matrix, weighted with the relevant angles of incidence, could then be used for H*(10) rate calculations. CONCLUSIONS In this paper the ML-EM method has been implemented for ambient dose equivalent measurements with spectrometric detectors. The results obtained during the verifications carried out with a simulated pulse-height spectrum and in reference photon beams are in excellent agreement with the expected values, considering the involved uncertainties. The incident gamma-ray fluence spectrum can be obtained in fractions of a second by unfolding the measured spectrum and the ambient dose equivalent can be accurately calculated using the relevant energy dependent coefficients, i.e. the mass energy transfer coefficient in air and the air kerma to ambient dose equivalent coefficient. A flat-energy response with respect to the quantity H*(10) is thus expected when using well-characterized spectrometric detectors. The application of the ML-EM method, like that of any other H*(10) calculation procedure using measured pulse-height spectra, requires the knowledge of the detector response matrix, as described above. Unfortunately, detector response matrixes are not provided by the manufacturers of spectrometric systems, and therefore they must be created by the users. The detector response matrix may be determined very accurately by Monte Carlo simulations, considering that a detailed characterization of the detection probe geometry and materials is possible. Although these simulations could be time demanding, they need to be carried out only once for a particular detection system. Well-characterized spectrometric detectors, with an isotropic response, can be used for ambient dose equivalent measurements without any knowledge of the incident photon field. Furthermore, the possibility of incorporating the detector resolution into the response matrix allows the improvement of the energy resolution in the calculated photon spectra, facilitating the identification of the involved radionuclides and the evaluation of their contributions to the measured ambient dose equivalent rate. Due to the diversity of currently available methods for ambient dose rate calculations from measured pulse-height spectra, the analysis of the advantages and disadvantages of the implemented ML-EM method in comparison with each of these existing procedures is beyond the scope of this work. The readers have been provided with the information required to evaluate the adequacy of the described approach for their particular applications or in their particular cases. Moreover, several comparison exercises that have been carried out in the framework of the MetroERM project among the implemented H*(10) calculation methods will be subject of further publications. ACKNOWLEDGEMENTS The author would like to thank CIEMAT for providing the working conditions required to work out this article. The author is also grateful to Dr José Antonio Suárez Navarro for his comments and suggestions. FUNDING This work was supported by EMRP within the JRP-Env57 (MetroERM) project. The EMRP is jointly funded by the EMRP participating countries within EURAMET and the European Union. REFERENCES 1 International Commission on Radiation Units and Measurements . Determination of dose equivalents resulting from external radiation sources. ICRU Report 39. ICRU Publications ( 1985 ). 2 European Council Directive 2013/59/EURATOM on basic safety standards for protection against the dangers arising from exposure to ionising radiation and repealing. Directives 89/618/Euratom, 90/641/Euratom, 96/29/Euratom, 97/43/Euratom and 2003/122/Euratom. OJ of the EU. L13; 57: 1–73 (2014). 3 Neumaier , S. , Dombrowski , H. and Kessler , P. Metrology for radiological early warning networks in Europe (‘METROERM’)—a join European metrology research project . Health Phys. 111 ( 2 ), 100 – 105 ( 2016 ). Google Scholar CrossRef Search ADS PubMed 4 Dombrowski , H. Area dose rate values derived from NaI or LaBr3 spectra . Radiat. Prot. Dosim. 160 ( 4 ), 269 – 276 ( 2014 ). Google Scholar CrossRef Search ADS 5 EUROPEAN ASSOCIATION OF NATIONAL METROLOGY INSTITUTES . Final Publishable Join Research Project Summary Report for ENV57 MetroERM—Metrology for Radiological Early Warning Networks in Europe. Available on https://www.euramet.org/research-innovation/search-research-projects. Issued in September ( 2017 ). 6 Shepp , L. A. and Vardi , Y. Maximum likelihood reconstruction for emission tomography . IEEE Trans. Med. Imaging MI-1 ( 2 ), 113 – 122 ( 1982 ). Google Scholar CrossRef Search ADS 7 Shepp , L. A. , Vardi , Y. and Kaufman , L. A statistical model for positron emission tomography . J. Am. Stat. Assoc. 80 ( 389 ), 8 – 20 ( 1985 ). Google Scholar CrossRef Search ADS 8 Meng , L. J. and Ramsden , D. An inter-comparison of three spectral-deconvolution algorithms for gamma-ray spectroscopy . IEEE Trans. Nucl. Sci. 47 ( 4 ), 1329 – 1336 ( 2000 ). Google Scholar CrossRef Search ADS 9 Hubbell , J. H. and Seltzer , S. M. Tables of X-Ray Mass Attenuation Coefficients and Mass Energy-Absorption Coefficients (version 1.4) ( Gaithersburg, MD : National Institute of Standards and Technology ) ( 2004 ) [Online] Available: http://physics.nist.gov/xaamdi. 10 International Commission on Radiation Units and Measurements . Conversion coefficients for use in radiological protection against external radiation. Bethesda, MN: ICRU Report 57 ( 1998 ). 11 International Commission on Radiation Units and Measurements . Fundamental quantities and units for ionizing radiation. ICRU Report 85. ICRU Publications ( 2011 ). 12 Camp , A. , Vargas , A. and Fernández-Varea , J. M. Determination of LaBr3(Ce) internal background using a HPGe detector and Monte Carlo simulations . Appl. Radiat. Isot. 109 , 512 – 517 ( 2016 ). Google Scholar CrossRef Search ADS PubMed 13 Box , G. E. P. and Müller , M. E. A Note on the generation of random normal deviates . Ann. Math. Stat. 29 ( 2 ), 610 – 611 ( 1958 ). Google Scholar CrossRef Search ADS 14 Casanovas , R. , Morant , J. J. and Salvadó , M. Energy and resolution calibration of NaI(Tl) and LaBr3(Ce) scintillators and validation of an EGS5 Monte Carlo user code for efficiency calculations . Nucl. Instrum. Methods Phys. Res. A 675 , 78 – 83 ( 2012 ). Google Scholar CrossRef Search ADS 15 Camp , A. and Vargas , A. Ambient dose estimation H*(10) from LaBr3(Ce) spectra . Radiat. Prot. Dosim. 160 ( 4 ), 264 – 268 ( 2014 ). Google Scholar CrossRef Search ADS 16 Hendricks , J. S. , McKinney , G. W. , Fensin , M. L. , James , M. R. , Johns , R. C. , Durkee , J. W. , Finch , J. P. , Pelowitz , D. B. , Waters , L. S. and Johnson , M. W. ‘MCNPX 2.6.0 Extensions’, Los Alamos National Laboratory. LA-UR-08-2216 ( 2008 ). © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Radiation Protection DosimetryOxford University Press

Published: May 4, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off