/[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.2 by pam-fi, Thu Jun 1 09:03:09 2006 UTC revision 1.5 by pam-ts, Wed Jun 4 07:57:04 2014 UTC
# Line 1  Line 1 
1    
2    
3  *************************************************************  *************************************************************
4  *  *
5  *     Routines to compute the NPOINT track intersection points  *     Routines to compute the NPOINT track intersection points
# Line 24  Line 25 
25  *     DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)    *     DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)  
26  *     - old routine  *     - old routine
27  *      *    
28  *     DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,AL_P,IFAIL)  *     DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
29  *     - as the old one, but:  *     - as the previous one, but:
30  *     -- the projected angles are given as output  *     -- the projected angles are given as output
31  *     -- the track lengths are given as output:  *     -- the track lengths are given as output:
32  *     ---- for planes above the reference plane, the length until  *     ---- for planes above the reference plane, the length until
# Line 33  Line 34 
34  *     ---- for planes below the reference plane, the length until  *     ---- for planes below the reference plane, the length until
35  *          the higher closest plane (reference plane included)  *          the higher closest plane (reference plane included)
36  *  *
37    *     March 2008 --> optimized by Paolo
38    *
39    *     DOTRACK3(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,ZINI,IFAIL)
40    *     - as the previous one, but with the coordinate ZINI of the reference
41    *       plane given as input parameter
42    *
43    *
44  **************************************************************  **************************************************************
45    
46    
47        SUBROUTINE        SUBROUTINE
48       $     DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)       $     DOTRACK3(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P
49         $     ,ZINI,IFAIL)
50    
51        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
52    
53        DIMENSION VECT(8),VECTINI(8),VOUT(8)  
54          DIMENSION VECT(7),VECTINI(7),VOUT(7)
55        DATA TOLL/1.d-8/                  DATA TOLL/1.d-8/          
56    
57  *     tolerance in reaching the next plane during the tracking procedure  *     tolerance in reaching the next plane during the tracking procedure
58  *     -----------------------------------------------  *     -----------------------------------------------
59  *     I/O parameters  *     I/O parameters
60        PARAMETER (NPOINT_MAX=100)        PARAMETER (NPOINT_MAX=500)
61        DIMENSION ZIN(NPOINT_MAX)        DIMENSION ZIN(NPOINT_MAX)
62        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
63        DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)        DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
64        DIMENSION TLOUT(NPOINT_MAX)        DIMENSION TLOUT(NPOINT_MAX)
65        DIMENSION AL_P(5)        DIMENSION AL_P(5)
66  *     -----------------------------------------------        *     -----------------------------------------------      
67        DATA ZINI/23.5/           !z coordinate of the reference plane  *      DATA ZINI/23.5/           !z coordinate of the reference plane
68                REAL*8 L
69    
70  *     ==================================================================  *     ==================================================================
71  *     divide the track in two parts: below and above the reference plane  *     divide the track in two parts: below and above the reference plane
72  *     ==================================================================        *     ==================================================================      
# Line 81  Line 94 
94        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
95        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))
96        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
97        DO I=MAX(IUPDOWN,1),NPOINT        DO I=MAX(IUPDOWN,1),NPOINT
98           step=vout(3)-zin(i)           L = 0.0
99           VOUT(8)= 0   10      DO J=1,7
 c$$$         print*,'DOWN ',i,' - Track from ',  
 c$$$     $        vout(3),' to ',zin(i),' step ',step  
  10      DO J=1,8  
100              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
101              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
102           ENDDO           ENDDO
103             IF(VOUT(6).GE.0.) THEN
104                IFAIL=1
105    c$$$            if(TRKVERBOSE)      
106    c$$$     $           PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
107                RETURN
108             ENDIF
109             step=(zin(i)-vect(3))/VOUT(6)
110   11      continue   11      continue
111           CALL GRKUTA2(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
112             L = L + STEP
113           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
114              IFAIL=1              IFAIL=1
115  c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'  c$$$            if(TRKVERBOSE)
116  c$$$            print*,'charge',charge  c$$$     $      PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
117  c$$$            print*,'vect',vect  c$$$            if(TRKVERBOSE)print*,'charge',charge
118  c$$$            print*,'vout',vout  c$$$            if(TRKVERBOSE)print*,'vect',vect
119  c$$$            print*,'step',step  c$$$            if(TRKVERBOSE)print*,'vout',vout
120    c$$$            if(TRKVERBOSE)print*,'step',step
121    c$$$            if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
122              RETURN              RETURN
123           ENDIF           ENDIF
124           Z=VOUT(3)           Z=VOUT(3)
125           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
126           IF(Z.GT.ZIN(I)+TOLL) GOTO 10                 IF(Z.GT.ZIN(I)+TOLL) GOTO 10      
127           IF(Z.LE.ZIN(I)-TOLL) THEN           IF(Z.LE.ZIN(I)-TOLL) THEN
128                L = L - STEP
129              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
130              DO J=1,8              DO J=1,7
131                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
132              ENDDO              ENDDO
133              GOTO 11              GOTO 11
# Line 118  c$$$            print*,'step',step Line 139  c$$$            print*,'step',step
139              THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)              THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
140              THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)              THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
141           ENDIF           ENDIF
142           TLOUT(I) = VOUT(8)           TLOUT(I) = L
143  c         print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP  c         print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP
144        ENDDO        ENDDO
145    
   
   
146  *     ==================================================================  *     ==================================================================
147  *     track from reference plane UP  *     track from reference plane UP
148  *     ==================================================================  *     ==================================================================
# Line 141  c         print*,'D ',VECT(3),VOUT(3),VO Line 160  c         print*,'D ',VECT(3),VOUT(3),VO
160        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))
161        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
162        DO I=MIN((IUPDOWN-1),NPOINT),1,-1        DO I=MIN((IUPDOWN-1),NPOINT),1,-1
163           step=vout(3)-zin(i)           L = 0.0
164           step = -step  
165           VOUT(8)=0   20      DO J=1,7
 c$$$         print*,'UP ',i,' - Track from ',  
 c$$$     $        vout(3),' to ',zin(i),' step ',step  
  20      DO J=1,8  
166              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
167              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
168           ENDDO           ENDDO
169             IF(VOUT(6).LE.0.) THEN
170                IFAIL=1
171    c$$$            if(TRKVERBOSE)      
172    c$$$     $           PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
173                RETURN
174             ENDIF
175             step=(zin(i)-vect(3))/VOUT(6)
176   22      continue   22      continue
177           CALL GRKUTA2(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
178             L = L + STEP
179           IF(VOUT(3).LT.VECT(3)) THEN           IF(VOUT(3).LT.VECT(3)) THEN
180              IFAIL=1              IFAIL=1
181  c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (UP)'  c$$$            if(TRKVERBOSE)
182  c$$$            print*,'charge',charge  c$$$     $      PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
183  c$$$            print*,'vect',vect  c$$$            if(TRKVERBOSE)print*,'charge',charge
184  c$$$            print*,'vout',vout  c$$$            if(TRKVERBOSE)print*,'vect',vect
185  c$$$            print*,'step',step  c$$$            if(TRKVERBOSE)print*,'vout',vout
186    c$$$            if(TRKVERBOSE)print*,'step',step
187    c$$$            if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
188              RETURN              RETURN
189           ENDIF           ENDIF
190           Z=VOUT(3)           Z=VOUT(3)
191           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
192           IF(Z.LT.ZIN(I)-TOLL) GOTO 20                 IF(Z.LT.ZIN(I)-TOLL) GOTO 20      
193           IF(Z.GE.ZIN(I)+TOLL) THEN           IF(Z.GE.ZIN(I)+TOLL) THEN
194                L = L - STEP
195              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))              STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
196              DO J=1,8              DO J=1,7
197                 VECT(J)=VECTINI(J)                 VECT(J)=VECTINI(J)
198              ENDDO              ENDDO
199              GOTO 22              GOTO 22
# Line 178  c$$$            print*,'step',step Line 205  c$$$            print*,'step',step
205              THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)              THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
206              THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)              THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
207           ENDIF           ENDIF
208           TLOUT(I) = VOUT(8)           TLOUT(I) = L
   
209        ENDDO        ENDDO
210    
211        RETURN        RETURN
212        END        END
   
213  ************************************************************************  ************************************************************************
214  *  *
215  *  *
216  *  *
217  *  *
218  ************************************************************************  ************************************************************************
219        SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)  
220          SUBROUTINE
221         $     DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
222    
223        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
224    
225        DIMENSION VECT(7),VECTINI(7),VOUT(7)  
226        DATA TOLL/1.d-8/            ccc      DIMENSION VECT(7),VECTINI(7),VOUT(7)
227    ccc      DATA TOLL/1.d-8/          
228    
229  *     tolerance in reaching the next plane during the tracking procedure  *     tolerance in reaching the next plane during the tracking procedure
230  *     -----------------------------------------------  *     -----------------------------------------------
231  *     I/O parameters  *     I/O parameters
232        PARAMETER (NPOINT_MAX=100)        PARAMETER (NPOINT_MAX=100)
233        DIMENSION ZIN(NPOINT_MAX)        DIMENSION ZIN(NPOINT_MAX)
234        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)        DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
235          DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
236          DIMENSION TLOUT(NPOINT_MAX)
237        DIMENSION AL_P(5)        DIMENSION AL_P(5)
238  *     -----------------------------------------------        *     -----------------------------------------------      
239        DATA ZINI/23.5/           !z coordinate of the reference plane        DATA ZINI/23.5/           !z coordinate of the reference plane
240          ccc      REAL*8 L
 *     ==================================================================  
 *     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  
241    
242  *     ==================================================================        CALL DOTRACK3(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P
243  *     track from reference plane DOWN       $     ,ZINI,IFAIL)
 *     ==================================================================  
 *     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  
 c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'  
 c$$$            print*,'charge',charge  
 c$$$            print*,'vect',vect  
 c$$$            print*,'vout',vout  
 c$$$            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  
244    
245    ccc
246    ccc   OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD
247    ccc
248    
249    c$$$*     ==================================================================
250    c$$$*     divide the track in two parts: below and above the reference plane
251    c$$$*     ==================================================================      
252    c$$$      IUPDOWN=0
253    c$$$      DO I=1,NPOINT  
254    c$$$         IF(ZIN(I).LT.ZINI)THEN
255    c$$$            if(i.ne.1)IUPDOWN=I
256    c$$$            GOTO 88
257    c$$$         ENDIF
258    c$$$         IUPDOWN=NPOINT+1
259    c$$$      ENDDO
260    c$$$ 88   CONTINUE
261    c$$$
262    c$$$*     ==================================================================
263    c$$$*     track from reference plane DOWN
264    c$$$*     ==================================================================
265    c$$$*     parameters for GRKUTA
266    c$$$      IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
267    c$$$      IF(AL_P(5).EQ.0) CHARGE=1.
268    c$$$      VOUT(1)=AL_P(1)
269    c$$$      VOUT(2)=AL_P(2)
270    c$$$      VOUT(3)=ZINI            
271    c$$$      VOUT(4)=AL_P(3)*DCOS(AL_P(4))
272    c$$$      VOUT(5)=AL_P(3)*DSIN(AL_P(4))
273    c$$$      VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
274    c$$$      IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
275    c$$$      IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
276    c$$$      DO I=MAX(IUPDOWN,1),NPOINT
277    c$$$         L = 0.0
278    c$$$ 10      DO J=1,7
279    c$$$            VECT(J)=VOUT(J)
280    c$$$            VECTINI(J)=VOUT(J)
281    c$$$         ENDDO
282    c$$$         IF(VOUT(6).GE.0.) THEN
283    c$$$            IFAIL=1
284    c$$$c$$$            if(TRKVERBOSE)      
285    c$$$c$$$     $           PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
286    c$$$            RETURN
287    c$$$         ENDIF
288    c$$$         step=(zin(i)-vect(3))/VOUT(6)
289    c$$$ 11      continue
290    c$$$         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
291    c$$$         L = L + STEP
292    c$$$         IF(VOUT(3).GT.VECT(3)) THEN
293    c$$$            IFAIL=1
294    c$$$c$$$            if(TRKVERBOSE)
295    c$$$c$$$     $      PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
296    c$$$c$$$            if(TRKVERBOSE)print*,'charge',charge
297    c$$$c$$$            if(TRKVERBOSE)print*,'vect',vect
298    c$$$c$$$            if(TRKVERBOSE)print*,'vout',vout
299    c$$$c$$$            if(TRKVERBOSE)print*,'step',step
300    c$$$c$$$            if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
301    c$$$            RETURN
302    c$$$         ENDIF
303    c$$$         Z=VOUT(3)
304    c$$$         IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
305    c$$$         IF(Z.GT.ZIN(I)+TOLL) GOTO 10      
306    c$$$         IF(Z.LE.ZIN(I)-TOLL) THEN
307    c$$$            L = L - STEP
308    c$$$            STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
309    c$$$            DO J=1,7
310    c$$$               VECT(J)=VECTINI(J)
311    c$$$            ENDDO
312    c$$$            GOTO 11
313    c$$$         ENDIF
314    c$$$ 100     XOUT(I)=VOUT(1)
315    c$$$         YOUT(I)=VOUT(2)
316    c$$$         ZIN(I)=VOUT(3)
317    c$$$         IF(VOUT(3).ne.0)THEN
318    c$$$            THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
319    c$$$            THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
320    c$$$         ENDIF
321    c$$$         TLOUT(I) = L
322    c$$$c         print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP
323    c$$$      ENDDO
324    c$$$
325    c$$$*     ==================================================================
326    c$$$*     track from reference plane UP
327    c$$$*     ==================================================================
328    c$$$*     parameters for GRKUTA:
329    c$$$*     -opposite charge
330    c$$$*     -opposite momentum direction
331    c$$$      IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5))
332    c$$$      IF(AL_P(5).EQ.0) CHARGE=-1.
333    c$$$      VOUT(1)=AL_P(1)
334    c$$$      VOUT(2)=AL_P(2)
335    c$$$      VOUT(3)=ZINI            
336    c$$$      VOUT(4)=-AL_P(3)*DCOS(AL_P(4))
337    c$$$      VOUT(5)=-AL_P(3)*DSIN(AL_P(4))
338    c$$$      VOUT(6)=1.*DSQRT(1.-AL_P(3)**2)
339    c$$$      IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
340    c$$$      IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
341    c$$$      DO I=MIN((IUPDOWN-1),NPOINT),1,-1
342    c$$$         L = 0.0
343    c$$$
344    c$$$ 20      DO J=1,7
345    c$$$            VECT(J)=VOUT(J)
346    c$$$            VECTINI(J)=VOUT(J)
347    c$$$         ENDDO
348    c$$$         IF(VOUT(6).LE.0.) THEN
349    c$$$            IFAIL=1
350    c$$$c$$$            if(TRKVERBOSE)      
351    c$$$c$$$     $           PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
352    c$$$            RETURN
353    c$$$         ENDIF
354    c$$$         step=(zin(i)-vect(3))/VOUT(6)
355    c$$$ 22      continue
356    c$$$         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
357    c$$$         L = L + STEP
358    c$$$         IF(VOUT(3).LT.VECT(3)) THEN
359    c$$$            IFAIL=1
360    c$$$c$$$            if(TRKVERBOSE)
361    c$$$c$$$     $      PRINT *,'dofit (grkuta): WARNING ===> backward track!!'
362    c$$$c$$$            if(TRKVERBOSE)print*,'charge',charge
363    c$$$c$$$            if(TRKVERBOSE)print*,'vect',vect
364    c$$$c$$$            if(TRKVERBOSE)print*,'vout',vout
365    c$$$c$$$            if(TRKVERBOSE)print*,'step',step
366    c$$$c$$$            if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1
367    c$$$            RETURN
368    c$$$         ENDIF
369    c$$$         Z=VOUT(3)
370    c$$$         IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
371    c$$$         IF(Z.LT.ZIN(I)-TOLL) GOTO 20      
372    c$$$         IF(Z.GE.ZIN(I)+TOLL) THEN
373    c$$$            L = L - STEP
374    c$$$            STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
375    c$$$            DO J=1,7
376    c$$$               VECT(J)=VECTINI(J)
377    c$$$            ENDDO
378    c$$$            GOTO 22
379    c$$$         ENDIF
380    c$$$ 200     XOUT(I)=VOUT(1)
381    c$$$         YOUT(I)=VOUT(2)
382    c$$$         ZIN(I)=VOUT(3)
383    c$$$         IF(VOUT(3).ne.0)THEN
384    c$$$            THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
385    c$$$            THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
386    c$$$         ENDIF
387    c$$$         TLOUT(I) = L
388    c$$$      ENDDO
389    c$$$
390          RETURN
391          END
392    
393    ************************************************************************
394    *
395    *
396    *
397    *
398    ************************************************************************
399          SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)
400    
401  *     ==================================================================        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 *     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  
 c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (UP)'  
 c$$$            print*,'charge',charge  
 c$$$            print*,'vect',vect  
 c$$$            print*,'vout',vout  
 c$$$            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)  
402    
403        ENDDO  *     -----------------------------------------------
404    *     I/O parameters
405          PARAMETER (NPOINT_MAX=100)
406          DIMENSION ZIN(NPOINT_MAX)
407          DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
408          DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
409          DIMENSION TLOUT(NPOINT_MAX)
410          DIMENSION AL_P(5)
411    *     -----------------------------------------------      
412    
413          CALL DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
414    
415    ccc
416    ccc   OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD
417    ccc
418    c$$$      DIMENSION VECT(7),VECTINI(7),VOUT(7)
419    c$$$      DATA TOLL/1.d-8/          
420    c$$$*     tolerance in reaching the next plane during the tracking procedure
421    c$$$*     -----------------------------------------------
422    c$$$*     I/O parameters
423    c$$$      PARAMETER (NPOINT_MAX=100)
424    c$$$      DIMENSION ZIN(NPOINT_MAX)
425    c$$$      DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
426    c$$$      DIMENSION AL_P(5)
427    c$$$*     -----------------------------------------------      
428    c$$$      DATA ZINI/23.5/           !z coordinate of the reference plane
429    c$$$      
430    c$$$*     ==================================================================
431    c$$$*     divide the track in two parts: below and above the reference plane
432    c$$$*     ==================================================================      
433    c$$$      IUPDOWN=0
434    c$$$      DO I=1,NPOINT  
435    c$$$         IF(ZIN(I).LT.ZINI)THEN
436    c$$$            if(i.ne.1)IUPDOWN=I
437    c$$$            GOTO 88
438    c$$$         ENDIF
439    c$$$         IUPDOWN=NPOINT+1
440    c$$$      ENDDO
441    c$$$ 88   CONTINUE
442    c$$$
443    c$$$*     ==================================================================
444    c$$$*     track from reference plane DOWN
445    c$$$*     ==================================================================
446    c$$$*     parameters for GRKUTA
447    c$$$      IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
448    c$$$      IF(AL_P(5).EQ.0) CHARGE=1.
449    c$$$      VOUT(1)=AL_P(1)
450    c$$$      VOUT(2)=AL_P(2)
451    c$$$      VOUT(3)=ZINI            
452    c$$$      VOUT(4)=AL_P(3)*DCOS(AL_P(4))
453    c$$$      VOUT(5)=AL_P(3)*DSIN(AL_P(4))
454    c$$$      VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
455    c$$$      IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
456    c$$$      IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
457    c$$$      DO I=MAX(IUPDOWN,1),NPOINT
458    c$$$         step=vout(3)-zin(i)
459    c$$$c$$$         print*,'DOWN ',i,' - Track from ',
460    c$$$c$$$     $        vout(3),' to ',zin(i),' step ',step
461    c$$$ 10      DO J=1,7
462    c$$$            VECT(J)=VOUT(J)
463    c$$$            VECTINI(J)=VOUT(J)
464    c$$$         ENDDO
465    c$$$ 11      continue
466    c$$$         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
467    c$$$         IF(VOUT(3).GT.VECT(3)) THEN
468    c$$$            IFAIL=1
469    c$$$c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)'
470    c$$$c$$$            print*,'charge',charge
471    c$$$c$$$            print*,'vect',vect
472    c$$$c$$$            print*,'vout',vout
473    c$$$c$$$            print*,'step',step
474    c$$$           RETURN
475    c$$$         ENDIF
476    c$$$         Z=VOUT(3)
477    c$$$         IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
478    c$$$         IF(Z.GT.ZIN(I)+TOLL) GOTO 10      
479    c$$$         IF(Z.LE.ZIN(I)-TOLL) THEN
480    c$$$            STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
481    c$$$            DO J=1,7
482    c$$$               VECT(J)=VECTINI(J)
483    c$$$            ENDDO
484    c$$$            GOTO 11
485    c$$$         ENDIF
486    c$$$ 100     XOUT(I)=VOUT(1)
487    c$$$         YOUT(I)=VOUT(2)
488    c$$$         ZIN(I)=VOUT(3)
489    c$$$*         THXOUT(I)=
490    c$$$*         THYOUT(I)=
491    c$$$      ENDDO
492    c$$$
493    c$$$
494    c$$$
495    c$$$*     ==================================================================
496    c$$$*     track from reference plane UP
497    c$$$*     ==================================================================
498    c$$$*     parameters for GRKUTA:
499    c$$$*     -opposite charge
500    c$$$*     -opposite momentum direction
501    c$$$      IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5))
502    c$$$      IF(AL_P(5).EQ.0) CHARGE=-1.
503    c$$$      VOUT(1)=AL_P(1)
504    c$$$      VOUT(2)=AL_P(2)
505    c$$$      VOUT(3)=ZINI            
506    c$$$      VOUT(4)=-AL_P(3)*DCOS(AL_P(4))
507    c$$$      VOUT(5)=-AL_P(3)*DSIN(AL_P(4))
508    c$$$      VOUT(6)=1.*DSQRT(1.-AL_P(3)**2)
509    c$$$      IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
510    c$$$      IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
511    c$$$      DO I=MIN((IUPDOWN-1),NPOINT),1,-1
512    c$$$         step=vout(3)-zin(i)
513    c$$$         step = -step
514    c$$$c$$$         print*,'UP ',i,' - Track from ',
515    c$$$c$$$     $        vout(3),' to ',zin(i),' step ',step
516    c$$$ 20      DO J=1,7
517    c$$$            VECT(J)=VOUT(J)
518    c$$$            VECTINI(J)=VOUT(J)
519    c$$$         ENDDO
520    c$$$ 22      continue
521    c$$$         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
522    c$$$         IF(VOUT(3).LT.VECT(3)) THEN
523    c$$$            IFAIL=1
524    c$$$c$$$            PRINT *,'=== WARNING ===> tracciamento invertito (UP)'
525    c$$$c$$$            print*,'charge',charge
526    c$$$c$$$            print*,'vect',vect
527    c$$$c$$$            print*,'vout',vout
528    c$$$c$$$            print*,'step',step
529    c$$$            RETURN
530    c$$$         ENDIF
531    c$$$         Z=VOUT(3)
532    c$$$         IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
533    c$$$         IF(Z.LT.ZIN(I)-TOLL) GOTO 20      
534    c$$$         IF(Z.GE.ZIN(I)+TOLL) THEN
535    c$$$            STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
536    c$$$            DO J=1,7
537    c$$$               VECT(J)=VECTINI(J)
538    c$$$            ENDDO
539    c$$$            GOTO 22
540    c$$$         ENDIF
541    c$$$ 200     XOUT(I)=VOUT(1)
542    c$$$         YOUT(I)=VOUT(2)
543    c$$$         ZIN(I)=VOUT(3)
544    c$$$
545    c$$$      ENDDO
546    
547        RETURN        RETURN
548        END        END
   
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23