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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (hide annotations) (download)
Fri Jun 6 11:24:03 2014 UTC (10 years, 5 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 mocchiut 1.1
2 pam-fi 1.3
3 mocchiut 1.1 *************************************************************
4     *
5 pam-fi 1.2 * Routines to compute the NPOINT track intersection points
6 mocchiut 1.1 * 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 pam-fi 1.2 * by solving the equations of motion with Runge-Kuta method
12 mocchiut 1.1 *
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 pam-fi 1.2 * Routines:
24 mocchiut 1.1 *
25 pam-fi 1.2 * DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)
26     * - old routine
27     *
28 pam-fi 1.3 * DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
29 pam-fi 1.4 * - as the previous one, but:
30 pam-fi 1.2 * -- 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 mocchiut 1.1 *
37 pam-fi 1.3 * March 2008 --> optimized by Paolo
38 pam-fi 1.4 *
39 pam-ts 1.5 * DOTRACK3(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,ZINI,IFAIL)
40 pam-fi 1.4 * - as the previous one, but with the coordinate ZINI of the reference
41     * plane given as input parameter
42     *
43     *
44 mocchiut 1.1 **************************************************************
45    
46 pam-fi 1.4
47 pam-fi 1.2 SUBROUTINE
48 pam-fi 1.4 $ DOTRACK3(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P
49     $ ,ZINI,IFAIL)
50 pam-fi 1.2
51     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
52    
53 pam-fi 1.3
54     DIMENSION VECT(7),VECTINI(7),VOUT(7)
55 pam-fi 1.2 DATA TOLL/1.d-8/
56 pam-fi 1.3
57 pam-fi 1.2 * tolerance in reaching the next plane during the tracking procedure
58     * -----------------------------------------------
59     * I/O parameters
60 pam-ts 1.5 PARAMETER (NPOINT_MAX=500)
61 pam-fi 1.2 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 pam-fi 1.4 * DATA ZINI/23.5/ !z coordinate of the reference plane
68 pam-fi 1.3 REAL*8 L
69    
70 pam-fi 1.2 * ==================================================================
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 pam-fi 1.3 DO I=MAX(IUPDOWN,1),NPOINT
98     L = 0.0
99     10 DO J=1,7
100 pam-fi 1.2 VECT(J)=VOUT(J)
101     VECTINI(J)=VOUT(J)
102     ENDDO
103 pam-fi 1.3 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 pam-fi 1.2 11 continue
111 pam-fi 1.3 CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
112     L = L + STEP
113 pam-fi 1.2 IF(VOUT(3).GT.VECT(3)) THEN
114     IFAIL=1
115 pam-fi 1.3 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 pam-fi 1.2 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 pam-fi 1.3 L = L - STEP
129 pam-fi 1.2 STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
130 pam-fi 1.3 DO J=1,7
131 pam-fi 1.2 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 pam-fi 1.3 TLOUT(I) = L
143 pam-fi 1.2 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 pam-fi 1.3 L = 0.0
164    
165     20 DO J=1,7
166 pam-fi 1.2 VECT(J)=VOUT(J)
167     VECTINI(J)=VOUT(J)
168     ENDDO
169 pam-fi 1.3 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 pam-fi 1.2 22 continue
177 pam-fi 1.3 CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
178     L = L + STEP
179 pam-fi 1.2 IF(VOUT(3).LT.VECT(3)) THEN
180     IFAIL=1
181 pam-fi 1.3 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 pam-fi 1.2 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 pam-fi 1.3 L = L - STEP
195 pam-fi 1.2 STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
196 pam-fi 1.3 DO J=1,7
197 pam-fi 1.2 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 pam-fi 1.3 TLOUT(I) = L
209 pam-fi 1.2 ENDDO
210    
211     RETURN
212     END
213 pam-fi 1.4 ************************************************************************
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 pam-ts 1.6 PARAMETER (NPOINT_MAX=500)
233 pam-fi 1.4 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 pam-fi 1.2
393     ************************************************************************
394     *
395     *
396     *
397     *
398     ************************************************************************
399 mocchiut 1.1 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 mocchiut 1.7 PARAMETER (NPOINT_MAX=500) ! EM it was =100
406 mocchiut 1.1 DIMENSION ZIN(NPOINT_MAX)
407     DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
408 pam-fi 1.3 DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
409     DIMENSION TLOUT(NPOINT_MAX)
410 mocchiut 1.1 DIMENSION AL_P(5)
411     * -----------------------------------------------
412    
413 pam-fi 1.3 CALL DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
414 mocchiut 1.1
415 pam-fi 1.3 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 mocchiut 1.1
547     RETURN
548     END

  ViewVC Help
Powered by ViewVC 1.1.23