Simulation of discrete Bremsstrahlung by electrons

+-------------+                                               +----------##
| Geant 3.16  |               GEANT User's Guide              | PHYS341  ##
+-------------+                                               +----------##
                                   

Author(s) : G.N. Patrick, L.Urban Submitted: 26.10.84 Origin : Same Revised: 19.12.92

Subroutines

                    +--------------------------------+
                    |CALL GBREME  |
                    +--------------------------------+
                                  

GBREME generates a photon from Bremsstrahlung of an electron by treating it as a discrete process. The secondary photon energy is sampled from a parameterization of the Bremsstrahlung calculation of Seltzer and Berger [bib-SEL1] for electron energies below 10 GeV, and from the screened Bethe-Heitler cross-section above 10 GeV. Midgal corrections are applied in both cases. For the angular distribution of the photon, is calls the function GBTETH.

Input : via common block /GCTRAK/

Output: via common block /GCKING/.

GBREME is called automatically from the tracking routine GTELEC if, and when, the parent electron reaches its radiation point during the tracking stage of GEANT.

           THETA = GBTETH(ENER,PARTM,EFRAC)
                                  
GBTETH calculates the angular + - distribution for the e e -pair in pair production and for the photon in Bremsstrahlung. In case of Bremsstrahlung it gives the scaled angle for a photon radiated by an electron (mass PARTM) of energy ENER. The energy of photon is EFRAC times the initial energy of the electron. GBTETH is called by GBREME.

Method

The Coulomb corrected Bethe-Heitler differential cross-section for production of a photon of energy epsilonE (epsilon= k/E, where k is the photon energy) by an electron of primary energy E is given as in [bib-WILL], [bib-BUTC], [bib-EGS3]:

                                            2                                           2
    ((dsigma(Z, E,epsilon))/(depsilon)) =((r alphaZ[Z+xi(Z)])/ (epsilon)){[1+(1-epsilon) ][Phi (delta)- F(Z)]-((2)/(3))(1- epsilon)[Phi (delta)-F(Z)]}(1)
                                            0                                                 1                                        2

where, Phi (delta) are the screening functions depending on the screening i variable delta

                          1/3
        delta =((136m)/ (Z   E))  ((epsilon)/((1-epsilon)))          m =electron mass
                                                     2
          Phi (delta)= 20.867-3. 242delta+ 0.625delta
      .      1                                          }             delta<=1
                                                    2
          Phi (delta)= 20.209-1. 930delta-0.086delta
        Phi (delta) =Phi  (delta)= 21.12-4. 184ln(delta+ 0.952)
           1            2                                             delta>1
        F(Z)   ={  4/ 3lnZ          .                                  E<0. 05GeV
                   4/ 3lnZ +4f (Z)                                     E> =0.05 GeV
                              c
                           2/3             1/3
        xi(Z) =((ln (1440/Z   ))/(ln (183/Z   )- f (Z)))
                                                  c
        f (Z)    the Coulomb correction function
         c

                                                         2        3
    f (Z)    =    a(1/ (1+ a)+0. 20206-0.0369a +0. 0083*a -0. 002a )
     c
                            2
        a    =    (alpha* Z)
    alpha    =    1. /137.

The kinematically allowed region for epsilon is

    epsilon ==((k )/ (E))<=epsilon< =1-((m)/ (E))==epsilon             (2)
           c     c                                        1

where, k is the photon cut-off energy below which the Bremsstrahlung is c treated as a continuous energy loss (BCUTE in the program). The cross-section (1) can be decomposed as (the normalization is not important !)

                             2
    ((dsigma)/ (depsilon))= i=1 alpha f (epsilon)g (epsilon)           (3)
                                     i i          i

  where
                                                                                                                          2        2
       alpha = F  ln((epsilon )/ (epsilon ))                                               alpha  =F  ((3)/ (8)) ((epsilon -epsilon )/(1- epsilon ))
            1   10           1           c                                                      2   20                                           c
                                                                                                                          1        c
                                                                                                                             2        2
       f (epsilon)= ((1)/(epsilon ln(epsilon / epsilon )))                                 f (epsilon) =((2epsilon)/ (epsilon -epsilon ))
        1                                   1         c                                     2                                1        c
       g (epsilon)= ((F )/(F  ))((1- epsilon)/(1-epsilon ))M(epsilon)                      g (epsilon) =((F )/ (F  )).C (epsilon)
        1              1    10                          c                                   2              2     20    M
  and
       F (delta)= 3Phi (delta)-Phi (delta)- 2F(Z)                                          F   =F (delta   )
        1             1           2                                                         10   1      min
       F (delta)= 2Phi (delta)-2F(Z)                                                       F   =F (delta   )
        2             1                                                                     20   2      min
                           1/3
       delta   = ((136m)/(Z   E))((epsilon )/ (1-epsilon ))
            min                           c             c

delta is the minimum value of the variable delta, and C (epsilon) is min M the Midgal correction factor [bib-MIGD]:

                                 2                 2
    C (epsilon) =((1 +C / epsilon )/(1 +C / epsilon ))
     M                 0         1       0

where

                   2
    C  =((nr  lamb-a)/ (pi))
     0      0

n
electron density in the medium
r
classical electron radius 0
lamb-
reduced Compton wavelength of the electron.

This correction decreases the cross-section for low photon energy. Using the decomposition (3) the sampling of epsilon can be done by (for a set of random numbers: 0<=r < =1) i

  1. selecting the integer i (1 or 2) if r <=BPAR =alpha / (alpha + alpha ), then i= 1 0 1 1 2 if r >BPAR, then i= 2 0
  2. sampling epsilon from f (epsilon) i r1 epsilon= epsilonc(epsilon /epsilon ) when i =1 1 c 2 2 2 epsilon= sqrt(epsilon + r (epsilonm -epsilon )) when i =2 1 c c
  3. calculating the rejection function g (epsilon) and i if r >g (epsilon) rejecting epsilon, starting again from 1 2 i if r <=g (epsilon) accepting epsilon 2 i

As in the case of the pair production after sampling epsilon from f (epsilon), we have to check that i

                      1/3
    delta =((136m)/ (Z   E))((epsilon)/(1-epsilon  ))<=delta    =exp (((21. 12-f(Z))/(4. 184)))-0.952
                                                 c          max

and, if this is not the case, we have to start again from step 1.

The decomposition (3) used here is simpler and more effective than that used in the earlier versions of GEANT and in EGS, and the method outlined above has no convergence problems.

After the successful sampling of epsilon, GBREME generates the polar angles of the radiated photon with respect to an axis defined along the parent electron's momentum. The energy-angle distribution is given by Tsai [bib-TSAI], [bib-TSAI-err] as following:

                                           2 2       4                    2                       4     2
    ((dsigma)/ (dkdOmega))    =    ((2alpha e )/(pikm )){[((2y- 2)/((1+ l) ))+ ((12l(1-y))/((1 +l) ))](Z  +Z).
                                               2         2                      4        2          2
                                   +.[((2-2y- y )/((1+ l) ))-((4l(1- y))/((1+ l) ))][X-2Z f((alphaZ)  )]}
                                                                                                            (4)

where k is the photon energy, E is the initial electron energy, y= k/E 2 2 2 and l= E theta /m . This distribution is quite complicated for sampling and, furthermore, for a variable u= Etheta/m, shows a very weak dependence on Z, E(k), y= k/E. Thus, the distribution can be approximated by a function

              -au     -3au
    f(u) =C(ue    +due    )                                            (5)

where

                  2
    C    =    ((9a )/ (9+ d))
    a    =    0. 625
    d    =    27. 0

The sampling of the function f(u) can be done in the following way (r i are uniformly distributed random numbers in [0,1]):

    -au -3au
  1. Choose between ue and due , weights: 9/(9+ d) and d/(9+ d). If r <9/(9+ d) then a= b, 1 else b= 3a. -bu
  2. Sample ue , u=-(ln r + lnr )/b. 2 3

    The azimuthal angle, Phi, is generated isotropically. This information is used to calculate the momentum vector of the radiated photon, to transform it to the GEANT coordinate system and to store the result into common block /GCKING/. Also, the momentum of the parent electron is updated.

    Restrictions

    1. Effects due to breakdown of the Born approximation at low energies are ignored.
    2. Target materials composed of compounds or mixtures are treated identically to elements using the effective atomic number Z calculated in GSMIXT (this is not the case when computing the mean free path!).