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