1

Monte-Carlo Simulation Con’t: •

Compton Effect Monte-Carlo:

•

If the interaction taking place is Compton, then, we first need to sample the KleinNishina cross-section (ignoring binding energies) to find the zenith angle or energy of either the recoil electron or scattered photon. Everything else will follow.

•

Of the two simple methods to get the angle of the scattered photon (rejection sampling and numerical inversion), we will concentrate on the first. One uses the latter when the cross-section cannot be integrated analytically or when other inefficiencies happen with rejection sampling. One numerically approximates the differential cross-section for numerical inversion.

•

Rejection Sampling:

•

Consider a “number dartboard,” a box (figure below) with a random number, R1 for the solid angle (all possibilities: 0 to 4π), and another, R2, for impact parameter of hitting the electron cross-section (!) from 0 to the size of the electron squared: squaring the classical electron radius, r0. Recall that r0 = e 2 / mc 2 = 2.818 x10 −13 cm .

Lecture 14 MP 501 Kissick 2016

2 •

From Rock’s old notes:

d2

d eσ (φ ) d Ωφ

∝ r02

d1

•

The two random numbers pick a point (throw a dart at the space) in the space that is scaled as follows: Solid angle: p1 ≡ 4πR1 , and the relationship between these and the scattered photon angle, φ is:

φ = cos −1 (1 − p1 / 2π ) = cos −1 (1 − 2 R1 ) Cross-section: p2 = r02 R2 , and Choose the Klein-Nishina cross-section:

d eσ (φ ) r02  hν '   hν hν '  =  + − sin 2 (φ )    dΩφ 2  hν   hν ' hν  2

{recall: dΩφ ≡ 2π sin dφ } •

The probability of the interaction if solid angle is p1 or visa-versa is d1/d2.

Lecture 14 MP 501 Kissick 2016

3 •

If the chosen point is 1. above the d eσ (φ ) / dΩ line , then it is rejected:

r 2  hν '   hν hν '  p2 > 0  + − sin 2 (φ )    2  hν   hν ' hν  2

2. If below, then it is accepted:

r 2  hν '   hν hν '  p2 ≤ 0  + − sin 2 (φ )    2  hν   hν ' hν  2

•

Once accepted, get the angle from φ = cos −1 (1 − p1 / 2π ) = cos −1 (1 − 2 R1 ) .

•

Angles of φ that correspond to smaller d eσ (φ ) / dΩφ are chosen less often than values of φ that have larger cross-section.

•

For another example of using rejection sampling to invert a Gaussian function, see the appendix for this lecture set.

•

Also in the appendix for this lecture set is the Kahn algorithm for calculating the angles for the Compton scattered photon. The codes are provided too.

•

Rejection Sampling Efficiency – an estimate:

•

The Klein-Nishina cross-section decreases with increasing photon energy, and the efficiency is equal to the area under the curve relative to the area in the whole box. Therefore, the efficiency is going to be highest for low energy photons where the Thomson differential cross-section is applicable – for this estimate, let’s assume it’s all this: Area under curve ~ e σ 0=

8π 2 r0 3

π

Area in total box ~ 2πr02 sin θdθ = 4πr02

∫ 0

Therefore, a rough estimate of efficiency is

efficiency =

(8π / 3)r02 2 = 4πr02 3

Lecture 14 MP 501 Kissick 2016

4

•

For low energies, the pair of random numbers chosen will be sufficient only about 66% of the time, and at 1MeV, it turns out that only 20% of the guesses will be accepted, so this method is not very efficient at all! Luckily, there are more efficient methods that you will hopefully learn in other classes for the KleinNishina cross-section sampling.

•

Pair Producton Monte-Carlo:

•

At this point still, we are under the assumption that the charged particles do not transport, but instead are absorbed in the same voxel as they are created. Therefore, all we need to do is subtract off the 1.022MeV that we had to pay in order to create the two particles:

E (i, j , k ) updated ⇐ E (i, j , k ) previous + w(hν − 1.022 MeV ) •

The annihilation photons are assumed to emerge from the center of the voxel in which the pair production interaction occurred. The initial direction (α1,β 1) of one of the photons is chosen from an isotropic dictribution in the same way the other photons were. We could then assume that the other photon will go in the opposite direction: (α2=α1+π, β2=β1+π).

•

Monte-Carlo Transport of Secondary Photons:

•

If they are below the cutoff, (hν)cut, then they deposit their energy locally in that voxel, otherwise it will transport like any other photon.

•

Note that we will need to find the secondary particle direction relative to the initial particle’s direction. That is done with tedious trigonometry, but (if I did it right), the answer will be the following (derived in the appendix for this lecture):

•

Here we define the angles as the “angle cosines:” Initial direction cosines:

α 0 = cos θ 0 , k = 1 − α 02 , β 0 = sin θ 0 sin ϕ 0 , γ 0 = sin θ 0 cos ϕ 0 .

Lecture 14 MP 501 Kissick 2016

5

Redirected direction cosines, with above definitions: Note – check my work!

α = α 0 cos θ − k sin θ sin ϕ β = γ 0 cos θ + (α 0γ 0 / k ) sin θ sin ϕ − ( β 0 / k ) sin θ cos ϕ γ = β 0 cos θ + (α 0 β 0 / k ) sin θ sin ϕ − (γ 0 / k ) sin θ cos ϕ

•

Post-Simulation Processing to Determine Dose Rate:

•

When the simulation is complete, the dose in each voxel per unit decay is determined by the following:

Dose(i, j , k ) / decay =

E (i, j , k ) / N tot M (i, j , k )

Where E (i, j , k ) is the energy deposited in the i,j,k voxel and M (i, j , k ) is that voxel’s mass. Notice that the dose per decay is the same as the dose rate per unit activity but just normalizing to the decay rate.

•

Monte-Carlo Simulation of Charged Particles:

•

Now let’s drop the big assumption of ignoring the charged particle transport, and consider it separately here.

•

This time the situation is much trickier. In the photon case, we had easy, discrete interactions, but with charged particles, there are many continuous interactions along the particle’s path.

•

We will need to approximate the charged particle case with “condensed histories.”

•

In this approach, the charged particle’s path is broken up into discrete steps. The kinetic energy loss and multiple scattering will get lumped together into these steps.

Lecture 14 MP 501 Kissick 2016

6

•

The following figure illustrates the results of simulations of the dose depth curves in water for 3 types of methods – in ascending order of realism: -- If the primary electrons do not scatter or produce secondary knock-on electrons (labeled “CSDA straight-ahead”) as in primitive “Class I” algorithms. -- if primary electrons scattered but did not produce knock-ons (labeled “CSDA with ang scatt”), as in more advanced “Class I” algorithms. -- if the electrons have energy-loss straggling, are allowed to scatter and the transport of secondary electrons is included (labeled “Straggling with ang scat and secondaries”) as in “Class II” algorithms.

Lecture 14 MP 501 Kissick 2016

7

•

Notes on the above important figure: -- The inclusion of angular scattering causes the position at the end-of-range to be moved to shallower depths and be blurred. -- The inclusion of energy-loss straggling causes the depth-dose curve to extend beyond the mean range of electrons, and -- The forward transport of knock-ons reduces the dose near the surface. -- This figure illustrates the amount of dose contributed by secondary knock-on electrons (∆=1 keV) and the amount of dose to water contributed by bremsstrahlung produced in water.

•

In “Class II” electron transport codes, collisional energy losses above a threshold energy ∆ are treated as discrete interactions, and the knock-ons are then treated as a separate particles. The restricted mass stopping power can be used to determine the pathlength travelled in each step.

•

For example, if the kinetic energy lost in a step is 10%, then the transport for a charged particle with no knock-ons created would look like the following:

T = T0

T = 0.8T0

T = 0.9T0

T = Tcutoff

l l'

•

The pathlength travelled in a step, l’ would be given by:

l' = •

0.1 ⋅ T ρ (dT / ρdx) ∆

Like hiking in the woods, l < l’

Lecture 14 MP 501 Kissick 2016

8 •

The relationship between l and l’ has been quantified below:

•

As the above figure shows, it pays to keep the step small: typically, energy losses of 0.5% to 10% is used.

•

To determine if a discrete knock-on or bremsstrahlung event has occurred, use either the Möller electron-electron cross-section, or use the the Bhabha positronelectron interactions in the following for the probability of producing a knock-on electron between energies T '− dT ' / 2 and T '+ dT ' / 2 :

f (T ' ) = ρl '

T ' + dT ' / 2

NA N dσ h (dσ h / dT ')dT ' ≈ ρl ' A δT ' ∫ A T '− dT '/ 2 A dT '

Where, dσ h / dT ' is either the Möller or the Bhabha relations, and δT ' is the bin width in the following histogram of energy transfers to the knock-on electron:

Lecture 14 MP 501 Kissick 2016

9

δT ' P (T ' )

f (T ' ) =

dP(T ' ) dT '

T'

∆ •

Following the standard Monte-Carlo trick, we can form a cumulative probability distribution, and invert it as follows: T 'max

only random numbers in this range will produce a δ-ray

T'

∆ 0

0.5

rmin

1.0

r •

When the random number, r, is less than rmin, then no knock-on is produced. If it is above the minimum, then the above equation can be used to determine the kinetic energy of the knock-on electron.

•

Bremsstrahlung interactions can be handled in a similar way.

Lecture 14 MP 501 Kissick 2016

10 •

A production quality Monte-Carlo code, EGS4, operates according to the following flow-chart:

•

All together, the amount of energy lost by the charged particle is given by:

 dT 

δT =   ρl '+{∑ Tknock −on + ∑ hν brem }  ρdx  ∆ •

The amount of energy deposited in the medium is:

 dT   ρl ' E (i, j , k ) updated ⇐ E (i, j , k ) previous + w  ρdx  ∆ •

We also need to determine the amount of scattering in each step.

Lecture 14 MP 501 Kissick 2016

11 •

IF we assume that the scattering angles are Gaussian distributed, then the following picture (left), and it has a cumulative probability distribution (the “error function”), when inverted (right):

(

f θ / θ2

)

(

Pθ / θ2

)

θ

r 0

0.5

1.0

θ / θ2 •

The following shows a simulation with 6 primary electrons on a simple water phantom: Note, the electrons lose 3% of their remaining kinetic energy each step, and one can see the steps in the figure!

•

The easiest (classic for materials engineering applications) simplest Monte-Carlo code can be obtained free at: http://www.srim.org/ {a very fun site to explore, but the physics is quite serious and good!}

•

Next lectures, (after EXAM), will be concerned with practical dose calculations/measurements.

Lecture 14 MP 501 Kissick 2016