source: trunk/examples/extended/electromagnetic/TestEm13/geant3/src/gtgama.F @ 1309

Last change on this file since 1309 was 807, checked in by garnier, 16 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.