/[PAMELA software]/eventviewer/ground/src/track.for
ViewVC logotype

Annotation of /eventviewer/ground/src/track.for

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Tue Dec 6 10:12:57 2005 UTC (18 years, 11 months ago) by mocchiut
Branch point for: MAIN, EventViewer
Initial revision

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 starting from the
21     * reference plane (ZINI)
22     *
23     **************************************************************
24    
25     SUBROUTINE TRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)
26    
27     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
28    
29     DIMENSION VECT(7),VECTINI(7),VOUT(7)
30     DATA TOLL/1.d-8/
31     * tolerance in reaching the next plane during the tracking procedure
32     * -----------------------------------------------
33     * I/O parameters
34     PARAMETER (NPOINT_MAX=100)
35     DIMENSION ZIN(NPOINT_MAX)
36     DIMENSION XOUT(NPOINT_MAX)
37     DIMENSION YOUT(NPOINT_MAX)
38     DIMENSION AL_P(5)
39     DOUBLE PRECISION Z
40     DOUBLE PRECISION STEP
41     DOUBLE PRECISION CHARGE
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     Z = 0.
49     IUPDOWN=0
50     DO I=1,NPOINT
51     IF(ZIN(I).LT.ZINI)THEN
52     if(i.ne.1)IUPDOWN=I
53     GOTO 88
54     ENDIF
55     IUPDOWN=NPOINT+1
56     ENDDO
57     88 CONTINUE
58     c DO I=1,NPOINT
59     c print *,' TRACK >>>>>>>>> ',I,' ZIN ',ZIN(i)
60     c enddo
61     c print *,' iupdown = ',iupdown
62    
63    
64     * ==================================================================
65     * track from reference plane DOWN
66     * ==================================================================
67     * parameters for GRKUTA
68     IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
69     IF(AL_P(5).EQ.0) CHARGE=1.
70     VOUT(1)=AL_P(1)
71     VOUT(2)=AL_P(2)
72     VOUT(3)=ZINI
73     VOUT(4)=AL_P(3)*DCOS(AL_P(4))
74     VOUT(5)=AL_P(3)*DSIN(AL_P(4))
75     VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
76     IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
77     IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
78     c
79     DO I=MAX(IUPDOWN,1),NPOINT
80     step=vout(3)-zin(i)
81     c print*,'DOWN ',i,' - Track from ',
82     c $ vout(3),' to ',zin(i),' step ',step
83     10 DO J=1,7
84     VECT(J)=VOUT(J)
85     VECTINI(J)=VOUT(J)
86     c print *,'vect ',j,' ',vect(j)
87     ENDDO
88     11 continue
89     CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
90     IF(VOUT(3).GT.VECT(3)) THEN
91     IFAIL=1
92     PRINT *,'=== WARNING ===> inverted tracking (DOWN)'
93     print*,'charge',charge
94     print*,'vect',vect
95     print*,'vout',vout
96     print*,'step',step
97     RETURN
98     ENDIF
99     Z=VOUT(3)
100     c if(Z.ne.Z) goto 11
101     IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
102     IF(Z.GT.ZIN(I)+TOLL) GOTO 10
103     IF(Z.LE.ZIN(I)-TOLL) THEN
104     STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
105     DO J=1,7
106     VECT(J)=VECTINI(J)
107     ENDDO
108     GOTO 11
109     ENDIF
110     100 XOUT(I)=VOUT(1)
111     YOUT(I)=VOUT(2)
112     ZIN(I)=VOUT(3)
113     c print *,'vout3 ',VOUT(3)
114     ENDDO
115    
116     * ==================================================================
117     * track from refernce plane UP
118     * ==================================================================
119     * parameters for GRKUTA:
120     * -opposite charge
121     * -opposite momentum direction
122     IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5))
123     IF(AL_P(5).EQ.0) CHARGE=-1.
124     VOUT(1)=AL_P(1)
125     VOUT(2)=AL_P(2)
126     VOUT(3)=ZINI
127     VOUT(4)=-AL_P(3)*DCOS(AL_P(4))
128     VOUT(5)=-AL_P(3)*DSIN(AL_P(4))
129     VOUT(6)=1.*DSQRT(1.-AL_P(3)**2)
130     IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
131     IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
132     DO I=MIN((IUPDOWN-1),NPOINT),1,-1
133     step=vout(3)-zin(i)
134     step = -step
135     c print*,'UP ',i,' - Track from ',
136     c $ vout(3),' to ',zin(i),' step ',step
137     20 DO J=1,7
138     VECT(J)=VOUT(J)
139     VECTINI(J)=VOUT(J)
140     ENDDO
141     22 continue
142     CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
143     IF(VOUT(3).LT.VECT(3)) THEN
144     IFAIL=1
145     PRINT *,'=== WARNING ===> inverted tracking (UP)'
146     print*,'charge',charge
147     print*,'vect',vect
148     print*,'vout',vout
149     print*,'step',step
150     RETURN
151     ENDIF
152     Z=VOUT(3)
153     IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
154     IF(Z.LT.ZIN(I)-TOLL) GOTO 20
155     IF(Z.GE.ZIN(I)+TOLL) THEN
156     STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
157     DO J=1,7
158     VECT(J)=VECTINI(J)
159     ENDDO
160     GOTO 22
161     ENDIF
162     200 XOUT(I)=VOUT(1)
163     YOUT(I)=VOUT(2)
164     ZIN(I)=VOUT(3)
165     ENDDO
166    
167     RETURN
168     END
169    
170     include './grkuta.for'
171     include './inter_B.for'

  ViewVC Help
Powered by ViewVC 1.1.23