This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / erdecks / ertrch.F
1 *CMZ :          24/11/95  16.28.16  by  S.Ravndal
2 *CMZ :  3.21/02 29/03/94  15.41.49  by  S.Giani
3 *-- Author :
4       SUBROUTINE ERTRCH
5 C.
6 C.    ******************************************************************
7 C.    *                                                                *
8 C.    *    Average charged track is extrapolated by one step           *
9 C.    *                                                                *
10 C.    *    ==>Called by : ERTRGO                                       *
11 C.    *       Original routine : GTHADR                                *
12 C.    *       Authors   M.Maire, E.Nagy  *********                     *
13 C.    *                                                                *
14 C.    ******************************************************************
15 C.
16 #include "geant321/gcbank.inc"
17 #include "geant321/gconst.inc"
18 #include "geant321/gccuts.inc"
19 #include "geant321/gcphys.inc"
20 #include "geant321/gcjloc.inc"
21 #include "geant321/gckine.inc"
22 #include "geant321/gcmate.inc"
23 #include "geant321/gcmulo.inc"
24 #include "geant321/gctmed.inc"
25 #include "geant321/gctrak.inc"
26 #include "geant321/gcunit.inc"
27 #include "geant321/ertrio.inc"
28 #include "geant321/erwork.inc"
29
30 #if defined(CERNLIB_SINGLE)
31       PARAMETER (EPSMAC=1.E-11)
32 #endif
33 #if !defined(CERNLIB_SINGLE)&&!defined(CERNLIB_IBM)
34       PARAMETER (EPSMAC=5.E-6)
35 #endif
36 #if !defined(CERNLIB_SINGLE)&&defined(CERNLIB_IBM)
37       PARAMETER (EPSMAC=5.E-5)
38 #endif
39 #if !defined(CERNLIB_SINGLE)
40       DOUBLE PRECISION GKR,DEMEAN,STOPP1,STOPP2,STOPMX,STOPRG,STOPC
41       DOUBLE PRECISION EKIPR
42 #endif
43       REAL VNEXT(6)
44       SAVE CFLD,CHARG2,RMASS,CUTPRO,IKCUT,STOPC
45 C.
46 C.    ------------------------------------------------------------------
47 *
48 *
49 * *** Update local pointers if medium has changed
50 *
51       IF (IUPD.EQ.0) THEN
52          IUPD  = 1
53          CHARG2 = CHARGE*CHARGE
54          IF (IPART.LE.3) THEN
55             CUTEK  = CUTELE
56             RMASS  = 1.
57             JRANG  = LQ(JMA-15)
58          ELSE IF (IPART.LE.6) THEN
59             CUTEK  = CUTMUO
60             RMASS  = 1.
61             JRANG  = LQ(JMA-16)
62          ELSE
63             CUTEK  = CUTHAD
64             RMASS  = PMASS/AMASS
65             JRANG  = LQ(JMA-16) + NEK1
66          ENDIF
67          CUTPRO = MAX(CUTEK*RMASS,ELOW(1))
68          IKCUT = GEKA*LOG10(CUTPRO) + GEKB
69          GKR   = (CUTPRO - ELOW(IKCUT))/(ELOW(IKCUT+1) - ELOW(IKCUT))
70          STOPC = (1.-GKR)*Q(JRANG+IKCUT) + GKR*Q(JRANG+IKCUT+1)
71          CFLD  = 0.
72          IF (FIELDM.NE.0.) CFLD = 3333.*DEGRAD*TMAXFD/ABS(FIELDM*CHARGE)
73       ENDIF
74 *
75 * *** Compute current step size
76 *
77       STEP   = BIG
78       GEKRT1 = 1. - GEKRAT
79 *
80 * *** Step limitation due to energy loss (stopping range) ?
81 *
82       IF (ILOSS*DEEMAX.GT.0.) THEN
83          STOPP1 = GEKRT1*Q(JRANG+IEKBIN) + GEKRAT*Q(JRANG+IEKBIN+1)
84          STOPMX = (STOPP1 - STOPC)/(RMASS*CHARG2)
85          EKF  = (1. - BACKTR*DEEMAX)*GEKIN*RMASS
86          IF (EKF.LT.ELOW(1))    EKF = ELOW(1)
87          IF (EKF.GE.ELOW(NEK1)) EKF = ELOW(NEK1)*0.99
88          IKF=GEKA*LOG10(EKF)+GEKB
89          GKR=(EKF-ELOW(IKF))/(ELOW(IKF+1)-ELOW(IKF))
90          STOPP2 = (1.-GKR)*Q(JRANG+IKF) + GKR*Q(JRANG+IKF+1)
91          SLOSP  = ABS (STOPP1 - STOPP2)
92          STEP   = SLOSP/(RMASS*CHARG2)
93       ENDIF
94 *
95 * *** Step limitation due to energy loss in magnetic field ?
96 *
97       IF (IFIELD*FIELDM.NE.0.) THEN
98          SFIELD = CFLD*VECT(7)
99          IF (SFIELD.LT.STEP) STEP = SFIELD
100       ENDIF
101 *
102 * *** Compute point where to store error matrix
103 *
104       LERST  = 0
105       STEPER = BIG
106       ASCL1  = BIG
107       DO 20 IPR = 1,NEPRED
108          STEPE  = BIG
109          IF (LELENG) STEPE = ERLENG(IPR) - SLENG
110          IF (LEPLAN) THEN
111             SCAL1 = 0.
112             SCAL2 = 0.
113             DO 18 I=1,3
114                SCAL1 = SCAL1 + ERPLO(I,4,IPR)*(ERPLO(I,3,IPR)-VECT(I))
115                SCAL2 = SCAL2 + ERPLO(I,4,IPR)*VECT(I+3)
116    18       CONTINUE
117             STEPE = SCAL1/SCAL2
118          ENDIF
119          IF (STEPE.LE.PREC) STEPE = BIG
120          IF (STEPE.LT.STEPER) THEN
121             STEPER = STEPE
122             INLIST = IPR
123             IF (LEPLAN) ASCL1  = ABS (SCAL1)
124          ENDIF
125    20 CONTINUE
126       IF (STEPER.LE.STEP)  THEN
127          STEP  = STEPER
128          LERST = 1
129       ENDIF
130 *
131 * *** Step limitation due to geometry ?
132 *
133       LNEXT  = 0
134       IF (STEP.GE.0.95*SAFETY) THEN
135          CALL GTNEXT
136          IF (IGNEXT.NE.0) THEN
137             STEP   = SNEXT + PREC
138             LNEXT = 1
139             IF ((STEPER-SNEXT).GT.(2*PREC)) LERST = 0
140          ENDIF
141       ENDIF
142 *
143 * *** Linear transport when no field or very short step
144 *
145       IF (IFIELD.EQ.0.OR.STEP.LE.2*PREC) THEN
146         IF (IGNEXT.NE.0) THEN
147           DO 25 I = 1,3
148             VECTMP = VECT(I) +STEP*VECT(I+3)
149             IF(VECTMP.EQ.VECT(I)) THEN
150 *
151 * *** Correct for machine precision
152 *
153                   IF(VECT(I+3).NE.0.) THEN
154                      VECTMP =
155      +               VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
156 #if defined(DEBUG)
157                      WRITE(CHMAIL, 10000)
158                      CALL GMAIL(0,0)
159                      WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
160                      CALL GMAIL(0,0)
161 10000 FORMAT(' Boundary correction in ERTRCH: ',
162      +       '    GEKIN      NUMED       STEP      SNEXT')
163 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
164 #endif
165                   ENDIF
166             ENDIF
167             VOUT(I) = VECTMP
168    25     CONTINUE
169             INWVOL = 2
170             NMEC = NMEC +1
171             LMEC(NMEC) = 1
172         ELSE
173             DO 30 I = 1,3
174                VOUT(I)  = VECT(I) +STEP*VECT(I+3)
175    30       CONTINUE
176         ENDIF
177         DO 35 I = 4,6
178           VOUT(I)  = VECT(I)
179    35   CONTINUE
180         GOTO 74
181       END IF
182 *
183 * *** otherwise, swim particle in magnetic field
184 *
185       NMEC = NMEC +1
186       LMEC(NMEC) = 4
187 *
188    50 LERST = 0
189       LNEXT = 0
190       CALL GUSWIM (CHTR , STEP, VECT, VOUT)
191 *
192 *     When near to boundary, take proper action (cut-step,crossing...)
193       IF (STEP.GE.SAFETY) THEN
194          INEAR = 0
195          IF (IGNEXT.NE.0) THEN
196            DO 51 I = 1,3
197                VNEXT(I+3) = VECT(I+3)
198                VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
199    51      CONTINUE
200            DO 52 I = 1,3
201              IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 55
202    52      CONTINUE
203            INEAR = 1
204          ENDIF
205 *
206    55    CALL GINVOL (VOUT,ISAME)
207          IF (ISAME.EQ.0) THEN
208            IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
209              INWVOL = 2
210              NMEC = NMEC +1
211              LMEC(NMEC) = 1
212              LNEXT = 1
213            ELSE
214 *              Cut step
215              STEP = 0.5*STEP
216              IF (LMEC(NMEC).NE.24) THEN
217                NMEC = NMEC +1
218                LMEC(NMEC) = 24
219              ENDIF
220              GOTO 50
221            ENDIF
222          ENDIF
223       ENDIF
224 *
225 *
226 *     preset plane reached  ?
227    74 CONTINUE
228       IF ((LEPLAN).AND.(STEP.GE.ASCL1)) THEN
229          SCAL3 = 0.
230          DO 28 I=1,3
231             SCAL3=SCAL3+ERPLO(I,4,INLIST)*(ERPLO(I,3,INLIST)-VOUT(I))
232    28    CONTINUE
233          ASCL3 = ABS(SCAL3)
234          SSCL1 = ASCL1/SCAL1
235          IF (SCAL3*SSCL1.LT. -PREC) THEN
236 *            Cut step
237              STEP  = STEP*(ASCL1/(ASCL1+ASCL3))
238              NMEC  = NMEC +1
239              LMEC(NMEC) = 24
240              GOTO 50
241          ELSE
242            IF(ASCL3.LE.PREC) LERST = 1
243          ENDIF
244       ENDIF
245 *
246       DO 75 I=1,6
247            VECT(I) = VOUT(I)
248    75 CONTINUE
249 *
250       IF (LELENG.AND.(STEP.GE.STEPER)) LERST = 1
251 *
252       SLENG = SLENG + STEP
253 *
254 * *** Now apply selected mechanisms
255 *
256       IF (LNEXT.EQ.1) THEN
257           INWVOL = 2
258           NMEC = NMEC + 1
259           LMEC(NMEC) = 1
260       ENDIF
261 *
262 * *** apply energy loss : find the kinetic energy corresponding
263 *      to the new stopping range = stopmx -/+ step
264 *      (take care of the back tracking !)
265 *
266       IF (ILOSS*DEEMAX.GT.0) THEN
267          NMEC = NMEC +1
268          LMEC(NMEC) = 3
269          CALL ERLAND (STEP,Z,A,DENS,VECT(7),GETOT,AMASS,DEDX2)
270          DEDX2  = DEDX2*CHARG2*CHARG2
271          STOPRG = STOPP1 - BACKTR*STEP*RMASS*CHARG2
272          IKF = IEKBIN
273          IF (BACKTR.LE.0.) THEN
274    95       IF (STOPRG.LT.Q(JRANG+IKF)) THEN
275                IKF = IKF - 1
276                IF (IKF.GT.1) GO TO 95
277             ENDIF
278          ELSE
279    96       IF (STOPRG.GE.Q(JRANG+IKF+1)) THEN
280                IKF = IKF + 1
281                IF (IKF.LT.NEK1) GO TO 96
282             ENDIF
283          ENDIF
284          GKR = (STOPRG - Q(JRANG+IKF)) / (Q(JRANG+IKF+1) - Q(JRANG+IKF))
285          EKIPR = (1. -GKR)*ELOW(IKF) + GKR*ELOW(IKF+1)
286          GEKINT =  EKIPR/RMASS
287          IF (GEKINT.GT.CUTEK) THEN
288             DESTEP = ABS (GEKIN - GEKINT)
289             GEKIN  = GEKINT
290             GETOT  = GEKIN + AMASS
291             VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
292             CALL GEKBIN
293          ELSE
294             DESTEP = GEKINT
295             GEKIN  = 0.
296             GETOT  = AMASS
297             VECT(7)= 0.
298             INWVOL = 0
299             ISTOP  = 2
300             NMEC = NMEC + 1
301             LMEC(NMEC) = 30
302          ENDIF
303       ENDIF
304 *
305 * *** Propagate error matrix
306 *
307       IF (.NOT. LEONLY) CALL ERPROP
308 *
309 * *** Store informations
310 *
311       IF(LERST.EQ.1) THEN
312          NMEC = NMEC + 1
313          LMEC(NMEC) = 27
314          CALL ERSTOR
315       ENDIF
316 *
317       END