1 |
|
2 |
|
3 |
************************************************************* |
4 |
* |
5 |
* Routines to compute the NPOINT track intersection points |
6 |
* with planes of z-coordinate given by ZIN |
7 |
* given the track parameters |
8 |
* |
9 |
* The routine is based on GRKUTA, which computes the |
10 |
* trajectory of a charged particle in a magnetic field |
11 |
* by solving the equations of motion with Runge-Kuta method |
12 |
* |
13 |
* Variables that have to be assigned when the subroutine |
14 |
* is called: |
15 |
* |
16 |
* ZIN(1-NPOINT) ----> z coordinates of the planes |
17 |
* AL_P(1-5) ----> track-parameter vector |
18 |
* |
19 |
* NB !!! |
20 |
* The routine works properly only if the |
21 |
* planes are numbered in descending order |
22 |
* |
23 |
* Routines: |
24 |
* |
25 |
* DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL) |
26 |
* - old routine |
27 |
* |
28 |
* DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL) |
29 |
* - as the previous one, but: |
30 |
* -- the projected angles are given as output |
31 |
* -- the track lengths are given as output: |
32 |
* ---- for planes above the reference plane, the length until |
33 |
* the lower closest plane (reference plane included) |
34 |
* ---- for planes below the reference plane, the length until |
35 |
* the higher closest plane (reference plane included) |
36 |
* |
37 |
* March 2008 --> optimized by Paolo |
38 |
* |
39 |
* DOTRACK3(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,ZINI,IFAIL) |
40 |
* - as the previous one, but with the coordinate ZINI of the reference |
41 |
* plane given as input parameter |
42 |
* |
43 |
* |
44 |
************************************************************** |
45 |
|
46 |
|
47 |
SUBROUTINE |
48 |
$ DOTRACK3(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P |
49 |
$ ,ZINI,IFAIL) |
50 |
|
51 |
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
52 |
|
53 |
|
54 |
DIMENSION VECT(7),VECTINI(7),VOUT(7) |
55 |
DATA TOLL/1.d-8/ |
56 |
|
57 |
* tolerance in reaching the next plane during the tracking procedure |
58 |
* ----------------------------------------------- |
59 |
* I/O parameters |
60 |
PARAMETER (NPOINT_MAX=500) |
61 |
DIMENSION ZIN(NPOINT_MAX) |
62 |
DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX) |
63 |
DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX) |
64 |
DIMENSION TLOUT(NPOINT_MAX) |
65 |
DIMENSION AL_P(5) |
66 |
* ----------------------------------------------- |
67 |
* DATA ZINI/23.5/ !z coordinate of the reference plane |
68 |
REAL*8 L |
69 |
|
70 |
* ================================================================== |
71 |
* divide the track in two parts: below and above the reference plane |
72 |
* ================================================================== |
73 |
IUPDOWN=0 |
74 |
DO I=1,NPOINT |
75 |
IF(ZIN(I).LT.ZINI)THEN |
76 |
if(i.ne.1)IUPDOWN=I |
77 |
GOTO 88 |
78 |
ENDIF |
79 |
IUPDOWN=NPOINT+1 |
80 |
ENDDO |
81 |
88 CONTINUE |
82 |
|
83 |
* ================================================================== |
84 |
* track from reference plane DOWN |
85 |
* ================================================================== |
86 |
* parameters for GRKUTA |
87 |
IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5)) |
88 |
IF(AL_P(5).EQ.0) CHARGE=1. |
89 |
VOUT(1)=AL_P(1) |
90 |
VOUT(2)=AL_P(2) |
91 |
VOUT(3)=ZINI |
92 |
VOUT(4)=AL_P(3)*DCOS(AL_P(4)) |
93 |
VOUT(5)=AL_P(3)*DSIN(AL_P(4)) |
94 |
VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2) |
95 |
IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) |
96 |
IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
97 |
DO I=MAX(IUPDOWN,1),NPOINT |
98 |
L = 0.0 |
99 |
10 DO J=1,7 |
100 |
VECT(J)=VOUT(J) |
101 |
VECTINI(J)=VOUT(J) |
102 |
ENDDO |
103 |
IF(VOUT(6).GE.0.) THEN |
104 |
IFAIL=1 |
105 |
c$$$ if(TRKVERBOSE) |
106 |
c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' |
107 |
RETURN |
108 |
ENDIF |
109 |
step=(zin(i)-vect(3))/VOUT(6) |
110 |
11 continue |
111 |
CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
112 |
L = L + STEP |
113 |
IF(VOUT(3).GT.VECT(3)) THEN |
114 |
IFAIL=1 |
115 |
c$$$ if(TRKVERBOSE) |
116 |
c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' |
117 |
c$$$ if(TRKVERBOSE)print*,'charge',charge |
118 |
c$$$ if(TRKVERBOSE)print*,'vect',vect |
119 |
c$$$ if(TRKVERBOSE)print*,'vout',vout |
120 |
c$$$ if(TRKVERBOSE)print*,'step',step |
121 |
c$$$ if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1 |
122 |
RETURN |
123 |
ENDIF |
124 |
Z=VOUT(3) |
125 |
IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100 |
126 |
IF(Z.GT.ZIN(I)+TOLL) GOTO 10 |
127 |
IF(Z.LE.ZIN(I)-TOLL) THEN |
128 |
L = L - STEP |
129 |
STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) |
130 |
DO J=1,7 |
131 |
VECT(J)=VECTINI(J) |
132 |
ENDDO |
133 |
GOTO 11 |
134 |
ENDIF |
135 |
100 XOUT(I)=VOUT(1) |
136 |
YOUT(I)=VOUT(2) |
137 |
ZIN(I)=VOUT(3) |
138 |
IF(VOUT(3).ne.0)THEN |
139 |
THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.) |
140 |
THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) |
141 |
ENDIF |
142 |
TLOUT(I) = L |
143 |
c print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP |
144 |
ENDDO |
145 |
|
146 |
* ================================================================== |
147 |
* track from reference plane UP |
148 |
* ================================================================== |
149 |
* parameters for GRKUTA: |
150 |
* -opposite charge |
151 |
* -opposite momentum direction |
152 |
IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5)) |
153 |
IF(AL_P(5).EQ.0) CHARGE=-1. |
154 |
VOUT(1)=AL_P(1) |
155 |
VOUT(2)=AL_P(2) |
156 |
VOUT(3)=ZINI |
157 |
VOUT(4)=-AL_P(3)*DCOS(AL_P(4)) |
158 |
VOUT(5)=-AL_P(3)*DSIN(AL_P(4)) |
159 |
VOUT(6)=1.*DSQRT(1.-AL_P(3)**2) |
160 |
IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) |
161 |
IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
162 |
DO I=MIN((IUPDOWN-1),NPOINT),1,-1 |
163 |
L = 0.0 |
164 |
|
165 |
20 DO J=1,7 |
166 |
VECT(J)=VOUT(J) |
167 |
VECTINI(J)=VOUT(J) |
168 |
ENDDO |
169 |
IF(VOUT(6).LE.0.) THEN |
170 |
IFAIL=1 |
171 |
c$$$ if(TRKVERBOSE) |
172 |
c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' |
173 |
RETURN |
174 |
ENDIF |
175 |
step=(zin(i)-vect(3))/VOUT(6) |
176 |
22 continue |
177 |
CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
178 |
L = L + STEP |
179 |
IF(VOUT(3).LT.VECT(3)) THEN |
180 |
IFAIL=1 |
181 |
c$$$ if(TRKVERBOSE) |
182 |
c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' |
183 |
c$$$ if(TRKVERBOSE)print*,'charge',charge |
184 |
c$$$ if(TRKVERBOSE)print*,'vect',vect |
185 |
c$$$ if(TRKVERBOSE)print*,'vout',vout |
186 |
c$$$ if(TRKVERBOSE)print*,'step',step |
187 |
c$$$ if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1 |
188 |
RETURN |
189 |
ENDIF |
190 |
Z=VOUT(3) |
191 |
IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200 |
192 |
IF(Z.LT.ZIN(I)-TOLL) GOTO 20 |
193 |
IF(Z.GE.ZIN(I)+TOLL) THEN |
194 |
L = L - STEP |
195 |
STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) |
196 |
DO J=1,7 |
197 |
VECT(J)=VECTINI(J) |
198 |
ENDDO |
199 |
GOTO 22 |
200 |
ENDIF |
201 |
200 XOUT(I)=VOUT(1) |
202 |
YOUT(I)=VOUT(2) |
203 |
ZIN(I)=VOUT(3) |
204 |
IF(VOUT(3).ne.0)THEN |
205 |
THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.) |
206 |
THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) |
207 |
ENDIF |
208 |
TLOUT(I) = L |
209 |
ENDDO |
210 |
|
211 |
RETURN |
212 |
END |
213 |
************************************************************************ |
214 |
* |
215 |
* |
216 |
* |
217 |
* |
218 |
************************************************************************ |
219 |
|
220 |
SUBROUTINE |
221 |
$ DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL) |
222 |
|
223 |
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
224 |
|
225 |
|
226 |
ccc DIMENSION VECT(7),VECTINI(7),VOUT(7) |
227 |
ccc DATA TOLL/1.d-8/ |
228 |
|
229 |
* tolerance in reaching the next plane during the tracking procedure |
230 |
* ----------------------------------------------- |
231 |
* I/O parameters |
232 |
PARAMETER (NPOINT_MAX=500) |
233 |
DIMENSION ZIN(NPOINT_MAX) |
234 |
DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX) |
235 |
DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX) |
236 |
DIMENSION TLOUT(NPOINT_MAX) |
237 |
DIMENSION AL_P(5) |
238 |
* ----------------------------------------------- |
239 |
DATA ZINI/23.5/ !z coordinate of the reference plane |
240 |
ccc REAL*8 L |
241 |
|
242 |
CALL DOTRACK3(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P |
243 |
$ ,ZINI,IFAIL) |
244 |
|
245 |
ccc |
246 |
ccc OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD |
247 |
ccc |
248 |
|
249 |
c$$$* ================================================================== |
250 |
c$$$* divide the track in two parts: below and above the reference plane |
251 |
c$$$* ================================================================== |
252 |
c$$$ IUPDOWN=0 |
253 |
c$$$ DO I=1,NPOINT |
254 |
c$$$ IF(ZIN(I).LT.ZINI)THEN |
255 |
c$$$ if(i.ne.1)IUPDOWN=I |
256 |
c$$$ GOTO 88 |
257 |
c$$$ ENDIF |
258 |
c$$$ IUPDOWN=NPOINT+1 |
259 |
c$$$ ENDDO |
260 |
c$$$ 88 CONTINUE |
261 |
c$$$ |
262 |
c$$$* ================================================================== |
263 |
c$$$* track from reference plane DOWN |
264 |
c$$$* ================================================================== |
265 |
c$$$* parameters for GRKUTA |
266 |
c$$$ IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5)) |
267 |
c$$$ IF(AL_P(5).EQ.0) CHARGE=1. |
268 |
c$$$ VOUT(1)=AL_P(1) |
269 |
c$$$ VOUT(2)=AL_P(2) |
270 |
c$$$ VOUT(3)=ZINI |
271 |
c$$$ VOUT(4)=AL_P(3)*DCOS(AL_P(4)) |
272 |
c$$$ VOUT(5)=AL_P(3)*DSIN(AL_P(4)) |
273 |
c$$$ VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2) |
274 |
c$$$ IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) |
275 |
c$$$ IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
276 |
c$$$ DO I=MAX(IUPDOWN,1),NPOINT |
277 |
c$$$ L = 0.0 |
278 |
c$$$ 10 DO J=1,7 |
279 |
c$$$ VECT(J)=VOUT(J) |
280 |
c$$$ VECTINI(J)=VOUT(J) |
281 |
c$$$ ENDDO |
282 |
c$$$ IF(VOUT(6).GE.0.) THEN |
283 |
c$$$ IFAIL=1 |
284 |
c$$$c$$$ if(TRKVERBOSE) |
285 |
c$$$c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' |
286 |
c$$$ RETURN |
287 |
c$$$ ENDIF |
288 |
c$$$ step=(zin(i)-vect(3))/VOUT(6) |
289 |
c$$$ 11 continue |
290 |
c$$$ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
291 |
c$$$ L = L + STEP |
292 |
c$$$ IF(VOUT(3).GT.VECT(3)) THEN |
293 |
c$$$ IFAIL=1 |
294 |
c$$$c$$$ if(TRKVERBOSE) |
295 |
c$$$c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' |
296 |
c$$$c$$$ if(TRKVERBOSE)print*,'charge',charge |
297 |
c$$$c$$$ if(TRKVERBOSE)print*,'vect',vect |
298 |
c$$$c$$$ if(TRKVERBOSE)print*,'vout',vout |
299 |
c$$$c$$$ if(TRKVERBOSE)print*,'step',step |
300 |
c$$$c$$$ if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1 |
301 |
c$$$ RETURN |
302 |
c$$$ ENDIF |
303 |
c$$$ Z=VOUT(3) |
304 |
c$$$ IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100 |
305 |
c$$$ IF(Z.GT.ZIN(I)+TOLL) GOTO 10 |
306 |
c$$$ IF(Z.LE.ZIN(I)-TOLL) THEN |
307 |
c$$$ L = L - STEP |
308 |
c$$$ STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) |
309 |
c$$$ DO J=1,7 |
310 |
c$$$ VECT(J)=VECTINI(J) |
311 |
c$$$ ENDDO |
312 |
c$$$ GOTO 11 |
313 |
c$$$ ENDIF |
314 |
c$$$ 100 XOUT(I)=VOUT(1) |
315 |
c$$$ YOUT(I)=VOUT(2) |
316 |
c$$$ ZIN(I)=VOUT(3) |
317 |
c$$$ IF(VOUT(3).ne.0)THEN |
318 |
c$$$ THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.) |
319 |
c$$$ THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) |
320 |
c$$$ ENDIF |
321 |
c$$$ TLOUT(I) = L |
322 |
c$$$c print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP |
323 |
c$$$ ENDDO |
324 |
c$$$ |
325 |
c$$$* ================================================================== |
326 |
c$$$* track from reference plane UP |
327 |
c$$$* ================================================================== |
328 |
c$$$* parameters for GRKUTA: |
329 |
c$$$* -opposite charge |
330 |
c$$$* -opposite momentum direction |
331 |
c$$$ IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5)) |
332 |
c$$$ IF(AL_P(5).EQ.0) CHARGE=-1. |
333 |
c$$$ VOUT(1)=AL_P(1) |
334 |
c$$$ VOUT(2)=AL_P(2) |
335 |
c$$$ VOUT(3)=ZINI |
336 |
c$$$ VOUT(4)=-AL_P(3)*DCOS(AL_P(4)) |
337 |
c$$$ VOUT(5)=-AL_P(3)*DSIN(AL_P(4)) |
338 |
c$$$ VOUT(6)=1.*DSQRT(1.-AL_P(3)**2) |
339 |
c$$$ IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) |
340 |
c$$$ IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
341 |
c$$$ DO I=MIN((IUPDOWN-1),NPOINT),1,-1 |
342 |
c$$$ L = 0.0 |
343 |
c$$$ |
344 |
c$$$ 20 DO J=1,7 |
345 |
c$$$ VECT(J)=VOUT(J) |
346 |
c$$$ VECTINI(J)=VOUT(J) |
347 |
c$$$ ENDDO |
348 |
c$$$ IF(VOUT(6).LE.0.) THEN |
349 |
c$$$ IFAIL=1 |
350 |
c$$$c$$$ if(TRKVERBOSE) |
351 |
c$$$c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' |
352 |
c$$$ RETURN |
353 |
c$$$ ENDIF |
354 |
c$$$ step=(zin(i)-vect(3))/VOUT(6) |
355 |
c$$$ 22 continue |
356 |
c$$$ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
357 |
c$$$ L = L + STEP |
358 |
c$$$ IF(VOUT(3).LT.VECT(3)) THEN |
359 |
c$$$ IFAIL=1 |
360 |
c$$$c$$$ if(TRKVERBOSE) |
361 |
c$$$c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' |
362 |
c$$$c$$$ if(TRKVERBOSE)print*,'charge',charge |
363 |
c$$$c$$$ if(TRKVERBOSE)print*,'vect',vect |
364 |
c$$$c$$$ if(TRKVERBOSE)print*,'vout',vout |
365 |
c$$$c$$$ if(TRKVERBOSE)print*,'step',step |
366 |
c$$$c$$$ if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1 |
367 |
c$$$ RETURN |
368 |
c$$$ ENDIF |
369 |
c$$$ Z=VOUT(3) |
370 |
c$$$ IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200 |
371 |
c$$$ IF(Z.LT.ZIN(I)-TOLL) GOTO 20 |
372 |
c$$$ IF(Z.GE.ZIN(I)+TOLL) THEN |
373 |
c$$$ L = L - STEP |
374 |
c$$$ STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) |
375 |
c$$$ DO J=1,7 |
376 |
c$$$ VECT(J)=VECTINI(J) |
377 |
c$$$ ENDDO |
378 |
c$$$ GOTO 22 |
379 |
c$$$ ENDIF |
380 |
c$$$ 200 XOUT(I)=VOUT(1) |
381 |
c$$$ YOUT(I)=VOUT(2) |
382 |
c$$$ ZIN(I)=VOUT(3) |
383 |
c$$$ IF(VOUT(3).ne.0)THEN |
384 |
c$$$ THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.) |
385 |
c$$$ THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) |
386 |
c$$$ ENDIF |
387 |
c$$$ TLOUT(I) = L |
388 |
c$$$ ENDDO |
389 |
c$$$ |
390 |
RETURN |
391 |
END |
392 |
|
393 |
************************************************************************ |
394 |
* |
395 |
* |
396 |
* |
397 |
* |
398 |
************************************************************************ |
399 |
SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL) |
400 |
|
401 |
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
402 |
|
403 |
* ----------------------------------------------- |
404 |
* I/O parameters |
405 |
PARAMETER (NPOINT_MAX=500) ! EM it was =100 |
406 |
DIMENSION ZIN(NPOINT_MAX) |
407 |
DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX) |
408 |
DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX) |
409 |
DIMENSION TLOUT(NPOINT_MAX) |
410 |
DIMENSION AL_P(5) |
411 |
* ----------------------------------------------- |
412 |
|
413 |
CALL DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL) |
414 |
|
415 |
ccc |
416 |
ccc OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD |
417 |
ccc |
418 |
c$$$ DIMENSION VECT(7),VECTINI(7),VOUT(7) |
419 |
c$$$ DATA TOLL/1.d-8/ |
420 |
c$$$* tolerance in reaching the next plane during the tracking procedure |
421 |
c$$$* ----------------------------------------------- |
422 |
c$$$* I/O parameters |
423 |
c$$$ PARAMETER (NPOINT_MAX=100) |
424 |
c$$$ DIMENSION ZIN(NPOINT_MAX) |
425 |
c$$$ DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX) |
426 |
c$$$ DIMENSION AL_P(5) |
427 |
c$$$* ----------------------------------------------- |
428 |
c$$$ DATA ZINI/23.5/ !z coordinate of the reference plane |
429 |
c$$$ |
430 |
c$$$* ================================================================== |
431 |
c$$$* divide the track in two parts: below and above the reference plane |
432 |
c$$$* ================================================================== |
433 |
c$$$ IUPDOWN=0 |
434 |
c$$$ DO I=1,NPOINT |
435 |
c$$$ IF(ZIN(I).LT.ZINI)THEN |
436 |
c$$$ if(i.ne.1)IUPDOWN=I |
437 |
c$$$ GOTO 88 |
438 |
c$$$ ENDIF |
439 |
c$$$ IUPDOWN=NPOINT+1 |
440 |
c$$$ ENDDO |
441 |
c$$$ 88 CONTINUE |
442 |
c$$$ |
443 |
c$$$* ================================================================== |
444 |
c$$$* track from reference plane DOWN |
445 |
c$$$* ================================================================== |
446 |
c$$$* parameters for GRKUTA |
447 |
c$$$ IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5)) |
448 |
c$$$ IF(AL_P(5).EQ.0) CHARGE=1. |
449 |
c$$$ VOUT(1)=AL_P(1) |
450 |
c$$$ VOUT(2)=AL_P(2) |
451 |
c$$$ VOUT(3)=ZINI |
452 |
c$$$ VOUT(4)=AL_P(3)*DCOS(AL_P(4)) |
453 |
c$$$ VOUT(5)=AL_P(3)*DSIN(AL_P(4)) |
454 |
c$$$ VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2) |
455 |
c$$$ IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) |
456 |
c$$$ IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
457 |
c$$$ DO I=MAX(IUPDOWN,1),NPOINT |
458 |
c$$$ step=vout(3)-zin(i) |
459 |
c$$$c$$$ print*,'DOWN ',i,' - Track from ', |
460 |
c$$$c$$$ $ vout(3),' to ',zin(i),' step ',step |
461 |
c$$$ 10 DO J=1,7 |
462 |
c$$$ VECT(J)=VOUT(J) |
463 |
c$$$ VECTINI(J)=VOUT(J) |
464 |
c$$$ ENDDO |
465 |
c$$$ 11 continue |
466 |
c$$$ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
467 |
c$$$ IF(VOUT(3).GT.VECT(3)) THEN |
468 |
c$$$ IFAIL=1 |
469 |
c$$$c$$$ PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)' |
470 |
c$$$c$$$ print*,'charge',charge |
471 |
c$$$c$$$ print*,'vect',vect |
472 |
c$$$c$$$ print*,'vout',vout |
473 |
c$$$c$$$ print*,'step',step |
474 |
c$$$ RETURN |
475 |
c$$$ ENDIF |
476 |
c$$$ Z=VOUT(3) |
477 |
c$$$ IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100 |
478 |
c$$$ IF(Z.GT.ZIN(I)+TOLL) GOTO 10 |
479 |
c$$$ IF(Z.LE.ZIN(I)-TOLL) THEN |
480 |
c$$$ STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) |
481 |
c$$$ DO J=1,7 |
482 |
c$$$ VECT(J)=VECTINI(J) |
483 |
c$$$ ENDDO |
484 |
c$$$ GOTO 11 |
485 |
c$$$ ENDIF |
486 |
c$$$ 100 XOUT(I)=VOUT(1) |
487 |
c$$$ YOUT(I)=VOUT(2) |
488 |
c$$$ ZIN(I)=VOUT(3) |
489 |
c$$$* THXOUT(I)= |
490 |
c$$$* THYOUT(I)= |
491 |
c$$$ ENDDO |
492 |
c$$$ |
493 |
c$$$ |
494 |
c$$$ |
495 |
c$$$* ================================================================== |
496 |
c$$$* track from reference plane UP |
497 |
c$$$* ================================================================== |
498 |
c$$$* parameters for GRKUTA: |
499 |
c$$$* -opposite charge |
500 |
c$$$* -opposite momentum direction |
501 |
c$$$ IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5)) |
502 |
c$$$ IF(AL_P(5).EQ.0) CHARGE=-1. |
503 |
c$$$ VOUT(1)=AL_P(1) |
504 |
c$$$ VOUT(2)=AL_P(2) |
505 |
c$$$ VOUT(3)=ZINI |
506 |
c$$$ VOUT(4)=-AL_P(3)*DCOS(AL_P(4)) |
507 |
c$$$ VOUT(5)=-AL_P(3)*DSIN(AL_P(4)) |
508 |
c$$$ VOUT(6)=1.*DSQRT(1.-AL_P(3)**2) |
509 |
c$$$ IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) |
510 |
c$$$ IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 |
511 |
c$$$ DO I=MIN((IUPDOWN-1),NPOINT),1,-1 |
512 |
c$$$ step=vout(3)-zin(i) |
513 |
c$$$ step = -step |
514 |
c$$$c$$$ print*,'UP ',i,' - Track from ', |
515 |
c$$$c$$$ $ vout(3),' to ',zin(i),' step ',step |
516 |
c$$$ 20 DO J=1,7 |
517 |
c$$$ VECT(J)=VOUT(J) |
518 |
c$$$ VECTINI(J)=VOUT(J) |
519 |
c$$$ ENDDO |
520 |
c$$$ 22 continue |
521 |
c$$$ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) |
522 |
c$$$ IF(VOUT(3).LT.VECT(3)) THEN |
523 |
c$$$ IFAIL=1 |
524 |
c$$$c$$$ PRINT *,'=== WARNING ===> tracciamento invertito (UP)' |
525 |
c$$$c$$$ print*,'charge',charge |
526 |
c$$$c$$$ print*,'vect',vect |
527 |
c$$$c$$$ print*,'vout',vout |
528 |
c$$$c$$$ print*,'step',step |
529 |
c$$$ RETURN |
530 |
c$$$ ENDIF |
531 |
c$$$ Z=VOUT(3) |
532 |
c$$$ IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200 |
533 |
c$$$ IF(Z.LT.ZIN(I)-TOLL) GOTO 20 |
534 |
c$$$ IF(Z.GE.ZIN(I)+TOLL) THEN |
535 |
c$$$ STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) |
536 |
c$$$ DO J=1,7 |
537 |
c$$$ VECT(J)=VECTINI(J) |
538 |
c$$$ ENDDO |
539 |
c$$$ GOTO 22 |
540 |
c$$$ ENDIF |
541 |
c$$$ 200 XOUT(I)=VOUT(1) |
542 |
c$$$ YOUT(I)=VOUT(2) |
543 |
c$$$ ZIN(I)=VOUT(3) |
544 |
c$$$ |
545 |
c$$$ ENDDO |
546 |
|
547 |
RETURN |
548 |
END |