This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gtrak / gtneut.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:44  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.24  by  S.Giani
11 *-- Author :
12       SUBROUTINE GTNEUT
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *   Neutral hadron type track. Computes step size and propagates *
17 C.    *    particle through step.                                      *
18 C.    *                                                                *
19 C.    *   ==>Called by : GTRACK                                        *
20 C.    *      Authors     R.Brun, F.Bruyant  *********                  *
21 C.    *                                                                *
22 C.    ******************************************************************
23 C.
24 #include "geant321/gcbank.inc"
25 #include "geant321/gccuts.inc"
26 #include "geant321/gckine.inc"
27 #include "geant321/gcking.inc"
28 #include "geant321/gconsp.inc"
29 #include "geant321/gcphys.inc"
30 #include "geant321/gcstak.inc"
31 #include "geant321/gctmed.inc"
32 #include "geant321/gctrak.inc"
33 #if defined(CERNLIB_USRJMP)
34 #include "geant321/gcjump.inc"
35 #endif
36 #if defined(CERNLIB_DEBUG)
37 #include "geant321/gcunit.inc"
38 #endif
39 #if !defined(CERNLIB_SINGLE)
40       PARAMETER (EPSMAC=1.E-6)
41 #endif
42 #if defined(CERNLIB_SINGLE)
43       PARAMETER (EPSMAC=1.E-11)
44 #endif
45 C.
46 C.    ------------------------------------------------------------------
47 *
48 * *** Particle below energy threshold ? Special treatment
49 *
50       IF (GEKIN.LE.CUTNEU) THEN
51          GEKIN  = 0.
52          GETOT  = AMASS
53          VECT(7)= 0.
54          ISTOP  = 2
55          NMEC = NMEC + 1
56          LMEC(NMEC) = 30
57          IF (IHADR.EQ.0.AND.AMASS.GT.0.) THEN
58             IF (IDCAY.NE.0) THEN
59                GAMMA  = GETOT/AMASS
60                TOFG = TOFG +GAMMA*SUMLIF/CLIGHT
61                SUMLIF = 0.
62                IF (TOFG.GE.TOFMAX) THEN
63                   ISTOP = 4
64                   NMEC = NMEC + 1
65                   LMEC(NMEC) = 22
66                   NGKINE = 0
67                ELSE
68                   NMEC = NMEC + 1
69                   LMEC(NMEC) = 5
70                   ISTOP =1
71                   CALL GDECAY
72                ENDIF
73             ENDIF
74             GO TO 999
75          ENDIF
76          IPROC = 12
77          GO TO 40
78       ENDIF
79 *
80 * *** Compute step size
81 *
82       IPROC  = 103
83       STEP   = STEMAX
84 *
85 *  **   Step limitation due to hadron interaction ?
86 *
87       IF (IHADR.GT.0) THEN
88 #if !defined(CERNLIB_USRJMP)
89          CALL GUPHAD
90 #endif
91 #if defined(CERNLIB_USRJMP)
92          CALL JUMPT0(JUPHAD)
93 #endif
94          IF (SHADR.LT.STEP) THEN
95             IF (SHADR.LE.0.) SHADR = PREC
96             STEP  = SHADR
97             IPROC = 12
98          ENDIF
99       ENDIF
100 *
101 *  **   Step limitation due to decay ?
102 *
103       IF ((IDCAY.NE.0).AND.(AMASS.GT.0.)) THEN
104          SDCAY = SUMLIF*VECT(7)/AMASS
105          IF (SDCAY.LT.STEP) THEN
106             STEP  = SDCAY
107             IPROC = 5
108          ENDIF
109       ENDIF
110 *
111       IF (STEP.LT.0.) STEP = 0.
112 *
113 *  **   Step limitation due to geometry ?
114 *
115       IF (STEP.GE.SAFETY) THEN
116          CALL GTNEXT
117          IF (IGNEXT.NE.0) THEN
118             STEP = SNEXT + PREC
119             IPROC = 0
120             INWVOL = 2
121             NMEC = NMEC + 1
122             LMEC(NMEC)  = 1
123          ENDIF
124 *
125 *        Update SAFETY in stack companions, if any
126          IF (IQ(JSTAK+3).NE.0) THEN
127             DO 10 IST = IQ(JSTAK+3),IQ(JSTAK+1)
128                JST = JSTAK +3 +(IST-1)*NWSTAK
129                Q(JST+11) = SAFETY
130    10       CONTINUE
131             IQ(JSTAK+3) = 0
132          ENDIF
133 *
134       ELSE
135          IQ(JSTAK+3) = 0
136       ENDIF
137 *
138 * *** Linear transport
139 *
140       IF (INWVOL.EQ.2) THEN
141          DO 20 I = 1,3
142             VECTMP  = VECT(I) +STEP*VECT(I+3)
143             IF(VECTMP.EQ.VECT(I)) THEN
144 *
145 * *** Correct for machine precision
146 *
147                IF(VECT(I+3).NE.0.) THEN
148                   VECTMP = VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*
149      +            EPSMAC
150                   IF(NMEC.GT.0) THEN
151                      IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
152                   ENDIF
153                   NMEC=NMEC+1
154                   LMEC(NMEC)=104
155 #if defined(CERNLIB_DEBUG)
156                   WRITE(CHMAIL, 10000)
157                   CALL GMAIL(0,0)
158                   WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
159                   CALL GMAIL(0,0)
160 10000 FORMAT(' Boundary correction in GTNEUT: ',
161      +       '    GEKIN      NUMED       STEP      SNEXT')
162 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
163 #endif
164                ENDIF
165             ENDIF
166             VECT(I) = VECTMP
167    20    CONTINUE
168       ELSE
169          DO 30 I = 1,3
170             VECT(I)  = VECT(I) +STEP*VECT(I+3)
171    30    CONTINUE
172       ENDIF
173 *
174       SLENG = SLENG +STEP
175 *
176 * *** Update time of flight
177 *
178       SUMLIF = SUMLIF -STEP*AMASS/VECT(7)
179       TOFG = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
180       IF (TOFG.GE.TOFMAX) THEN
181          ISTOP = 4
182          NMEC  = NMEC +1
183          LMEC(NMEC) = 22
184          GO TO 999
185       ENDIF
186 *
187 * *** Update interaction probabilities
188 *
189       IF (IHADR.GT.0) ZINTHA = ZINTHA*(1.-STEP/SHADR)
190 *
191 * *** apply the selected process if any
192 *
193    40 IF (IPROC.EQ.0) GO TO 999
194       NMEC = NMEC +1
195       LMEC(NMEC) = IPROC
196 *
197 *  **   Hadron interaction ?
198 *
199       IF (IPROC.EQ.12) THEN
200 #if !defined(CERNLIB_USRJMP)
201          CALL GUHADR
202 #endif
203 #if defined(CERNLIB_USRJMP)
204          CALL JUMPT0(JUHADR)
205 #endif
206 *          Check time cut-off for decays at rest
207          IF (LMEC(NMEC).EQ.5) THEN
208             TOFG   = TOFG +SUMLIF/CLIGHT
209             SUMLIF = 0.
210             IF (TOFG.GE.TOFMAX) THEN
211                ISTOP = 4
212                LMEC(NMEC) = 22
213                NGKINE = 0
214             ENDIF
215          ENDIF
216 *
217 *  **   Decay ?
218 *
219       ELSE IF (IPROC.EQ.5) THEN
220          ISTOP = 1
221          CALL GDECAY
222       ENDIF
223 *                                                             END GTNEUT
224   999 END