/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/track.f
ViewVC logotype

Contents of /DarthVader/TrackerLevel2/src/F77/track.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (show annotations) (download)
Fri Jun 6 11:24:03 2014 UTC (10 years, 6 months ago) by mocchiut
Branch: MAIN
CVS Tags: v10RED, v10REDr01, HEAD
Changes since 1.6: +1 -1 lines
New version, STDOUT under _debug flag, GetDeflection and GetdEdx added to ExtTrack, compilation warnings fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23