GEANT321/gphys/gbreme.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.2  1997/11/13 09:25:40  gunter
   6 * Correction By Laszlo Urban; protect against divide by 0. if AL becomes 0.
   7 * around GEKIN 2.93E-5.
   8 * Reported by harald@psiclu.cern.ch (signed by keller@biomed.ee.ethz.ch)
   9 *
  10 * Revision 1.1.1.1  1995/10/24 10:21:22  cernlib
  11 * Geant
  12 *
  13 *
  14 #include "geant321/pilot.h"
  15 *CMZ :  3.21/03 02/08/94  19.37.36  by  S.Ravndal
  16 *-- Author :
  17       SUBROUTINE GBREME
  18 C.
  19 C.    ******************************************************************
  20 C.    *                                                                *
  21 C.    *  Simulates discrete hard BREMSSTRAHLUNG by electrons.          *
  22 C.    *                                                                *
  23 C.    *  The secondary photon energy is sampled from a                 *
  24 C.    *  parametrization of the bremsstrahlung calculation of          *
  25 C.    *  Seltzer and Berger (NIM B12,p.95(1985)) for electron          *
  26 C.    *  energies 1 keV - 10 GeV . For higher energies the             *
  27 C.    *  parametrization agrees with the screened Bethe-Heitler        *
  28 C.    *  bremsstrahlung spectrum.Migdal corrections are applied        *
  29 C.    *      by default. The routine works ( together                  *
  30 C.    *  with the routines GBRELE and GBRSGE ) without the Migdal      *
  31 C.    *  corrections using the Patchy switch +USE,BETHE.               *
  32 C.    *                                                                *
  33 C.    *  NOTE :                                                        *
  34 C.    *      BCUTE is the cut-off energy above which the photon energy *
  35 C.    *      spectrum is generated.                                    *
  36 C.    *                                                                *
  37 C.    *    ==>Called by : GTELEC                                       *
  38 C.    *       Authors    L.Urban    29/10/93 *********                 *
  39 C.    *                                                                *
  40 C.    *  13/11/97 bug corrected by L.Urban                             *
  41 C.    ******************************************************************
  42 C.
  43 #include "geant321/gcbank.inc"
  44 #include "geant321/gcphys.inc"
  45 #include "geant321/gcjloc.inc"
  46 #include "geant321/gconsp.inc"
  47 #include "geant321/gctrak.inc"
  48 #include "geant321/gcking.inc"
  49 #include "geant321/gcmate.inc"
  50 #include "geant321/gccuts.inc"
  51       DIMENSION RNDM(2)
  52       PARAMETER (CORFAC=0.805485E-10)
  53       PARAMETER (TLIM=0.000999)
  54       PARAMETER  (
  55      +AL00=-0.205398E+01,AL01= 0.238815E-01,AL02= 0.525483E-03,
  56      +AL10=-0.769748E-01,AL11=-0.691499E-01,AL12= 0.222453E-02,
  57      +AL20= 0.406463E-01,AL21=-0.101281E-01,AL22= 0.340919E-03,
  58      +A10 = 0.467733E+01,A11 =-0.619012E+00,A12 = 0.202225E-01,
  59      +A20 =-0.734101E+01,A21 = 0.100462E+01,A22 =-0.320985E-01,
  60      +A30 = 0.293119E+01,A31 =-0.403761E+00,A32 = 0.125153E-01,
  61      +BL00= 0.104133E+01,BL01=-0.943291E-02,BL02=-0.454758E-03,
  62      +BL10= 0.119253E+00,BL11= 0.407467E-01,BL12=-0.130718E-02,
  63      +BL20=-0.159391E-01,BL21= 0.727752E-02,BL22=-0.194405E-03,
  64      +B10 = 0.423071E+01,B11 =-0.610995E+00,B12 = 0.195531E-01,
  65      +B20 =-0.712527E+01,B21 = 0.969160E+00,B22 =-0.274255E-01,
  66      +B30 = 0.269925E+01,B31 =-0.363283E+00,B32 = 0.955316E-02)
  67 C.    ------------------------------------------------------------------
  68 C.
  69 C             Ensure cut-off avoids infra-red catastrophe.
  70 C
  71       IF (GEKIN.LE.BCUTE) GO TO 30
  72       KCASE=NAMEC(9)
  73 C
  74 *******************************>  Z3=Q(JPROB+2)
  75       Z3=(Z*(Z+1.))**0.3333333
  76       IF(Z3.LE.0.)GO TO 30
  77       Z32=Z3**2
  78       EEL1   = GETOT
  79       XC=BCUTE/GEKIN
  80       ALXC=LOG(XC)
  81       U=LOG(GEKIN/EMASS)
  82       U2=U**2
  83       V=LOG(Z)
  84 *
  85       IF(GEKIN.LE.TLIM) THEN
  86          AL0=AL00+AL01*Z3+AL02*Z32
  87          AL1=AL10+AL11*Z3+AL12*Z32
  88          AL2=AL20+AL21*Z3+AL22*Z32
  89          AL=AL0+AL1*U+AL2*U2
  90          BL0=BL00+BL01*Z3+BL02*Z32
  91          BL1=BL10+BL11*Z3+BL12*Z32
  92          BL2=BL20+BL21*Z3+BL22*Z32
  93          BL=BL0+BL1*U+BL2*U2
  94          GMAX=1.+AL*XC+BL*XC**2
  95          IF(GEKIN.LT.0.0001) THEN
  96             G1=1.+AL+BL
  97             IF(G1.GT.GMAX) GMAX=G1
  98 *
  99             IF(ABS(AL).GT.1.e-6) THEN
 100               X0=-BL/(2.*AL)
 101               IF((XC.LT.X0).AND.(X0.LT.1.)) THEN
 102                  G0=1.+AL*X0+BL*X0**2
 103                  IF(G0.GT.GMAX) GMAX=G0
 104               ENDIF
 105             ENDIF
 106 *
 107          ENDIF
 108       ELSE
 109          U3=U2*U
 110          A1=A10+A11*Z3+A12*Z32
 111          A2=A20+A21*Z3+A22*Z32
 112          A3=A30+A31*Z3+A32*Z32
 113          AH=1.+A1/U+A2/U2+A3/U3
 114          B1=B10+B11*Z3+B12*Z32
 115          B2=B20+B21*Z3+B22*Z32
 116          B3=B30+B31*Z3+B32*Z32
 117          BH=0.75+B1/U+B2/U2+B3/U3
 118 *
 119          F=4*V-0.55*V**2
 120          DEL0=136.*EMASS/(Z3*EEL1)
 121          EPC=XC*GEKIN/EEL1
 122          DC=DEL0*EPC/(1.-EPC)
 123          CC=42.392-F
 124 *
 125          IF(DC.LE.1.) THEN
 126             DC2=DC**2
 127             F1=(42.392-7.796*DC+1.961*DC2-F)/CC
 128             F2=(41.734-6.484*DC+1.250*DC2-F)/CC
 129          ELSE
 130             F1=(42.24-8.368*LOG(DC+0.952)-F)/CC
 131             IF(F1.LT.0.) F1=0.
 132             F2=F1
 133          ENDIF
 134 *
 135          GMAX=(1.-AH*EPC)*F1+BH*EPC**2*F2
 136       ENDIF
 137 *
 138       CORR0=CORFAC*DENS*Z/A
 139       EPM=GEKIN/EEL1
 140       SC0=1.+CORR0/(EPM*EPM)
 141 *
 142 *     sample photon energy  ( according to 1/Ephoton)
 143 *
 144    10 CALL GRNDM(RNDM,2)
 145 *
 146       X=EXP(RNDM(1)*ALXC)
 147       EP=X*GEKIN/EEL1
 148 *
 149 *     Migdal correction for Ephoton->0. or no correction (Bethe)
 150 *
 151 #if !defined(CERNLIB_BETHE)
 152       CORR=SC0/(1.+CORR0/(EP*EP))
 153 #endif
 154 #if defined(CERNLIB_BETHE)
 155       CORR=1.
 156 #endif
 157 *
 158 *     calculate rejection function g(x)
 159 *
 160       IF(GEKIN.LE.TLIM) THEN
 161          G=1.+AL*X+BL*X**2
 162       ELSE
 163          D=DEL0*EP/(1.-EP)
 164          IF(D.LE.1.) THEN
 165             D2=D**2
 166             F1=(42.392-7.796*D+1.961*D2-F)/CC
 167             F2=(41.734-6.484*D+1.250*D2-F)/CC
 168          ELSE
 169             F1=(42.24-8.368*LOG(D+0.952)-F)/CC
 170             IF(F1.LT.0.) F1=0.
 171             F2=F1
 172          ENDIF
 173          G=(1.-AH*EP)*F1+BH*EP**2*F2
 174       ENDIF
 175       G=G*CORR/GMAX
 176       IF(RNDM(2).GT.G) GOTO 10
 177 *
 178 *     photon energy is sampled according to the Seltzer-Berger spectrum
 179 *
 180       EGAMMA=EEL1*EP
 181 C
 182 C        CUT ON ENERGY THRESHOLD ?
 183 C
 184       IF((IBREM.NE.1).OR.(EGAMMA.LE.CUTGAM)) THEN
 185          DESTEP = DESTEP + EGAMMA
 186          GO TO 20
 187       ENDIF
 188 C
 189 C             Generate emitted photon angles with respect to a Z-axis
 190 C             defined along parent track. PHI is generated isotropically
 191 C             and THETA is assigned a universal angular distribution
 192 C
 193       EMASS1 = EMASS
 194       THETA  = GBTETH(EEL1, EMASS1, EP)*EMASS/EEL1
 195       SINTH  = SIN(THETA)
 196       COSTH  = COS(THETA)
 197       CALL GRNDM(RNDM,1)
 198       PHI    = TWOPI*RNDM(1)
 199       COSPHI = COS(PHI)
 200       SINPHI = SIN(PHI)
 201 C
 202 C             Polar co-ordinates to momentum components.
 203 C
 204       NGKINE = NGKINE + 1
 205       GKIN(1,1)=EGAMMA*SINTH*COSPHI
 206       GKIN(2,1)=EGAMMA*SINTH*SINPHI
 207       GKIN(3,1)=EGAMMA*COSTH
 208       GKIN(4,1)=EGAMMA
 209       GKIN(5,1)=1.
 210       TOFD(1)  =0.
 211       GPOS(1,1) = VECT(1)
 212       GPOS(2,1) = VECT(2)
 213       GPOS(3,1) = VECT(3)
 214 C
 215 C             Rotate photon into GEANT system
 216 C
 217       CALL GVROT(VECT(4),GKIN)
 218 C
 219 C             Correct track for lost energy
 220 C
 221    20 CONTINUE
 222       GEKIN = GEKIN - EGAMMA
 223       GETOT = GEKIN + EMASS
 224       VECT(7)=SQRT (ABS((GETOT+EMASS)*GEKIN))
 225       CALL GEKBIN
 226 C
 227 C             Update probabilities
 228 C
 229    30 CALL GRNDM(RNDM,1)
 230       ZINTBR=-LOG(RNDM(1))
 231 C
 232       END