/[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.3 by pam-fi, Wed Mar 5 17:00:20 2008 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  *     modified on 05/2006 to give as output also the angles and  *     DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)  
26  *     theh track-lengths  *     - old routine
27    *    
28    *     DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
29    *     - 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  *  *
37    *     March 2008 --> optimized by Paolo
38  **************************************************************  **************************************************************
39    
40        SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)        SUBROUTINE
41         $     DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
42    
43        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
44    
45    
46        DIMENSION VECT(7),VECTINI(7),VOUT(7)        DIMENSION VECT(7),VECTINI(7),VOUT(7)
47        DATA TOLL/1.d-8/                  DATA TOLL/1.d-8/          
48    
49  *     tolerance in reaching the next plane during the tracking procedure  *     tolerance in reaching the next plane during the tracking procedure
50  *     -----------------------------------------------  *     -----------------------------------------------
51  *     I/O parameters  *     I/O parameters
# Line 38  Line 53 
53        DIMENSION ZIN(NPOINT_MAX)        DIMENSION ZIN(NPOINT_MAX)
54        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
55        DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)        DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
56          DIMENSION TLOUT(NPOINT_MAX)
57        DIMENSION AL_P(5)        DIMENSION AL_P(5)
58  *     -----------------------------------------------        *     -----------------------------------------------      
59        DATA ZINI/23.5/           !z coordinate of the reference plane        DATA ZINI/23.5/           !z coordinate of the reference plane
60                REAL*8 L
61    
62  *     ==================================================================  *     ==================================================================
63  *     divide the track in two parts: below and above the reference plane  *     divide the track in two parts: below and above the reference plane
64  *     ==================================================================        *     ==================================================================      
# Line 69  Line 86 
86        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
87        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))
88        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
89        DO I=MAX(IUPDOWN,1),NPOINT        DO I=MAX(IUPDOWN,1),NPOINT
90           step=vout(3)-zin(i)           L = 0.0
 c$$$         print*,'DOWN ',i,' - Track from ',  
 c$$$     $        vout(3),' to ',zin(i),' step ',step  
91   10      DO J=1,7   10      DO J=1,7
92              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
93              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
94           ENDDO           ENDDO
95             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   11      continue   11      continue
103           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
104             L = L + STEP
105           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
106              IFAIL=1              IFAIL=1
107              PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'  c$$$            if(TRKVERBOSE)
108              print*,'charge',charge  c$$$     $      PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
109              print*,'vect',vect  c$$$            if(TRKVERBOSE)print*,'charge',charge
110              print*,'vout',vout  c$$$            if(TRKVERBOSE)print*,'vect',vect
111              print*,'step',step  c$$$            if(TRKVERBOSE)print*,'vout',vout
112    c$$$            if(TRKVERBOSE)print*,'step',step
113    c$$$            if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
114              RETURN              RETURN
115           ENDIF           ENDIF
116           Z=VOUT(3)           Z=VOUT(3)
117           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
118           IF(Z.GT.ZIN(I)+TOLL) GOTO 10                 IF(Z.GT.ZIN(I)+TOLL) GOTO 10      
119           IF(Z.LE.ZIN(I)-TOLL) THEN           IF(Z.LE.ZIN(I)-TOLL) THEN
120                L = L - STEP
121              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
122              DO J=1,7              DO J=1,7
123                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
# Line 101  c$$$     $        vout(3),' to ',zin(i), Line 127  c$$$     $        vout(3),' to ',zin(i),
127   100     XOUT(I)=VOUT(1)   100     XOUT(I)=VOUT(1)
128           YOUT(I)=VOUT(2)           YOUT(I)=VOUT(2)
129           ZIN(I)=VOUT(3)           ZIN(I)=VOUT(3)
130  *         THXOUT(I)=           IF(VOUT(3).ne.0)THEN
131  *         THYOUT(I)=              THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
132                THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
133             ENDIF
134             TLOUT(I) = L
135    c         print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP
136        ENDDO        ENDDO
137    
   
   
138  *     ==================================================================  *     ==================================================================
139  *     track from reference plane UP  *     track from reference plane UP
140  *     ==================================================================  *     ==================================================================
# Line 124  c$$$     $        vout(3),' to ',zin(i), Line 152  c$$$     $        vout(3),' to ',zin(i),
152        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))
153        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
154        DO I=MIN((IUPDOWN-1),NPOINT),1,-1        DO I=MIN((IUPDOWN-1),NPOINT),1,-1
155           step=vout(3)-zin(i)           L = 0.0
156           step = -step  
 c$$$         print*,'UP ',i,' - Track from ',  
 c$$$     $        vout(3),' to ',zin(i),' step ',step  
157   20      DO J=1,7   20      DO J=1,7
158              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
159              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
160           ENDDO           ENDDO
161             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   22      continue   22      continue
169           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
170             L = L + STEP
171           IF(VOUT(3).LT.VECT(3)) THEN           IF(VOUT(3).LT.VECT(3)) THEN
172              IFAIL=1              IFAIL=1
173              PRINT *,'=== WARNING ===> tracciamento invertito (UP)'  c$$$            if(TRKVERBOSE)
174              print*,'charge',charge  c$$$     $      PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
175              print*,'vect',vect  c$$$            if(TRKVERBOSE)print*,'charge',charge
176              print*,'vout',vout  c$$$            if(TRKVERBOSE)print*,'vect',vect
177              print*,'step',step  c$$$            if(TRKVERBOSE)print*,'vout',vout
178    c$$$            if(TRKVERBOSE)print*,'step',step
179    c$$$            if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
180              RETURN              RETURN
181           ENDIF           ENDIF
182           Z=VOUT(3)           Z=VOUT(3)
183           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
184           IF(Z.LT.ZIN(I)-TOLL) GOTO 20                 IF(Z.LT.ZIN(I)-TOLL) GOTO 20      
185           IF(Z.GE.ZIN(I)+TOLL) THEN           IF(Z.GE.ZIN(I)+TOLL) THEN
186                L = L - STEP
187              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
188              DO J=1,7              DO J=1,7
189                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
# Line 156  c$$$     $        vout(3),' to ',zin(i), Line 193  c$$$     $        vout(3),' to ',zin(i),
193   200     XOUT(I)=VOUT(1)   200     XOUT(I)=VOUT(1)
194           YOUT(I)=VOUT(2)           YOUT(I)=VOUT(2)
195           ZIN(I)=VOUT(3)           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             TLOUT(I) = L
201        ENDDO        ENDDO
202    
203        RETURN        RETURN
204        END        END
205    
206    ************************************************************************
207    *
208    *
209    *
210    *
211    ************************************************************************
212          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          DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
222          DIMENSION TLOUT(NPOINT_MAX)
223          DIMENSION AL_P(5)
224    *     -----------------------------------------------      
225    
226          CALL DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
227    
228    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    
360          RETURN
361          END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23