*************************************************************
*
* 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,TLOUT,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)
*
* March 2008 --> optimized by Paolo
**************************************************************
SUBROUTINE
$ DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,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 THXOUT(NPOINT_MAX),THYOUT(NPOINT_MAX)
DIMENSION TLOUT(NPOINT_MAX)
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
* ==================================================================
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
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 GRKUTA(CHARGE,STEP,VECT,VOUT)
L = L + STEP
IF(VOUT(3).GT.VECT(3)) THEN
IFAIL=1
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,7
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) = L
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
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 GRKUTA(CHARGE,STEP,VECT,VOUT)
L = L + STEP
IF(VOUT(3).LT.VECT(3)) THEN
IFAIL=1
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,7
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) = L
ENDDO
RETURN
END
************************************************************************
*
*
*
*
************************************************************************
SUBROUTINE DOTRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
* -----------------------------------------------
* 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)
* -----------------------------------------------
CALL DOTRACK2(NPOINT,ZIN,XOUT,YOUT,THXOUT,THYOUT,TLOUT,AL_P,IFAIL)
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