- Email: [email protected]
Multiple scattering of proton via stochastic differential equations
Nuclear Instruments and Methods in Physics Research B 356–357 (2015) 160–163
Contents lists available at ScienceDirect
Nuclear Instruments and Methods in Physics Research B journal homepage: www.elsevier.com/locate/nimb
Multiple scattering of proton via stochastic differential equations M.R. Kia, Houshyar Noshad ⇑ Department of Energy Engineering and Physics, Amirkabir University of Technology (Tehran Polytechnic), P.O. Box 15875-4413, Hafez Avenue, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 19 February 2015 Received in revised form 27 April 2015 Accepted 29 April 2015 Available online 15 May 2015 Keywords: Landau probability density function Moliere probability density function Stochastic differential equation Inelastic collision Elastic scattering
a b s t r a c t Multiple scattering of protons through a target is explained by a set of coupled stochastic differential equations. The motion of protons in matter is calculated by analytical random sampling from Moliere and Landau probability density functions (PDF). To satisfy the Vavilov theory, the moments for energy distribution of a 49.1 MeV proton beam in aluminum target are obtained. The skewness for the PDF of energy demonstrates that the energy distribution of protons in thin thickness becomes a Landau function, whereas, by increasing the thickness of the target it does not follow a Gaussian function completely. Afterwards, the depth-dose distributions are calculated for a 60 MeV proton beam traversing soft tissue and for a 160 MeV proton beam travelling through water. The results prove that when elastic scattering is taken into account, the Bragg-peak position is decreased, while the dose deposited in the Bragg region is increased. The results obtained in this article are benchmarked by comparison of our results with the experimental data reported in the literature. Ó 2015 Elsevier B.V. All rights reserved.
1. Introduction Development of accurate algorithms for heavy charged particle transport through matter is important due to its applications in therapy and high RBE compared to conventional radiotherapy with photons and electrons. In general, inelastic collision and elastic scattering are two principle features characterizing the passage of heavy charged particles through matter. Inelastic collision was first calculated by Bohr using the classical argument and later by Bethe-Bloch taking into account quantum mechanics [1,2]. This process for a beam of particles is subject to energy fluctuation, because of the stochastic nature of the energy loss process. Landau demonstrated that the energy probability distribution in thin thickness is highly asymmetric with a pronounced high energy tail [3]. Later Vavilov predicted that the PDF of energy has a shape between the Gaussian and Landau function based on the parameter j [4]. This theory was experimentally verified by Maccabee and Tschalar [5,6]. In multiple scattering where the energy loss is negligible, Moliere expressed the polar angle distribution as series for investigation on charged particle deviation [7]. Moreover, Bethe presented Moliere distribution function for small angles in a simple mathematical approach [8]. Moliere function can be simplified in Gaussian form if Dx=Lr > 103 , where Dx and Lr are the thickness and radiation length of absorber, respectively [9]. In this case, Noshad employed the results of the TRIM ⇑ Corresponding author. Tel.: +98 21 64545254; fax: +98 21 66495519. E-mail address: [email protected] (H. Noshad). http://dx.doi.org/10.1016/j.nimb.2015.04.078 0168-583X/Ó 2015 Elsevier B.V. All rights reserved.
computer code and studied beam deviation in matter [10]. Later Mertens et al., derived an analytical method for flux of particles with the Gaussian PDF to analyze the beam deviation [11]. One can find a semi-analytical method based on the Fermi–Eyges equation for calculation of dose distribution [11,12] developed by Hollmark et al. [13]. The results of considering inelastic scattering for dose distribution are investigated in [14–16]. In this article, we derived a set of coupled stochastic differential equations in order to investigate on inelastic collision and elastic scattering by analytical random sampling from the Moliere and Landau PDFs. Besides, the evolution for the PDF of energy as well as the depth-dose distribution are investigated in different materials. We also studied the effect of elastic scattering on depth-dose distribution via the moments of the energy PDF. The results obtained from our stochastic model are in good agreement with the experimental data. 2. Theory and methods 2.1. Inelastic collision With due attention to the stochastic nature of energy loss for passage of a proton beam through matter, the following formalism has been considered for the stochastic differential equation for evolution of energy T i ðxÞ in depth x and thickness interval Dx,
T i ðx þ DxÞ ¼ T i ðxÞ þ n k þ ln
n2mc2 bi ðxÞ2 ð1 bi ðxÞ2 ÞI
! bi ðxÞ2 C þ 1 :
ð1Þ
M.R. Kia, H. Noshad / Nuclear Instruments and Methods in Physics Research B 356–357 (2015) 160–163
161
where k is a constant. Our results are in good agreement with the experimental data with k = 0.007 for water and tissue [21]; whereas, k = 0.00648 and 0.00640 correspond to aluminum and silicon, respectively [22,23]. To obtain k, two random numbers f1 2 ½0; 0:9931 and f2 2 ½0; 1 are generated. If 0 6 f1 < 0:0141, then k is obtained as follows 14 13 12 k ¼ 106 ð0:15948f15 2 1:277f2 þ 4:672f2 10:35f2 10 9 8 7 þ 15:50f11 2 16:60f2 þ 13:10f2 7:749f2 þ 3:452f2
1:156f62 Þ þ 28745f52 5202f42 þ 6621f32 563f22 þ 30:50f2 3:4:
ð3Þ
Moreover, if 0:0141 6 f1 < 0:1136, the procedure can be followed as
k ¼ 3:052f82 þ 14:42f72 29:31f62 þ 33:77f52 24:67f42 þ 12:33f32 4:725f22 þ 2:226f2 2:
ð4Þ
If 0:1136 6 f1 < 0:6642, the algorithm can be continued as below
k ¼ 0:5071f82 þ 4:565f72 11:78f62 þ 16:26f52 13:35f42 þ 8:298f32 3:103f22 þ 3:622f2 1;
ð5Þ
and if 0:6642 6 f1 < 0:9435, then the value for k can be calculated from the following prescription Fig. 1. The energy probability density function for a 19.68 MeV proton beam in an aluminum foil at x = 0.0051 g/cm2 in (a) and x = 0.099 g/cm2 in (b).
In the relation, C ¼ 0:5772156 is the Euler–Mascheroni constant and bi ðxÞ is the velocity coefficient of the ith charged particle at step x. The parameter n and k denote the mean energy loss in thickness Dx and Landau parameter, respectively. The mean excitation potential of target material I can be obtained from theoretical formalism and experimental data [19,20]. In order to employ the Landau PDF, the value of thickness Dx is not constant and it is calculated as follows
Dx ¼ k
1:022b4 A ð1 b2 Þ0:1535qZ
;
ð2Þ
k ¼ 293:86f82 968:52f72 þ 1340:7f62 988:77f52 8 þ 419:7f42 97:852f32 þ 14:704f22 þ 3:1669f2 þ 3:
ð6Þ
Finally, if 0:9435 6 f1 6 0:9934, the value for k can be determined as follows 9 8 7 k ¼ 16956:88f10 2 73891:79f2 þ 139407:7f2 148468:2f2
þ 97854:35f62 41107:26f52 þ 10945:64f42 1745:207f32 þ 166:5348f22 þ 11:37218f1 þ 20:
ð7Þ
According to the above Eqs. (3)–(7), k is obtained via the aforementioned analytical functions. The associated error for random sampling of parameter k in the interval 3:4 6 k 6 150 is less than 1%
Fig. 2. The moments for PDF of the energy for a 49.1 MeV proton beam in an aluminum target.
162
M.R. Kia, H. Noshad / Nuclear Instruments and Methods in Physics Research B 356–357 (2015) 160–163
in comparison with that obtained from Landau PDF. In order to derive the Eqs. (3)–(7) we can indicate the following stages. First, the Landau PDF was fitted by a polynomial function for the ith interval using a least squared method. Afterwards, the integral Rk gðkÞ ¼ ki /i ðkÞ dk is analytically solved. The weight of each interval is considered by solving the aforementioned integral in the interval defined by ½kmin ; kmax . Therefore, the equation gðkÞ n ¼ 0 is numerically solved via Newton’s method for 106 values of n for each interval. At the end, the data were appropriately fitted by a polynomial function. 2.2. Elastic scattering In a homogenous target, there is no difference for energy loss in various directions, thus, the coordinate system is located on the velocity vector of the particle. When the particle traversing t, it rotates by an angle h. Thus, the total angle of the particle becomes hi ðxÞ þ h with respect to the initial beam direction. This process has stochastic nature, so that h should be obtained by random sampling from probability density function of angle in each step for any particle. By defining the parameter c as the following expression
c¼
8831Z 1:333 qt ; 2 b2 A 1:13 þ b3:76Z 2 1372
ð8Þ
from the Moliere series taking into account two terms of 1/Bn. Therefore, the value of the energy for each particle in step x can be obtained by using the one dimensional finite element interpolation method. Finally, by generating a random number f 2 ½0; 1, if f 6 0:5 then h becomes h, and otherwise, the sign of h is not changed. 3. Results and discussion To investigateon energy distribution of a proton beam, the PDF of energy can be obtained in each step using the method presented in Section 2. Fig. 1 shows the PDF for a 19.68 MeV proton beam in an aluminum foil at x = 0.0051 g/cm2 and x = 0.099 g/cm2. It is evident that the PDFs of the energy are similar to the Landau-shaped with a long tail at the early steps. These results are favorably consistent with the experimental data reported in [6]. Furthermore, as the particle passages through the target, the PDF of energy fairly follows a Gaussian function. This behavior for the PDF of energy can be studied with the moments of the probability density function. These moments for a 49.1 MeV proton beam passing through an aluminum target are depicted in Fig. 2. In accordance with Fig. 2a, the value of kurtosis after a few steps has fluctuations around 3, which is equal to the value of kurtosis for a Gaussian function. Whereas, in Fig. 2b, the value of skewness has
where the parameters q, A and Z stand for target density, effective atomic mass and atomic number, respectively. The thickness t is equal to the thickness of particle deviation and it should be chosen so that it satisfies the condition c P 112. The Moliere parameter B with the associated error less than 3% for c lies in the interval 102 and 109 can be determined as follows [26]
B ¼ 1:153 þ 2:583 ln c:
ð9Þ
According to our condition c P 112, the value of parameter B is larger than 13. One can see that, its value is increased around 20, when the particle energy is reduced to a few MeV. In this case, the terms, which are dependent on 1/Bn in the Moliere series [1], can be ignored. Consequently, the value of X i ðx þ DxÞ can be determined by analytical random sampling of h from Moliere PDF as follows
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X i ðx þ 1Þ ¼ X i ðxÞ þ t cos hi ðxÞ þ 0:4 Z 2 Bqt ln n=ðp2 b2 AÞ ;
ð10Þ
while for y direction it takes the following form
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Y i ðx þ 1Þ ¼ Y i ðxÞ þ t sin hi ðxÞ þ 0:4 Z 2 Bqt ln n=ðp2 b2 AÞ :
ð11Þ
In the relation, the parameter n is a random number between 0 and 1. The associated error due to random sampling of h in our analysis is evaluated to be less than 3% as compared to the values obtained
Fig. 3. The trajectories of protons with an initial energy of 60 MeV for tissue.
Fig. 4. The depth-dose distributions for a 60 MeV proton beam in soft tissue (a), water (b) and 160 MeV in water (c).
M.R. Kia, H. Noshad / Nuclear Instruments and Methods in Physics Research B 356–357 (2015) 160–163
163
account. In accordance with Fig. 5a, the value of the mean energy reduces when elastic scattering is considered. Moreover, according to Fig. 5b, the variance of energy increases due to elastic scattering effect. It is concluded that, when the elastic scattering is taken into account, the position of the Bragg-peak decreases, and consequently, the energy deposited in the target increases due to conservation of energy. 4. Conclusion
Fig. 5. The effect of elastic scattering on the mean energy and variance of the proton energy probability density function is shown. For this purpose, the two PDFs of the proton energy are obtained for each interval with and without considering elastic scattering. (a) The difference between the mean energies of the two PDFs is illustrated as a function of x. (b) The difference between the variances of the two PDFs is shown as a function of x.
fluctuations around 0.2, and it never becomes zero unlike the Gaussian function, and consequently, the PDF of energy is asymmetric. These asymmetry values are obtained as 0.18 for water and soft tissue, while the value of 0.12 is found for silicon target. As shown in Fig. 2c and d, the values of the variances at 2.605 g/cm2 and 2.675 g/cm2 are calculated to be 1.115 (MeV)2 and 1.625 (MeV)2, respectively. Moreover, the mean values for the PDF of energy are 11.66 MeV at 2.605 g/cm2 and 9.33 MeV at 2.675 g/cm2. The results are in better agreement with the experimental data reported in the literature [6] than those obtained from other computer codes [27]. The trajectories of protons with initial energy of 60 MeV in soft tissue are depicted in Fig. 3. One can see that this figure is in good agreement with the Monte Carlo computer code [10]. The dose distributions with cut-off energy of 0.1 MeV are illustrated in Fig. 4. As shown in Fig. 4a, the position of the Bragg-peak for soft tissue is 35.53 mm, which is the same for the data obtained from the TRIM computer code reported by Noshad [10]. The Bragg-peak positions for water with a 60 MeV proton beam is determined to be as 30.44 mm, which is obtained from Fig. 4b. This value is in favorably agreement with the data reported in the literature [16]. Furthermore, Fig. 4c is in excellent agreement with the work carried out by Paganetti for a 160 MeV proton beam [14]. To investigate the elastic scattering effects on dose distribution, the moments of the PDF of energy for inelastic collision is subtracted from the moments when elastic scattering is taken into
In this article, coupled stochastic differential equations were presented to investigate the motion of protons through matter. These equations were solved by analytical random sampling from Landau and Moliere probability density functions. The results for the passage of a proton beam for water, soft tissue, aluminum and silicon satisfy the Vavilov theory, and they are in excellent agreement with the experimental data. Computation of the skewness for the PDF of energy demonstrated that the energy distribution does not follow a Gaussian function completely. Moreover, the depth-dose distribution with and without elastic scattering were studied. Our results for a 160 MeV proton beam in water show that by considering elastic scattering, the Bragg-peak decreases by 7.3 mm, while the value of mean energy in the Bragg region increases by 0.6 MeV. We drew a conclusion that, the Bragg-peak position decreases and the dose deposited in the Bragg region increases while the elastic scattering is taken into account. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.nimb.2015.04. 078. References [1] W.R. Leo, Techniques for Nuclear and Particle Physics Experiments, SpringerVerlag, Berlin Hiedelberg New York, 1987. [2] Brain J. McParland, Nuclear Medicine Radiation Dosimetry, Springer-Verlag, Berlin Hiedelberg New York, 2010. [3] L. Landau, J. Phys. U.S.S.R. 8 (1944) 201. [4] P.V. Vavilov, Sov. Phys. JETP 5 (1957) 4. [5] H.D. Maccabee, M.R. Rajv, C.A. Tobias, Phys. Rev 165 (1968) 2. [6] C. Tschalxer, H.D. Maccabeel, Phys. Rev. A 1 (1970) 7. [7] G. Moliere, Z. Naturforsch. 2a, 133 (1947); 3a, 78 (1948). [8] H.A. Bethe, Phys. Rev 89 (1953) 6. [9] V.L. Highland, Nucl. Instr. Meth 129 (1975) 497–499. [10] H. Noshad, N. Givechi, Radiat. Meas. 39 (2005) 521–524. [11] C.J. Mertens, J.W. Wilson, S.A. Walker, J. Tweed, Adv. Space Res. 40 (2007) 1357–1367; M. Hollmark, J. Uhrdin, D. Belkic, I. Gudowska, A. Brahme, Med. Biol. 49 (2004) 3247–3265. [12] V.Y. Kuperman, Nucl. Instr. Meth. B 187 (2002) 15–19. [13] M. Hollmark, I. Gudowska, D. Belkic, A. Brahme, N. Sobolevsky, Med. Biol. 53 (2008) 3477–3491. [14] H. Paganetti, Phys. Med. Biol. 47 (2002) 747–764. [15] U. Titt, B. Bednarz, H. Paganetti, Phys. Phys. Med. Biol. 57 (2012) 6381–6393. [16] A. Rosenfeld, A.J. Wroe, I. Cornelius, M. Reinhard, D. Alexiev, IEEE Trans. Nucl. Sci. 51 (2004) 6. [19] S.M. Seltzer, M.J. Berger, Int. J. Appl. Radiat. Isot 33 (1982) 1189–1218. [20] C.A. Jayachandran, Phys. Med. Biol. 16 (1971) 4. [21] D.I. Thwaites, Radiat. Prot. Dosimetry. 13 (1985) 4. [22] D.O. Caldwell, J.R. Richardson, Phys. Rev. 94 (1954) 1. [23] J.J. Khaya, T.M. Amos, H. Bichsel, Phys. Rev. 176 (1968) 2. [26] T. Mukoyama, Y. Watanabe, Bull. Inst. Chem. Res 56 (1978) 1. [27] I. Evseev et al., INAC 2011, Belo Horizonte, MG, Brazil.
Contents lists available at ScienceDirect
Nuclear Instruments and Methods in Physics Research B journal homepage: www.elsevier.com/locate/nimb
Multiple scattering of proton via stochastic differential equations M.R. Kia, Houshyar Noshad ⇑ Department of Energy Engineering and Physics, Amirkabir University of Technology (Tehran Polytechnic), P.O. Box 15875-4413, Hafez Avenue, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 19 February 2015 Received in revised form 27 April 2015 Accepted 29 April 2015 Available online 15 May 2015 Keywords: Landau probability density function Moliere probability density function Stochastic differential equation Inelastic collision Elastic scattering
a b s t r a c t Multiple scattering of protons through a target is explained by a set of coupled stochastic differential equations. The motion of protons in matter is calculated by analytical random sampling from Moliere and Landau probability density functions (PDF). To satisfy the Vavilov theory, the moments for energy distribution of a 49.1 MeV proton beam in aluminum target are obtained. The skewness for the PDF of energy demonstrates that the energy distribution of protons in thin thickness becomes a Landau function, whereas, by increasing the thickness of the target it does not follow a Gaussian function completely. Afterwards, the depth-dose distributions are calculated for a 60 MeV proton beam traversing soft tissue and for a 160 MeV proton beam travelling through water. The results prove that when elastic scattering is taken into account, the Bragg-peak position is decreased, while the dose deposited in the Bragg region is increased. The results obtained in this article are benchmarked by comparison of our results with the experimental data reported in the literature. Ó 2015 Elsevier B.V. All rights reserved.
1. Introduction Development of accurate algorithms for heavy charged particle transport through matter is important due to its applications in therapy and high RBE compared to conventional radiotherapy with photons and electrons. In general, inelastic collision and elastic scattering are two principle features characterizing the passage of heavy charged particles through matter. Inelastic collision was first calculated by Bohr using the classical argument and later by Bethe-Bloch taking into account quantum mechanics [1,2]. This process for a beam of particles is subject to energy fluctuation, because of the stochastic nature of the energy loss process. Landau demonstrated that the energy probability distribution in thin thickness is highly asymmetric with a pronounced high energy tail [3]. Later Vavilov predicted that the PDF of energy has a shape between the Gaussian and Landau function based on the parameter j [4]. This theory was experimentally verified by Maccabee and Tschalar [5,6]. In multiple scattering where the energy loss is negligible, Moliere expressed the polar angle distribution as series for investigation on charged particle deviation [7]. Moreover, Bethe presented Moliere distribution function for small angles in a simple mathematical approach [8]. Moliere function can be simplified in Gaussian form if Dx=Lr > 103 , where Dx and Lr are the thickness and radiation length of absorber, respectively [9]. In this case, Noshad employed the results of the TRIM ⇑ Corresponding author. Tel.: +98 21 64545254; fax: +98 21 66495519. E-mail address: [email protected] (H. Noshad). http://dx.doi.org/10.1016/j.nimb.2015.04.078 0168-583X/Ó 2015 Elsevier B.V. All rights reserved.
computer code and studied beam deviation in matter [10]. Later Mertens et al., derived an analytical method for flux of particles with the Gaussian PDF to analyze the beam deviation [11]. One can find a semi-analytical method based on the Fermi–Eyges equation for calculation of dose distribution [11,12] developed by Hollmark et al. [13]. The results of considering inelastic scattering for dose distribution are investigated in [14–16]. In this article, we derived a set of coupled stochastic differential equations in order to investigate on inelastic collision and elastic scattering by analytical random sampling from the Moliere and Landau PDFs. Besides, the evolution for the PDF of energy as well as the depth-dose distribution are investigated in different materials. We also studied the effect of elastic scattering on depth-dose distribution via the moments of the energy PDF. The results obtained from our stochastic model are in good agreement with the experimental data. 2. Theory and methods 2.1. Inelastic collision With due attention to the stochastic nature of energy loss for passage of a proton beam through matter, the following formalism has been considered for the stochastic differential equation for evolution of energy T i ðxÞ in depth x and thickness interval Dx,
T i ðx þ DxÞ ¼ T i ðxÞ þ n k þ ln
n2mc2 bi ðxÞ2 ð1 bi ðxÞ2 ÞI
! bi ðxÞ2 C þ 1 :
ð1Þ
M.R. Kia, H. Noshad / Nuclear Instruments and Methods in Physics Research B 356–357 (2015) 160–163
161
where k is a constant. Our results are in good agreement with the experimental data with k = 0.007 for water and tissue [21]; whereas, k = 0.00648 and 0.00640 correspond to aluminum and silicon, respectively [22,23]. To obtain k, two random numbers f1 2 ½0; 0:9931 and f2 2 ½0; 1 are generated. If 0 6 f1 < 0:0141, then k is obtained as follows 14 13 12 k ¼ 106 ð0:15948f15 2 1:277f2 þ 4:672f2 10:35f2 10 9 8 7 þ 15:50f11 2 16:60f2 þ 13:10f2 7:749f2 þ 3:452f2
1:156f62 Þ þ 28745f52 5202f42 þ 6621f32 563f22 þ 30:50f2 3:4:
ð3Þ
Moreover, if 0:0141 6 f1 < 0:1136, the procedure can be followed as
k ¼ 3:052f82 þ 14:42f72 29:31f62 þ 33:77f52 24:67f42 þ 12:33f32 4:725f22 þ 2:226f2 2:
ð4Þ
If 0:1136 6 f1 < 0:6642, the algorithm can be continued as below
k ¼ 0:5071f82 þ 4:565f72 11:78f62 þ 16:26f52 13:35f42 þ 8:298f32 3:103f22 þ 3:622f2 1;
ð5Þ
and if 0:6642 6 f1 < 0:9435, then the value for k can be calculated from the following prescription Fig. 1. The energy probability density function for a 19.68 MeV proton beam in an aluminum foil at x = 0.0051 g/cm2 in (a) and x = 0.099 g/cm2 in (b).
In the relation, C ¼ 0:5772156 is the Euler–Mascheroni constant and bi ðxÞ is the velocity coefficient of the ith charged particle at step x. The parameter n and k denote the mean energy loss in thickness Dx and Landau parameter, respectively. The mean excitation potential of target material I can be obtained from theoretical formalism and experimental data [19,20]. In order to employ the Landau PDF, the value of thickness Dx is not constant and it is calculated as follows
Dx ¼ k
1:022b4 A ð1 b2 Þ0:1535qZ
;
ð2Þ
k ¼ 293:86f82 968:52f72 þ 1340:7f62 988:77f52 8 þ 419:7f42 97:852f32 þ 14:704f22 þ 3:1669f2 þ 3:
ð6Þ
Finally, if 0:9435 6 f1 6 0:9934, the value for k can be determined as follows 9 8 7 k ¼ 16956:88f10 2 73891:79f2 þ 139407:7f2 148468:2f2
þ 97854:35f62 41107:26f52 þ 10945:64f42 1745:207f32 þ 166:5348f22 þ 11:37218f1 þ 20:
ð7Þ
According to the above Eqs. (3)–(7), k is obtained via the aforementioned analytical functions. The associated error for random sampling of parameter k in the interval 3:4 6 k 6 150 is less than 1%
Fig. 2. The moments for PDF of the energy for a 49.1 MeV proton beam in an aluminum target.
162
M.R. Kia, H. Noshad / Nuclear Instruments and Methods in Physics Research B 356–357 (2015) 160–163
in comparison with that obtained from Landau PDF. In order to derive the Eqs. (3)–(7) we can indicate the following stages. First, the Landau PDF was fitted by a polynomial function for the ith interval using a least squared method. Afterwards, the integral Rk gðkÞ ¼ ki /i ðkÞ dk is analytically solved. The weight of each interval is considered by solving the aforementioned integral in the interval defined by ½kmin ; kmax . Therefore, the equation gðkÞ n ¼ 0 is numerically solved via Newton’s method for 106 values of n for each interval. At the end, the data were appropriately fitted by a polynomial function. 2.2. Elastic scattering In a homogenous target, there is no difference for energy loss in various directions, thus, the coordinate system is located on the velocity vector of the particle. When the particle traversing t, it rotates by an angle h. Thus, the total angle of the particle becomes hi ðxÞ þ h with respect to the initial beam direction. This process has stochastic nature, so that h should be obtained by random sampling from probability density function of angle in each step for any particle. By defining the parameter c as the following expression
c¼
8831Z 1:333 qt ; 2 b2 A 1:13 þ b3:76Z 2 1372
ð8Þ
from the Moliere series taking into account two terms of 1/Bn. Therefore, the value of the energy for each particle in step x can be obtained by using the one dimensional finite element interpolation method. Finally, by generating a random number f 2 ½0; 1, if f 6 0:5 then h becomes h, and otherwise, the sign of h is not changed. 3. Results and discussion To investigateon energy distribution of a proton beam, the PDF of energy can be obtained in each step using the method presented in Section 2. Fig. 1 shows the PDF for a 19.68 MeV proton beam in an aluminum foil at x = 0.0051 g/cm2 and x = 0.099 g/cm2. It is evident that the PDFs of the energy are similar to the Landau-shaped with a long tail at the early steps. These results are favorably consistent with the experimental data reported in [6]. Furthermore, as the particle passages through the target, the PDF of energy fairly follows a Gaussian function. This behavior for the PDF of energy can be studied with the moments of the probability density function. These moments for a 49.1 MeV proton beam passing through an aluminum target are depicted in Fig. 2. In accordance with Fig. 2a, the value of kurtosis after a few steps has fluctuations around 3, which is equal to the value of kurtosis for a Gaussian function. Whereas, in Fig. 2b, the value of skewness has
where the parameters q, A and Z stand for target density, effective atomic mass and atomic number, respectively. The thickness t is equal to the thickness of particle deviation and it should be chosen so that it satisfies the condition c P 112. The Moliere parameter B with the associated error less than 3% for c lies in the interval 102 and 109 can be determined as follows [26]
B ¼ 1:153 þ 2:583 ln c:
ð9Þ
According to our condition c P 112, the value of parameter B is larger than 13. One can see that, its value is increased around 20, when the particle energy is reduced to a few MeV. In this case, the terms, which are dependent on 1/Bn in the Moliere series [1], can be ignored. Consequently, the value of X i ðx þ DxÞ can be determined by analytical random sampling of h from Moliere PDF as follows
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X i ðx þ 1Þ ¼ X i ðxÞ þ t cos hi ðxÞ þ 0:4 Z 2 Bqt ln n=ðp2 b2 AÞ ;
ð10Þ
while for y direction it takes the following form
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Y i ðx þ 1Þ ¼ Y i ðxÞ þ t sin hi ðxÞ þ 0:4 Z 2 Bqt ln n=ðp2 b2 AÞ :
ð11Þ
In the relation, the parameter n is a random number between 0 and 1. The associated error due to random sampling of h in our analysis is evaluated to be less than 3% as compared to the values obtained
Fig. 3. The trajectories of protons with an initial energy of 60 MeV for tissue.
Fig. 4. The depth-dose distributions for a 60 MeV proton beam in soft tissue (a), water (b) and 160 MeV in water (c).
M.R. Kia, H. Noshad / Nuclear Instruments and Methods in Physics Research B 356–357 (2015) 160–163
163
account. In accordance with Fig. 5a, the value of the mean energy reduces when elastic scattering is considered. Moreover, according to Fig. 5b, the variance of energy increases due to elastic scattering effect. It is concluded that, when the elastic scattering is taken into account, the position of the Bragg-peak decreases, and consequently, the energy deposited in the target increases due to conservation of energy. 4. Conclusion
Fig. 5. The effect of elastic scattering on the mean energy and variance of the proton energy probability density function is shown. For this purpose, the two PDFs of the proton energy are obtained for each interval with and without considering elastic scattering. (a) The difference between the mean energies of the two PDFs is illustrated as a function of x. (b) The difference between the variances of the two PDFs is shown as a function of x.
fluctuations around 0.2, and it never becomes zero unlike the Gaussian function, and consequently, the PDF of energy is asymmetric. These asymmetry values are obtained as 0.18 for water and soft tissue, while the value of 0.12 is found for silicon target. As shown in Fig. 2c and d, the values of the variances at 2.605 g/cm2 and 2.675 g/cm2 are calculated to be 1.115 (MeV)2 and 1.625 (MeV)2, respectively. Moreover, the mean values for the PDF of energy are 11.66 MeV at 2.605 g/cm2 and 9.33 MeV at 2.675 g/cm2. The results are in better agreement with the experimental data reported in the literature [6] than those obtained from other computer codes [27]. The trajectories of protons with initial energy of 60 MeV in soft tissue are depicted in Fig. 3. One can see that this figure is in good agreement with the Monte Carlo computer code [10]. The dose distributions with cut-off energy of 0.1 MeV are illustrated in Fig. 4. As shown in Fig. 4a, the position of the Bragg-peak for soft tissue is 35.53 mm, which is the same for the data obtained from the TRIM computer code reported by Noshad [10]. The Bragg-peak positions for water with a 60 MeV proton beam is determined to be as 30.44 mm, which is obtained from Fig. 4b. This value is in favorably agreement with the data reported in the literature [16]. Furthermore, Fig. 4c is in excellent agreement with the work carried out by Paganetti for a 160 MeV proton beam [14]. To investigate the elastic scattering effects on dose distribution, the moments of the PDF of energy for inelastic collision is subtracted from the moments when elastic scattering is taken into
In this article, coupled stochastic differential equations were presented to investigate the motion of protons through matter. These equations were solved by analytical random sampling from Landau and Moliere probability density functions. The results for the passage of a proton beam for water, soft tissue, aluminum and silicon satisfy the Vavilov theory, and they are in excellent agreement with the experimental data. Computation of the skewness for the PDF of energy demonstrated that the energy distribution does not follow a Gaussian function completely. Moreover, the depth-dose distribution with and without elastic scattering were studied. Our results for a 160 MeV proton beam in water show that by considering elastic scattering, the Bragg-peak decreases by 7.3 mm, while the value of mean energy in the Bragg region increases by 0.6 MeV. We drew a conclusion that, the Bragg-peak position decreases and the dose deposited in the Bragg region increases while the elastic scattering is taken into account. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.nimb.2015.04. 078. References [1] W.R. Leo, Techniques for Nuclear and Particle Physics Experiments, SpringerVerlag, Berlin Hiedelberg New York, 1987. [2] Brain J. McParland, Nuclear Medicine Radiation Dosimetry, Springer-Verlag, Berlin Hiedelberg New York, 2010. [3] L. Landau, J. Phys. U.S.S.R. 8 (1944) 201. [4] P.V. Vavilov, Sov. Phys. JETP 5 (1957) 4. [5] H.D. Maccabee, M.R. Rajv, C.A. Tobias, Phys. Rev 165 (1968) 2. [6] C. Tschalxer, H.D. Maccabeel, Phys. Rev. A 1 (1970) 7. [7] G. Moliere, Z. Naturforsch. 2a, 133 (1947); 3a, 78 (1948). [8] H.A. Bethe, Phys. Rev 89 (1953) 6. [9] V.L. Highland, Nucl. Instr. Meth 129 (1975) 497–499. [10] H. Noshad, N. Givechi, Radiat. Meas. 39 (2005) 521–524. [11] C.J. Mertens, J.W. Wilson, S.A. Walker, J. Tweed, Adv. Space Res. 40 (2007) 1357–1367; M. Hollmark, J. Uhrdin, D. Belkic, I. Gudowska, A. Brahme, Med. Biol. 49 (2004) 3247–3265. [12] V.Y. Kuperman, Nucl. Instr. Meth. B 187 (2002) 15–19. [13] M. Hollmark, I. Gudowska, D. Belkic, A. Brahme, N. Sobolevsky, Med. Biol. 53 (2008) 3477–3491. [14] H. Paganetti, Phys. Med. Biol. 47 (2002) 747–764. [15] U. Titt, B. Bednarz, H. Paganetti, Phys. Phys. Med. Biol. 57 (2012) 6381–6393. [16] A. Rosenfeld, A.J. Wroe, I. Cornelius, M. Reinhard, D. Alexiev, IEEE Trans. Nucl. Sci. 51 (2004) 6. [19] S.M. Seltzer, M.J. Berger, Int. J. Appl. Radiat. Isot 33 (1982) 1189–1218. [20] C.A. Jayachandran, Phys. Med. Biol. 16 (1971) 4. [21] D.I. Thwaites, Radiat. Prot. Dosimetry. 13 (1985) 4. [22] D.O. Caldwell, J.R. Richardson, Phys. Rev. 94 (1954) 1. [23] J.J. Khaya, T.M. Amos, H. Bichsel, Phys. Rev. 176 (1968) 2. [26] T. Mukoyama, Y. Watanabe, Bull. Inst. Chem. Res 56 (1978) 1. [27] I. Evseev et al., INAC 2011, Belo Horizonte, MG, Brazil.