************************************************************* * * 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 equations of motion with Runge-Kuta method * * Variables that have to be assigned when the subroutine * is called: * * ZIN(1-NPOINT) ----> z coordinates of the planes * AL_P(1-5) ----> track-parameter vector * * NB !!! * The routine works properly only if the * planes are numbered in descending order * * Routines: * * 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) 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 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) ENDDO RETURN END