]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/erdecks/ertrch.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix needed on Sun and Alpha
[u/mrichter/AliRoot.git] / GEANT321 / erdecks / ertrch.F
CommitLineData
fe4da5cc 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
5C.
6C. ******************************************************************
7C. * *
8C. * Average charged track is extrapolated by one step *
9C. * *
10C. * ==>Called by : ERTRGO *
11C. * Original routine : GTHADR *
12C. * Authors M.Maire, E.Nagy ********* *
13C. * *
14C. ******************************************************************
15C.
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
45C.
46C. ------------------------------------------------------------------
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)
16110000 FORMAT(' Boundary correction in ERTRCH: ',
162 + ' GEKIN NUMED STEP SNEXT')
16310100 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