--- DarthVader/TrackerLevel2/src/F77/track.f 2006/05/19 13:15:56 1.1 +++ DarthVader/TrackerLevel2/src/F77/track.f 2006/06/01 09:03:09 1.2 @@ -1,13 +1,13 @@ ************************************************************* * -* Routine to compute the NPOINT track intersection points +* Routines to compute the NPOINT track intersection points * with planes of z-coordinate given by ZIN * given the track parameters * * The routine is based on GRKUTA, which computes the * trajectory of a charged particle in a magnetic field -* by solving the equatoins of motion with Runge-Kuta method +* by solving the equations of motion with Runge-Kuta method * * Variables that have to be assigned when the subroutine * is called: @@ -19,12 +19,178 @@ * The routine works properly only if the * planes are numbered in descending order * +* Routines: * -* modified on 05/2006 to give as output also the angles and -* theh track-lengths +* DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL) +* - old routine +* +* DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,AL_P,IFAIL) +* - as the old one, but: +* -- the projected angles are given as output +* -- the track lengths are given as output: +* ---- for planes above the reference plane, the length until +* the lower closest plane (reference plane included) +* ---- for planes below the reference plane, the length until +* the higher closest plane (reference plane included) * ************************************************************** + SUBROUTINE + $ DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL) + + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + DIMENSION VECT(8),VECTINI(8),VOUT(8) + 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) + VOUT(8)= 0 +c$$$ print*,'DOWN ',i,' - Track from ', +c$$$ $ vout(3),' to ',zin(i),' step ',step + 10 DO J=1,8 + VECT(J)=VOUT(J) + VECTINI(J)=VOUT(J) + ENDDO + 11 continue + CALL GRKUTA2(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,8 + VECT(J)=VECTINI(J) + ENDDO + GOTO 11 + ENDIF + 100 XOUT(I)=VOUT(1) + YOUT(I)=VOUT(2) + ZIN(I)=VOUT(3) + IF(VOUT(3).ne.0)THEN + 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) +c print*,'D ',VECT(3),VOUT(3),VOUT(8),STEP + 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 + VOUT(8)=0 +c$$$ print*,'UP ',i,' - Track from ', +c$$$ $ vout(3),' to ',zin(i),' step ',step + 20 DO J=1,8 + VECT(J)=VOUT(J) + VECTINI(J)=VOUT(J) + ENDDO + 22 continue + CALL GRKUTA2(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,8 + VECT(J)=VECTINI(J) + ENDDO + GOTO 22 + ENDIF + 200 XOUT(I)=VOUT(1) + YOUT(I)=VOUT(2) + ZIN(I)=VOUT(3) + IF(VOUT(3).ne.0)THEN + 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) + + ENDDO + + RETURN + END + +************************************************************************ +* +* +* +* +************************************************************************ SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) @@ -37,7 +203,6 @@ PARAMETER (NPOINT_MAX=100) DIMENSION ZIN(NPOINT_MAX) DIMENSION XOUT(NPOINT_MAX),YOUT(NPOINT_MAX) - DIMENSION THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX) DIMENSION AL_P(5) * ----------------------------------------------- DATA ZINI/23.5/ !z coordinate of the reference plane @@ -81,12 +246,12 @@ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) IF(VOUT(3).GT.VECT(3)) THEN IFAIL=1 - PRINT *,'=== WARNING ===> tracciamento invertito (DOWN)' - print*,'charge',charge - print*,'vect',vect - print*,'vout',vout - print*,'step',step - RETURN +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 @@ -136,11 +301,11 @@ CALL GRKUTA(CHARGE,STEP,VECT,VOUT) IF(VOUT(3).LT.VECT(3)) THEN IFAIL=1 - PRINT *,'=== WARNING ===> tracciamento invertito (UP)' - print*,'charge',charge - print*,'vect',vect - print*,'vout',vout - print*,'step',step +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)