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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by mocchiut, Fri May 19 13:15:56 2006 UTC revision 1.7 by mocchiut, Fri Jun 6 11:24:03 2014 UTC
# Line 1  Line 1 
1    
2    
3  *************************************************************  *************************************************************
4  *  *
5  *     Routine to compute the NPOINT track intersection points  *     Routines to compute the NPOINT track intersection points
6  *     with planes of z-coordinate given by ZIN  *     with planes of z-coordinate given by ZIN
7  *     given the track parameters  *     given the track parameters
8  *  *
9  *     The routine is based on GRKUTA, which computes the  *     The routine is based on GRKUTA, which computes the
10  *     trajectory of a charged particle in a magnetic field  *     trajectory of a charged particle in a magnetic field
11  *     by solving the equatoins of motion with Runge-Kuta method  *     by solving the equations of motion with Runge-Kuta method
12  *  *
13  *     Variables that have to be assigned when the subroutine  *     Variables that have to be assigned when the subroutine
14  *     is called:  *     is called:
# Line 19  Line 20 
20  *     The routine works properly only if the  *     The routine works properly only if the
21  *     planes are numbered in descending order  *     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  *  *
 *     modified on 05/2006 to give as output also the angles and  
 *     theh track-lengths  
43  *  *
44  **************************************************************  **************************************************************
45    
46        SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)  
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)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
52    
53    
54        DIMENSION VECT(7),VECTINI(7),VOUT(7)        DIMENSION VECT(7),VECTINI(7),VOUT(7)
55        DATA TOLL/1.d-8/                  DATA TOLL/1.d-8/          
56    
57  *     tolerance in reaching the next plane during the tracking procedure  *     tolerance in reaching the next plane during the tracking procedure
58  *     -----------------------------------------------  *     -----------------------------------------------
59  *     I/O parameters  *     I/O parameters
60        PARAMETER (NPOINT_MAX=100)        PARAMETER (NPOINT_MAX=500)
61        DIMENSION ZIN(NPOINT_MAX)        DIMENSION ZIN(NPOINT_MAX)
62        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
63        DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)        DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
64          DIMENSION TLOUT(NPOINT_MAX)
65        DIMENSION AL_P(5)        DIMENSION AL_P(5)
66  *     -----------------------------------------------        *     -----------------------------------------------      
67        DATA ZINI/23.5/           !z coordinate of the reference plane  *      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  *     divide the track in two parts: below and above the reference plane
72  *     ==================================================================        *     ==================================================================      
# Line 69  Line 94 
94        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
95        IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))        IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
96        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
97        DO I=MAX(IUPDOWN,1),NPOINT        DO I=MAX(IUPDOWN,1),NPOINT
98           step=vout(3)-zin(i)           L = 0.0
 c$$$         print*,'DOWN ',i,' - Track from ',  
 c$$$     $        vout(3),' to ',zin(i),' step ',step  
99   10      DO J=1,7   10      DO J=1,7
100              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
101              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
102           ENDDO           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   11      continue
111           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
112             L = L + STEP
113           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
114              IFAIL=1              IFAIL=1
115              PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'  c$$$            if(TRKVERBOSE)
116              print*,'charge',charge  c$$$     $      PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
117              print*,'vect',vect  c$$$            if(TRKVERBOSE)print*,'charge',charge
118              print*,'vout',vout  c$$$            if(TRKVERBOSE)print*,'vect',vect
119              print*,'step',step  c$$$            if(TRKVERBOSE)print*,'vout',vout
120    c$$$            if(TRKVERBOSE)print*,'step',step
121    c$$$            if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
122              RETURN              RETURN
123           ENDIF           ENDIF
124           Z=VOUT(3)           Z=VOUT(3)
125           IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100           IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
126           IF(Z.GT.ZIN(I)+TOLL) GOTO 10                 IF(Z.GT.ZIN(I)+TOLL) GOTO 10      
127           IF(Z.LE.ZIN(I)-TOLL) THEN           IF(Z.LE.ZIN(I)-TOLL) THEN
128                L = L - STEP
129              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
130              DO J=1,7              DO J=1,7
131                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
# Line 101  c$$$     $        vout(3),' to ',zin(i), Line 135  c$$$     $        vout(3),' to ',zin(i),
135   100     XOUT(I)=VOUT(1)   100     XOUT(I)=VOUT(1)
136           YOUT(I)=VOUT(2)           YOUT(I)=VOUT(2)
137           ZIN(I)=VOUT(3)           ZIN(I)=VOUT(3)
138  *         THXOUT(I)=           IF(VOUT(3).ne.0)THEN
139  *         THYOUT(I)=              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        ENDDO
145    
   
   
146  *     ==================================================================  *     ==================================================================
147  *     track from reference plane UP  *     track from reference plane UP
148  *     ==================================================================  *     ==================================================================
# Line 124  c$$$     $        vout(3),' to ',zin(i), Line 160  c$$$     $        vout(3),' to ',zin(i),
160        IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))        IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
161        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
162        DO I=MIN((IUPDOWN-1),NPOINT),1,-1        DO I=MIN((IUPDOWN-1),NPOINT),1,-1
163           step=vout(3)-zin(i)           L = 0.0
164           step = -step  
 c$$$         print*,'UP ',i,' - Track from ',  
 c$$$     $        vout(3),' to ',zin(i),' step ',step  
165   20      DO J=1,7   20      DO J=1,7
166              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
167              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
168           ENDDO           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   22      continue
177           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
178             L = L + STEP
179           IF(VOUT(3).LT.VECT(3)) THEN           IF(VOUT(3).LT.VECT(3)) THEN
180              IFAIL=1              IFAIL=1
181              PRINT *,'=== WARNING ===> tracciamento invertito (UP)'  c$$$            if(TRKVERBOSE)
182              print*,'charge',charge  c$$$     $      PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
183              print*,'vect',vect  c$$$            if(TRKVERBOSE)print*,'charge',charge
184              print*,'vout',vout  c$$$            if(TRKVERBOSE)print*,'vect',vect
185              print*,'step',step  c$$$            if(TRKVERBOSE)print*,'vout',vout
186    c$$$            if(TRKVERBOSE)print*,'step',step
187    c$$$            if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
188              RETURN              RETURN
189           ENDIF           ENDIF
190           Z=VOUT(3)           Z=VOUT(3)
191           IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200           IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
192           IF(Z.LT.ZIN(I)-TOLL) GOTO 20                 IF(Z.LT.ZIN(I)-TOLL) GOTO 20      
193           IF(Z.GE.ZIN(I)+TOLL) THEN           IF(Z.GE.ZIN(I)+TOLL) THEN
194                L = L - STEP
195              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
196              DO J=1,7              DO J=1,7
197                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
# Line 156  c$$$     $        vout(3),' to ',zin(i), Line 201  c$$$     $        vout(3),' to ',zin(i),
201   200     XOUT(I)=VOUT(1)   200     XOUT(I)=VOUT(1)
202           YOUT(I)=VOUT(2)           YOUT(I)=VOUT(2)
203           ZIN(I)=VOUT(3)           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        ENDDO
210    
211        RETURN        RETURN
212        END        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
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23