Radiation Transport Calculation by Monte Carlo Method H
Radiation Transport Calculation by Monte Carlo Method H. Hirayama, Y. Namito KEK, High Energy Accelerator Research Organization 2009. 8. 8
Monte Carlo Method • A mathematical method to solve problem using random number is called as “Monte Carlo Method” – Named by J. von Neumann S. M. Ulam • Accordingly, generation of random number is the most important technique in Monte Carlo method.
Method to generate random number • Utilize dice or roulette – Very slow • Utilize table of random number – Statistical nature is well studied. – Need to hold total number necessary as data. – Not most speedy way of generating random number. • Utilize physical phenomena such as disintegration of radioisotope. – Cumbersome in converting to numbers. – Have problem in stability and reproducibility.
Pseudorandom number • Select initial seed random number R 0 , adequately. Generate next random number by recursion equation Rn+1= f(Rn). • The rest after divided by m is treated as next random number. • The total number of integer less than m is m. Then, pseudorandom number has finite periodic length. • Good pseudorandom number, – – Can be generated QUICKLY. Have long periodic length. Have reproducibility. Have good statistical nature. • A pseudorandom number can be obtained by dividing the generated pseudorandom number by m.
Pseudorandom number by linear congruential method • Linear congruential method introduced by D. H. Lehmer is most widely used; Rn+1=mod(a. Rn+b, m) – mod(a. Rn+b, m) is a rest when a. Rn+b is divided by m. • a, b and m are positive integer. m is a maximum integer which compiler can handle. Name RANDU SLAC RAN 1 SLAC RAN 6 a 65539 b 0 m 231 69069 663608491 0 0 231
Generation of pseudorandom number by pocket calculator • Generate 10 pseudorandom number by setting R 0=3, a=5, b=0 , and m=16. • Certify that the same pseudorandom number is produced. • What is a cycle length? • Generate pseudorandom number again from a different R 0. • Excel is NOT recommended to use.
n 0 1 Rn R 0=3 Rn*5 Rn+1=mod(Rn*a, m) 15 15/16=0・・・Rest 15 R 1=15 2 3 mod(Rn*a, m): The rest when Rn*a is divided by m.
n Rn Rn*5 Rn+1=mod(Rn*a, m) 0 R 0=3 1 R 1=15 15 15/16=0・・・Rest 15 R 1=15 2 R 2=11 3 R 3=7 4 R 4=3 55 55/16=3・・・Rest 7 R 3=7 75 75/16=4・・・Rest 11 R 2=11 35 35/16=2・・・Rest 3 R 4=3
Calculation p using random number • Choose ten pairs of random numbers ( , ) from left to right at an arbitrary place in Table 1 (made using SLAC RAN 6). • Count number of pair which satisfy following condition. AREA= p /4 p= 4 x (number of pairs which satisfy the condition. )/(Total number of pairs)
How to use Table 1
Trial No ξ η R 1 2 3 0. 896 0. 759 0. 251 0. 618 0. 690 0. 094 1. 088 1. 026 0. 268 4 5 6 7 8 9 10 0. 371 0. 492 0. 789 0. 397 0. 576 0. 517 0. 909 0. 148 0. 519 0. 567 0. 179 0. 341 0. 583 0. 380 0. 399 0. 715 0. 972 0. 435 0669 0. 779 0. 985 A/10=0. 8 (A/10)*4=3. 2 R≦ 1 ○ ○ ○ ○ A=8
Other way of generating random number • Marasaglia-Zaman random number – G. Masaglia and A. Zaman, “A New Class of Random Number Generator”, Annals of Applied Probability 1(1991)462 -480. – Long period length – 2144 ~1043 – A little bit cumbersome for controlling random number – Run on any 32 -bit computer • RANLUX random number – F. James, “A Fortran implementation of the high-quality pseudorandom number generators”, Comp. Phys. Comm. 79 (1994) 111 -114. – Period length is 10 171 – By using seed number of 1 -231, independent series of random number can be produced. No overlap of them is expected. – Utilized in egs 5 as a default random number generator.
Sampling from a discrete probability distribution • Let’s write the physical processes which are independent and rebellion each other as x 1, x 2, . . . , xn. And let’s assume that their probability as p 1, p 2, . . . . , pn. (For example, photoelectric effect, Compton scattering, and pair production in interaction of photon with matter. • Let’s write a random number which distributes between 0 and 1 uniformly as . We choose an event xi when following condition is satisfied. F(xi) is called as Cumulative distribution function.
Introduction of sampling from a discrete probability distribution (1) Example) Sample interaction using a random number from a distribution of, Photoelectric effect:30%, Compton scattering : 50%, Pair production: 20%
Introduction of sampling from a disrete probability distribution (2) Random Probability of Compton scattering ”Cumulative distribution function” or “Adding up calculation”.
Example of sampling from a discrete probability distribution. • Let’s write probability of photoelectric effect, Compton scattering, and pair production as Pphoto, PCompt , Ppair. Pphoto +PCompt + Ppair =1 • Photoelectric effect is sampled when • Compton scattering is sampled when • Pair production is sampled when
Sampling from a continuous probability distribution • Let’s write a physical process emerge in a region of x and x+dx with a probability of f(x)dx [a x b]. This f(x) is called as probability density function (PDF). • Cumulative distribution function (CDF: F(x)) • Let’s write a random number which distribute between 0 and 1 uniformly as . Sampling procedure is written as, x is obtained as , l This x can be obtained by direct calculation if this equation is analytically solved. This is called as “Direct sampling method”.
F(x) 1 0 a x b x
Example of direct sampling – Calculation of flight path length • Let’s write an interaction probability of one incident particle per unit distance as St. The probability that first interaction occur between l and l+dl is, l: Flight path length l: Mean free path : Random number (1 - and are equivalent. )
Infinite medium Photon Incident condition:Energy, position, direction e 0, x 0, y 0, z 0, u 0, v 0, w 0
Sample distance l toward interaction point l=-ln( )/ Coordinate after movement x=x 0+u 0 l, y=y 0+v 0 l, z=z 0+w 0 l l Initial condition:Energy, position, and direction e 0, x 0, y 0, z 0, u 0, v 0, w 0
Sample kind of interaction Photoelectric effect:a, Compton scattering : b, Pair production :c <=a/(a+b+c): Photoelectric effect a/(a+b+c)< <=(a+b)/(a+b+c): Compton scattering >(a+b)/(a+b+c): Pair production x=x 0+u 0 l, y=y 0+v 0 l, z=z 0+w 0 l l Initial Condition: Energy, position, and direction e 0, x 0, y 0, z 0, u 0, v 0, w 0
Electron Trace of generated particle Sample energy and direction of each particle. Photon l x=x 0+u 0 l, y=y 0+v 0 l, z=z 0+w 0 l Initial Condition: Energy, position, and direction e 0, x 0, y 0, z 0, u 0, v 0, w 0
d>l: Move toward interaction point d<=l:Move by a distance of d Same material: Additional move by a distance of l-d Region boundary Different material :Sample interaction point again Interaction point x=x 0+u 0 l, y=y 0+v 0 l, z=z 0+w 0 l l d Calculate straight path length toward boundary (d). Initial Condition: Energy, position, and direction e 0, x 0, y 0, z 0, u 0, v 0, w 0
Record information Particle moves:Energy deposition Path length Region boundary Boundary crossing Photoelectric effect:Energy deposition Below cut off energy Stop tracing (Ex: Out of boudary) Calculate straight distance (d) toward boundary l d Initial Condition: Energy, position, and direction e 0, x 0, y 0, z 0, u 0, v 0, w 0
As the mean free path of elastic scattering of electron and positron is nm to m range, direct treatment of this process is unrealistic from the point of calculation efficiency. Sample distance toward interaction point Condensed History Technique Divide a distance toward major interaction into many fine steps, evaluate changing of direction and position and route distance using multiple scattering model. l=-ln( )/ l x x 11 x l 2 x l 3 l 4 l 5 x x l 6 l 7 x l 8 Charged particle looses part of its energy while moving via ionization or excitation. Energy deposition at each step is, Route distance x stopping power (d. E/dx) Initial condition :Energy, position, and direction Electron or positron e 0, x 0, y 0, z 0, u 0, v 0, w 0
Transport calculation by hand calculation • Use random numbers in Table 1 (Made by SLAC RAN 6) – Can start at arbitrary point and proceed to arbitrary direction in a series. (No jump, no change direction) when use random number. – Effective digit is 3 in calculation results
How to use random number in Table 1
Photon transport calculation 1 by hand calculation (Fig. 1) • As is shown in Fig. 1, there is a material A of 50 cm thickness. – Let’s assume that a photon of 0. 5 Me. V incident to material A vertically from left side – Let’s assume that mean free path is 20 cm. – Let’s assume that ratio of photoelectric effect and Compton scattering as 1: 1. – Let’s assume that photon energy and direction are not changed after Compton scattering. • Example 1 – Initial random number: 0. 234 -- l=-20. 0 x ln(0. 234)=29. 0 – 29. 0(cm)<50. 0(cm) – Next random number: 0. 208 (<0. 5) – Photo electric effect (Terminate) • Example 2 – – – Next random number: 0. 906 -- l=-20. 0 x ln(0. 906)=1. 97(cm)<50. 0(cm) Next random number : 0. 716 (>0. 5) – Compton scattering Next random number: 0. 996 -- l=-20. 0 x ln(0. 996)=0. 0802(cm)<50. 0 -1. 97(cm)
Single layer No. d(cm) Random number l(cm) d>l d l Random Photo. Compt. number Exp. 1 50. 0 0. 234 29. 0 * 0. 208 Exp. 2 50. 0 0. 906 1. 97 * 0. 716 * 48. 03 0. 996 0. 0802 * 0. 600 * 47. 95 0. 183 34. 0 0. 868 * 13. 95 0. 351 20. 9 * *
Photon transport calculation by hand 2 (Fig. 2) • As is shown in Fig. 2, there are material B of 20 cm thickness after material A of 30 cm thickness. – Let’s assume that a photon of 0. 5 Me. V incident onto material A vertically from left side. – The mean free path and ratio of photoelectric effect to Compton in material A is the same as previous problem. – Let’s assume that mean free path in material B to be 3 cm. – Let’s assume that Photo to Compton ration in material B to be 3: 1. – Like previous example, photon energy and direction are assumed to be un-changed after Compton scattering.
Example of calculation • Used random number in table. 1 – Initial random number : 0. 329 -- l=-20. 0 x ln(0. 329)=22. 2 – 22. 2(cm)<30. 0(cm) – Next random number 0. 612 (>0. 5) – Compton scatt. – Next random number: 0. 234 --l=-20. 0 x ln(0. 234)=29. 0 – 29. 0(cm)>30. 0 -22. 2(cm) – Move toward boundary between A and B (30. 0 cm) – Next random number : 0. 281 --l=-3. 0 x ln(0. 281)=3. 80 – 3. 80(cm)<20. 0(cm)
Medium A d(cm) Random number 30. 0 0. 329 22. 2 7. 8 0. 234 29. 0 No. Exp. 1 l(cm) d>l d l Random Photo. Compt number * 0. 612 * * Medium B d(cm) Random number l(cm) d>l d l Random Photo. Compt number 20. 00 0. 281 3. 80 * 0. 906 * 16. 20 0. 716 1. 00 * 0. 996 * 15. 20 0. 600 1. 53 * 0. 183 *
Complex but realistic calculation of photon • Geometry is Al slab of 10 cm thickness. Trace photon trajectory in following assumptions. • Incident photon energy is 0. 5 Me. V. • Photon scattering angle in Compton scattering is in the unit of 90 degree and scattering probability toward each angle is the same. This is independent of photon energy. • Scattered photon energy is calculated as,
Complex but realistic calculation of photon • Scattering azimuth angle is 0 degree and 180 degree with the same probability (1: 1). (Compton scattering occurs in X-Z plane. Left side from propagation direction is defined as 0 degree. ) • Read mfp and branching ratio of interaction from Figs. 4 and 5, respectively. • Photon cut off energy is 0. 05 Me. V.
Electron trajectory without secondary particle • Geometry is Al slab of 1 mm thickness as Fig. 7. • 1. 0 Me. V electron incident vertically from left side. • Ignore correction of route distance due to multiple scattering. • Ignore production of secondary particles, such as ray and bremsstrahlung photon. • Electron step size is 0. 01 cm regardless of electron energy.
Electron trajectory without secondary particle • Changing of electron direction due to multiple scattering is in unit of 90 degree with the same probability. This is independent of electron energy. – 0°-- 1/3, 90°-- 1/3, 180°-- 1/3 • Azimuth angle after multiple scattering is either 0 degree and 180 degree with the same probability. 0º and 180º are left side and right side of propagation. • Energy loss due to inelastic scattering is 0. 04 Me. V per 0. 01 cm. This is independent of electron energy. • Electron cut off energy is 0. 01 Me. V.
Record of modification • English version was made on 16 Jan 2013 from Japanese version of 2009. 8. 8. • Subscript number in inequality in Page 14 was reduced by 1. 2013. 1. 17* • Slides of “Reference Web page” and “Example of random walk” were deleted* 2015. 7. 27
- Slides: 46