--- DarthVader/TrackerLevel2/src/F77/track.f 2006/06/01 09:03:09 1.2 +++ DarthVader/TrackerLevel2/src/F77/track.f 2008/03/05 17:00:20 1.3 @@ -1,4 +1,5 @@ + ************************************************************* * * Routines to compute the NPOINT track intersection points @@ -24,7 +25,7 @@ * DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL) * - old routine * -* DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,AL_P,IFAIL) +* DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL) * - as the old one, but: * -- the projected angles are given as output * -- the track lengths are given as output: @@ -33,6 +34,7 @@ * ---- for planes below the reference plane, the length until * the higher closest plane (reference plane included) * +* March 2008 --> optimized by Paolo ************************************************************** SUBROUTINE @@ -40,8 +42,10 @@ IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION VECT(8),VECTINI(8),VOUT(8) + + DIMENSION VECT(7),VECTINI(7),VOUT(7) DATA TOLL/1.d-8/ + * tolerance in reaching the next plane during the tracking procedure * ----------------------------------------------- * I/O parameters @@ -53,7 +57,8 @@ DIMENSION AL_P(5) * ----------------------------------------------- DATA ZINI/23.5/ !z coordinate of the reference plane - + REAL*8 L + * ================================================================== * divide the track in two parts: below and above the reference plane * ================================================================== @@ -81,32 +86,40 @@ 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) - VOUT(8)= 0 -c$$$ print*,'DOWN ',i,' - Track from ', -c$$$ $ vout(3),' to ',zin(i),' step ',step - 10 DO J=1,8 + DO I=MAX(IUPDOWN,1),NPOINT + L = 0.0 + 10 DO J=1,7 VECT(J)=VOUT(J) VECTINI(J)=VOUT(J) ENDDO + IF(VOUT(6).GE.0.) THEN + IFAIL=1 +c$$$ if(TRKVERBOSE) +c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' + RETURN + ENDIF + step=(zin(i)-vect(3))/VOUT(6) 11 continue - CALL GRKUTA2(CHARGE,STEP,VECT,VOUT) + CALL GRKUTA(CHARGE,STEP,VECT,VOUT) + L = L + STEP 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 +c$$$ if(TRKVERBOSE) +c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' +c$$$ if(TRKVERBOSE)print*,'charge',charge +c$$$ if(TRKVERBOSE)print*,'vect',vect +c$$$ if(TRKVERBOSE)print*,'vout',vout +c$$$ if(TRKVERBOSE)print*,'step',step +c$$$ if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1 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 + L = L - STEP STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) - DO J=1,8 + DO J=1,7 VECT(J)=VECTINI(J) ENDDO GOTO 11 @@ -118,12 +131,10 @@ THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.) THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) ENDIF - TLOUT(I) = VOUT(8) + TLOUT(I) = L c print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP ENDDO - - * ================================================================== * track from reference plane UP * ================================================================== @@ -141,32 +152,40 @@ 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 - VOUT(8)=0 -c$$$ print*,'UP ',i,' - Track from ', -c$$$ $ vout(3),' to ',zin(i),' step ',step - 20 DO J=1,8 + L = 0.0 + + 20 DO J=1,7 VECT(J)=VOUT(J) VECTINI(J)=VOUT(J) ENDDO + IF(VOUT(6).LE.0.) THEN + IFAIL=1 +c$$$ if(TRKVERBOSE) +c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' + RETURN + ENDIF + step=(zin(i)-vect(3))/VOUT(6) 22 continue - CALL GRKUTA2(CHARGE,STEP,VECT,VOUT) + CALL GRKUTA(CHARGE,STEP,VECT,VOUT) + L = L + STEP 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 +c$$$ if(TRKVERBOSE) +c$$$ $ PRINT *,'dofit (grkuta): WARNING ===> backward track!!' +c$$$ if(TRKVERBOSE)print*,'charge',charge +c$$$ if(TRKVERBOSE)print*,'vect',vect +c$$$ if(TRKVERBOSE)print*,'vout',vout +c$$$ if(TRKVERBOSE)print*,'step',step +c$$$ if(TRKVERBOSE)print*,'DeltaB',DELTA0,DELTA1 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 + L = L - STEP STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) - DO J=1,8 + DO J=1,7 VECT(J)=VECTINI(J) ENDDO GOTO 22 @@ -178,8 +197,7 @@ THXOUT(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.) THYOUT(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.) ENDIF - TLOUT(I) = VOUT(8) - + TLOUT(I) = L ENDDO RETURN @@ -195,135 +213,149 @@ IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION VECT(7),VECTINI(7),VOUT(7) - DATA TOLL/1.d-8/ -* tolerance in reaching the next plane during the tracking procedure * ----------------------------------------------- * I/O parameters PARAMETER (NPOINT_MAX=100) DIMENSION ZIN(NPOINT_MAX) DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX) + DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX) + DIMENSION TLOUT(NPOINT_MAX) DIMENSION AL_P(5) * ----------------------------------------------- - DATA ZINI/23.5/ !z coordinate of the reference plane - -* ================================================================== -* 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 -* ================================================================== -* track from reference plane DOWN -* ================================================================== -* 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 - - - -* ================================================================== -* 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) + CALL DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL) - ENDDO +ccc +ccc OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD-OLD +ccc +c$$$ DIMENSION VECT(7),VECTINI(7),VOUT(7) +c$$$ DATA TOLL/1.d-8/ +c$$$* tolerance in reaching the next plane during the tracking procedure +c$$$* ----------------------------------------------- +c$$$* I/O parameters +c$$$ PARAMETER (NPOINT_MAX=100) +c$$$ DIMENSION ZIN(NPOINT_MAX) +c$$$ DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX) +c$$$ DIMENSION AL_P(5) +c$$$* ----------------------------------------------- +c$$$ DATA ZINI/23.5/ !z coordinate of the reference plane +c$$$ +c$$$* ================================================================== +c$$$* divide the track in two parts: below and above the reference plane +c$$$* ================================================================== +c$$$ IUPDOWN=0 +c$$$ DO I=1,NPOINT +c$$$ IF(ZIN(I).LT.ZINI)THEN +c$$$ if(i.ne.1)IUPDOWN=I +c$$$ GOTO 88 +c$$$ ENDIF +c$$$ IUPDOWN=NPOINT+1 +c$$$ ENDDO +c$$$ 88 CONTINUE +c$$$ +c$$$* ================================================================== +c$$$* track from reference plane DOWN +c$$$* ================================================================== +c$$$* parameters for GRKUTA +c$$$ IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5)) +c$$$ IF(AL_P(5).EQ.0) CHARGE=1. +c$$$ VOUT(1)=AL_P(1) +c$$$ VOUT(2)=AL_P(2) +c$$$ VOUT(3)=ZINI +c$$$ VOUT(4)=AL_P(3)*DCOS(AL_P(4)) +c$$$ VOUT(5)=AL_P(3)*DSIN(AL_P(4)) +c$$$ VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2) +c$$$ IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) +c$$$ IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 +c$$$ DO I=MAX(IUPDOWN,1),NPOINT +c$$$ step=vout(3)-zin(i) +c$$$c$$$ print*,'DOWN ',i,' - Track from ', +c$$$c$$$ $ vout(3),' to ',zin(i),' step ',step +c$$$ 10 DO J=1,7 +c$$$ VECT(J)=VOUT(J) +c$$$ VECTINI(J)=VOUT(J) +c$$$ ENDDO +c$$$ 11 continue +c$$$ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) +c$$$ IF(VOUT(3).GT.VECT(3)) THEN +c$$$ IFAIL=1 +c$$$c$$$ PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)' +c$$$c$$$ print*,'charge',charge +c$$$c$$$ print*,'vect',vect +c$$$c$$$ print*,'vout',vout +c$$$c$$$ print*,'step',step +c$$$ RETURN +c$$$ ENDIF +c$$$ Z=VOUT(3) +c$$$ IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100 +c$$$ IF(Z.GT.ZIN(I)+TOLL) GOTO 10 +c$$$ IF(Z.LE.ZIN(I)-TOLL) THEN +c$$$ STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) +c$$$ DO J=1,7 +c$$$ VECT(J)=VECTINI(J) +c$$$ ENDDO +c$$$ GOTO 11 +c$$$ ENDIF +c$$$ 100 XOUT(I)=VOUT(1) +c$$$ YOUT(I)=VOUT(2) +c$$$ ZIN(I)=VOUT(3) +c$$$* THXOUT(I)= +c$$$* THYOUT(I)= +c$$$ ENDDO +c$$$ +c$$$ +c$$$ +c$$$* ================================================================== +c$$$* track from reference plane UP +c$$$* ================================================================== +c$$$* parameters for GRKUTA: +c$$$* -opposite charge +c$$$* -opposite momentum direction +c$$$ IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5)) +c$$$ IF(AL_P(5).EQ.0) CHARGE=-1. +c$$$ VOUT(1)=AL_P(1) +c$$$ VOUT(2)=AL_P(2) +c$$$ VOUT(3)=ZINI +c$$$ VOUT(4)=-AL_P(3)*DCOS(AL_P(4)) +c$$$ VOUT(5)=-AL_P(3)*DSIN(AL_P(4)) +c$$$ VOUT(6)=1.*DSQRT(1.-AL_P(3)**2) +c$$$ IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5)) +c$$$ IF(AL_P(5).EQ.0.) VOUT(7)=1.E8 +c$$$ DO I=MIN((IUPDOWN-1),NPOINT),1,-1 +c$$$ step=vout(3)-zin(i) +c$$$ step = -step +c$$$c$$$ print*,'UP ',i,' - Track from ', +c$$$c$$$ $ vout(3),' to ',zin(i),' step ',step +c$$$ 20 DO J=1,7 +c$$$ VECT(J)=VOUT(J) +c$$$ VECTINI(J)=VOUT(J) +c$$$ ENDDO +c$$$ 22 continue +c$$$ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) +c$$$ IF(VOUT(3).LT.VECT(3)) THEN +c$$$ IFAIL=1 +c$$$c$$$ PRINT *,'=== WARNING ===> tracciamento invertito (UP)' +c$$$c$$$ print*,'charge',charge +c$$$c$$$ print*,'vect',vect +c$$$c$$$ print*,'vout',vout +c$$$c$$$ print*,'step',step +c$$$ RETURN +c$$$ ENDIF +c$$$ Z=VOUT(3) +c$$$ IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200 +c$$$ IF(Z.LT.ZIN(I)-TOLL) GOTO 20 +c$$$ IF(Z.GE.ZIN(I)+TOLL) THEN +c$$$ STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3)) +c$$$ DO J=1,7 +c$$$ VECT(J)=VECTINI(J) +c$$$ ENDDO +c$$$ GOTO 22 +c$$$ ENDIF +c$$$ 200 XOUT(I)=VOUT(1) +c$$$ YOUT(I)=VOUT(2) +c$$$ ZIN(I)=VOUT(3) +c$$$ +c$$$ ENDDO RETURN END -