An algorithm for Monte Carlo simulation of bremsstrahlung emission by electrons

An algorithm for Monte Carlo simulation of bremsstrahlung emission by electrons 1 2 3 first-author 2 3 2 3 corresp Corresponding Author, Email address: peng.he@fds.org.cn Corresponding Author, Email address: peng.he@fds.org.cn peng.he@fds.org.cn School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230027, China. School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230027, China. Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academy of Sciences, Hefei 230031, China. Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academy of Sciences, Hefei 230031, China. Collaborative Innovation Center of Radiation Medicine of Jiangsu Higher Education Institutions, Suzhou 215006, China. Collaborative Innovation Center of Radiation Medicine of Jiangsu Higher Education Institutions, Suzhou 215006, China. An algorithm for Monte Carlo simulation of bremsstrahlung emission by electrons based on the framework of SuperMC is presented in this paper. With efficient and accurate methods to sample the angular distribution and energy of bremsstrahlung photons. The photon energy is sampled according to scaled energy-loss differential cross sections tabulated by Seltzer and Berger. A novel hybrid model for photon angular distribution by low and high energy incident electrons are developed. The model uses Tsai’s full form of angular distribution function with atomic form factors for high energy incident electrons. For electrons of < 500 keV, a simple efficient and accurate analytical distribution function is developed, using adjustable parameters determined from the fitting of numerical values of the shape functions tabulated by Kissel et al . The efficiency of sampling photon energy is 80%. Our angular sampling algorithm for high energy electron bremsstrahlung based on Tsai distribution function is very efficient (sampling efficiency~70%) in the useful photon energy range. Bremsstrahlung Monte Carlo SuperMC Differential cross section Angular distribution 1 Introduction 1 As a result of the acceleration or especially the deceleration of charged particles caused by the electrostatic field of atoms, electrons (or positrons) emit bremsstrahlung. Bremsstrahlung emission is of basic interest in many fields. Examples of applications where the bremsstrahlung interaction plays a vital role are, X-ray generators, radiation shielding, medical physics, and electron probe microanalysis [ 1 1 ]. In medical physics, the bremsstrahlung phenomenon is very important in both diagnostics radiology and radiotherapy because the fundamental characteristics of X-ray machine is described by the spectral distribution of photon beam originating from the target. Bremsstrahlung photons have much higher penetration than other charged particles and thus have large contribution in dose in deep regions of the irradiated samples, therefore an accurate Monte Carlo simulation of bremsstrahlung is required in a general purpose Monte Carlo code such as SuperMC [ 2 2 ]. SuperMC is a general purpose, easy-to-use Monte Carlo simulation program, developed by FDS Team [ 2 2 - 5 5 2-5 ], for modeling and simulating nuclear facilities. The released version of SuperMC can perform the simulation of neutrons, photons and coupled neutron and photon transport while Monte Carlo simulation of electron transport has been under development and will be released soon, this research work is part of the development of electron transport. The main technical features of SuperMC include hybrid MC-deterministic methods and the adoption of advanced information technologies, while the main usability features include automatic modeling of geometry and physics, visualization and virtual simulation and cloud computing services. SuperMC has been applied to fusion reactor FDS I[ 6 6 ], FDS II[ 7 7 ], FDS III[ 8 8 ], ITER [ 9 9 - 10 10 9-10 ] and other reactor studies[ 11 11 ] and is also designed for MC-based dose calculation engine of ARTS (Accurate/Advanced Radiotherapy System) [ 12 12 ]. Bremsstrahlung process is described by the differential cross section (DCS) in energy and direction of emitted photons and the direction of outgoing electrons. Atomic DCSs of bremsstrahlung emission by electrons are derived from conventional quantum electrodynamics methods and comprehensive review on the DCSs is described by Koch and Motz [ 13 13 ] and Tsai [ 14 14 ]. Efforts have been made to use accurate methods for sampling the photon energy and, especially, the angular distribution of bremsstrahlung photons, but still, there is need to use simple, efficient and accurate methods to sample energy and angular distribution. The 2BS formula of Koch and Motz was implemented by Alex F. Bielajew et al . [ 15 15 ] in EGS4 [ 16 16 ] for the angular distribution of bremsstrahlung photons. Acosta et al . [ 17 17 ] developed a parameterization with adjustable parameters determined from the fitting of shape functions provided by Kissel et al . [ 18 18 ] based on the calculations of partial wave methods by Tseng et al . [ 19 19 ]. Geant4 [ 20 20 ] uses a parameterization of Tsai [ 14 14 ] DDCS for sampling photon angular distribution, while Rodrigues et al . [ 21 21 ] implemented 2BN formula of Koch and Motz, in Geant4 as more accurate option for low energy region. EGSnrc [ 22 22 ] incorporated modified version of 2BS formula which converges to 2BN for the treatment of low energy region. In this article, an algorithm is described for Monte Carlo simulation of electron bremsstrahlung with kinetic energy range of 1 keV–10 GeV and for Z = 1–100. We also describe the efficient and accurate methods by which the angular distribution and energy of bremsstrahlung photons are sampled. Photon energy is sampled according to most reliable Seltzer and Berger [ 23 23 ] scaled DCSs and for angular distribution of bremsstrahlung photons, a novel hybrid model is developed. In Section 2, DCSs used for sampling energy and angular distribution of emitted photons, and corresponding random sampling algorithms, are presented. Results of the newly developed MC simulation algorithm are presented in Section 3. 2 Monte Carlo simulation of electron bremsstrahlung 1 Bremsstrahlung process is described by the atomic differential cross section, in the direction θ and energy K of the emitted photons and the direction of the outgoing electrons. We consider the process in which electron of energy E (kinetic energy T ) is accelerated or decelerated by the screened coulomb field of the target atom and emit bremsstrahlung photons of energy K in the interval [0, E ] in the direction of polar angle θ . Angular deflections of the incident electron due to bremsstrahlung are already accounted for by electric scattering in SuperMC and assumed that the direction of electron is not altered by the emission of bremsstrahlung photon. In conventional MC simulations, bremsstrahlung emission is modeled using analytical DCSs obtained from simple approximations and this was done to minimize the computer memory requirements. With the advancement in computer technology, however, it is now convenient to use the combination of numerical DCSs and analytical formulae, which are more accurate and reliable at present. In this study we used the models and numerical data which are considered as most reliable at present and their form suit for Monte Carlo simulation. 2.1. Sampling photon energy 2 In this section we describe the algorithm to sample photon energy, K , from scaled energy-loss total DCSs (electron-electron, electron-nucleus), χ ( Z , T , w = K / T ), tabulated by Seltzer and Berger [ 23 23 ] which are based on the combination of numerical partial wave results and high energy theory. <math display="block" id="M1"><mrow><mfrac><mrow><msup><mi>β</mi><mn>2</mn></msup></mrow><mrow><msup><mi>Z</mi><mn>2</mn></msup></mrow></mfrac><mi>k</mi><mfrac><mrow><mtext>d</mtext><mi>σ</mi></mrow><mrow><mtext>d</mtext><mi>K</mi></mrow></mfrac><mo> </mo><mo>=</mo><mi>χ</mi><mrow><mo>(</mo><mrow><mi>Z</mi><mo>,</mo><mi>T</mi><mo>,</mo><mi>w</mi><mo>=</mo><mi>K</mi><mo>/</mo><mi>T</mi></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where, w is reduced photon energy, <math id="M2"><mrow><mo> </mo><mi>β</mi></mrow></math> is incident electron velocity in units of the velocity of light and T is the kinetic energy of incident electron. The set includes results for 57 values of electron kinetic energies from 1 keV to 10 GeV versus 30 values of reduced photon energy points span on the interval [0, 1], and for atomic numbers Z = 1–100. The values of electron kinetic energies are suitably spaced (spacing is nearly logarithmic) to allow accurate linear interpolation in ln T . The continuous scaled DCSs are, therefore, obtained by cubic spline interpolation in w = K / T and linear interpolation in ln T . For sampling photon energy, a simple and efficient composite sampling procedure is implemented. The distribution function <math id="M3"><mrow><mi>p</mi><mrow><mo>(</mo><mi>w</mi><mo>)</mo></mrow></mrow></math> for the reduced photon energy, w = K / T , and for given Z and T over the interval [ K c / T , 1] is given as, <math display="block" id="M4"><mrow><mi>p</mi><mrow><mo>(</mo><mi>w</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mi>w</mi></mfrac><mo>×</mo><mi>χ</mi><mrow><mo>(</mo><mrow><mi>Z</mi><mo>,</mo><mi>T</mi><mo>,</mo><mi>w</mi></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where K c is threshold photon energy. As can be seen from Fig. 1 Fig. 1 , χ ( Z , T , w ) is relatively flat over the interval [ K c / T , 1] and can be used for rejection sampling while direct (inverse-transform) sampling method is used to sample w from the distribution <math id="M5"><mrow><msup><mi>w</mi><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> in the interval [ K c / T , 1]. The algorithm is then i. Calculate χ max (find maximum value of χ in K / T grid for given T and Z ). ii. Generate two random numbers r 1 and r 2 over the interval [0, 1]. iii. Sample w from the distribution w −1 as w = exp[ln( K c r 1 / T )] iv. If r 2 ≤ χ ( Z , T , w )/ χ max , accept w ; else, go to Step ii. Efficiency and accuracy of this sampling algorithm is discussed in Section 3.1. A number of comparisons between measured and calculated scaled differential cross section in photon energy have been done [ 23 23 ]. These comparisons show that Seltzer and Berger’ scaled differential cross sections are the most reliable representation of bremsstrahlung at present. 2.2. Hybrid model for photon angular distribution 2 A new hybrid model for sampling the direction of bremsstrahlung photons has been developed which characterizes the explicit treatment of photon angular distribution for low and high energy incident electrons. High energy analytical formula of DDCS (differential in energy and angle of bremsstrahlung photon) described by Tsai (1974) [ 14 14 ] based on high energy theory of Davies, Bethe, Maximon and Olsen (DBMO) [ 24 24 , 25 25 ] with emphasis on form factor screening corrections is used to sample the direction of emitted photons. However at low energies, 2BS formula of Koch and Motz and Tsai’s DDCS give rise to deviations which could reach 50 Degrees on the most probable angle [ 21 21 ]. This is of particular importance when simulating bremsstrahlung emission in very thin targets and addressing topics like medical imaging which works with X rays of about a few keV to a few hundred keV. Therefore an accurate treatment of the angular distribution of the bremsstrahlung photons, needs to be done for low energy incident electrons. Thus a new, accurate and fairly simple analytical distribution function has been developed for photon angular sampling with adjustable parameters determined from the fitting of the numerical values of shape functions tabulated by Kissel et al . [ 18 18 ]. 2.2.1 High energy bremsstrahlung angular distribution 3 In this study, for high energy electrons, the angular distribution of emitted photons has been sampled accurately, according to the “exact form” of Tsai’s double differential cross section (DDCS) which includes atomic form factors. Thus the sampling is fully consistent with theory. It should be noted that the parameterization of the angular distribution instead of “exact form” may produce uncertainties and the accuracy of the parameterization needs to check against, many, possible electron and photon energies and materials. The double differential cross section, energy angle distribution of electron bremsstrahlung, described by Tsai [ 14 14 ] is <math display="block" id="M6"><mrow><mfrac><mrow><mtext>d</mtext><mi>σ</mi></mrow><mrow><mtext>d</mtext><mi>Ω</mi><mtext>d</mtext><mi>K</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><mn>2</mn><msup><mi>α</mi><mn>2</mn></msup></mrow><mrow><mi>π</mi><msubsup><mi>m</mi><mi>e</mi><mn>4</mn></msubsup></mrow></mfrac><mfrac><mrow><msup><mi>E</mi><mn>2</mn></msup></mrow><mi>K</mi></mfrac><mrow><mo>{</mo><mrow><mrow><mo>[</mo><mrow><mfrac><mrow><mn>2</mn><mi>y</mi><mo>−</mo><mn>2</mn></mrow><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><mn>12</mn><msup><mi>u</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mn>4</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow><mrow><mo>(</mo><mrow><msup><mi>Z</mi><mn>2</mn></msup><mo>+</mo><mi>Z</mi></mrow><mo>)</mo></mrow><mo>+</mo><mrow><mo>[</mo><mrow><mfrac><mrow><mn>2</mn><mo>−</mo><mn>2</mn><mi>y</mi><mo>+</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></mfrac><mo>−</mo><mfrac><mrow><mn>4</mn><msup><mi>u</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mn>4</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow><mrow><mo>(</mo><mrow><mi>X</mi><mo>−</mo><mn>2</mn><msup><mi>Z</mi><mn>2</mn></msup><msub><mi>f</mi><mi>c</mi></msub><mrow><mo>(</mo><mrow><msup><mrow><mrow><mo>(</mo><mrow><mi>α</mi><mi>Z</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow><mo>}</mo></mrow><mo>,</mo></mrow></math> where y = K / E and u = θE / m e ( m e = 0.511 MeV), α=1/137, Z is atomic number, and f c is coulomb correction. For the calculation of the function X, Tsai [ 14 14 ] provided simple expressions: <math display="block" id="M7"><mrow><mi>X</mi><mo>=</mo><msub><mi>X</mi><mrow><mtext>el</mtext></mrow></msub><mo>+</mo><msub><mi>X</mi><mrow><mtext>inel</mtext></mrow></msub><mo>=</mo><msup><mi>Z</mi><mn>2</mn></msup><mrow><mo>[</mo><mrow><mi>l</mi><mi>n</mi><mfrac><mrow><msup><mi>a</mi><mn>2</mn></msup><msubsup><mi>m</mi><mi>e</mi><mn>2</mn></msubsup><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mrow><msup><mi>a</mi><mn>2</mn></msup><msub><mi>t</mi><mrow><mtext>min</mtext></mrow></msub><mo>+</mo><mn>1</mn></mrow></mfrac><mo>−</mo><mn>1</mn></mrow><mo>]</mo></mrow><mo>+</mo><mi>Z</mi><mrow><mo>[</mo><mrow><mi>l</mi><mi>n</mi><mfrac><mrow><msup><mi>a</mi><mrow><mo>'</mo><mn>2</mn></mrow></msup><msubsup><mi>m</mi><mi>e</mi><mn>2</mn></msubsup><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mrow><msup><mi>a</mi><mrow><mo>'</mo><mn>2</mn></mrow></msup><msub><mi>t</mi><mrow><mtext>min</mtext></mrow></msub><mo>+</mo><mn>1</mn></mrow></mfrac><mo>−</mo><mn>1</mn></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> where a =184.15 e −1/2 Z −1/3 / m e , and a '=1194 e −1/2 Z −2/3 / m e are atomic parameters. The minimum momentum transfer used to calculate X is defined as t min = { ym e 2 (1+ u 2 )/[2E(1− y )]} 2 . To sample the direction of bremsstrahlung photon from Tsai’s DDCS, the angular distribution function is defined as <math display="block" id="M8"><mrow><mi>f</mi><mrow><mo>(</mo><mi>u</mi><mo>)</mo></mrow><mo>=</mo><mi>h</mi><mrow><mo>(</mo><mi>u</mi><mo>)</mo></mrow><mo>×</mo><mi>R</mi><mrow><mo>(</mo><mi>u</mi><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math display="block" id="M9"><mrow><mi>R</mi><mrow><mo>(</mo><mi>u</mi><mo>)</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><mfrac><mrow><mn>2</mn><mi>y</mi><mo>−</mo><mn>2</mn></mrow><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mrow/></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><mn>12</mn><msup><mi>u</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mn>3</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow><mrow><mo>(</mo><mrow><msup><mi>Z</mi><mn>2</mn></msup><mo>+</mo><mi>Z</mi></mrow><mo>)</mo></mrow><mo>+</mo><mrow><mo>[</mo><mrow><mfrac><mrow><mn>2</mn><mo>−</mo><mn>2</mn><mi>y</mi><mo>+</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mrow/></msup></mrow></mfrac><mo>−</mo><mfrac><mrow><mn>4</mn><msup><mi>u</mi><mn>2</mn></msup><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>y</mi></mrow><mo>)</mo></mrow></mrow><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mn>3</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow><mrow><mo>(</mo><mrow><mi>X</mi><mo>−</mo><mn>2</mn><msup><mi>Z</mi><mn>2</mn></msup><mi>f</mi><mrow><mo>(</mo><mrow><msup><mrow><mrow><mo>(</mo><mrow><mi>α</mi><mi>Z</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M10"><mrow><mi>h</mi><mrow><mo>(</mo><mi>u</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> As can be seen from Fig. 2 Fig. 2 , R ( u ) can be employed as rejection function and therefore a composite approach of direct sampling and rejection sampling is used to sample u from Eq. (5). Here u is directly sampled from Eq. (7) and then R ( u ) is used as rejection function. In order to employ rejection sampling, the location of the maximum, <math id="M11"><mrow><msub><mi>u</mi><mrow><mi>m</mi><mi>a</mi><mi>x</mi></mrow></msub></mrow></math> , of rejection function must be known so that it allows the efficient determination of the rejection function normalization, <math id="M12"><mrow><msup><mrow><mrow><mo>[</mo><mrow><mi>R</mi><mrow><mo>(</mo><mrow><msub><mi>u</mi><mrow><mi>m</mi><mi>a</mi><mi>x</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> . Unfortunately, finding rejection function normalization in complete generality proved to be very difficult. However some observations have been made to estimate the location of the maximum of the rejection function in an efficient way. We calculated <math id="M13"><mrow><msub><mi>u</mi><mrow><mi>m</mi><mi>a</mi><mi>x</mi></mrow></msub></mrow></math> using numerical methods for Z = 13, 47 and 82 and for electron kinetic energies of 0.01, 0.1, 0.3, 0.5, 0.7, 1, 10, 100, 500 and 1000 MeV, and for 38 values of reduced photon energies equally spaced between 0.05 to 1 (Where 0.05 is set as threshold reduced photon energy) and it has been observed that the location of the maximum of the rejection function exist in a very small interval of <math id="M14"><mrow><mi>u</mi><mo>∈</mo><mrow><mo>[</mo><mrow><mn>0.7</mn><mo>,</mo><mo> </mo><mn>2.4</mn></mrow><mo>]</mo></mrow></mrow></math> as shown in Fig. 3. Fig. 3. From these calculations we found that the true maxima is 3% higher (at the most) than estimated maxima at one of the values of <math id="M15"><mrow><mi>u</mi><mo>=</mo><mn>1.0</mn><mo>,</mo><mo> </mo><mi>u</mi><mo>=</mo><mn>1.5</mn></mrow></math> or <math id="M16"><mrow><mi>u</mi><mo>=</mo><mn>2.5</mn></mrow></math> , therefore the maxima of the rejection function can be calculated as <math display="block" id="M17"><mrow><msub><mi>R</mi><mrow><mtext>max</mtext></mrow></msub><mo>=</mo><mn>1.03</mn><mo>×</mo><mi>m</mi><mi>a</mi><mi>x</mi><mrow><mo>[</mo><mrow><mi>R</mi><mrow><mo>(</mo><mrow><mn>1.0</mn></mrow><mo>)</mo></mrow><mo>,</mo><mi>R</mi><mrow><mo>(</mo><mrow><mn>1.5</mn></mrow><mo>)</mo></mrow><mo>,</mo><mi>R</mi><mrow><mo>(</mo><mrow><mn>2.5</mn></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>.</mo></mrow></math> Since true maxima is 3% (at the most) higher than the estimated maxima, therefore, this small error can be corrected by introducing an additional multiplicative factor to the estimated maxima, which is 1.03 (3%) for entire energy range. The algorithm is then i. Calculate R max from Eq. (8) for given Z , E and K . ii. Generate two random numbers r 1 and r 2 over the interval [0, 1]. iii. Calculate normalization constant ( N u ) of h ( u ) over [0,π E / m e ] as N u = 1/arctan(π E / m e ) iv. Sample u from h ( u ) as u = tan( r 1 / N u ) v. If r 2 ≤ R ( u )/ R max , accept u ; else, go to step ii . The proposed sampling algorithm requires several logarithmic calculations in step (i) and in step (v) which makes it slightly slow, but holds high sampling efficiency. 2.2.2. Low energy bremsstrahlung angular distribution 3 Numerical values of the shape functions based on the partial wave methods [ 19 19 ], have been published by Kissel et al . [ 18 18 ] for 144 benchmark values of electron kinetic energy (1–500 keV), atomic number ( Z =1–92) and reduced photon energy. For electron kinetic energies of < 500 keV, the angular distribution is sampled from an analytical distribution function with adjustable parameters determined from the fitting of Kissel et al . [ 18 18 ] data. Kissel and coworkers have provided the analytical form of the shape functions in terms of Legendre polynomial, but their analytical expression does not suit for the sampling of polar angle θ . Therefore, in this study, a fairly simple parameterization of the partial wave double differential cross section with adjustable parameters is proposed as: <math display="block" id="M18"><mrow><mi>ρ</mi><mrow><mo>(</mo><mi>q</mi><mo>)</mo></mrow><mo>=</mo><mi>A</mi><msup><mi>e</mi><mrow><mn>1</mn><mo>+</mo><mi>q</mi><mo>−</mo><msup><mi>e</mi><mi>q</mi></msup></mrow></msup><mo>,</mo><mi>q</mi><mo>=</mo><mfrac><mrow><mi>B</mi><mo>−</mo><mi>θ</mi></mrow><mi>C</mi></mfrac><mo>,</mo></mrow></math> where θ = 0–π is the polar angle of bremsstrahlung photon; and A , B and C are adjustable parameters, which are determined by fitting to the numerical values of shape functions for 144 benchmark cases and stored in tables for rapid retrieves of their values. Some differences were found between the data and fits for lower photon energies, but still the fits are accurate enough for sampling angular distribution of bremsstrahlung photons. Comparisons between analytical function determined by Eq. (9) with the original shape functions tabulated by Kissel et al . [ 18 18 ] are shown in Fig. 4. Fig. 4. The best side of the proposed analytical form is that ρ has analytical integral and the corresponding cumulative probability distribution function has analytical inverse, therefore direct sampling procedure (inverse transform method) was applied to sample θ. The normalization constant N A is constant chosen to normalize the integral of ρ ( θ ) in the overall range of θ to unity. <math display="block" id="M19"><mrow><msub><mi>N</mi><mi>A</mi></msub><mo>=</mo><mi>A</mi><mo>×</mo><mi>C</mi><mo>×</mo><mi>exp</mi><mrow><mo>(</mo><mrow><msup><mi>e</mi><mrow><mo>−</mo><msup><mi>e</mi><mrow><mfrac><mrow><mi>B</mi><mo>+</mo><mi>π</mi></mrow><mi>C</mi></mfrac></mrow></msup></mrow></msup><mo>−</mo><msup><mi>e</mi><mrow><mo>−</mo><msup><mi>e</mi><mrow><mfrac><mi>B</mi><mi>C</mi></mfrac></mrow></msup></mrow></msup></mrow><mo>)</mo></mrow></mrow></math> Since the calculation of N A involves five exponential calculations therefore in order to make algorithm efficient, we calculated N A for 144 benchmarked values and stored in tables along with adjustable parameters for their rapid retrieves, for particular T , Z , and K / T . The polar angle is then sampled as <math display="block" id="M20"><mrow><mi>θ</mi><mo>=</mo><mi>B</mi><mo>−</mo><mi>C</mi><mo>×</mo><mi>ln</mi><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>ln</mi><mrow><mo>(</mo><mrow><mfrac><mi>r</mi><mrow><msub><mi>N</mi><mi>A</mi></msub><mo>×</mo><mi>A</mi><mo>×</mo><mi>C</mi></mrow></mfrac><mo>+</mo><msup><mi>e</mi><mrow><mn>1</mn><mo>−</mo><msup><mi>e</mi><mrow><mfrac><mi>B</mi><mi>C</mi></mfrac></mrow></msup></mrow></msup></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where r is random number uniformly distributed over [0, 1]. This algorithm requires one log(log) calculation and two exponential calculations in Eq. (11). As it uses direct sampling method, the algorithm is more efficient than the sampling method proposed by Acosta et al . [ 17 17 ] which uses four random number, on average, to sample each value of polar angle. 2.2.3. Comparison between Kissel shape functions, Tsai DDCS and fitting function 3 Figure 4 Figure 4 shows comparisons of the polar angle distribution of bremsstrahlung photons using Tsai DDCS, Kissel et al . data and the analytical distribution function, Eq. (9), with reduced photon energy, K / T = 0.6, for 50 keV, 100 keV and K / T = 0.8 for 500 keV electrons incident on Silver and Gold. For low energy about 50 keV and 100 keV there is a significant deviation in the most probable angle between Tsai’s DDCS and Kissel’s data/Analytical function (equation 9). With an increase in the electron energy all three approaches tend to overlap and present a good agreement. Since Kissel et al . data is very accurate for electron energies below 500 keV, therefore we recommended to use analytical distribution function given by Eq. (9) in this energy range. For energies of > 500 keV, Tsai DDCS is a good choice as it retains the high sampling efficiency and physics accuracy 3 Results and Discussion 1 3.1. Photon Energy, sampling efficiency and accuracy 2 3.1.1 Sampling Accuracy 3 The accuracy of the sampling algorithm is determined by comparing sampled photon energy distribution with interpolated values of scaled DCSs (Seltzer and Berger) [ 23 23 ]. Two examples are presented. In figure 5a figure 5a , photon energy distribution is plotted for 100 keV electrons incident on Silver with K c =10 keV. Similar comparison is plotted in Fig.5(b) Fig.5(b) except for 100 MeV electrons. The smooth curves are from equation (2) and histograms are sampled values of photon energies from MC simulation. Fig. 5 Fig. 5 verified that the photon energy sampling algorithm works as expected. All plots in Fig. 5 Fig. 5 are normalized to unity. 3.1.2. Sampling efficiency 3 The sampling efficiency, P , of the photon energy sampling algorithm is given in Table 1 Table 1 . The sampling efficiency is defined as the ratio of the accepted sampled values to the total number of trial. We have calculated the efficiency in terms of incident electron energy. Although the sampling efficiency severely depends on the threshold photon energy and also less on Z , the qualitative behavior remains the same for other values of K c . We are therefore more interested in calculating P for different electron energies while keeping K c = 1keV and Z = 43 (Ag). 1 1 E (MeV) 1 1 0.05 1 1 0.1 1 1 0.5 1 1 1 1 1 10 1 1 50 1 1 100 1 1 500 1 1 1000 1 1 10000 1 1 P (%) 1 1 90 1 1 84 1 1 79 1 1 80 1 1 87 1 1 91 1 1 92 1 1 94 1 1 95 1 1 95 From Table 1 Table 1 , the sampling efficiency in energy region of E < 100 keV and E > 1 MeV is higher than that in 100 keV–1 MeV, which is obvious from the shape of the energy spectrum, as shown in Fig. 1. Fig. 1. Even for higher incident electron energies, the energy spectrum is flat except at end points of K/T . 3.2. Photon Angular Distribution, sampling accuracy and efficiency 2 This section concerns the verification that the angular sampling algorithms work as expected. The verification is performed by comparing the sampled angular distribution with the theoretical expressions. In Fig. 6 Fig. 6 , the sampled angular distribution from Tsai’ DDCS are plotted with the theoretical expression given in Eq. (5). The accuracy of angular sampling using analytical function in Eq. (9) is determined by fitting of Kissel’s data, as plotted in Fig. 7. Fig. 7. The smooth curves are from theoretical expressions, and the histograms are sampled angular distribution from Monte Carlo simulation. All plots are normalized to unity. Error bars in Figs. 6 Figs. 6 and 7 7 are uncertainty (±3σ) in each bin. The sampling efficiency of our proposed angular sampling algorithm based on Tsai’s DDCS and angular sampling algorithm proposed by Bielajew et al . (1989) implemented in EGS5 [ 26 26 ] and Geant4 [ 20 20 ] are plotted in Fig. 8. Fig. 8. The sampling efficiency is plotted for incident electron energies from 500 keV to 1000 MeV for silver and gold, and K / T = {0,0.5,1}. One can realize that the sampling efficiency based on Tsai’s DDCS is generally very good for both Z = 47 and Z = 79 and about 70% of average for K/T ~ 0,0.5 and reduces by few percent when K/T =1 . The sampling efficiency for Z = 79 decreases, at most, to 55% at 1 MeV, only when K/T = 1 it retains 62% up to 1000 MeV. So, the sampling technique employed here is more efficient in the useful photon energy range (lower reduced photon energies) than the sampling technique used to sample the angular distribution of photons, reported by Alex F. Bielajew et al . (1989) [ 15 15 ] using 2BS formula of Koch and Motz [ 13 13 ] which is implemented in EGS5 [ 26 26 ] and Geant4 [ 20 20 ]. 4 Conclusion 1 It can be concluded that the proposed algorithm provides accurate and efficient methods for sampling energy and direction of bremsstrahlung photons. The efficiency of sampling photon energy is more than 80% with this new algorithm. Hybrid model provides very accurate treatment of photon angular distribution by low as well as high energy incident electrons. Angular sampling which is based on Tsai DDCS is very efficient for useful photon energy range, compared to other available angular sampling techniques in high energy region (Alex F. Bielajew et al . 1989). Distribution function used for photon angular sampling using Tsai DDCS has contribution of atomic form factors and hence explains physics phenomena very well especially in high energy region. The parameterization of the shape functions given by Eq.(9) provides a very accurate fit to the Kissel’s benchmark data ( Fig. 4 Fig. 4 ) therefore simulated angular distribution is fully consistent with the most accurate and reliable theory at present and the sampling algorithm is more efficient than the previously developed method by Acosta et al . [ 24 24 ], which uses four random numbers, on average, to sample one photon energy. It is worth mentioning that the approaches used for sampling energy and direction of bremsstrahlung photons can be used for a wide range of incident electron energies and materials and hence suitable for general purpose Monte Carlo simulations. References F. Salvat , J.M. Fernández-Varea , J. Sempau , et al . Monte Carlo simulation of bremsstrahlung emission by electrons . Radiat Phys Chem , 2006 , 75 : 1201 - 1219 . DOI: 10.1016/j.radphyschem.2005.05.008 http://doi.org/10.1016/j.radphyschem.2005.05.008 . Monte Carlo simulation of bremsstrahlung emission by electrons Y. Wu , J. Song , H. Zheng , et al . CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC . Ann Nucl Energy , 2015 , 82 : 161 - 168 . DOI: 10.1016/j.anucene.2014.08.058 http://doi.org/10.1016/j.anucene.2014.08.058 . CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC Y. Wu . CAD-based interface programs for fusion neutron transport simulation . Fusion Eng Des , 2009 , 84 : 1987 - 1992 . DOI: 10.1016/j.fusengdes.2008.12.041 http://doi.org/10.1016/j.fusengdes.2008.12.041 . CAD-based interface programs for fusion neutron transport simulation Y. Li , L. Lu , A. Ding , et al . Benchmarking of MCAM 4.0 with the ITER 3D model . Fusion Eng Des , 2007 , 82 : 2861 - 2866 . DOI: 10.1016/j.fusengdes.2008.12.041 http://doi.org/10.1016/j.fusengdes.2008.12.041 . Benchmarking of MCAM 4.0 with the ITER 3D model Y. Wu , Z. Xie , U. Fischer . A Discrete ordinates nodal method for one-dimensional neutron transport calculation in curvilinear geometries . Nucl Sci Eng , 1999 , 133 ( 3 ): 350 - 357 . A Discrete ordinates nodal method for one-dimensional neutron transport calculation in curvilinear geometries Y. Wu , S. Zheng , X. Zhu , et al . Conceptual design of the fusion-driven subcritical system FDS-I . Fusion Eng Des , 2006 , 81 , Part B : 1305 - 1311 . DOI: 10.1016/j.fusengdes.2005.10.015 http://doi.org/10.1016/j.fusengdes.2005.10.015 . Conceptual design of the fusion-driven subcritical system FDS-I Y. Wu , FDS Team . Conceptual design of the China fusion power plant FDS-II . Fusion Eng Des , 2008 , 83 ( 10-12 ): 1683 - 1689 . DOI: 10.1016/j.fusengdes.2008.06.048 http://doi.org/10.1016/j.fusengdes.2008.06.048 . Conceptual design of the China fusion power plant FDS-II Y. Wu , FDS Team , Fusion-based hydrogen production reactor and its material selection . J Nucl Mater , 2009 , 386-388 : 122 - 126 . DOI: 10.1016/j.jnucmat.2008.12.075 http://doi.org/10.1016/j.jnucmat.2008.12.075 . Fusion-based hydrogen production reactor and its material selection Y. Wu , FDS Team . Conceptual design and testing strategy of a dual functional Lithium-Lead test blanket module in ITER and EAST . Nucl Fusion , 2007 , 47 ( 11 ): 1533 - 1539 . Conceptual design and testing strategy of a dual functional Lithium-Lead test blanket module in ITER and EAST Y. Wu , FDS Team . Design analysis of the China dual-functional lithium lead (DFLL) test blanket module in ITER . Fusion Eng Des , 2007 , 82 : 1893 - 1903 . DOI: 10.1016/j.fusengdes.2007.08.012 http://doi.org/10.1016/j.fusengdes.2007.08.012 . Design analysis of the China dual-functional lithium lead (DFLL) test blanket module in ITER Y. Wu . Progress in fusion-driven hybrid system studies in China . Fusion Eng Des , 2002 , 63-64 : 73 - 80 . DOI: 10.1016/S0920-3796(02)00239-9 http://doi.org/10.1016/S0920-3796(02)00239-9 . Progress in fusion-driven hybrid system studies in China Y. Wu , G. Song , R.F. Cao , et al . Development of accurate/advanced radiotherapy treatment planning and quality assurance system (ARTS) . Chinese Phys C , 2008 , 32 ( Suppl. II ): 177 - 182 . Development of accurate/advanced radiotherapy treatment planning and quality assurance system (ARTS) H.W. Koch , J.W. Motz . Bremsstrahlung cross section formulas and related data . Rev Mod Phys , 1959 , 31 : 920 - 955 . DOI: 10.1103/RevModPhys.31.920 http://doi.org/10.1103/RevModPhys.31.920 . Bremsstrahlung cross section formulas and related data Y.S. Tsai . Pair production and bremsstrahlung of charged leptons . Rev Mod Phys , 1974 , 46 : 815 - 851 . DOI: 10.1103/RevModPhys.46.815 http://doi.org/10.1103/RevModPhys.46.815 . Pair production and bremsstrahlung of charged leptons A.F. Bielajew , M. Rahde , S.C. Chen . Improved bremsstrahlung photon angular sampling in EGS4 code system . NRC Ottawa , Canada , PIRS-0203R , 1989 . Improved bremsstrahlung photon angular sampling in EGS4 code system W.R. Nelson , H. Hirayama , D.W.O. Rogers . The EGS4 code system . Stanford Linear Accelerator Centre , Stanford. CA . Report SLAC-265 , 1985 . The EGS4 code system E. Acosta , X. Llovet , F. Salvat . Monte Carlo simulation of bresmsstrahlung emission by electrons . Appl Phys Lett , 2002 , 80 : 3228 - 3230 . DOI: 10.1063/1.1473684 http://doi.org/10.1063/1.1473684 . Monte Carlo simulation of bresmsstrahlung emission by electrons L. Kissel , C.A. Quarles , R.H. Pratt . Shape functions for atomic field bremsstrahlung from electrons of kinetic energy 1-500 keV on selected neutral atoms 1≤Z≤92 . At Dta Nucl Data Tables , 1983 , 28 : 381 - 460 . DOI: 10.1016/0092-640X(83)90001-3 http://doi.org/10.1016/0092-640X(83)90001-3 . Shape functions for atomic field bremsstrahlung from electrons of kinetic energy 1-500 keV on selected neutral atoms 1≤Z≤92 H.K. Tseng , R.H. Pratt , C.M. Lee . Electron bremsstrahlung angular distribution in the 1-500 keV energy range . Phys Rev A , 1979 , 19 : 187 - 195 . DOI: 10.1103/PhysRevA.19.187 http://doi.org/10.1103/PhysRevA.19.187 . Electron bremsstrahlung angular distribution in the 1-500 keV energy range S. Agostinelli , J. Allison , K. Amako , et al . Geant4- a simulation toolkit . Nucl Instrum Meth A , 2003 , 506 : 250 - 303 . DOI: 10.1016/S0168-9002(03)01368-8 http://doi.org/10.1016/S0168-9002(03)01368-8 . Geant4- a simulation toolkit P. Rodrigues , R. Moura , C. Ortigao , et al . Geant 4 applications and developments for medical physics experiments . IEEE Trans Nucl Sci , 2004 , 51 : 1412 - 1419 . DOI: 10.1109/TNS.2004.832314 http://doi.org/10.1109/TNS.2004.832314 . Geant 4 applications and developments for medical physics experiments I. Kawrakow , D.W.O. Rogers . The EGSnrc code system. Monte Carlo simulation of electron and photon transport. PIRS 701-2001 The EGSnrc code system. Monte Carlo simulation of electron and photon transport. PIRS 701-2001 S.M. Seltzer , M.J. Berger . Bremsstrahlung energy spectra from electrons with kinetic energy 1 keV-10 GeV incident on screened nuclei and orbital electrons of neutral atoms with Z = 1-100 . At. Data Nucl. Data Tables , 1986 , 35 : 345 - 418 . DOI: 10.1016/0092-640X(86)90014-8 http://doi.org/10.1016/0092-640X(86)90014-8 . Bremsstrahlung energy spectra from electrons with kinetic energy 1 keV-10 GeV incident on screened nuclei and orbital electrons of neutral atoms with Z = 1-100 H. Davies , H.A. Bethe , L.C. Maximon . Theory of bremsstrahlung and pair production. II. Integral cross section for pair production . Phys Rev , 1954 , 93 : 788 - 795 . DOI: 10.1103/PhysRev.93.788 http://doi.org/10.1103/PhysRev.93.788 . Theory of bremsstrahlung and pair production. II. Integral cross section for pair production H. Olsen , L.C. Maximon . Photon and electron polarization in high-energy bremsstrahlung and pair production with screening . Phys Rev , 1989 , 114 : 887 - 904 . DOI: 10.1103/PhysRev.114.887 http://doi.org/10.1103/PhysRev.114.887 . Photon and electron polarization in high-energy bremsstrahlung and pair production with screening H. Hirayama , Y. Namito , A.F. Bielajew , et al . The EGS5 code system . Stanford Linear Accelerator Centre , Stanford. CA . Report SLAC-R-730 , 2015 . The EGS5 code system The authors would like to show their great appreciation to other members of FDS Team for their support and contribution to this research. 10.1007/s41365-017-0218-7 This work was supported by the National Natural Science Foundation of China (No. 81101132 and 11305203), the Strategic Priority Research Program of Chinese Academy of Sciences (No. XDA03040000) and the Natural Science Foundation of Anhui Province of China (No. 1508085QH180). 2016-05-12 2016-09-27 2016-10-07 2017-03-31 2017-05-01 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2017 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-017-0218-7 http://dx.doi.org/10.1007/s41365-017-0218-7 . 2017 Nuclear Science and Techniques 05 28