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 |
---|