C===============================================================================
      SUBROUTINE CTVMCF (PRINT, NV,  VXI)
C===============================================================================
C  CTVMCF calculates "conversion constraint" contributions to the derivative
C         vector and matrix for vertex Nv (08/29/94 JPB and WJA)
C--Input parameters
C        Input data in the include file CTVMFT (COMMON /CTVMFC/)
C  A physically reasonable constraint to impose in fitting a vertex at which
C  a photon conversion is believed to have occurred is that the invariant 
C  mass of the two tracks be as small as possible, i.e. that the opening angle
C  of the two tracks at the vertex be zero.  Where the tracks are copunctual,
C  they should be collinear.  
C  In r-phi, a convenient way to express this constraint is that at the radius
C  of conversion, the sum of the azimuth and the turning angle is equal for
C  both tracks, i.e. Phi1+2*Crv1*S1(r) = Phi2+2*Crv2*S2(r), where S_i(r) is
C  the arc length (projected into r-phi) of track i at radius r, and Phi_i 
C  and Crv_i are the azimuth and half-curvature of track i.  In other words,
C  we equate the r-phi momentum directions of the two tracks at the vertex.
C  In r-z, the constraint is clearly that Ctg1 = Ctg2.
C-------------------------------------------------------------------------------
        IMPLICIT NONE
#include <Stntuple/fit/ctvmft.inc>
      INTEGER PRINT
      INTEGER Nv
      DOUBLE PRECISION VXI(MaxDim)
      INTEGER Trk1,Trk2,NvF,NcF,I,J,IT1F,IT2F
      REAL XV,YV,Crv1,Crv2,Phi1,Phi2,Ctg1,Ctg2,Phs1,Phs2,S1,S2
      REAL DS1DX,DS2DX,DS1DY,DS2DY,DSDCrv1,DSDCrv2,DSDPhi1,DSDPhi2
      REAL CPhi1,CPhi2,SPhi1,SPhi2,CPhs1,CPhs2,SPhs1,SPhs2
      REAL SHOUT,MINUS/-1/,ZERO/0/
      REAL PI
      PARAMETER (PI=3.141592653589793238)
      CHARACTER*80 STRING
C-------------------------------------------------------------------------------
C  Find the first two tracks for this vertex
      Trk1 = 0
      Trk2 = 0
      DO I=1,NTRACK
        IF (TrkVtx(I,Nv)) THEN
          IF (Trk1.EQ.0) THEN
            Trk1 = I
          ELSE
            Trk2 = I
            GO TO 100
          END IF
        END IF
      END DO
 100  CONTINUE
C  Save some possibly useful offsets into VMAT
      IT1F = TOFF(Trk1)
      IT2F = TOFF(Trk2)
      NvF = VOFF(Nv)
      NcF = COFF(Nv)
C  Vertex position
      XV = XYZVRT(1,Nv)
      YV = XYZVRT(2,Nv)
C  Track parameters
      Crv1 = PAR0(1,Trk1)+PARDIF(1,Trk1)
      Crv2 = PAR0(1,Trk2)+PARDIF(1,Trk2)
      Phi1 = PAR0(2,Trk1)+PARDIF(2,Trk1)
      Phi2 = PAR0(2,Trk2)+PARDIF(2,Trk2)
      Ctg1 = PAR0(3,Trk1)+PARDIF(3,Trk1)
      Ctg2 = PAR0(3,Trk2)+PARDIF(3,Trk2)
C  Simple functions of track parameters
      SPhi1 = SIN(Phi1)
      SPhi2 = SIN(Phi2)
      CPhi1 = COS(Phi1)
      CPhi2 = COS(Phi2)
C  Track turning angles (to get momentum direction)
      SPhs1 = 2*Crv1*(XV*CPhi1+YV*SPhi1)
      SPhs2 = 2*Crv2*(XV*CPhi2+YV*SPhi2)
      CPhs1 = SQRT((1+SPhs1)*(1-SPhs1))
      CPhs2 = SQRT((1+SPhs2)*(1-SPhs2))
      Phs1 = ASIN(SPhs1)
      Phs2 = ASIN(SPhs2)
C  Arc length in XY plane
      S1 = 0.5*Phs1/Crv1
      S2 = 0.5*Phs2/Crv2
C  d(arc length)/d(vertex position)
      DS1DX = CPhi1/CPhs1
      DS2DX = CPhi2/CPhs2
      DS1DY = SPhi1/CPhs1
      DS2DY = SPhi2/CPhs2
C  d(arc length)/d(track parameters)
      IF (ABS(CPhs1).GT.0.005) THEN
        DSDCrv1 = 0.5*(SPhs1/CPhs1-Phs1)/(Crv1*Crv1)
      ELSE
        DSDCrv1 = 0.33333333*SPhs1**3*(0.9*SPhs1*SPhs1-1)/(Crv1*Crv1)
      END IF
      IF (ABS(CPhs2).GT.0.005) THEN
        DSDCrv2 = 0.5*(SPhs2/CPhs2-Phs2)/(Crv2*Crv2)
      ELSE
        DSDCrv2 = 0.33333333*SPhs2**3*(0.9*SPhs2*SPhs2-1)/(Crv2*Crv2)
      END IF
      DSDPhi1 = (-XV*SPhi1+YV*CPhi1)/CPhs1
      DSDPhi2 = (-XV*SPhi2+YV*CPhi2)/CPhs2
C  Constraint in r-phi, C_rphi = (Phi1+2*Crv1*S1)-(Phi2+2*Crv2*S2)
      VXI(NcF+1) = (Phi2+Phs2)-(Phi1+Phs1)
      IF (VXI(NcF+1).LT.-PI) VXI(NcF+1) = VXI(NcF+1)+2*PI
      IF (VXI(NcF+1).GT.+PI) VXI(NcF+1) = VXI(NcF+1)-2*PI
C  Optional constraint in r-z, C_rz = Ctg1-Ctg2
      IF (Cvtx(Nv).EQ.2) THEN
        VXI(NcF+2) = Ctg2-Ctg1
      END IF
c-----------------------------------------------------------------------
C  Derivatives of the r-phi constraint C_rphi with respect to the fit 
C  parameters (X,Y,Z,Crv1,Phi1,Ctg1,Crv2,Phi2,Ctg2)
C d/dXv
c-----------------------------------------------------------------------
      VMAT(NVF+1,NCF+1) = 2*(Crv1*DS1DX-Crv2*DS2DX)  
C d/dYv
      VMAT(NVF+2,NCF+1) = 2*(Crv1*DS1DY-Crv2*DS2DY)
C d/dCrv1
      VMAT(IT1F+1,NCF+1) = +(1+2*(S1+Crv1*DSDCrv1))  
C d/dCrv2
      VMAT(IT2F+1,NCF+1) = -(1+2*(S2+Crv2*DSDCrv2))  
C d/dPhi1
      VMAT(IT1F+2,NCF+1) = +(1+2*Crv1*DSDPhi1)       
C d/dPhi2
      VMAT(IT2F+2,NCF+1) = -(1+2*Crv2*DSDPhi2)
c-----------------------------------------------------------------------  
C  Optional derivatives of the r-z constraint C_rz with respect to
C  the fit parameters
c-----------------------------------------------------------------------
      IF (Cvtx(Nv).EQ.2) THEN
C d/dCtg1
        VMAT(IT1F+3,NcF+2) = +1                      
C d/dCtg2
        VMAT(IT2F+3,NcF+2) = -1                      
      END IF
C  Symmetrize
      DO I=1,MATDIM-1
        DO J=I+1,MATDIM
          VMAT(J,I) = VMAT(I,J)
        END DO
      END DO
C  All done
      RETURN
      END