/[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.3 - (hide annotations) (download)
Wed Mar 5 17:00:20 2008 UTC (16 years, 9 months ago) by pam-fi
Branch: MAIN
CVS Tags: v5r00, v6r01, v6r00
Changes since 1.2: +186 -154 lines
modified TrkSinglet, optimized DoTrack2, fixed bug in evaluation of effective angle

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.2 * - as the old 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 mocchiut 1.1 *
37 pam-fi 1.3 * March 2008 --> optimized by Paolo
38 mocchiut 1.1 **************************************************************
39    
40 pam-fi 1.2 SUBROUTINE
41     $ DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
42    
43     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
44    
45 pam-fi 1.3
46     DIMENSION VECT(7),VECTINI(7),VOUT(7)
47 pam-fi 1.2 DATA TOLL/1.d-8/
48 pam-fi 1.3
49 pam-fi 1.2 * tolerance in reaching the next plane during the tracking procedure
50     * -----------------------------------------------
51     * I/O parameters
52     PARAMETER (NPOINT_MAX=100)
53     DIMENSION ZIN(NPOINT_MAX)
54     DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
55     DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
56     DIMENSION TLOUT(NPOINT_MAX)
57     DIMENSION AL_P(5)
58     * -----------------------------------------------
59     DATA ZINI/23.5/ !z coordinate of the reference plane
60 pam-fi 1.3 REAL*8 L
61    
62 pam-fi 1.2 * ==================================================================
63     * divide the track in two parts: below and above the reference plane
64     * ==================================================================
65     IUPDOWN=0
66     DO I=1,NPOINT
67     IF(ZIN(I).LT.ZINI)THEN
68     if(i.ne.1)IUPDOWN=I
69     GOTO 88
70     ENDIF
71     IUPDOWN=NPOINT+1
72     ENDDO
73     88 CONTINUE
74    
75     * ==================================================================
76     * track from reference plane DOWN
77     * ==================================================================
78     * parameters for GRKUTA
79     IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
80     IF(AL_P(5).EQ.0) CHARGE=1.
81     VOUT(1)=AL_P(1)
82     VOUT(2)=AL_P(2)
83     VOUT(3)=ZINI
84     VOUT(4)=AL_P(3)*DCOS(AL_P(4))
85     VOUT(5)=AL_P(3)*DSIN(AL_P(4))
86     VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
87     IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
88     IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
89 pam-fi 1.3 DO I=MAX(IUPDOWN,1),NPOINT
90     L = 0.0
91     10 DO J=1,7
92 pam-fi 1.2 VECT(J)=VOUT(J)
93     VECTINI(J)=VOUT(J)
94     ENDDO
95 pam-fi 1.3 IF(VOUT(6).GE.0.) THEN
96     IFAIL=1
97     c$$$ if(TRKVERBOSE)
98     c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
99     RETURN
100     ENDIF
101     step=(zin(i)-vect(3))/VOUT(6)
102 pam-fi 1.2 11 continue
103 pam-fi 1.3 CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
104     L = L + STEP
105 pam-fi 1.2 IF(VOUT(3).GT.VECT(3)) THEN
106     IFAIL=1
107 pam-fi 1.3 c$$$ if(TRKVERBOSE)
108     c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
109     c$$$ if(TRKVERBOSE)print*,'charge',charge
110     c$$$ if(TRKVERBOSE)print*,'vect',vect
111     c$$$ if(TRKVERBOSE)print*,'vout',vout
112     c$$$ if(TRKVERBOSE)print*,'step',step
113     c$$$ if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
114 pam-fi 1.2 RETURN
115     ENDIF
116     Z=VOUT(3)
117     IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
118     IF(Z.GT.ZIN(I)+TOLL) GOTO 10
119     IF(Z.LE.ZIN(I)-TOLL) THEN
120 pam-fi 1.3 L = L - STEP
121 pam-fi 1.2 STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
122 pam-fi 1.3 DO J=1,7
123 pam-fi 1.2 VECT(J)=VECTINI(J)
124     ENDDO
125     GOTO 11
126     ENDIF
127     100 XOUT(I)=VOUT(1)
128     YOUT(I)=VOUT(2)
129     ZIN(I)=VOUT(3)
130     IF(VOUT(3).ne.0)THEN
131     THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
132     THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
133     ENDIF
134 pam-fi 1.3 TLOUT(I) = L
135 pam-fi 1.2 c print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP
136     ENDDO
137    
138     * ==================================================================
139     * track from reference plane UP
140     * ==================================================================
141     * parameters for GRKUTA:
142     * -opposite charge
143     * -opposite momentum direction
144     IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5))
145     IF(AL_P(5).EQ.0) CHARGE=-1.
146     VOUT(1)=AL_P(1)
147     VOUT(2)=AL_P(2)
148     VOUT(3)=ZINI
149     VOUT(4)=-AL_P(3)*DCOS(AL_P(4))
150     VOUT(5)=-AL_P(3)*DSIN(AL_P(4))
151     VOUT(6)=1.*DSQRT(1.-AL_P(3)**2)
152     IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
153     IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
154     DO I=MIN((IUPDOWN-1),NPOINT),1,-1
155 pam-fi 1.3 L = 0.0
156    
157     20 DO J=1,7
158 pam-fi 1.2 VECT(J)=VOUT(J)
159     VECTINI(J)=VOUT(J)
160     ENDDO
161 pam-fi 1.3 IF(VOUT(6).LE.0.) THEN
162     IFAIL=1
163     c$$$ if(TRKVERBOSE)
164     c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
165     RETURN
166     ENDIF
167     step=(zin(i)-vect(3))/VOUT(6)
168 pam-fi 1.2 22 continue
169 pam-fi 1.3 CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
170     L = L + STEP
171 pam-fi 1.2 IF(VOUT(3).LT.VECT(3)) THEN
172     IFAIL=1
173 pam-fi 1.3 c$$$ if(TRKVERBOSE)
174     c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
175     c$$$ if(TRKVERBOSE)print*,'charge',charge
176     c$$$ if(TRKVERBOSE)print*,'vect',vect
177     c$$$ if(TRKVERBOSE)print*,'vout',vout
178     c$$$ if(TRKVERBOSE)print*,'step',step
179     c$$$ if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
180 pam-fi 1.2 RETURN
181     ENDIF
182     Z=VOUT(3)
183     IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
184     IF(Z.LT.ZIN(I)-TOLL) GOTO 20
185     IF(Z.GE.ZIN(I)+TOLL) THEN
186 pam-fi 1.3 L = L - STEP
187 pam-fi 1.2 STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
188 pam-fi 1.3 DO J=1,7
189 pam-fi 1.2 VECT(J)=VECTINI(J)
190     ENDDO
191     GOTO 22
192     ENDIF
193     200 XOUT(I)=VOUT(1)
194     YOUT(I)=VOUT(2)
195     ZIN(I)=VOUT(3)
196     IF(VOUT(3).ne.0)THEN
197     THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
198     THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
199     ENDIF
200 pam-fi 1.3 TLOUT(I) = L
201 pam-fi 1.2 ENDDO
202    
203     RETURN
204     END
205    
206     ************************************************************************
207     *
208     *
209     *
210     *
211     ************************************************************************
212 mocchiut 1.1 SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)
213    
214     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
215    
216     * -----------------------------------------------
217     * I/O parameters
218     PARAMETER (NPOINT_MAX=100)
219     DIMENSION ZIN(NPOINT_MAX)
220     DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
221 pam-fi 1.3 DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
222     DIMENSION TLOUT(NPOINT_MAX)
223 mocchiut 1.1 DIMENSION AL_P(5)
224     * -----------------------------------------------
225    
226 pam-fi 1.3 CALL DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
227 mocchiut 1.1
228 pam-fi 1.3 ccc
229     ccc OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD
230     ccc
231     c$$$ DIMENSION VECT(7),VECTINI(7),VOUT(7)
232     c$$$ DATA TOLL/1.d-8/
233     c$$$* tolerance in reaching the next plane during the tracking procedure
234     c$$$* -----------------------------------------------
235     c$$$* I/O parameters
236     c$$$ PARAMETER (NPOINT_MAX=100)
237     c$$$ DIMENSION ZIN(NPOINT_MAX)
238     c$$$ DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
239     c$$$ DIMENSION AL_P(5)
240     c$$$* -----------------------------------------------
241     c$$$ DATA ZINI/23.5/ !z coordinate of the reference plane
242     c$$$
243     c$$$* ==================================================================
244     c$$$* divide the track in two parts: below and above the reference plane
245     c$$$* ==================================================================
246     c$$$ IUPDOWN=0
247     c$$$ DO I=1,NPOINT
248     c$$$ IF(ZIN(I).LT.ZINI)THEN
249     c$$$ if(i.ne.1)IUPDOWN=I
250     c$$$ GOTO 88
251     c$$$ ENDIF
252     c$$$ IUPDOWN=NPOINT+1
253     c$$$ ENDDO
254     c$$$ 88 CONTINUE
255     c$$$
256     c$$$* ==================================================================
257     c$$$* track from reference plane DOWN
258     c$$$* ==================================================================
259     c$$$* parameters for GRKUTA
260     c$$$ IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
261     c$$$ IF(AL_P(5).EQ.0) CHARGE=1.
262     c$$$ VOUT(1)=AL_P(1)
263     c$$$ VOUT(2)=AL_P(2)
264     c$$$ VOUT(3)=ZINI
265     c$$$ VOUT(4)=AL_P(3)*DCOS(AL_P(4))
266     c$$$ VOUT(5)=AL_P(3)*DSIN(AL_P(4))
267     c$$$ VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
268     c$$$ IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
269     c$$$ IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
270     c$$$ DO I=MAX(IUPDOWN,1),NPOINT
271     c$$$ step=vout(3)-zin(i)
272     c$$$c$$$ print*,'DOWN ',i,' - Track from ',
273     c$$$c$$$ $ vout(3),' to ',zin(i),' step ',step
274     c$$$ 10 DO J=1,7
275     c$$$ VECT(J)=VOUT(J)
276     c$$$ VECTINI(J)=VOUT(J)
277     c$$$ ENDDO
278     c$$$ 11 continue
279     c$$$ CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
280     c$$$ IF(VOUT(3).GT.VECT(3)) THEN
281     c$$$ IFAIL=1
282     c$$$c$$$ PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'
283     c$$$c$$$ print*,'charge',charge
284     c$$$c$$$ print*,'vect',vect
285     c$$$c$$$ print*,'vout',vout
286     c$$$c$$$ print*,'step',step
287     c$$$ RETURN
288     c$$$ ENDIF
289     c$$$ Z=VOUT(3)
290     c$$$ IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
291     c$$$ IF(Z.GT.ZIN(I)+TOLL) GOTO 10
292     c$$$ IF(Z.LE.ZIN(I)-TOLL) THEN
293     c$$$ STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
294     c$$$ DO J=1,7
295     c$$$ VECT(J)=VECTINI(J)
296     c$$$ ENDDO
297     c$$$ GOTO 11
298     c$$$ ENDIF
299     c$$$ 100 XOUT(I)=VOUT(1)
300     c$$$ YOUT(I)=VOUT(2)
301     c$$$ ZIN(I)=VOUT(3)
302     c$$$* THXOUT(I)=
303     c$$$* THYOUT(I)=
304     c$$$ ENDDO
305     c$$$
306     c$$$
307     c$$$
308     c$$$* ==================================================================
309     c$$$* track from reference plane UP
310     c$$$* ==================================================================
311     c$$$* parameters for GRKUTA:
312     c$$$* -opposite charge
313     c$$$* -opposite momentum direction
314     c$$$ IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5))
315     c$$$ IF(AL_P(5).EQ.0) CHARGE=-1.
316     c$$$ VOUT(1)=AL_P(1)
317     c$$$ VOUT(2)=AL_P(2)
318     c$$$ VOUT(3)=ZINI
319     c$$$ VOUT(4)=-AL_P(3)*DCOS(AL_P(4))
320     c$$$ VOUT(5)=-AL_P(3)*DSIN(AL_P(4))
321     c$$$ VOUT(6)=1.*DSQRT(1.-AL_P(3)**2)
322     c$$$ IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
323     c$$$ IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
324     c$$$ DO I=MIN((IUPDOWN-1),NPOINT),1,-1
325     c$$$ step=vout(3)-zin(i)
326     c$$$ step = -step
327     c$$$c$$$ print*,'UP ',i,' - Track from ',
328     c$$$c$$$ $ vout(3),' to ',zin(i),' step ',step
329     c$$$ 20 DO J=1,7
330     c$$$ VECT(J)=VOUT(J)
331     c$$$ VECTINI(J)=VOUT(J)
332     c$$$ ENDDO
333     c$$$ 22 continue
334     c$$$ CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
335     c$$$ IF(VOUT(3).LT.VECT(3)) THEN
336     c$$$ IFAIL=1
337     c$$$c$$$ PRINT *,'=== WARNING ===> tracciamento invertito (UP)'
338     c$$$c$$$ print*,'charge',charge
339     c$$$c$$$ print*,'vect',vect
340     c$$$c$$$ print*,'vout',vout
341     c$$$c$$$ print*,'step',step
342     c$$$ RETURN
343     c$$$ ENDIF
344     c$$$ Z=VOUT(3)
345     c$$$ IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
346     c$$$ IF(Z.LT.ZIN(I)-TOLL) GOTO 20
347     c$$$ IF(Z.GE.ZIN(I)+TOLL) THEN
348     c$$$ STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
349     c$$$ DO J=1,7
350     c$$$ VECT(J)=VECTINI(J)
351     c$$$ ENDDO
352     c$$$ GOTO 22
353     c$$$ ENDIF
354     c$$$ 200 XOUT(I)=VOUT(1)
355     c$$$ YOUT(I)=VOUT(2)
356     c$$$ ZIN(I)=VOUT(3)
357     c$$$
358     c$$$ ENDDO
359 mocchiut 1.1
360     RETURN
361     END

  ViewVC Help
Powered by ViewVC 1.1.23