]>
Commit | Line | Data |
---|---|---|
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 | ||
9 | Cmodif: SUBROUTINE GHELIX (CHARGE, STEP, VECT, VOUT) changed into: | |
10 | SUBROUTINE EXTRAP_ONESTEP_HELIX (CHARGE, STEP, VECT, VOUT) | |
11 | C. | |
12 | C. ****************************************************************** | |
13 | C. * * | |
14 | C. * Performs the tracking of one step in a magnetic field * | |
15 | C. * The trajectory is assumed to be a helix in a constant field * | |
16 | C. * taken at the mid point of the step. * | |
17 | C. * Parameters: * | |
18 | C. * input * | |
19 | C. * STEP =arc length of the step asked * | |
20 | C. * VECT =input vector (position,direction cos and momentum) * | |
21 | C. * CHARGE= electric charge of the particle * | |
22 | C. * output * | |
23 | C. * VOUT = same as VECT after completion of the step * | |
24 | C. * * | |
25 | C. * ==>Called by : <USER>, GUSWIM * | |
26 | C. * Author M.Hansroul ********* * | |
27 | C. * Modified S.Egli, S.V.Levonian * | |
28 | C. * Modified V.Perevoztchikov | |
29 | C. * * | |
30 | C. ****************************************************************** | |
31 | C. | |
32 | ||
33 | Cmodif: 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) | |
41 | C. | |
42 | C. ------------------------------------------------------------------ | |
43 | C. | |
44 | C units are kgauss,centimeters,gev/c | |
45 | C | |
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) | |
51 | C | |
52 | Cmodif: 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 | |
59 | Cmodif: 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 | |
111 | C | |
112 | 999 END | |
113 | ||
114 | Cmodif: SUBROUTINE GHELX3 (FIELD, STEP, VECT, VOUT) changed into: | |
115 | SUBROUTINE EXTRAP_ONESTEP_HELIX3 (FIELD, STEP, VECT, VOUT) | |
116 | C. | |
117 | C. ****************************************************************** | |
118 | C. * * | |
119 | C. * Tracking routine in a constant field oriented * | |
120 | C. * along axis 3 * | |
121 | C. * Tracking is performed with a conventional * | |
122 | C. * helix step method * | |
123 | C. * * | |
124 | C. * ==>Called by : <USER>, GUSWIM * | |
125 | C. * Authors R.Brun, M.Hansroul ********* * | |
126 | C * Rewritten V.Perevoztchikov | |
127 | C. * * | |
128 | C. ****************************************************************** | |
129 | C. | |
130 | ||
131 | Cmodif: 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) | |
138 | C. | |
139 | C. ------------------------------------------------------------------ | |
140 | C. | |
141 | C units are kgauss,centimeters,gev/c | |
142 | C | |
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 | ||
180 | C | |
181 | 999 END | |
182 | ||
183 | Cmodif: SUBROUTINE GRKUTA (CHARGE,STEP,VECT,VOUT) changed into: | |
184 | SUBROUTINE EXTRAP_ONESTEP_RUNGEKUTTA (CHARGE,STEP,VECT,VOUT) | |
185 | C. | |
186 | C. ****************************************************************** | |
187 | C. * * | |
188 | C. * Runge-Kutta method for tracking a particle through a magnetic * | |
189 | C. * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of * | |
190 | C. * Standards, procedure 25.5.20) * | |
191 | C. * * | |
192 | C. * Input parameters * | |
193 | C. * CHARGE Particle charge * | |
194 | C. * STEP Step size * | |
195 | C. * VECT Initial co-ords,direction cosines,momentum * | |
196 | C. * Output parameters * | |
197 | C. * VOUT Output co-ords,direction cosines,momentum * | |
198 | C. * User routine called * | |
199 | C. * CALL GUFLD(X,F) * | |
200 | C. * * | |
201 | C. * ==>Called by : <USER>, GUSWIM * | |
202 | C. * Authors R.Brun, M.Hansroul ********* * | |
203 | C. * V.Perevoztchikov (CUT STEP implementation) * | |
204 | C. * * | |
205 | C. * * | |
206 | C. ****************************************************************** | |
207 | C. | |
208 | Cmodif: no condition from CERNLIB_SINGLE for double precision | |
209 | IMPLICIT DOUBLE PRECISION(A-H,O-Z) | |
210 | Cmodif: 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 | |
240 | Cmodif: 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 | ||
273 | Cmodif: 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 | ||
301 | Cmodif: 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 |