
==== Front
Z Med Phys
Z Med Phys
Zeitschrift für Medizinische Physik
0939-3889
1876-4436
Elsevier

S0939-3889(22)00133-7
10.1016/j.zemedi.2022.11.012
Short Communication
Note on uncertainty in Monte Carlo dose calculations and its relation to microdosimetry
Hartmann Günther H. g.hartmann@dkfz.de
a⁎
Menzel Hans G. b
a German Cancer Research Center (DKFZ), Heidelberg, Germany
b International Commission on Radiation Units and Measurements (ICRU), Germany
⁎ Corresponding author: Günther H. Hartmann, German Cancer Research Center (DKFZ), Heidelberg, Germany. g.hartmann@dkfz.de
26 12 2022
8 2024
26 12 2022
34 3 468476
18 8 2022
29 11 2022
© 2022 Published by Elsevier GmbH on behalf of DGMP, ÖGMP and SSRMP.
2022

https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Purpose

The Type A standard uncertainty in Monte Carlo (MC) dose calculations is usually determined using the “history by history” method. Its applicability is based on the assumption that the central limit theorem (CLT) can be applied such that the dispersion of repeated calculations can be modeled by a Normal distribution. The justification for this assumption, however, is not obvious. The concept of stochastic quantities used in the field of microdosimetry offers an alternative approach to assess uncertainty. This leads to a new and simple expression.

Methods

The value of the MC determined absorbed dose is considered a random variable which is comparable to the stochastic quantity specific energy, z. This quantity plays an important role in microdosimetry and in the definition of the quantity absorbed dose, D. One of the main features of z is that it is itself the product of two other random variables, specifically of the mean dose contribution in a ‘single event’ and of the mean number of such events. The term ‘single event’ signifies the sum of energies imparted by all correlated particles to the matter in a given volume. The similarity between the MC calculated absorbed dose and the specific energy is used to establish the ‘event by event’ method for the determination of the uncertainty. MC dose calculations were performed to test and compare both methods.

Results

It is shown that the dispersion of values obtained by MC dose calculations indeed depend on the product of the mean absorbed dose per event, and the number of events. Applying methods to obtain the variance of a product of two random variables, a simple formula for the assessment of uncertainties is obtained which is slightly different from the ‘history by history’ method. Interestingly, both formulas yield indistinguishable results. This finding is attributed to the large number of histories used in MC simulations. Due to the fact that the values of a MC calculated absorbed dose are the product of two approximately Normal distributions it can be demonstrated that the resulting product is also approximately normally distributed.

Conclusions

The event by event approach appears to be more suitable than the history by history approach because it takes into account the randomness of the number of events involved in MC dose calculations. Under the condition of large numbers of histories, however, both approaches lead to the same simple expression for the determination of uncertainty in MC dose calculations. It is suggested to replace the formula currently used by the new expression. Finally, it turned out that the concept and ideas that were developed in the field of microdosimetry already 50 years ago can be usefully applied also in MC calculations.

Keywords

Monte Carlo calculation
Absorbed dose
Uncertainty
Microdosimetry
Specific energy
==== Body
pmc1 Introduction

1.1 Background

It is in the nature of MC simulations applied to the calculation of absorbed dose that the results are subject to a random distribution. This becomes obvious when the absorbed dose is determined on the basis of series of repeated calculations under the condition that each time a different seed for the initialization of the random generator is used. Due to this randomness the true value of the absorbed dose cannot be exactly determined even by a sample of MC simulations. Consequently, the finally reported result is the mean value of a sample, and it is complete only when accompanied by a statement of its associated uncertainty. This situation is very similar to measurements of a quantity of interest. The document “Guide to the expression of uncertainty in measurement” [1] establishes general rules to deal with the subject of uncertainties. Recommendations in this document are also applicable for MC calculations.

State of the art for the determination of uncertainty in MC calculations is the so-called ‘history by history’ method. It was suggested by Sempau et al. [2] and in-depth explained by Walter et al. where its implementation in the MC codes of the EGSnrc system [3] is described. Further applications are addressed in the literature [4], [5].

The validity of the history by history method is based on two assumptions:(1) The total energy absorbed by the matter in a volume of interest is the mean of a stochastic quantity called energy imparted. This quantity is the sum of all energy depositions which take place during a given history. The term ‘history’ comprises the simulation of the energy depositions by a single primary particle and all its secondaries. It is assumed that the single values of the quantity energy imparted which belong to different histories (and thus to different primary particles) are not correlated.

(2) At a given number of histories performed during a MC simulation, the values of energy imparted present a sample of their underlying distribution. Consequently, the sample mean is also a random variable. As this sample mean is obtained by summation, it is assumed that the CLT can be applied. The CLT [6] states that when a large number of sample values of a random variable are summed up, the resulting sum approaches a Normal distribution even if this random variable itself is not normally distributed.

Under these assumptions, the variance and hence the standard uncertainty of the absorbed dose in the VOI per number of histories can be determined using the sample variance of a Normal distribution. Moreover, this approach can be combined with a ‘clever trick’ of programming [2], [3] such that the standard uncertainty can be easily updated after each new history (‘on the fly’).

1.2 Is the CLT really applicable?

The CLT can be applied at many situations. It therefore appears reasonable to apply the CLT at the sum of the energy imparted. It is, however, not obvious whether this application is justified. It is a fact that the distribution of the energy imparted per history includes a very large number of zero contributions. As an example, for a spherical air-filled cavity of 5 mm in diameter placed at 5 cm depth in a water phantom and irradiated with a 10 × 10 cm2 6 MV photon beam, the ratio between zero and non-zero values of energy imparted is in the order of 10000 to 1. This characteristic can be mathematically expressed by saying that the density distribution of the energy imparted per history includes to a large part a Dirac delta function for the probability of no energy imparted. The influence of this characteristic on the applicability of the CLT therefore needs a closer look.

1.3 Modified approach based on stochastic quantities for microdosimetry

A possible approach to answer the question whether the distribution of the sum of the energy imparted per history complies with a Normal distribution simply consists of a statistical test procedure. For that, however, a large number of repeated MC calculations must be carried out. In this work a further approach for the determination of uncertainty is provided which is based on properties of stochastic quantities as introduced in the field of microdosimetry. The aim of this paper is to present this alternative approach and to discuss consequences for the determination of uncertainty.

2 Methods

2.1 Definitions and properties of the stochastic radiological quantity of specific energy

The formalism described below is based on known properties of stochastic radiological quantities as given in the ICRU Report 85a [7]. For this work the specific energy, z is of interest. This quantity was already defined in 1968 [8], and then included in the ICRU report 36 “Microdosimetry” [9] and again in ICRU Report 85a [7]. In order to facilitate further reading, the following definitions are repeated:(a) The energy deposit, εi, is the energy deposited in a single interaction, i.

(b) The energy imparted, ε, to the matter in a given volume is the sum of all energy deposits in the volume.

(c) The specific energy, z, is the quotient of ε by m, where ε is the energy imparted by ionizing radiation to matter in a volume of mass m (subsequently denoted as volume of interest and abbreviated as VOI).

(d) Since the energy deposit in a single interaction is a random variable, the energy imparted and the specific energy are random variables too.

The term “energy-deposition event” (short: event) also needs a clear definition. It represents the imparting of energy to matter in a volume by statistically correlated particles only. Statistically correlated means that the energy imparted to the matter in that volume is due to a single primary particle and to all secondaries produced by it. With this definition, the following distinction can be made: The entire energy imparted in a VOI can belong to one energy-deposition event or to a number of several energy deposition events; for example, the single energy deposits might belong to one or several independent particle trajectories. This approach is comparable to the introduction of the energy imparted per history as described above. However, an event is analog to a history only in the case that a dose contribution really occurs in the VOI. Correspondingly, a history with zero contribution in the VOI is not counted as an event. Therefore, at a given number of histories, the number of events is usually much smaller than the number of histories.

Main properties of the specific energy according to ICRU Report 36 [10] and ICRU Report 85a [7] are:• The distribution of z is characterized by the probability density, f(z,D), which depends on the absorbed dose D in the VOI.

• The specific energy z can be due to one or more energy deposition events. The number of events, Nevent which contributes to a specific energy is, in general, distributed at random and can be described by a Poisson distribution.

• The distribution of the specific energy deposited in single events is termed a single event distribution described by the probability density, fs(z).

• The mean value of the specific energy with respect to the single event distribution fs(z):(1) z¯s=∫zfszdz

is called the frequency-mean of the specific energy per event. It is a non-stochastic quantity.

It was A.M. Kellerer who has demonstrated in the early days of microdosimetry in 1968 that there exist two fundamental relations between the probability density of the specific energy, f(z,D) and the probability density of the specific energy deposited in a single event, fs(z)dz. The first relation refers to the mean of the specific energy, z¯ in a VOI, which can be expressed as the product of the mean number of events, N¯event and the corresponding frequency-mean of the specific energy in a single event:(2) z¯=∫zf(z,D)dz=N¯event·z¯s

The second relation refers to the relative variance of the specific energy, Vrelz:(3) Vrelz=Vzz¯2=VNeventE2Nevent+1ENeventVszEs2z

where VNevent and ENevent are the variance and the expectation value of Nevent occurring in the VOI, and Vsz and Esz are those of the specific energy in a single event being calculated using the single event distribution fs(z). [7], [8], [9].

2.2 Variance and relative uncertainty of the MC calculated absorbed dose

As already stated by G.A. Carlsson [10], the definition and the property of the specific energy z are directly transferable to the absorbed dose in a dosimeter and thus to the VOI as obtained from a MC calculation, dMC. In MC simulations, the equivalent of an energy deposit, εi is a single energy loss along a charged particle track, dei which is obtained either by a single step or by a condensed history step [11]. The dei from statistically correlated particles can be lumped together, yielding the absorbed dose qj in a single event j, when divided by the mass of the VOI:(4) qj=∑ideij/m

Its mean is given by:(5) q¯=1nevent∑i=1neventqi

The total MC calculated absorbed dose, dMC then becomes:(6) dMC=nevent·q¯

Note: The values of nevent and q¯ and in consequence also that of dMC are single realizations obtained in a MC calculation which are taken out from a probability sample of Nevent and Q¯ at a given number of histories. These probability samples are random variables denoted below as Nevent, Q¯ and DMC. Therefore, the uncertainty of a single value of dMC is given by the positive root of the variance of DMC.

It is of interest to make use of eq. (3) to determine the relative variance of DMC according to:(7) VrelDMC=VNeventE2Nevent+1ENeventVsQEs2Q

Unfortunately, the distribution of the absorbed dose in a single event, fs(q) cannot be simply described by a known mathematical expression such as a Normal distribution. As an example, Fig. 1 shows the MC calculated frequency distribution of the absorbed dose qj in a single event as obtained by a MC test calculation as described below. In this case it refers to the spherical air-filled cavity with a diameter of 5 mm positioned in a 10 × 10 cm2 field size at a depth of 5 cm in a water phantom. A Normal distribution having the same mean value is also shown for comparison. The comparison clearly confirms that the considerably asymmetric distribution of the qj in a single event is not well represented by a Normal distribution.Figure 1 Relative frequency distribution of the absorbed dose qj in a single event and that of a Normal distribution with identical mean value.

However, one can apply the CLT to the sum of the qj such that for the distribution of Q¯ a Normal distribution can now be assumed. This allows the variance of DMC to be calculated by the variance of the product of Nevent and Q¯, where both have approximatively a Normal distribution and known estimates for the mean and the variance.

There are several approaches to obtain the variance of the product of the two random variables Nevent and Q¯ (here the relative variance of DMC). It is generally assumed that nevent≫1.

(1) The most general approach for the variance of the product is given by [12]VrelDMC=VNevent·Q¯N¯event2·Q¯2=

(8) =VNeventN¯event2·VQ¯Q¯2+VNeventN¯event2+VQ¯Q¯2

(2) Following Ware and Lad [13], the product of two Normal distributions is not necessarily a Normal distribution. However, it approximates a Normal distribution if the following condition is fulfilled:(9) VNeventN¯event2·VQ¯Q¯2≪VNeventN¯event2+VQ¯Q¯2

This leads to:(10) VrelDMC=VNeventN¯event2+VQ¯Q¯2

Using the estimates for the variance of Nevent and Q¯ as given below, it is easy to demonstrate that the condition (9) is indeed fulfilled.

(3) Following GUM [1] the law of propagation for the combined uncertainty of a product can be applied which again leads to:(11) VrelDMC=VNeventN¯event2+VQ¯Q¯2

At all three approaches, the variance VNevent and VQ¯ in eq. (8), (10), or (11) are replaced by their known estimates:(12) V^Nevent=nevent

(13) V^Q¯=∑j=1neventqj-q¯2nevent2=∑j=1neventqj2nevent2-∑j=1neventqj2nevent3

Also the expectation values involved in D¯MC=ENevent·EQ¯ are replaced by the corresponding estimates:(14) E^Nevent=nevent

(15) E^Q¯=1nevent∑j=1neventqj

This leads for all three approaches to the following common expression for the relative Type A standard uncertainty uDMC/D¯MC:(16) uDMC/D¯MC=+∑j=1neventqj2∑j=1neventqj

Also note: The summation of the qj and qj2 leads to identical values regardless whether performed over the number of histories or over the number of events. Therefore, eq. (16) can also be written as:(17) uDMC/D¯MC=+∑j=1nhistqj2∑j=1nhistqj

The method to calculate the uncertainty based on eq. (17) is subsequently referred to as the ‘event by event’ method.

2.3 MC Test calculations

Absorbed dose values were calculated in a spherical air-filled cavity with a diameter of 5 mm positioned in a water phantom at a depth of 5 cm. SSD was 95 cm, a point source simulating the radiation quality of a 6 MV beam was used. Calculations were carried out with two different field sizes: 10 × 10 cm2 and 2 × 2 cm2. Absorbed dose in the cavity was calculated with a modified version of DOSXYZ of the EGSnrc system [14], [15].

In order to perform a statistical analysis, a number of 200 dose calculations (called batches) were carried out, each starting with a different seed for the random generator and performed with 108 histories.

The relative Type A standard uncertainty of absorbed dose per history was determined:(a) for each of the 200 batch calculations with 108 histories

(b) for a combination of an increasing number of batches each with 108 histories

(c) for single batches with an increasing number of histories

The methods used were the history-by-history method according to Walter et al. [3], the event by event method of this work, and the batch method as described in [3].

3 Results

3.1 Distribution of number of events and of the absorbed dose per event

Results refer to the 10 × 10 cm2 field size. In Fig. 2 the frequency distributions of the two random variables, the number of events and the mean absorbed dose per event as obtained from the 200 batch calculations with 108 histories are compared with the distribution functions associated with these random variables (Poisson distribution for Nevent, Normal distribution for Q¯). Distributions are shown as per cent deviation from their mean values which were: N¯event=3417 and Q¯= 0.165 μGy. These results well confirm the comparability of the stochastic characteristic between the specific energy and the MC calculated absorbed dose.Figure 2 Normalized frequency distributions of the nevent (left) and q¯ (right) obtained from 200 batch calculations; the associated distribution functions (see text) are shown as dashed curves.

The two random variables Nevent and Q¯ are not correlated. This is exemplified by the scatter plot in Fig. 3 which shows the number of events, nevent against the mean absorbed dose per event, q¯ as obtained in each of the 200 test calculations of absorbed dose.Figure 3 Scatter plot of the number of events, nevent against the mean absorbed dose per event resulting, q¯ from 200 MC calculations of absorbed dose. Dashed lines represent mean values.

3.2 Distribution of the MC calculated absorbed dose

The frequency distribution of dMC obtained as the product of nevent and q¯ is shown in

Fig. 4, again as per cent deviation from its mean value. The frequency distribution is compared with a Normal distribution with a standard deviation equal to ∑j=1nhistqj2/∑j=1nhistqj2.Figure 4 Normalized frequency distribution of the MC calculated absorbed dose as per cent deviation from its mean value. The associated Normal distribution function is shown as dashed curve.

3.3 Comparison of uncertainty results

In Fig. 5, the relative standard uncertainty is shown in dependence of the number of histories or the combined number of batches. Values obtained by the history by history method, by the event by event method of this work, and also the batch method were compared. Since the result of the batch method depends on which sample of batches is selected, the mean uncertainty and its variation at a given number of batches are shown. These values were obtained in a number of 100 random selections. Two observations are made:1. There is no numerical difference at all between the uncertainties evaluated by the history method and that by the event method.

2. The uncertainty determined with the batch method may considerably differ from that of the event by event method. In this test the difference becomes particularly obvious at batch numbers less than 10, such that a slight underestimation may result on average.

Figure 5 Comparison between the relative standard uncertainty obtained with the history by history method (solid line), the event by event method (open circles) and the batch method (triangles) as a function of the number of histories or the number of batches. The standard variation of the uncertainty as obtained with the batch method is marked as grey area. Left: 10 × 10 cm2 field size; right: 2 × 2 cm2 field size.

4 Discussion

Given that the distribution of the MC calculated absorbed dose values can indeed be well approximated by a Normal distribution, the approach of this work may appear somewhat sophisticated compared to the well-established history by history method. However, both methods are equivalent with respect to the results. This becomes particularly obvious when the expression for the relative standard uncertainty according to the history by history method is written as:(18) uDMC/D¯MC=+1Nhistory-1Nhistory·∑i=1Nhistoryqi2∑i=1Nhistoryqi2-1

For a large number of histories, equation (18) directly leads to the same simple expression of eq. (17) which was derived for the event by event method.

One should also mention, that the use of the event method was already addressed before by Ma et al. in the description of a MC dose calculation tool for radiotherapy treatment planning [16]. They stated that the relative statistical uncertainty can be estimated by an expression equivalent to eq. (16). They have, however, not mentioned that in this expression the summation over the number of events (which is not directly known in standard MC codes) can be simply replaced by the summation over the known number of histories.

One may argue that the approach in this work apparently does not lead to a change for the value of the Type A uncertainty in MC dose calculations. However, its derivation offers a better understanding of the additional influence of the randomness of the number of events. For instance, it can easily explain the well-known increase of uncertainty with field size at a given number of histories. It was also demonstrated that the event method is more reliable than the batch method as already stated for the history method by Walter et al. [3] and Sempau et al. [2].

Finally, the event-by-event method must be discussed with reference to the use of phase-space files as radiation source. Sempau et al. [4] as well as Walter et al. [2] have addressed the inevitability of a ‘latent variance’ contained in phase-space files. This term was introduced to distinguish between the uncertainty in a MC dose calculation due to the random nature of the particle transport to that due to the statistical fluctuations in the phase-space file. Another cause for uncertainty can be the recycling of particles in a phase-space source before moving on to the next one. In order to assess the overall uncertainty, the event by event method can therefore be applied in the simple form of the expression (17) provided that phase-space files are not used as radiation source.

5 Conclusion

It was shown that the concept and ideas that were developed in the field of microdosimetry to characterize the stochastic character of energy deposition on a microscopic scale can be usefully applied also at MC dose calculations. The main reason is that the stochastic character of radiation becomes visible not only in small volumes, but also under the condition of low doses, a condition which well applies to the MC determination of absorbed dose. For example, MC calculated values of absorbed dose in photon beams are typically in the order of μGy or less. Application of the microdosimetric approach demonstrates that DMC is the product of two random variables, specifically of the number of events and of the mean absorbed dose per event. This fact has two consequences: (1) It leads to an alternative and simple formula for the determination of uncertainty, and (2) it can provide a proof that for a large number of histories the results of repeated absorbed dose calculations approach a Normal distribution.

It is suggested, that the formula (18) of the history by history method which is implemented in MC programs for dose calculations such as in the EGSnrc system should be replaced by the new and simple expression of eq. (17) for two reasons: (1) the event by event method takes into account the randomness of the number of events and therefore appears to be more suitable, and (2) the expression is related to the usually large number of histories involved in MC calculations of absorbed dose

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

We gratefully acknowledge the comment by Frédéric Tessier from the National Research Council, Canada.
==== Refs
References

1 Evaluation of measurement data—Guide to the expression of uncertainty in measurement. JCGM 2008.
2 Sempau J. Wilderman S.J. Bielajew A.F. DPM, a fast, accurate Monte Carlo code optimized for photon and electron radiotherapy treatment planning dose calculations Phys Med Biol 45 2000 2263 2291 10958194
3 Walters B.R.B. Kawrakow I. Rogers D.W.O. History by history statistical estimators in the BEAM code system Med Phys 29 2002 2745 10.1118/1.1517611 12512706
4 Sempau J. Sánchez-Reyes A. Salvat F. Oulad ben Tahar H. Jiang S.B. Fernández-Varea J.M. Monte Carlo simulation of electron beams from an accelerator head using PENELOPE Phys Med Biol 46 2001 1163 1186 11324958
5 Andreo P. Monte Carlo techniques in medical radiation physics Phys Med Biol 36 1991 861 920 1886926
6 Fischer H. A history of the central limit theorem 2021 Springer New York, Dordrecht, Heidelberg, London
7 ICRU. Fundamental Quantities and Units for Ionizing Radiation (Revised), ICRU Report 85a. Bethesda: International Commission on Radiation Units and Measurements; 2011.
8 Kellerer AM. Mikrodosimetrie, Grundlagen einer Theorie der Strahlenqualität. GSF B-1 (1968)
9 ICRU. Microdosimetry, ICRU Report 36. Bethesda: International Commission on Radiation Units and Measurements; 1983.
10 Carlsson G.A. Theoretical basis for dosimetry Kase K.R. Bjärngard B.E. Attix F.H. The dosimetry of ionizing radiation 1985 Academic Press
11 Kawrakow I. Mainegra-Hing E. Rogers D.W.O. Tessier F. Walters B.R.B. The EGSnrc Code System: Monte Carlo Simulation of Electron and Photon Transport NRCC Report PIRS-701 2017
12 Goodman L.A. On the exact variance of products J Am Stat Assoc 55 1960 708 713 10.2307/2281592
13 Ware R. Lad F. Approximating the Distribution for Sums of Products of Normal Variables 2003 Department of Mathematics and Statistics, University of Canterbury, New Zealand Christchurch, New Zealand
14 Kawrakow I. Mainegra-Hing E. Rogers D.W.O. NRCC Report PIRS-877; EGSnrcMP: the Multi-Platform Environment for EGSnrc 2006 National Research Council of Canada Ottawa
15 Hartmann G.H. Andreo P. Kapsch R.-P. Zink K. Cema-based formalism for the determination of absorbed dose for high-energy photon beams Med Phys 48 11 2021 7461 7475 10.1002/mp.15266 34613620
16 Ma C.-A. Li J.S. Pawlicki T. Jiang S.B. Deng J. Lee M.C. Koumrian T. Luxton M. Brain S. A Monte Carlo dose calculation tool for radiotherapy treatment planning Phys Med Biol 47 2002 1671 1689 12069086
