GEANT321/gtrak/gtgama.F

   1 #include "geant321/pilot.h"
   2 *CMZ :  3.21/02 29/03/94  15.41.24  by  S.Giani
   3 *-- Author :
   4       SUBROUTINE GTGAMA
   5 C.
   6 C.    ******************************************************************
   7 C.    *                                                                *
   8 C.    *   Photon track. Computes step size and propagates particle     *
   9 C.    *    through step.                                               *
  10 C.    *                                                                *
  11 C.    *   ==>Called by : GTRACK                                        *
  12 C.    *      Authors    R.Brun, F.Bruyant L.Urban ********             *
  13 C.    *                                                                *
  14 C.    ******************************************************************
  15 C.
  16 #include "geant321/gcbank.inc"
  17 #include "geant321/gccuts.inc"
  18 #include "geant321/gcjloc.inc"
  19 #include "geant321/gconsp.inc"
  20 #include "geant321/gcphys.inc"
  21 #include "geant321/gcstak.inc"
  22 #include "geant321/gctmed.inc"
  23 #include "geant321/gcmulo.inc"
  24 #include "geant321/gctrak.inc"
  25 #if defined(CERNLIB_DEBUG)
  26 #include "geant321/gcunit.inc"
  27 #endif
  28 #if !defined(CERNLIB_SINGLE)
  29       PARAMETER (EPSMAC=1.E-6)
  30       DOUBLE PRECISION ONE,XCOEF1,XCOEF2,XCOEF3,ZERO
  31 #endif
  32 #if defined(CERNLIB_SINGLE)
  33       PARAMETER (EPSMAC=1.E-11)
  34 #endif
  35       PARAMETER (ONE=1,ZERO=0)
  36       PARAMETER (EPCUT=1.022E-3)
  37 C.
  38 C.    ------------------------------------------------------------------
  39 *
  40 * *** Particle below energy threshold ?  Short circuit
  41 *
  42 *
  43       IF (GEKIN.LE.CUTGAM) GOTO 998
  44 *
  45 * *** Update local pointers if medium has changed
  46 *
  47       IF(IUPD.EQ.0)THEN
  48          IUPD  = 1
  49          JPHOT = LQ(JMA-6)
  50          JCOMP = LQ(JMA-8)
  51          JPAIR = LQ(JMA-10)
  52          JPFIS = LQ(JMA-12)
  53          JRAYL = LQ(JMA-13)
  54       ENDIF
  55 *
  56 * *** Compute current step size
  57 *
  58       IPROC  = 103
  59       STEP   = STEMAX
  60       GEKRT1 = 1 .-GEKRAT
  61 *
  62 *  **   Step limitation due to pair production ?
  63 *
  64       IF (GETOT.GT.EPCUT) THEN
  65          IF (IPAIR.GT.0) THEN
  66             STEPPA = GEKRT1*Q(JPAIR+IEKBIN) +GEKRAT*Q(JPAIR+IEKBIN+1)
  67             SPAIR  = STEPPA*ZINTPA
  68             IF (SPAIR.LT.STEP) THEN
  69                STEP  = SPAIR
  70                IPROC = 6
  71             ENDIF
  72          ENDIF
  73       ENDIF
  74 *
  75 *  **   Step limitation due to Compton scattering ?
  76 *
  77       IF (ICOMP.GT.0) THEN
  78          STEPCO = GEKRT1*Q(JCOMP+IEKBIN) +GEKRAT*Q(JCOMP+IEKBIN+1)
  79          SCOMP  = STEPCO*ZINTCO
  80          IF (SCOMP.LT.STEP) THEN
  81             STEP  = SCOMP
  82             IPROC = 7
  83          ENDIF
  84       ENDIF
  85 *
  86 *  **   Step limitation due to photo-electric effect ?
  87 *
  88       IF (GEKIN.LT.0.4) THEN
  89          IF (IPHOT.GT.0) THEN
  90             STEPPH = GEKRT1*Q(JPHOT+IEKBIN) +GEKRAT*Q(JPHOT+IEKBIN+1)
  91             SPHOT  = STEPPH*ZINTPH
  92             IF (SPHOT.LT.STEP) THEN
  93                STEP  = SPHOT
  94                IPROC = 8
  95             ENDIF
  96          ENDIF
  97       ENDIF
  98 *
  99 *  **   Step limitation due to photo-fission ?
 100 *
 101       IF (JPFIS.GT.0) THEN
 102          STEPPF = GEKRT1*Q(JPFIS+IEKBIN) +GEKRAT*Q(JPFIS+IEKBIN+1)
 103          SPFIS  = STEPPF*ZINTPF
 104          IF (SPFIS.LT.STEP) THEN
 105             STEP  = SPFIS
 106             IPROC = 23
 107          ENDIF
 108       ENDIF
 109 *
 110 *  **   Step limitation due to Rayleigh scattering ?
 111 *
 112       IF (IRAYL.GT.0) THEN
 113          IF (GEKIN.LT.0.01) THEN
 114             STEPRA = GEKRT1*Q(JRAYL+IEKBIN) +GEKRAT*Q(JRAYL+IEKBIN+1)
 115             SRAYL  = STEPRA*ZINTRA
 116             IF (SRAYL.LT.STEP) THEN
 117                STEP  = SRAYL
 118                IPROC = 25
 119             ENDIF
 120          ENDIF
 121       ENDIF
 122 *
 123       IF (STEP.LT.0.) STEP = 0.
 124 *
 125 *  **   Step limitation due to geometry ?
 126 *
 127       IF (STEP.GE.SAFETY) THEN
 128          CALL GTNEXT
 129          IF (IGNEXT.NE.0) THEN
 130             STEP   = SNEXT + PREC
 131             INWVOL= 2
 132             IPROC = 0
 133             NMEC =  1
 134             LMEC(1)=1
 135          ENDIF
 136 *
 137 *        Update SAFETY in stack companions, if any
 138          IF (IQ(JSTAK+3).NE.0) THEN
 139             DO 10 IST = IQ(JSTAK+3),IQ(JSTAK+1)
 140                JST = JSTAK +3 +(IST-1)*NWSTAK
 141                Q(JST+11) = SAFETY
 142    10       CONTINUE
 143             IQ(JSTAK+3) = 0
 144          ENDIF
 145 *
 146       ELSE
 147          IQ(JSTAK+3) = 0
 148       ENDIF
 149 *
 150 * *** Linear transport
 151 *
 152       IF (INWVOL.EQ.2) THEN
 153          DO 20 I = 1,3
 154             VECTMP  = VECT(I) +STEP*VECT(I+3)
 155             IF(VECTMP.EQ.VECT(I)) THEN
 156 *
 157 * *** Correct for machine precision
 158 *
 159                IF(VECT(I+3).NE.0.) THEN
 160                   VECTMP = VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*
 161      +            EPSMAC
 162                   IF(NMEC.GT.0) THEN
 163                      IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
 164                   ENDIF
 165                   NMEC=NMEC+1
 166                   LMEC(NMEC)=104
 167 #if defined(CERNLIB_DEBUG)
 168                   WRITE(CHMAIL, 10000)
 169                   CALL GMAIL(0,0)
 170                   WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
 171                   CALL GMAIL(0,0)
 172 10000 FORMAT(' Boundary correction in GTGAMA: ',
 173      +       '    GEKIN      NUMED       STEP      SNEXT')
 174 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
 175 #endif
 176                ENDIF
 177             ENDIF
 178             VECT(I) = VECTMP
 179    20    CONTINUE
 180       ELSE
 181          DO 30 I = 1,3
 182             VECT(I)  = VECT(I) +STEP*VECT(I+3)
 183    30    CONTINUE
 184       ENDIF
 185 *
 186       SLENG = SLENG +STEP
 187 *
 188 * *** Update time of flight
 189 *
 190       TOFG = TOFG +STEP/CLIGHT
 191 *
 192 * *** Update interaction probabilities
 193 *
 194       IF (GETOT.GT.EPCUT) THEN
 195          IF (IPAIR.GT.0) ZINTPA = ZINTPA -STEP/STEPPA
 196       ENDIF
 197       IF (ICOMP.GT.0)    ZINTCO = ZINTCO -STEP/STEPCO
 198       IF (GEKIN.LT.0.4) THEN
 199          IF (IPHOT.GT.0) ZINTPH = ZINTPH -STEP/STEPPH
 200       ENDIF
 201       IF (JPFIS.GT.0)    ZINTPF = ZINTPF -STEP/STEPPF
 202       IF (IRAYL.GT.0) THEN
 203          IF (GEKIN.LT.0.01) ZINTRA = ZINTRA -STEP/STEPRA
 204       ENDIF
 205 *
 206       IF (IPROC.EQ.0) GO TO 999
 207       NMEC = 1
 208       LMEC(1) = IPROC
 209 *
 210 *  ** Pair production ?
 211 *
 212       IF (IPROC.EQ.6) THEN
 213          CALL GPAIRG
 214 *
 215 *  ** Compton scattering ?
 216 *
 217       ELSE IF (IPROC.EQ.7) THEN
 218          CALL GCOMP
 219 *
 220 *  ** Photo-electric effect ?
 221 *
 222       ELSE IF (IPROC.EQ.8) THEN
 223 *       Calculate range of the photoelectron ( with kin. energy Ephot)
 224 *
 225         IF(GEKIN.LE.0.001)  THEN
 226             JCOEF = LQ(JMA-17)
 227          IF(GEKRAT.LT.0.7) THEN
 228             I1 = MAX(IEKBIN-1,1)
 229          ELSE
 230             I1 = MIN(IEKBIN,NEKBIN-1)
 231          ENDIF
 232          I1 = 3*(I1-1)+1
 233          XCOEF1 = Q(JCOEF+I1)
 234          XCOEF2 = Q(JCOEF+I1+1)
 235          XCOEF3 = Q(JCOEF+I1+2)
 236          IF(XCOEF1.NE.0.) THEN
 237             STOPMX = -XCOEF2+SIGN(ONE,XCOEF1)*SQRT(XCOEF2**2 - (XCOEF3-
 238      +      GEKIN/XCOEF1))
 239          ELSE
 240             STOPMX = - (XCOEF3-GEKIN)/XCOEF2
 241          ENDIF
 242 *
 243 *        DO NOT call GPHOT if this (overestimated) range is smaller
 244 *                than SAFETY
 245 *
 246          IF (STOPMX.LE.SAFETY) GOTO 998
 247         ENDIF
 248
 249          CALL GPHOT
 250 *
 251 *  ** Rayleigh effect ?
 252 *
 253       ELSE IF (IPROC.EQ.25) THEN
 254          CALL GRAYL
 255 *
 256 *  ** Photo-fission ?
 257 *
 258       ELSE IF (IPROC.EQ.23) THEN
 259          CALL GPFIS
 260 *
 261       ENDIF
 262 *
 263          GOTO 999
 264 998      DESTEP = GEKIN
 265          GEKIN  = 0.
 266          GETOT  = 0.
 267          VECT(7)= 0.
 268          ISTOP  = 2
 269          NMEC   = 1
 270          LMEC(1)= 30
 271   999 END