Monte Carlo radiation transport codes How do they

Monte Carlo radiation transport codes How do they work ? Michel Maire (Lapp/Annecy) 9/29/2020 introduction to Monte Carlo radiation transport codes 1

Decay in flight (1) l An unstable particle have a time of life t initial momentum p ( velocity v ) distance to travel before decay AB = d 1 = t v v A (non relativist) • Geometry : the particle is inside a box. compute distance to boundary AC = d 2 • Transport the particle s = min (d 1, d 2) • if C < B : do nothing in C, but compute the time spent in flight : Dt = AC/v if B < C : decay the particle 9/29/2020 d 2 A introduction to Monte Carlo radiation transport codes B d 1 B C 2

Geometry (1) l l The apparatus is described as an assembly of volumes made of homogene Volumes can be embeded or assembled with boulean operations A A B B C C D E • when travelling inside the apparatus, the particle must know : • where I am ? locate the current volume • where I am going ? compute distance to next boundary 9/29/2020 introduction to Monte Carlo radiation transport codes 3

Geometry (2) l remember : a computer program in blind … A A . p B B C C D E • where I am ? locate the current volume • where I am going ? compute distance to next boundary example : a point P in a box p 9/29/2020 compute intersections with 6 planes introduction to Monte Carlo radiation transport codes 4

Decay in flight (2) • The time of life, t, is a random variable with probability density function : f(t) t • It can been demonstrated in a general way that the cumulative distribution function is itself a random variable with uniform probability on [0, 1] therefore : 1 - choose r uniformly random on [0, 1] 2 - compute t = F-1(r) • For the exponential law, this gives : t = -t ln(1 - r) = -t ln(r’) 9/29/2020 introduction to Monte Carlo radiation transport codes 5

Decay in flight (3) v l When the particle travel on a distance d, one must update thedelapsed. B time A t (t – d/v) • When t = 0, one must trigger the decay of the particle • • • for instance p 0 g g (~ 99%) g e+ e- (~1%) 0 0. 99 Select a channel according the branching ratio choose r uniformly on [0, 1] p q Generate the final state • • in the rest fram of the p 0 : d. W = sinq dq df apply Lorentz transform 9/29/2020 introduction to Monte Carlo radiation transport codes 6 1

Decay in flight : comments l the generation of the whole process needs at least 4 random numbers l the decay is the simplest but general scheme of the so called analogue Monte Carlo transport simulation 9/29/2020 introduction to Monte Carlo radiation transport codes 7

Compton scattering (1) g+ e- g + e. The distance before interaction, L, is a random variable g g L e- • η(E, ρ) is the probability of Compton interaction per cm • λ(E, ρ) = η-1 is the mean free path associated to the process (Compton) Sample l = -λ ln(r) with r uniform in [0, 1] 9/29/2020 introduction to Monte Carlo radiation transport codes 8

Compton scattering (2) l l l(E, r), and l, are dependent of the material one define the number of mean free path : l l 1 l 2 l 3 • nl is independent of the material and is a random variable with distribution : f(nl) = exp(-nl) • sample nl at origin of the track : nl = -ln(r) • update elapsed nl along the track : nl (nl – dli / li) • generate Compton scattering when nl = 0 9/29/2020 introduction to Monte Carlo radiation transport codes 9

Compton scattering (3) g • Let define : E 0 q e- • the differential cross section is : • sample e with the ‘acceptation-rejection’ method remark : the generation of the whole Compton scattering process needs at least 5 random numbers 9/29/2020 introduction to Monte Carlo radiation transport codes 10

MC : acceptation-rejection method (1) l • l l • • let f(x) a probability distribution. S 1 the surface under f assume we can enclose f(x) in a box ABCD, of surface S 0 choose a point P(x 1, y 1) uniformly random within S 0 accept P only if P belong to S 1 x will be sample according to the probability distribution f the envelope can be a distribution function e(x) simple enough to be sampled with inversion technique In this case x in sampled with e(x) and rejected with f(x) 9/29/2020 y D C f(x) S 0 y 1 P S 1 A x 1 x B y e(x) f(x) y 1 P A S 1 x 1 introduction to Monte Carlo radiation transport codes x B 11

MC : acceptation-rejection method (2) l assume that we can factorize : P(x) = K f(x) g(x) • • • f(x) : probability distribution simple enough to be inverted g(x) : ‘weight’ function with values in [0, 1] K > 0 : constant to assure proper normalization of f(x) and g(x) step 1 : choose x from f(x) by inversion method step 2 : accept-reject x with g(x) • even : P(x) = K 1 f 1(x) g 1(x) + K 2 f 2(x) g 2(x) + … step 0 : choose term i with probability Ki 9/29/2020 r 2 1 g(x) r 2 P A 0 B x x r 1 F(x) 1 r 1 introduction to Monte Carlo radiation transport codes 12

g 10 Me. V in Aluminium 9/29/2020 introduction to Monte Carlo radiation transport codes 13

Simulation of charged particles (e-/+) l l Deflection of charged particles in the Coulomb field of nuclei. • small deviation; pratically no energy loss In finite thickness, particles suffer many repeated elastic Coulomb scattering • > 106 interactions / mm The cumulative effect is a net deflection from the original particle direction Individual elastic collisions are grouped together to form 1 multiple scattering single atomic deviation macroscopic view condensed history technique (class 1 algorithms) 9/29/2020 introduction to Monte Carlo radiation transport codes 14

Multiple Coulomb scattering (1) l l longitudinal displacement : z (or geometrical path length) lateral displacement : r, F true (or corrected) path length : t angular deflection : q, f 9/29/2020 introduction to Monte Carlo radiation transport codes 15

Multiple Coulomb scattering (2) 9/29/2020 introduction to Monte Carlo radiation transport codes 16

Multiple Coulomb scattering (3) 50 mm Tungsten 10 e 600 ke. V 9/29/2020 introduction to Monte Carlo radiation transport codes 17

Ionization (1) A charged particle hits a quasi-free electron (d-ray) p E 0 e- T e- cut p Tmax 9/29/2020 introduction to Monte Carlo radiation transport codes T 18

Ionization (2) accounted in the condensed history of the incident particle d. E/dx is called (restricted) stopping power or linear energy transfered explicit creation of an e- : analogue simulation 9/29/2020 introduction to Monte Carlo radiation transport codes 19

Ionization (3) ‘hard’ inelasic collisions d-rays emission e- 200 Me. V proton 200 Me. V a 200 Me. V 1 cm Aluminium 9/29/2020 introduction to Monte Carlo radiation transport codes 20

Ionization (4) straggling : DE = [DE] + fluctuations e- 16 Me. V in water ( muls off ) 9/29/2020 introduction to Monte Carlo radiation transport codes 21

Condensed history algorithms group many charged particles track segments into one single ‘condensed’ step discrete collisions grouped collisions • • • elastic scattering on nucleus • multiple Coulomb scattering soft inelastic collisions • collision stopping power (restricted) soft bremsstrahlung emission • • ‘hard’ d-ray production • energy > cut ‘hard’ bremstrahlung emission • energy > cut positron annihilation radiative stopping power (restricted) 9/29/2020 introduction to Monte Carlo radiation transport codes 22

Principle of Monte Carlo dose computation • Simulate a large number of particle histories until all primary and secondary p • Calculate and store the amount of absorbed energy of each particle in each r • The statistical accuracy of the dose is determined by the number of particle h 9/29/2020 introduction to Monte Carlo radiation transport codes 23

A non exhaustive list of MC codes (1) l l ETRAN (Berger, Seltzer; NIST 1978) EGS 4 (Nelson, Hirayama, Rogers; SLAC 1985) www. slac. stanford. edu/egs l EGS 5 (Hirayama et al; KEK-SLAC 2005) rcwww. kek. jp/research/egs 5. html l EGSnrc (Kawrakow and Rogers; NRCC 2000) www. irs. inms. nrc. ca/inms/irs. html l Penelope (Salvat et al; U. Barcelona 1999) www. nea. fr/lists/penelope. html 9/29/2020 introduction to Monte Carlo radiation transport codes 24

A non exhaustive list of MC codes (2) l Fluka (Ferrari et al; CERN-INFN 2005) www. fluka. org l Geant 3 (Brun et al; CERN 1986) www. cern. ch l Geant 4 (Apostolakis et al; CERN++ 1999) geant 4. web. cern. ch/geant 4 l MARS (James and Mokhov; FNAL) www-ap. fnal. gov/MARS l MCNPX/MCNP 5 (LANL 1990) mcnpx. lanl. gov 9/29/2020 introduction to Monte Carlo radiation transport codes 25