source: trunk/examples/extended/electromagnetic/TestEm1/geant3/src/gtgama.F@ 1230

Last change on this file since 1230 was 807, checked in by garnier, 17 years ago

update
File size: 7.0 KB
Line 
1
2 SUBROUTINE GTGAMA
3C.
4C. ******************************************************************
5C. * *
6C. * Photon track. Computes step size and propagates particle *
7C. * through step. *
8C. * *
9C. * ==>Called by : GTRACK *
10C. * Authors R.Brun, F.Bruyant L.Urban ******** *
11C. * *
12C. ******************************************************************
13C.
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)
30C.
31*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
32 IABAN = NINT(DPHYS1)
33*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
34C. ------------------------------------------------------------------
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)
16810000 FORMAT(' Boundary correction in GTGAMA: ',
169 + ' GEKIN NUMED STEP SNEXT')
17010100 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
260998 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
Note: See TracBrowser for help on using the repository browser.