src/F77/track.f


*************************************************************
*
*     Routine to compute the NPOINT track intersection points
*     with planes of z-coordinate given by ZIN
*     given the track parameters
*
*     The routine is based on GRKUTA, which computes the 
*     trajectory of a charged particle in a magnetic field
*     by solving the equatoins of motion with Runge-Kuta method
*
*     Variables that have to be assigned when the subroutine 
*     is called:
*
*     ZIN(1-NPOINT) ----> z coordinates of the planes
*     AL_P(1-5)     ----> track-parameter vector
*
*     NB !!! 
*     The routine works properly only if the 
*     planes are numbered in descending order 
*     
*
*     modified on 05/2006 to give as output also the angles and
*     theh track-lengths 
*
**************************************************************

      SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      DIMENSION VECT(7),VECTINI(7),VOUT(7)
      DATA TOLL/1.d-8/          
*     tolerance in reaching the next plane during the tracking procedure
*     -----------------------------------------------
*     I/O parameters
      PARAMETER (NPOINT_MAX=100)
      DIMENSION ZIN(NPOINT_MAX)
      DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
      DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
      DIMENSION AL_P(5)
*     -----------------------------------------------      
      DATA ZINI/23.5/           !z coordinate of the reference plane
      
*     ==================================================================
*     divide the track in two parts: below and above the reference plane
*     ==================================================================      
      IUPDOWN=0
      DO I=1,NPOINT   
         IF(ZIN(I).LT.ZINI)THEN
            if(i.ne.1)IUPDOWN=I
            GOTO 88
         ENDIF
         IUPDOWN=NPOINT+1
      ENDDO
 88   CONTINUE

*     ==================================================================
*     track from reference plane DOWN
*     ==================================================================
*     parameters for GRKUTA
      IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
      IF(AL_P(5).EQ.0) CHARGE=1.
      VOUT(1)=AL_P(1)
      VOUT(2)=AL_P(2)
      VOUT(3)=ZINI             
      VOUT(4)=AL_P(3)*DCOS(AL_P(4))
      VOUT(5)=AL_P(3)*DSIN(AL_P(4))
      VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
      IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
      IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
      DO I=MAX(IUPDOWN,1),NPOINT 
         step=vout(3)-zin(i)
c$$$         print*,'DOWN ',i,' - Track from ',
c$$$     $        vout(3),' to ',zin(i),' step ',step
 10      DO J=1,7
            VECT(J)=VOUT(J)
            VECTINI(J)=VOUT(J)
         ENDDO
 11      continue
         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
         IF(VOUT(3).GT.VECT(3)) THEN
            IFAIL=1
            PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'
            print*,'charge',charge
            print*,'vect',vect
            print*,'vout',vout
            print*,'step',step
            RETURN
         ENDIF
         Z=VOUT(3)
         IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
         IF(Z.GT.ZIN(I)+TOLL) GOTO 10      
         IF(Z.LE.ZIN(I)-TOLL) THEN
            STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
            DO J=1,7
               VECT(J)=VECTINI(J)
            ENDDO
            GOTO 11
         ENDIF
 100     XOUT(I)=VOUT(1)
         YOUT(I)=VOUT(2)
         ZIN(I)=VOUT(3)
*         THXOUT(I)=
*         THYOUT(I)=
      ENDDO


*     ==================================================================
*     track from reference plane UP
*     ==================================================================
*     parameters for GRKUTA:
*     -opposite charge
*     -opposite momentum direction
      IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5))
      IF(AL_P(5).EQ.0) CHARGE=-1.
      VOUT(1)=AL_P(1)
      VOUT(2)=AL_P(2)
      VOUT(3)=ZINI             
      VOUT(4)=-AL_P(3)*DCOS(AL_P(4))
      VOUT(5)=-AL_P(3)*DSIN(AL_P(4))
      VOUT(6)=1.*DSQRT(1.-AL_P(3)**2)
      IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
      IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
      DO I=MIN((IUPDOWN-1),NPOINT),1,-1
         step=vout(3)-zin(i)
         step = -step
c$$$         print*,'UP ',i,' - Track from ',
c$$$     $        vout(3),' to ',zin(i),' step ',step
 20      DO J=1,7
            VECT(J)=VOUT(J)
            VECTINI(J)=VOUT(J)
         ENDDO
 22      continue
         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
         IF(VOUT(3).LT.VECT(3)) THEN
            IFAIL=1
            PRINT *,'=== WARNING ===> tracciamento invertito (UP)'
            print*,'charge',charge
            print*,'vect',vect
            print*,'vout',vout
            print*,'step',step
            RETURN
         ENDIF
         Z=VOUT(3)
         IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
         IF(Z.LT.ZIN(I)-TOLL) GOTO 20      
         IF(Z.GE.ZIN(I)+TOLL) THEN
            STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
            DO J=1,7
               VECT(J)=VECTINI(J)
            ENDDO
            GOTO 22
         ENDIF
 200     XOUT(I)=VOUT(1)
         YOUT(I)=VOUT(2)
         ZIN(I)=VOUT(3)

      ENDDO

      RETURN
      END

1	mocchiut	1.1
2			*************************************************************
3			*
4			* Routine to compute the NPOINT track intersection points
5			* with planes of z-coordinate given by ZIN
6			* given the track parameters
7			*
8			* The routine is based on GRKUTA, which computes the
9			* trajectory of a charged particle in a magnetic field
10			* by solving the equatoins of motion with Runge-Kuta method
11			*
12			* Variables that have to be assigned when the subroutine
13			* is called:
14			*
15			* ZIN(1-NPOINT) ----> z coordinates of the planes
16			* AL_P(1-5) ----> track-parameter vector
17			*
18			* NB !!!
19			* The routine works properly only if the
20			* planes are numbered in descending order
21			*
22			*
23			* modified on 05/2006 to give as output also the angles and
24			* theh track-lengths
25			*
26			**************************************************************
27
28			SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)
29
30			IMPLICIT DOUBLE PRECISION (A-H,O-Z)
31
32			DIMENSION VECT(7),VECTINI(7),VOUT(7)
33			DATA TOLL/1.d-8/
34			* tolerance in reaching the next plane during the tracking procedure
35			* -----------------------------------------------
36			* I/O parameters
37			PARAMETER (NPOINT_MAX=100)
38			DIMENSION ZIN(NPOINT_MAX)
39			DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
40			DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
41			DIMENSION AL_P(5)
42			* -----------------------------------------------
43			DATA ZINI/23.5/ !z coordinate of the reference plane
44
45			* ==================================================================
46			* divide the track in two parts: below and above the reference plane
47			* ==================================================================
48			IUPDOWN=0
49			DO I=1,NPOINT
50			IF(ZIN(I).LT.ZINI)THEN
51			if(i.ne.1)IUPDOWN=I
52			GOTO 88
53			ENDIF
54			IUPDOWN=NPOINT+1
55			ENDDO
56			88 CONTINUE
57
58			* ==================================================================
59			* track from reference plane DOWN
60			* ==================================================================
61			* parameters for GRKUTA
62			IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
63			IF(AL_P(5).EQ.0) CHARGE=1.
64			VOUT(1)=AL_P(1)
65			VOUT(2)=AL_P(2)
66			VOUT(3)=ZINI
67			VOUT(4)=AL_P(3)*DCOS(AL_P(4))
68			VOUT(5)=AL_P(3)*DSIN(AL_P(4))
69			VOUT(6)=-1.DSQRT(1.-AL_P(3)*2)
70			IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
71			IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
72			DO I=MAX(IUPDOWN,1),NPOINT
73			step=vout(3)-zin(i)
74			c$$$ print*,'DOWN ',i,' - Track from ',
75			c$$$ $ vout(3),' to ',zin(i),' step ',step
76			10 DO J=1,7
77			VECT(J)=VOUT(J)
78			VECTINI(J)=VOUT(J)
79			ENDDO
80			11 continue
81			CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
82			IF(VOUT(3).GT.VECT(3)) THEN
83			IFAIL=1
84			PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'
85			print*,'charge',charge
86			print*,'vect',vect
87			print*,'vout',vout
88			print*,'step',step
89			RETURN
90			ENDIF
91			Z=VOUT(3)
92			IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
93			IF(Z.GT.ZIN(I)+TOLL) GOTO 10
94			IF(Z.LE.ZIN(I)-TOLL) THEN
95			STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
96			DO J=1,7
97			VECT(J)=VECTINI(J)
98			ENDDO
99			GOTO 11
100			ENDIF
101			100 XOUT(I)=VOUT(1)
102			YOUT(I)=VOUT(2)
103			ZIN(I)=VOUT(3)
104			* THXOUT(I)=
105			* THYOUT(I)=
106			ENDDO
107
108
109
110			* ==================================================================
111			* track from reference plane UP
112			* ==================================================================
113			* parameters for GRKUTA:
114			* -opposite charge
115			* -opposite momentum direction
116			IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5))
117			IF(AL_P(5).EQ.0) CHARGE=-1.
118			VOUT(1)=AL_P(1)
119			VOUT(2)=AL_P(2)
120			VOUT(3)=ZINI
121			VOUT(4)=-AL_P(3)*DCOS(AL_P(4))
122			VOUT(5)=-AL_P(3)*DSIN(AL_P(4))
123			VOUT(6)=1.DSQRT(1.-AL_P(3)*2)
124			IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
125			IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
126			DO I=MIN((IUPDOWN-1),NPOINT),1,-1
127			step=vout(3)-zin(i)
128			step = -step
129			c$$$ print*,'UP ',i,' - Track from ',
130			c$$$ $ vout(3),' to ',zin(i),' step ',step
131			20 DO J=1,7
132			VECT(J)=VOUT(J)
133			VECTINI(J)=VOUT(J)
134			ENDDO
135			22 continue
136			CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
137			IF(VOUT(3).LT.VECT(3)) THEN
138			IFAIL=1
139			PRINT *,'=== WARNING ===> tracciamento invertito (UP)'
140			print*,'charge',charge
141			print*,'vect',vect
142			print*,'vout',vout
143			print*,'step',step
144			RETURN
145			ENDIF
146			Z=VOUT(3)
147			IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
148			IF(Z.LT.ZIN(I)-TOLL) GOTO 20
149			IF(Z.GE.ZIN(I)+TOLL) THEN
150			STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
151			DO J=1,7
152			VECT(J)=VECTINI(J)
153			ENDDO
154			GOTO 22
155			ENDIF
156			200 XOUT(I)=VOUT(1)
157			YOUT(I)=VOUT(2)
158			ZIN(I)=VOUT(3)
159
160			ENDDO
161
162			RETURN
163			END
164