/[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.2 by pam-fi, Thu Jun 1 09:03:09 2006 UTC
# Line 1  Line 1 
1    
2  *************************************************************  *************************************************************
3  *  *
4  *     Routine to compute the NPOINT track intersection points  *     Routines to compute the NPOINT track intersection points
5  *     with planes of z-coordinate given by ZIN  *     with planes of z-coordinate given by ZIN
6  *     given the track parameters  *     given the track parameters
7  *  *
8  *     The routine is based on GRKUTA, which computes the  *     The routine is based on GRKUTA, which computes the
9  *     trajectory of a charged particle in a magnetic field  *     trajectory of a charged particle in a magnetic field
10  *     by solving the equatoins of motion with Runge-Kuta method  *     by solving the equations of motion with Runge-Kuta method
11  *  *
12  *     Variables that have to be assigned when the subroutine  *     Variables that have to be assigned when the subroutine
13  *     is called:  *     is called:
# Line 19  Line 19 
19  *     The routine works properly only if the  *     The routine works properly only if the
20  *     planes are numbered in descending order  *     planes are numbered in descending order
21  *      *    
22    *     Routines:
23  *  *
24  *     modified on 05/2006 to give as output also the angles and  *     DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)  
25  *     theh track-lengths  *     - old routine
26    *    
27    *     DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,AL_P,IFAIL)
28    *     - as the old one, but:
29    *     -- the projected angles are given as output
30    *     -- the track lengths are given as output:
31    *     ---- for planes above the reference plane, the length until
32    *          the lower closest plane (reference plane included)
33    *     ---- for planes below the reference plane, the length until
34    *          the higher closest plane (reference plane included)
35  *  *
36  **************************************************************  **************************************************************
37    
38          SUBROUTINE
39         $     DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
40    
41          IMPLICIT DOUBLE PRECISION (A-H,O-Z)
42    
43          DIMENSION VECT(8),VECTINI(8),VOUT(8)
44          DATA TOLL/1.d-8/          
45    *     tolerance in reaching the next plane during the tracking procedure
46    *     -----------------------------------------------
47    *     I/O parameters
48          PARAMETER (NPOINT_MAX=100)
49          DIMENSION ZIN(NPOINT_MAX)
50          DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
51          DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
52          DIMENSION TLOUT(NPOINT_MAX)
53          DIMENSION AL_P(5)
54    *     -----------------------------------------------      
55          DATA ZINI/23.5/           !z coordinate of the reference plane
56          
57    *     ==================================================================
58    *     divide the track in two parts: below and above the reference plane
59    *     ==================================================================      
60          IUPDOWN=0
61          DO I=1,NPOINT  
62             IF(ZIN(I).LT.ZINI)THEN
63                if(i.ne.1)IUPDOWN=I
64                GOTO 88
65             ENDIF
66             IUPDOWN=NPOINT+1
67          ENDDO
68     88   CONTINUE
69    
70    *     ==================================================================
71    *     track from reference plane DOWN
72    *     ==================================================================
73    *     parameters for GRKUTA
74          IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
75          IF(AL_P(5).EQ.0) CHARGE=1.
76          VOUT(1)=AL_P(1)
77          VOUT(2)=AL_P(2)
78          VOUT(3)=ZINI            
79          VOUT(4)=AL_P(3)*DCOS(AL_P(4))
80          VOUT(5)=AL_P(3)*DSIN(AL_P(4))
81          VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
82          IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
83          IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
84          DO I=MAX(IUPDOWN,1),NPOINT
85             step=vout(3)-zin(i)
86             VOUT(8)= 0
87    c$$$         print*,'DOWN ',i,' - Track from ',
88    c$$$     $        vout(3),' to ',zin(i),' step ',step
89     10      DO J=1,8
90                VECT(J)=VOUT(J)
91                VECTINI(J)=VOUT(J)
92             ENDDO
93     11      continue
94             CALL GRKUTA2(CHARGE,STEP,VECT,VOUT)
95             IF(VOUT(3).GT.VECT(3)) THEN
96                IFAIL=1
97    c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'
98    c$$$            print*,'charge',charge
99    c$$$            print*,'vect',vect
100    c$$$            print*,'vout',vout
101    c$$$            print*,'step',step
102                RETURN
103             ENDIF
104             Z=VOUT(3)
105             IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
106             IF(Z.GT.ZIN(I)+TOLL) GOTO 10      
107             IF(Z.LE.ZIN(I)-TOLL) THEN
108                STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
109                DO J=1,8
110                   VECT(J)=VECTINI(J)
111                ENDDO
112                GOTO 11
113             ENDIF
114     100     XOUT(I)=VOUT(1)
115             YOUT(I)=VOUT(2)
116             ZIN(I)=VOUT(3)
117             IF(VOUT(3).ne.0)THEN
118                THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
119                THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
120             ENDIF
121             TLOUT(I) = VOUT(8)
122    c         print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP
123          ENDDO
124    
125    
126    
127    *     ==================================================================
128    *     track from reference plane UP
129    *     ==================================================================
130    *     parameters for GRKUTA:
131    *     -opposite charge
132    *     -opposite momentum direction
133          IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5))
134          IF(AL_P(5).EQ.0) CHARGE=-1.
135          VOUT(1)=AL_P(1)
136          VOUT(2)=AL_P(2)
137          VOUT(3)=ZINI            
138          VOUT(4)=-AL_P(3)*DCOS(AL_P(4))
139          VOUT(5)=-AL_P(3)*DSIN(AL_P(4))
140          VOUT(6)=1.*DSQRT(1.-AL_P(3)**2)
141          IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
142          IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
143          DO I=MIN((IUPDOWN-1),NPOINT),1,-1
144             step=vout(3)-zin(i)
145             step = -step
146             VOUT(8)=0
147    c$$$         print*,'UP ',i,' - Track from ',
148    c$$$     $        vout(3),' to ',zin(i),' step ',step
149     20      DO J=1,8
150                VECT(J)=VOUT(J)
151                VECTINI(J)=VOUT(J)
152             ENDDO
153     22      continue
154             CALL GRKUTA2(CHARGE,STEP,VECT,VOUT)
155             IF(VOUT(3).LT.VECT(3)) THEN
156                IFAIL=1
157    c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (UP)'
158    c$$$            print*,'charge',charge
159    c$$$            print*,'vect',vect
160    c$$$            print*,'vout',vout
161    c$$$            print*,'step',step
162                RETURN
163             ENDIF
164             Z=VOUT(3)
165             IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
166             IF(Z.LT.ZIN(I)-TOLL) GOTO 20      
167             IF(Z.GE.ZIN(I)+TOLL) THEN
168                STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
169                DO J=1,8
170                   VECT(J)=VECTINI(J)
171                ENDDO
172                GOTO 22
173             ENDIF
174     200     XOUT(I)=VOUT(1)
175             YOUT(I)=VOUT(2)
176             ZIN(I)=VOUT(3)
177             IF(VOUT(3).ne.0)THEN
178                THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
179                THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
180             ENDIF
181             TLOUT(I) = VOUT(8)
182    
183          ENDDO
184    
185          RETURN
186          END
187    
188    ************************************************************************
189    *
190    *
191    *
192    *
193    ************************************************************************
194        SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)        SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)
195    
196        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
# Line 37  Line 203 
203        PARAMETER (NPOINT_MAX=100)        PARAMETER (NPOINT_MAX=100)
204        DIMENSION ZIN(NPOINT_MAX)        DIMENSION ZIN(NPOINT_MAX)
205        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
       DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)  
206        DIMENSION AL_P(5)        DIMENSION AL_P(5)
207  *     -----------------------------------------------        *     -----------------------------------------------      
208        DATA ZINI/23.5/           !z coordinate of the reference plane        DATA ZINI/23.5/           !z coordinate of the reference plane
# Line 81  c$$$     $        vout(3),' to ',zin(i), Line 246  c$$$     $        vout(3),' to ',zin(i),
246           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
247           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
248              IFAIL=1              IFAIL=1
249              PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'  c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'
250              print*,'charge',charge  c$$$            print*,'charge',charge
251              print*,'vect',vect  c$$$            print*,'vect',vect
252              print*,'vout',vout  c$$$            print*,'vout',vout
253              print*,'step',step  c$$$            print*,'step',step
254              RETURN             RETURN
255           ENDIF           ENDIF
256           Z=VOUT(3)           Z=VOUT(3)
257           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
# Line 136  c$$$     $        vout(3),' to ',zin(i), Line 301  c$$$     $        vout(3),' to ',zin(i),
301           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
302           IF(VOUT(3).LT.VECT(3)) THEN           IF(VOUT(3).LT.VECT(3)) THEN
303              IFAIL=1              IFAIL=1
304              PRINT *,'=== WARNING ===> tracciamento invertito (UP)'  c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (UP)'
305              print*,'charge',charge  c$$$            print*,'charge',charge
306              print*,'vect',vect  c$$$            print*,'vect',vect
307              print*,'vout',vout  c$$$            print*,'vout',vout
308              print*,'step',step  c$$$            print*,'step',step
309              RETURN              RETURN
310           ENDIF           ENDIF
311           Z=VOUT(3)           Z=VOUT(3)
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.23