]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/extrap.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Cosmetics (Chrsitian)
[u/mrichter/AliRoot.git] / MUON / extrap.F
CommitLineData
0df6418f 1* One step extrapolation routines in Fortran:
2* EXTRAP_ONESTEP_HELIX, EXTRAP_ONESTEP_HELIX3, EXTRAP_ONESTEP_RUNGEKUTTA.
3* Taken respectively from Geant (gtrak) routines:
4* GHELIX, GHELIX3, GRKUTA.
5* Everything in double precision,
6* in order that the track fit with Minuit is converging.
7* Modifications from Geant are indicated with "Cmodif".
8
9Cmodif: SUBROUTINE GHELIX (CHARGE, STEP, VECT, VOUT) changed into:
10 SUBROUTINE EXTRAP_ONESTEP_HELIX (CHARGE, STEP, VECT, VOUT)
11C.
12C. ******************************************************************
13C. * *
14C. * Performs the tracking of one step in a magnetic field *
15C. * The trajectory is assumed to be a helix in a constant field *
16C. * taken at the mid point of the step. *
17C. * Parameters: *
18C. * input *
19C. * STEP =arc length of the step asked *
20C. * VECT =input vector (position,direction cos and momentum) *
21C. * CHARGE= electric charge of the particle *
22C. * output *
23C. * VOUT = same as VECT after completion of the step *
24C. * *
25C. * ==>Called by : <USER>, GUSWIM *
26C. * Author M.Hansroul ********* *
27C. * Modified S.Egli, S.V.Levonian *
28C. * Modified V.Perevoztchikov
29C. * *
30C. ******************************************************************
31C.
32
33Cmodif: everything in double precision
34 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
35
36 DIMENSION VECT(7),VOUT(7)
37 DIMENSION XYZ(3),H(4),HXP(3)
38 PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6,IPP=7)
39 PARAMETER (SIXTH = 1./6.)
40 PARAMETER (EC=2.9979251E-4)
41C.
42C. ------------------------------------------------------------------
43C.
44C units are kgauss,centimeters,gev/c
45C
46 VOUT(IPP) = VECT(IPP)
47 IF (CHARGE.EQ.0.) GO TO 10
48 XYZ(1) = VECT(IX) + 0.5 * STEP * VECT(IPX)
49 XYZ(2) = VECT(IY) + 0.5 * STEP * VECT(IPY)
50 XYZ(3) = VECT(IZ) + 0.5 * STEP * VECT(IPZ)
51C
52Cmodif: CALL GUFLD (XYZ, H) changed into:
53 CALL GUFLD_DOUBLE (XYZ, H)
54
55 H2XY = H(1)**2 + H(2)**2
56 H(4) = H(3)**2 + H2XY
57 IF (H(4).LE.1.E-12) GO TO 10
58 IF (H2XY.LE.1.E-12*H(4)) THEN
59Cmodif: CALL GHELX3 (CHARGE*H(3), STEP, VECT, VOUT) changed into:
60 CALL EXTRAP_ONESTEP_HELIX3 (CHARGE*H(3), STEP, VECT, VOUT)
61 GO TO 999
62 ENDIF
63 H(4) = SQRT(H(4))
64 H(1) = H(1)/H(4)
65 H(2) = H(2)/H(4)
66 H(3) = H(3)/H(4)
67 H(4) = H(4)*EC
68*
69 HXP(1) = H(2)*VECT(IPZ) - H(3)*VECT(IPY)
70 HXP(2) = H(3)*VECT(IPX) - H(1)*VECT(IPZ)
71 HXP(3) = H(1)*VECT(IPY) - H(2)*VECT(IPX)
72
73 HP = H(1)*VECT(IPX) + H(2)*VECT(IPY) + H(3)*VECT(IPZ)
74*
75 RHO = -CHARGE*H(4)/VECT(IPP)
76 TET = RHO * STEP
77 IF (ABS(TET).GT.0.15) THEN
78 SINT = SIN(TET)
79 SINTT = (SINT/TET)
80 TSINT = (TET-SINT)/TET
81 COS1T = 2.*(SIN(0.5*TET))**2/TET
82 ELSE
83 TSINT = SIXTH*TET**2
84 SINTT = (1. - TSINT)
85 SINT = TET*SINTT
86 COS1T = 0.5*TET
87 ENDIF
88*
89 F1 = STEP * SINTT
90 F2 = STEP * COS1T
91 F3 = STEP * TSINT * HP
92 F4 = -TET*COS1T
93 F5 = SINT
94 F6 = TET * COS1T * HP
95
96 VOUT(IX) = VECT(IX) + (F1*VECT(IPX) + F2*HXP(1) + F3*H(1))
97 VOUT(IY) = VECT(IY) + (F1*VECT(IPY) + F2*HXP(2) + F3*H(2))
98 VOUT(IZ) = VECT(IZ) + (F1*VECT(IPZ) + F2*HXP(3) + F3*H(3))
99
100 VOUT(IPX) = VECT(IPX) + (F4*VECT(IPX) + F5*HXP(1) + F6*H(1))
101 VOUT(IPY) = VECT(IPY) + (F4*VECT(IPY) + F5*HXP(2) + F6*H(2))
102 VOUT(IPZ) = VECT(IPZ) + (F4*VECT(IPZ) + F5*HXP(3) + F6*H(3))
103
104 GO TO 999
105
106 10 CONTINUE
107 DO 20 I = 1,3
108 VOUT(I) = VECT(I) + STEP * VECT(I+3)
109 VOUT(I+3) = VECT(I+3)
110 20 CONTINUE
111C
112 999 END
113
114Cmodif: SUBROUTINE GHELX3 (FIELD, STEP, VECT, VOUT) changed into:
115 SUBROUTINE EXTRAP_ONESTEP_HELIX3 (FIELD, STEP, VECT, VOUT)
116C.
117C. ******************************************************************
118C. * *
119C. * Tracking routine in a constant field oriented *
120C. * along axis 3 *
121C. * Tracking is performed with a conventional *
122C. * helix step method *
123C. * *
124C. * ==>Called by : <USER>, GUSWIM *
125C. * Authors R.Brun, M.Hansroul ********* *
126C * Rewritten V.Perevoztchikov
127C. * *
128C. ******************************************************************
129C.
130
131Cmodif: everything in double precision
132 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
133
134 DIMENSION VECT(7),VOUT(7),HXP(3)
135 PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6,IPP=7)
136 PARAMETER (SIXTH = 1./6.)
137 PARAMETER (EC=2.9979251E-4)
138C.
139C. ------------------------------------------------------------------
140C.
141C units are kgauss,centimeters,gev/c
142C
143 VOUT(IPP) = VECT(IPP)
144 H4 = FIELD * EC
145*
146 HXP(1) = - VECT(IPY)
147 HXP(2) = + VECT(IPX)
148
149 HP = VECT(IPZ)
150*
151 RHO = -H4/VECT(IPP)
152 TET = RHO * STEP
153 IF (ABS(TET).GT.0.15) THEN
154 SINT = SIN(TET)
155 SINTT = (SINT/TET)
156 TSINT = (TET-SINT)/TET
157 COS1T = 2.*(SIN(0.5*TET))**2/TET
158 ELSE
159 TSINT = SIXTH*TET**2
160 SINTT = (1. - TSINT)
161 SINT = TET*SINTT
162 COS1T = 0.5*TET
163 ENDIF
164*
165 F1 = STEP * SINTT
166 F2 = STEP * COS1T
167 F3 = STEP * TSINT * HP
168 F4 = -TET*COS1T
169 F5 = SINT
170 F6 = TET * COS1T * HP
171
172 VOUT(IX) = VECT(IX) + (F1*VECT(IPX) + F2*HXP(1))
173 VOUT(IY) = VECT(IY) + (F1*VECT(IPY) + F2*HXP(2))
174 VOUT(IZ) = VECT(IZ) + (F1*VECT(IPZ) + F3)
175
176 VOUT(IPX) = VECT(IPX) + (F4*VECT(IPX) + F5*HXP(1))
177 VOUT(IPY) = VECT(IPY) + (F4*VECT(IPY) + F5*HXP(2))
178 VOUT(IPZ) = VECT(IPZ) + (F4*VECT(IPZ) + F6)
179
180C
181 999 END
182
183Cmodif: SUBROUTINE GRKUTA (CHARGE,STEP,VECT,VOUT) changed into:
184 SUBROUTINE EXTRAP_ONESTEP_RUNGEKUTTA (CHARGE,STEP,VECT,VOUT)
185C.
186C. ******************************************************************
187C. * *
188C. * Runge-Kutta method for tracking a particle through a magnetic *
189C. * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
190C. * Standards, procedure 25.5.20) *
191C. * *
192C. * Input parameters *
193C. * CHARGE Particle charge *
194C. * STEP Step size *
195C. * VECT Initial co-ords,direction cosines,momentum *
196C. * Output parameters *
197C. * VOUT Output co-ords,direction cosines,momentum *
198C. * User routine called *
199C. * CALL GUFLD(X,F) *
200C. * *
201C. * ==>Called by : <USER>, GUSWIM *
202C. * Authors R.Brun, M.Hansroul ********* *
203C. * V.Perevoztchikov (CUT STEP implementation) *
204C. * *
205C. * *
206C. ******************************************************************
207C.
208Cmodif: no condition from CERNLIB_SINGLE for double precision
209 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
210Cmodif: REAL changed into DOUBLE PRECISION in following 2 lines
211 DOUBLE PRECISION CHARGE, STEP, VECT(*), VOUT(*), F(4)
212 DOUBLE PRECISION XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT
213 DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3)
214 EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
215 + (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
216*
217 PARAMETER (MAXIT = 1992, MAXCUT = 11)
218 PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32)
219 PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
220 PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
221 PARAMETER (PISQUA=.986960440109D+01)
222 PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
223*.
224*. ------------------------------------------------------------------
225*.
226* This constant is for units CM,GEV/C and KGAUSS
227*
228 ITER = 0
229 NCUT = 0
230 DO 10 J=1,7
231 VOUT(J)=VECT(J)
232 10 CONTINUE
233 PINV = EC * CHARGE / VECT(7)
234 TL = 0.
235 H = STEP
236*
237*
238 20 REST = STEP-TL
239 IF (ABS(H).GT.ABS(REST)) H = REST
240Cmodif: CALL GUFLD(VOUT,F) changed into:
241 CALL GUFLD_DOUBLE(VOUT,F)
242*
243* Start of integration
244*
245 X = VOUT(1)
246 Y = VOUT(2)
247 Z = VOUT(3)
248 A = VOUT(4)
249 B = VOUT(5)
250 C = VOUT(6)
251*
252 H2 = HALF * H
253 H4 = HALF * H2
254 PH = PINV * H
255 PH2 = HALF * PH
256 SECXS(1) = (B * F(3) - C * F(2)) * PH2
257 SECYS(1) = (C * F(1) - A * F(3)) * PH2
258 SECZS(1) = (A * F(2) - B * F(1)) * PH2
259 ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2)
260 IF (ANG2.GT.PISQUA) GO TO 40
261 DXT = H2 * A + H4 * SECXS(1)
262 DYT = H2 * B + H4 * SECYS(1)
263 DZT = H2 * C + H4 * SECZS(1)
264 XT = X + DXT
265 YT = Y + DYT
266 ZT = Z + DZT
267*
268* Second intermediate point
269*
270 EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
271 IF (EST.GT.H) GO TO 30
272
273Cmodif: CALL GUFLD(XYZT,F) changed into:
274 CALL GUFLD_DOUBLE(XYZT,F)
275 AT = A + SECXS(1)
276 BT = B + SECYS(1)
277 CT = C + SECZS(1)
278*
279 SECXS(2) = (BT * F(3) - CT * F(2)) * PH2
280 SECYS(2) = (CT * F(1) - AT * F(3)) * PH2
281 SECZS(2) = (AT * F(2) - BT * F(1)) * PH2
282 AT = A + SECXS(2)
283 BT = B + SECYS(2)
284 CT = C + SECZS(2)
285 SECXS(3) = (BT * F(3) - CT * F(2)) * PH2
286 SECYS(3) = (CT * F(1) - AT * F(3)) * PH2
287 SECZS(3) = (AT * F(2) - BT * F(1)) * PH2
288 DXT = H * (A + SECXS(3))
289 DYT = H * (B + SECYS(3))
290 DZT = H * (C + SECZS(3))
291 XT = X + DXT
292 YT = Y + DYT
293 ZT = Z + DZT
294 AT = A + TWO*SECXS(3)
295 BT = B + TWO*SECYS(3)
296 CT = C + TWO*SECZS(3)
297*
298 EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
299 IF (EST.GT.2.*ABS(H)) GO TO 30
300
301Cmodif: CALL GUFLD(XYZT,F) changed into:
302 CALL GUFLD_DOUBLE(XYZT,F)
303*
304 Z = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
305 Y = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H
306 X = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H
307*
308 SECXS(4) = (BT*F(3) - CT*F(2))* PH2
309 SECYS(4) = (CT*F(1) - AT*F(3))* PH2
310 SECZS(4) = (AT*F(2) - BT*F(1))* PH2
311 A = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD
312 B = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD
313 C = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD
314*
315 EST = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3)))
316 ++ ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3)))
317 ++ ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3)))
318*
319 IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30
320 ITER = ITER + 1
321 NCUT = 0
322* If too many iterations, go to HELIX
323 IF (ITER.GT.MAXIT) GO TO 40
324*
325 TL = TL + H
326 IF (EST.LT.(DLT32)) THEN
327 H = H*TWO
328 ENDIF
329 CBA = ONE/ SQRT(A*A + B*B + C*C)
330 VOUT(1) = X
331 VOUT(2) = Y
332 VOUT(3) = Z
333 VOUT(4) = CBA*A
334 VOUT(5) = CBA*B
335 VOUT(6) = CBA*C
336 REST = STEP - TL
337 IF (STEP.LT.0.) REST = -REST
338 IF (REST .GT. 1.E-5*ABS(STEP)) GO TO 20
339*
340 GO TO 999
341*
342** CUT STEP
343 30 NCUT = NCUT + 1
344* If too many cuts , go to HELIX
345 IF (NCUT.GT.MAXCUT) GO TO 40
346 H = H*HALF
347 GO TO 20
348*
349** ANGLE TOO BIG, USE HELIX
350 40 F1 = F(1)
351 F2 = F(2)
352 F3 = F(3)
353 F4 = SQRT(F1**2+F2**2+F3**2)
354 RHO = -F4*PINV
355 TET = RHO * STEP
356 IF(TET.NE.0.) THEN
357 HNORM = ONE/F4
358 F1 = F1*HNORM
359 F2 = F2*HNORM
360 F3 = F3*HNORM
361*
362 HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY)
363 HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ)
364 HXP(3) = F1*VECT(IPY) - F2*VECT(IPX)
365
366 HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ)
367*
368 RHO1 = ONE/RHO
369 SINT = SIN(TET)
370 COST = TWO*SIN(HALF*TET)**2
371*
372 G1 = SINT*RHO1
373 G2 = COST*RHO1
374 G3 = (TET-SINT) * HP*RHO1
375 G4 = -COST
376 G5 = SINT
377 G6 = COST * HP
378
379 VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1)
380 VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2)
381 VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3)
382
383 VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1)
384 VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2)
385 VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3)
386*
387 ELSE
388 VOUT(IX) = VECT(IX) + STEP*VECT(IPX)
389 VOUT(IY) = VECT(IY) + STEP*VECT(IPY)
390 VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ)
391*
392 ENDIF
393*
394 999 END