C===============================================================================
      SUBROUTINE DCALC (Nv,Mv, DXYZ, Dr,Dz, Dl)
C===============================================================================
C      Auxialliary subroutine for use with CTVMFT.
C      Calculates the distance between the two vertices Nv, Mv and its error.
C      Mostly useful only if Nv "points at" Mv.
C      Input parameters:  Nv, Mv  The vertex indices.
C          The rest of the Input data comes from the common block /CTVMFT/
C          which must contain a valid fit
C      Output parameters:
C      DXYZ = (unrotated) x, y, and z components of vertex displacement
C      Dr   = (signed) radial vertex separation (dot product of Dxy * VtxPxy)
C      Dz   = (signed) vertex separation (dot product of Dz * VtxPz)
C      Dl   = uncertainty of the three-vector (Dr,Dphi,Dz) (note Dphi~0)
C-------------------------------------------------------------------------------
        IMPLICIT NONE
#include <Stntuple/fit/ctvmft.inc>
      INTEGER Nv,Mv
      REAL Dxyz(3), Dr,Dz, Dl(3), FL(3), U
      INTEGER I,J,K, IERR, NvF,MvF
      DOUBLE PRECISION   S, ROT(3,3), E(3,3),F(3,3)
      CHARACTER*80 STRING
C-------------------------------------------------------------------------------
C  Setup vertex offset pointers, check for primary
      IF (Nv.LT. 1 .OR. Nv.GT.NVERTX) THEN
        IERR = 99
        call print_v('DCALC: vertex Nv %2i is invalid\n',nv)
        return
      ELSE IF (Mv.LT.0 .OR. Mv.GE.Nv) THEN
        IERR = 98
        call print_v('DCALC: vertex Mv %2i is invalid\n',mv)
        return
      END IF
      DO I=1,9
        ROT(I,1) = 0.0
      END DO
C            Pointers to information in VMAT for the two vertices
      NvF = VOFF(Nv)
      IF (Mv.GT.0) THEN
        MvF = VOFF(Mv)
      ELSE
        MvF =-1
        DO I=1,NVERTX
C the primary participated in this fit
          IF (VtxPnt(I,1).EQ.0) THEN  
            MvF = 0
          END IF
        END DO
      END IF
C covariance matrix on vertex separation
      DO I=1,3                         
      DO J=I,3
        E(I,J) = VMAT(NvF+I,NvF+J)
        IF (MvF.GE.0) THEN
          E(I,J) = E(I,J) + VMAT(MvF+I,MvF+J)
          E(I,J) = E(I,J) - VMAT(NvF+J,MvF+I) * 2.0
        ELSE
          E(I,J) = E(I,J) + EXYZPV(I,J)
        END IF
        E(J,I) = E(I,J)
      END DO
      END DO
C vertex separation vector
      IF (MvF.GE.0) THEN              
        DO I=1,3
        DXYZ(I) = XYZVRT(I,Nv) - XYZVRT(I,Mv)
        END DO
      ELSE
        DO I=1,3
        DXYZ(I) = XYZVRT(I,Nv) - XYZPV0(I)
        END DO
      END IF
      S = SQRT(DXYZ(1)**2 + DXYZ(2)**2)
      ROT(1,1) = DXYZ(1)/S
      ROT(1,2) = DXYZ(2)/S
      ROT(2,1) =-ROT(1,2)
      ROT(2,2) = ROT(1,1)
      ROT(3,3) = 1.0
      do i=1,2
        s = 0.0
        do j=1,2
          s = s + rot(i,j)*DXYZ(J)
        end do
        fl(i) = s
      end do
C pre and post multiply xy part
      DO I=1,2                         
      DO J=1,2
        S = 0.0
        DO K=1,2
          S = S + ROT(I,K) * E(K,J)
        END DO
        F(I,J) = S
      END DO
      END DO
      DO I=1,2
      DO J=1,2
        S = 0.0
        DO K=1,2
          S = S + F(I,K) * ROT(J,K)
        END DO
        E(I,J) = S
      END DO
      END DO
      U  = DXYZ(1)*VtxP4(1,Nv) + DXYZ(2)*VtxP4(2,Nv)
C (signed) xy vertex separation
      Dr = SIGN(SQRT(DXYZ(1)**2+DXYZ(2)**2),U) 
      U  = DXYZ(3)*VtxP4(3,Nv)
C (signed)  z vertex separation
      Dz = SIGN(ABS(DXYZ(3)),U)                
      DO I=1,3
        IF (E(I,I).GT.0.0) THEN
          Dl(I) = SQRT(E(I,I))
        ELSE
          Dl(I) = 0.0
        END IF
      END DO
      RETURN
      END