Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gtgama.F
Go to the documentation of this file.
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