# Appendix K - Monte Carlo Models

Many of the problems encountered in making a realistic detailed model of
a detector are so complex as to not be susceptible to analytic
techniques. Nevertheless, we can break the problem of modeling a very
complex system down into a series of choices for the relevant dynamical
variables. The Monte Carlo method allows us to choose those variables
and hence to construct such a model.

For example, suppose we want to model ![](media/image1.emf) ray
production in a medium by incident muons. The choices involved are first
at what depth is the ![](media/image2.emf) ray made. We choose x out of
a distribution ![](media/image3.emf) where ![](media/image4.emf) is the
cross section for ![](media/image5.emf) production. If there is indeed
![](media/image6.emf) production within the active volume, the dynamics
depends on a single variable, as has been mentioned previously. We can,
for example, pick the recoil electron kinetic energy from
![](media/image7.emf) from a distribution ![](media/image8.emf). Having
the energy, the recoil angle follows, as derived in Appendix A, from the
2-body kinematics. We now have the position and momentum vectors of the
electron at the point of production. We can find the final state muon as
well if we are to trace it further through the system.

If it is desired that the ![](media/image9.emf) ray be followed in its
path, other choices must be made. Suppose these are no fields, for
simplicity, so that the path is, at least locally, a straight line. We
pick a distance which is short with respect to both the energy loss,
dE~1~/dx, and for which the multiple scattering angle is small. Over
this distance the energy and angles of the electron can be considered to
be constant. We extrapolate the electron path as a straight line. We
pick a new energy ![](media/image10.emf) by removing the ionization
energy from the electron. We also pick new angles using a Gaussian
distribution with $\text{rms}\  = \ \theta_{o}$ about the initial
direction due to multiple scattering in the medium. Clearly, we can make
a rather complicated model out of individual, fairly simple, choices. In
fact, almost all the formulae appearing in the body of the text would,
at some point, be used in a Monte Carlo simulation of the complex
detectors used today. Therefore, we need to be conversant with al the
dynamics discussed in this text if we want to be able to make a
realistic model of a detector.

How do we actually go about choosing the individual dynamical quantities
which define the evolution of the system? First, we need a "random
number generator" which produces a uniformly distributed number r in the
interval (0, 1). These generators are widely available on almost all
computer platforms.

![](media/image11.emf) K.1

Suppose we have a probability distribution d[P]{.underline}(x)/dx
defined such that d[P]{.underline}(x) is the probability for x to occur
between x and x + dx. If x is constrained to the interval (x~min~,
x~max~) then the integral probability is normalized to one, so that the
cumulative probability from x~min~ to x should be uniformly distributed.

$\int_{x_{\text{min}}}^{x}{d\underline{P}(t\frac{)}{\int_{x_{\text{min}}}^{x_{\text{max}}}{\ d\underline{P}(t)\  = \ r}}}\  = \text{}\int_{x_{\text{min}}}^{x}{d\underline{P}}$
K.2

If r = 0 then x = x~min~, while if r = 1, x = x~max~.

A visual interpretation of $r = \int_{}^{}{d\underline{P}}$ is given in
Fig. K.1. The discrete case is shown when we can perform the integral to
find *N* intervals of *x* with equal probability,
$\Delta\underline{P} = \int_{}^{}{d\underline{P} \equiv g(x)}$. If the
inversion, $\Delta x = g^{- 1}(\Delta\underline{P})$ is possible then x
can be chosen by choosing index 1 \< *i* \< *N*, with *x* between x~min~
and x~max~.

![](media/image12.wmf)

Fig. K1: Visualization of intervals in x, $\Delta x$, over which the
integral $\Delta\underline{P}$ is the same.

If this expression can be analytically inverted, we can "pick x out of a
distribution d[P]{.underline}(x)/dx ." Let us examine some
possibilities. If d[P]{.underline}(x)/dx is uniform, we can find a
solution. An example arises in azimuthally uniform scattering, for
example unpolarized Compton scattering. The aximuthal angle *φ* is
chosen as,

![](media/image13.emf) K.3

A more complex behavior is power law dynamics. For example, we can pick
out of ![](media/image14.emf) for Rutherford scattering, or from 1/T^2^
for recoil ![](media/image15.emf) ray kinetic energies. By the way, the
student is strongly urged to try a few of these examples - on the
simplest computer that is accessible and sufficient.

![](media/image16.emf) K.4

The special case *α* = 0 reproduces Eq. K.3.

We have mentioned particle lifetimes several times. A model for unstable
particle decays would require picking a decay time out of an exponential
distribution. Another example is picking a free path L out of a mean
free path, exp(-L/\<L\>).

![](media/image17.emf) K.5

In the energy rather than the time domain any unstable particle is
described by a Lorentzian energy spectrum. The Fourier transform of a
decaying exponential with lifetime ![](media/image18.emf) is a
Lorentzian with energy full width
$\Gamma \  = \ \frac{\hslash}{\tau}$. The uncertainty relation
$\Delta E\Delta t = \hslash$ becomes $\Gamma\tau = \hslash$ in this
special case. This type of behavior also occurs in any damped resonant
system.

$\begin{array}{r}
\frac{\int_{x_{\text{min}}}^{x}\left\lbrack \frac{1}{\left( t - x_{o} \right)^{2} + (\frac{\Gamma}{2})^{2}} \right\rbrack \text{dt}}{\int_{x_{\text{min}}}^{x_{\text{max}}}\left\lbrack \frac{1}{\left( t - x_{o} \right)^{2} + (\frac{\Gamma}{2})^{2}} \right\rbrack \text{dt}} \\
\varphi_{\text{min}}\  = \ 2(x_{\text{min}} - x_{o}\frac{)}{\Gamma} \\
\varphi_{\text{max}}\  = \ 2(x_{\text{max}} - x_{o}\frac{)}{\Gamma} \\
x\  = \ x_{o} + \frac{\Gamma}{2}\text{tan}\left\{ \text{tan}^{- 1}\left( \varphi_{\text{min}} \right)\  + \ r\left\lbrack \text{tan}^{- 1}\left( \varphi_{\text{max}} \right) - \text{tan}^{- 1}\left( \varphi_{\text{min}} \right) \right\rbrack \right\}
\end{array}$ K.6

The most common distribution that is used is, perhaps, the Gaussian
distribution. The integral leads to a non-analytic error function, so it
seems that the technique fails. However, a trick avails; using the joint
probability and displacing the mean to zero, ![](media/image19.emf).

$\begin{array}{r}
d\underline{P}(x)\ d\underline{P}(y) \\
 = \text{}\ e^{\frac{- x^{2}}{2}\sigma^{2}}\text{}e^{\frac{- y^{2}}{2}\sigma^{2}}\text{dxdy} \\
 = \ e^{\frac{- r^{2}}{2}\sigma^{2}}\text{rdrd}\varphi \\
 = \ e^{\frac{- u}{2}\sigma^{2}}\frac{\text{du}d\varphi}{2},\ \text{}u\  = \ r^{2}
\end{array}$ K.7

Therefore, we can pick ![](media/image20.emf) from a uniform
distribution as in Eq. K.3 and r^2^ from an exponential distribution as
in Eq. K.5. The result is, ![](media/image21.emf),
![](media/image22.emf), two Gaussianly distributed variables x and y
with zero mean and $\text{rms}\  = \ \sigma$. The mean can be restored
by addition, $x^{'}\  = \ x\  + \ \left\langle x^{'} \right\rangle$,
$y^{'}\  = \ y\  + \ \left\langle y^{'} \right\rangle$.

Let us now consider angular distributions. In the case of isotropic
dynamics, e.g. decay of a spinless particle, we can use the previous
results for uniform distributions for ![](media/image23.emf) and
![](media/image24.emf). The solid angle is then a constant.

$\begin{array}{r}
\frac{ds}{d}\Omega = \text{const} = \frac{1}{4}p \\
d\Omega = d(\text{cos}q)df \\
f = 2pr \\
\text{cos}q = \text{cos}q_{\text{min}} + r(\text{cos}q_{\text{max}} - \text{cos}q_{\text{min}})
\end{array}$ K.8

Suppose we have a more complicated angular distribution. If, for
example, ![](media/image25.emf), then we can use the power law behavior,
Eq. K.4.

![](media/image26.emf) K.9

What about the dipole behavior which is seen in non-relativistic
radiation patterns? In that case the integral can be done, but it cannot
be analytically inverted. In fact, for most distributions, an analytic
inversion is not possible.

$\begin{array}{r}
\frac{d\sigma}{d}\Omega\  = \ \text{sin}^{2}\theta \\
\int_{}^{}{\text{sin}^{2}\theta}d(\text{cos}\theta)\  = \ \text{cos}\theta\  - \ \frac{\text{cos}^{3}\theta}{3}
\end{array}$ K.10

How do we proceed? Clearly, by performing the integral given in Eq. K.2
numerically. Assume that x is contained in the interval (x~min~,
x~max~). Assume also that within that interval there is a maximum value
of (d[P]{.underline}(x)/dx) called (d[P]{.underline}(x)/dx)~max~ We use
the "rejection" method, see Fig. K.2.

$\begin{array}{r}
\text{pick}\text{}r_{1}:\  = \ x_{\text{min}} + \ r_{1}(x_{\text{max}} - x_{\text{min}}) \\
\text{pick}\text{}r_{2}:\ \text{if}\ r_{2} < \ (d\underline{P}(x\frac{)}{d}\underline{P}(x)\text{dx})_{\text{max}}\ \text{then}\ \text{}\text{accept}\ x \\
\text{}\ \text{}\text{if}\ r_{2} > \ (d\underline{P}(x\frac{)}{d}\underline{P}(x)\text{dx})_{\text{max}}\text{}\ \text{then}\ \text{reject}\ x
\end{array}$ K.11

As can been easily seen in Fig. K.2, this procedure weights x by
(d[P]{.underline}(x)/dx). Clearly if (d[P]{.underline}(x)/dx) = 0, then
x is never accepted while if (d[P]{.underline}(x)/dx) \~
(d[P]{.underline}(x)/dx)~max~ then x is almost always accepted. For
example, in the case ![](media/image27.emf) we pick
![](media/image28.emf) uniformly between -1 and 1. If
![](media/image29.emf) is \< r~2~ then we accept that choice of
![](media/image30.emf) and continue. If not we repeat. The student
should explicitly try this procedure and verify that it works.

We can also combine methods by using the fact that joint probabilities
are multiplicative. Suppose we factor (d[P]{.underline}(x)/dx) into g(x)
and h(x) where ![](media/image31.emf).

$d\underline{P}(x\frac{)}{\text{dx}}\ \text{} = \ g(x)h(x)$ K.12

We might do this because h(x) is analytically invertable, for example.
We draw x from h(x) and then accept x if r \< g(x). Clearly the
probability to accept x is g(x) by the rejection technique.

![](media/image32.wmf)

Fig. K.2: Plot of distribution function (d[P]{.underline}(x)/dx)
indicating how one approximates it by choosing uniform width strips in x
and then accepting the strip x if (d[P]{.underline}(x)/dx)
/(d[P]{.underline}(x)/dx)~max~ \< a random number.

Several "processes" with probability P~i~ can be used in a further
generalization. For example, the total probability might be the total
cross section while the individual probabilities would refer to the
possible reactions making up the total reaction rate.

$d\underline{P}(x\frac{)}{\text{dx}}\ \text{} = \text{}\sum_{i}^{N}P_{i}\text{}g_{i}(x)\text{}h_{i}(x)$
K.13

First choose i by picking a process using a weight
P~i~/![](media/image33.emf). Then draw x from h~i~(x) either directly or
by rejection. Then accept x if r\< g~i~(x). Clearly very complex
distribution functions for x can be built up this way.