The Monte Carlo demonstrator widget

The Monte Carlo demonstrator widget

Interactive - have a play now. 

Monte Carlo calculations for radiation transport

Dr Chris Robbins, Grallator

The Monte Carlo demonstrator

Many radiation transport related calculations make use of the so-called Monte Carlo method. With this method individual photons (particles) are tracked from birth to death in a system. The paths taken and interactions made on the way provide a robust method of calculating physical values such as radiation intensity, photon/particle energy, dose rate etc.

The interactive calculator performs the Monte Carlo calculations for the simplest possible system to demonstrate the inverse-square law of radiation intensity. This can be performed with no photon/matter interaction to obtain the pure inverse-square law, or with photon absorption to obtain an inverse-square law modified by an exponential reduction in intensity. A typical starting screen is shown below.
 

The screen gives sliders to change the number of Monte Carlo particles tracked and the absorption coefficient of the material through which the photons are travelling. For the moment we'll stick to the setting where no absorption takes place. Values of the Monte Carlo predicted dose rate along with a 95% confidence interval for the error in the value are given at distances of 1m and 2m. The screen also gives the actual value as calculated from the well know inverse-square equation \[ D=\frac{100}{4\pi r^2}, \] where 100 is a normalising value to obtain the specified/measured dose rate at 1m (in unspecified relative units) and \(r\) is the distance from the source. Using these values the calculated dose rates at 1m and 2m are \[ \begin{align} D \left(r=1\right)&=7.958 \\ D \left(r=2\right)&=1.989 \end{align} \] It is observed that the dose rate at 1m is 4 times larger than the dose rate at 2m.

How does the Monte Carlo calculation perform? With 100 tracks the predictions are ball-park but the error ranges are quite large: of the order of 100%.

Recalculating with the same parameters produces a different result based on a new set of randomly tracked photons, as demonstrated below. Note, the changes to the mean are observed to be larger than the changes to the size of the error range.

It should be noted that in the lower image a value of zero is returned for the dose rate at 2m (and the error). The reason for this is that only 100 photons are being tracked. With such a low number there is a reasonable probability that none of the tracked photons actually passed through the Monte Carlo tally area.

Moving the number of tracks pointer to the right increases the number of tracks by an order of magnitude. After each change the "Calculate" button should be clicked to update the values. The examples below shows a couple of cases.

It is observed that as the number of tracks increases, the error estimate of the dose rate decreases. Note also that sometimes the actual value is outside of the confidence interval shown, as in the lower figure. Recall, the error bar is the 95% interval, so it would be expected that 5% of the time the actual value would be outside of this error interval.

The table below shows how the absolute and relative error falls as the number of tracks increases. Looking at the absolute values it is observed that the error reduces by one order of magnitude for every two orders of magnitude increase in the number of tracks. In fact the results show the error scales as \(\frac{1}{\sqrt{N}}\), where \(N\) is the number of tracks (and there is a good mathematical reason for this!).

It is also noted that the relative error for the 2m estimate is about twice the size of the 1m estimate. The reason for this is that there is a lower probability of tracks crossing a unit area at 2m than there is at 1m. In Monte Carlo, if you've got a small region of interest or one that's far from a source, you need to ensure there are enough tracks to sample the region.

Adding absorption

The second slider adds photon absorption to the model. This modifies the actual dose rate by multiplying by the usual exponential reduction function \[ D=\frac{100}{4\pi r^2}e^{-\mu r}, \] For a value of \(\mu = 0.5\), the value of \(e^{-\mu r}\) is 0.6065 and0.3679 for \(r=1\)m and \(r=2\)m respectively (rounded to 4 d.p.). \[ \begin{align} D \left(r=1\right)&=7.958 \times 0.6065 = 4.827 \\ D \left(r=2\right)&=1.989 \times 0.3679 = 0.732 \end{align} \]

The image below shows a Monte Carlo calculation using 10,000,000 tracks and an absorption coefficient of \(\mu = 0.5\) and giving a reasonably accurate estimate.

[Ionactive comment: What this indicates to us is like in many other circumstances,  the output is only as good as the input. Monte Carlo might give you the answer you desire or expect, but you have to know the statistics of the process, the sampling (tally) strategy and similar.]

What's under the bonnet

The following text is a very brief summary of the mathematics of Monte Carlo method and is by no means a rigorous treatment - please consult a textbook to fill in the places where statements are made without support or proof!

The Monte Carlo method is a numerical approach to solving problems using random sampling. The archetypal example that is often used is the calculation of \(\pi\) by picking random \( \left( x,y\right) \) points in a unit square and counting the fraction that satisfies \( x^2+y^2 \le 1 \) which, it can be shown, approaches \( \frac{\pi}{4}\) as the number of sample points increases. However, this example doesn't really address a concept in radiation transport so we won't use it! And we'll leave the "it can be shown" bit as an exercise for the reader...

A radiation problem in a two dimensional universe

One of the most fundamental concepts in radiation transport is the fall off of radiation intensity with distance, which is often referred to as the inverse square law. Let's consider a two-dimensional version of this problem where we have a source at the origin and we wish to calculate the number of photons crossing a unit length of arc at a distance \(r\), as shown below. Note, in this problem the photons are confined to the plane.

In the above diagram the red arcs have a length of 1 unit. A photon emitted by the source with a random angle \(\theta\) in radians as measured from the positive \(x\)-axis will cross a red arc at a distance \(r\) if (using the express \(s = r\theta\)) \[ 0 \leqslant \theta \leqslant \frac{1}{r}=\theta_{arc}. \] The red arcs will be our tally counting area. As photons can be emitted at any angle \( 0 \leqslant \theta \lt 2 \pi\) it should be reasonably obvious that the ratio of the number of photons that cross a red arc to the total number of photons emitted is equal to the ratio of the length of the arc to the circumference of the circle, \[ \frac{N_{arc}}{N_{emitted}}=\frac{l_{arc}}{c_{circle}}=\frac{1}{2\pi r}=\frac{r \theta_{arc}}{2\pi r}=\frac{\theta_{arc}}{2\pi}. \]

A Monte Carlo recipe to evaluate the photon intensity per unit length, \( I_{arc} \), at a distance \( r\) can be constructed for this set up as: 

• Set intersection counter to zero, \(N_{arc}=0\)
• loop N times:
• Select a uniform random number over the interval \(0 \leqslant u \lt1\)
• Calculate the associated random direction of photon emission as \(\theta = 2\pi u\)
• Add 1 to \(N_{arc}\) if \( 0 \leqslant \theta \leqslant \theta_{arc} \)

At the end of the calculation the value of \( I_{arc} \) is approximated by \[ I_{arc} \approx \frac{N_{arc}}{N} \]

NOTE!

As the direct calculation above shows the exact value is \( I_{arc} =\frac{1}{2\pi r} \), the radiation intensity in this two-dimensional universe falls as \(\frac{1}{r}\) and not as \(\frac{1}{r^2}\)(the inverse square law).For this reason radiation transport calculations have to be performed in three dimensions, and this will be addressed below! However, it is useful to continue in two dimensions a little longer to draw out more features of this simple, one variable Monte Carlo calculation.

How many photons should I track?

The obvious question to ask with the above algorithm is how accurate is the approximation given \(N\) paths, or turned around, how big should \(N\) be to get a desired accuracy? To look at this problem we need to consider two statistical measures: the expectation value and the variance.

Define the function \(f\left( \theta \right)\) on the domain \( 0 \le \theta \lt 2 \pi\) such that \[ f\left( \theta \right) = \begin{cases} 1 &\text{if } 0 \leqslant \theta \leqslant \theta_{arc} \\ 0 &\text{if } \theta \gt \theta_{arc} \end{cases} \] The function \(f\) is identified as a counting function and is the mathematical representation of the Monte Carlo recipe line:
 

• Add 1 to \(N_{arc}\) if \( 0 \leqslant \theta \leqslant \theta_{arc} \)

The expectation (mean) value of \(f\) over \( 0 \le \theta \lt 2 \pi\) is given by \[ \mu=\frac{\int_{0}^{2\pi}f (\theta)d\theta} {\int_{0}^{2\pi}d\theta}, \] where the denominator provides the appropriate normalisation, which is easily evauated as \(2\pi\). The numerator is evaluated as \[ \int_{0}^{2\pi}f (\theta)d\theta= \int_{0}^{\theta_{arc}} 1 \cdot d\theta + \int_{\theta_{arc}}^{2\pi} 0 \cdot d\theta = \theta_{arc} \] So that \[ \mu = \frac{\theta_{arc}}{2\pi}. \] This looks familiar from the earlier text above! The function \(f\) will also have a variance \(\sigma^2\)

Monte Carlo uses random sampling to generate a set of values of \(\theta\), \( \{ \theta_1, \theta_2, \ldots ,\theta_N\} \), from which is constructed the set of count function values \(f\left( \theta \right)\) \( F = \{ f_1, f_2, \ldots f_N \} \). The sample mean value is calculated as \[ \overline{f}=\frac{1}{N}\sum_{i=1}^{N}f_i \] and \(\overline{f} \rightarrow \mu\) as \(N\rightarrow \infty \).

So far an approximation of the mean has be constructed from a set of Monte Carlo trials but we still don't know how good the approximation is. One way to get an idea would be to rerun the Monte Carlo again and again with different sets of random angles. Each run would typically give a different sample mean value, but over a large number of cases the population of sample means converges to a normal distribution centred on the real mean, \(\mu\), and have some variance, \(\sigma_{\overline{f}}^2\), which is related to the variance of \(f\), \(\sigma^2\), by \[ Var\left[\overline{f}\right] = \sigma_{\overline{f}}^2 = \frac{\sigma^2}{N} \] This is the central limit theorem.

Although the actual value of \(\sigma^2\) may not be known for \(f\) (and it is not calculated here), it can be estimated from the sample variance as \[ s^2=\frac{1}{N-1}\sum_{i=1}^{N} \left( f_i - \overline{f}\right)^2 \approx \sigma^2, \] which gives \[ \sigma_{\overline{f}} \approx \sqrt{\frac{s^2}{N}}=\frac{s}{\sqrt{N}} \] This can be used to give a confidence interval estimate of the sample mean in the form \(\overline{f} \pm z \sigma_{\overline{f}} \), where \(z\) is determined by the required confidence interal. For a two-side 95% confidence interval \(z=1.96\) and Monte Carlo results are presented as \(\overline{f} \pm 1.96\sigma_{\overline{f}}\). The implication of the this result is that the error on the mean estimate scales as \( \frac{1}{\sqrt{N}} \), leading to the statement:

If you want to reduce the error on the caluclated estimate of the mean by a factor 10, you have to increase the number of Monte Carlo tracks by a factor of 100. 

Aside In practice when calculating the mean and variance you do not need to store all the \(f_i\) values and first loop to calculate the mean and then loop again to calculate the variance. Instead, as the counting values are either zero or non-zero, it is computationally more efficient to keep a running total of non-zero values from which the mean, \(\overline{f} \), is calculated by dividing this total by the number of tracks after all the tracks have been followed. To estimate the variance you can use the expectation value identity \[s^2 = \mathbb{E} \left[ \left(\ f_i - \overline{f} \right)^2 \right] =\mathbb{E} \left[ f_i^2 -2f_i\overline{f} + \left(\overline{f}\right)^2 \right] = \mathbb{E} \left[ f_i^2 \right] - 2\overline{f} \mathbb{E} \left[f_i \right] + \left(\overline{f}\right)^2 = \mathbb{E} \left[ f_i^2 \right] - \left(\overline{f}\right)^2 \] In this expression, \( \mathbb{E} \left[ f_i^2 \right] \) is the mean of the squared values of \(f_i\), which are computed as a running total. 

The above gives a method for error estimation given a number of tracks \(N\). However, with the assumption of a normally distrubuted error of the mean and an approximation for \(\sigma^2\), the expression can be rearranged to give an (approximate) number of tracks required to obtain a desired error of the mean as \[ N=\frac{s^2}{\sigma_{\overline{f}}^2} \]

Three dimensions

The above two-dimensional world was useful for describing how Monte Carlo works but doesn't result in the inverse-square law. To obtain a \(\frac{1}{r^2}\) fall off of intensity requires three dimensions. As with the previous case we'll use a point source and calculate the intensity per unit surface area on the surface of a sphere of radius \(r\) centred on the source.

Instead of one angle to define a path direction we now need two. The first angle, \(\phi\), is the azimuthal angle and has possible values \(0 \leqslant \phi \lt 2\pi\). The second angle, \(\theta\), is the polar angle and has possible values \(0 \leqslant \phi \lt \pi\).

The tally counting surface area is similar to the two dimensional calculation where a count is made if \( 0 \leqslant \phi \leqslant \phi_{A}\) and \( 0 \leqslant \theta \leqslant \theta_{A}\), for values of \(\phi_{A}\) and \(\theta_{A}\) that give a unit surface area. Such angles can be determined by assuming \(\phi_{A}=\theta_{A}=\alpha\) and solving \( r^2\alpha \left(1-\cos \alpha \right)-1=0 \), where \(r\) is the radius. (Derivation not given but this comes from the definition of a surface are element in polar coordinates being \( dA = r^2\sin\theta \space d\phi \space d\theta \).)

A naive way to run a Monte Carlo calculation on this geometry would be: 

• Set intersection counter to zero, \(N_{A}=0\)
• loop N times:
• Select a uniform random number over the interval \(0 \leqslant u \lt1\)
• Select a uniform random number over the interval \(0 \leqslant v \lt1\)
• Calculate the associated random direction of photon emission as \(\phi = 2\pi u\), \(\theta = \pi v\)
• Add 1 to \(N_{A}\) if \( 0 \leqslant \phi \leqslant \phi_{A}\) and \( 0 \leqslant \theta \leqslant \theta_{A}\)

However, this would give the wrong answer! The reason is that because of the definition of polar coordinates, the sampling is not uniformly distributed over the surface of a sphere but is instead biased towards the poles. The diagram below attempts to depict this. Eight sample azimuthal points are more closely packed when the polar angle is near the pole than they are when the polar angle is at the equator.

The correct way to run this Monte Carlo is to cosine correct the polar angle so that:

Set intersection counter to zero, \(N_{A}=0\)

loop N times:

Select a uniform random number over the interval \(0 \leqslant u \lt1\)

Select a uniform random number over the interval \(0 \leqslant v \lt1\)

Calculate the associated random direction of photon emission as \(\phi = 2\pi u\), \(\theta = \cos^{-1} \left( 2v-1\right) \)

Add 1 to \(N_{A}\) if \( 0 \leqslant \phi \leqslant \phi_{A}\) and \( 0 \leqslant \theta \leqslant \theta_{A}\)

Revisiting the confidence interval

The text above describes how the 95% confidence interval is estimated from the sample mean and sample variance by using the results of the central limit theorem. To demonstrate the validity of this (and now that we have a three-dimensional description) we can run the Monte Carlo simulation a number of times and plot a histogram of the distribution of mean values predicted. In the graph below the Monte Carlo simulator has been run 10,000 times with each case having 100,000 tracks, and the frequency of the predicted mean plotted.

The shape is Gaussian in appearance and has the following characteristics: \[ \begin{align} \text{mean of means} \space\space\space \overline{\overline{f}}&=7.957 \\ \\ \text{variance of means} \space\space\space \sigma_{\overline{f}}^2&=0.0071825 \end{align} \] which gives a 95% confidence interval of 0.166. This confidence interval based on the distribution of 100,000 means is very close to the value of 0.168 calculated from a single Monte Carlo sample as provided by the Monte Carlo simulator, see below.

As can be seen, the central limit theorem shows how useful and powerful it is!

Random numbers

The Monte Carlo method makes use of random numbers so you need a good random number generator! A good random number generator should have a long period (number of values produced before the sequence repeats) as well as satisfying a number of statistical tests to show there aren't any significant correlations in results. In the early days of the C language the minimum period for the inbuilt pseudo random number generator was 32767, i.e. a series of 32767 values would be produced before the sequence repeated. This is useless for Monte Carlo calculations! If you wanted to run 1,000,000 tracks you won't actually sample 1,000,000 unique tracks. Instead you'll be sampling the same 32767 tracks about 30 times and your errors won't get smaller. More modern generators such as the Mersenne Twister (MT19937) have much longer periods, in this case \(2^{19937}-1\), and have better statistical properties.