C======================================================================= =======
      SUBROUTINE MCALC (NTRKQ,MTRACK, FMASS,DMASS, P4)
C======================================================================= =======
C
C=======Description
C
C  Input parameters:
C  Requires the common in the CTVMFT include file to have relevant information
C  stored from a CTVMFT entrance.
C  NTRKQ = Number of tracks comprising the invariant mass.
C  MTRACK(i) selects the track specified in LIST(i) in the CTVMFT call
C  Output parameters:
C  FMASS = invariant mass of track combination
C  DMASS = computed error on the mass (will be returned as negative if
C           DMASS**2 was computed to be negative. This should not happen
C           if the matrix VMAT is valid)
C=======Author
C  John Marriner
C  6 April 1992
C----------------------------------------------------------------------- -------
        IMPLICIT NONE
#include <Stntuple/fit/ctvmft.inc>
      INTEGER NTRKQ
      INTEGER MTRACK(MAXTRK)
      REAL FMASS, DMASS
      INTEGER I,J, IERR, Nv,NvF, Nt,NtF, NDIM, Nvx,NVRTX(MaxVtx)
      REAL DMASSQ
      DOUBLE PRECISION SUM,QMASS, P4(4), DM(MAXDIM)
      CHARACTER*80 STRING
      CHARACTER*6  NAME
      DATA    NAME /'MCALC '/
C===============================================================================
      DO I=1,4
        P4(I) = 0.0D0
      END DO
      DO I=1,MAXDIM
        DM(I) = 0.0D0
      END DO
      FMASS = 0.0
      DMASS = 0.0
C complain
      IF (NTRKQ.GT.MAXTRK) THEN         
        IERR = 91
        call print_v('MCALC: too many tracks %2i\n',ntrkq)
        FMASS =-1.0
        GO TO 100
      END IF
      DO J=1,NTRKQ
        Nt = MTRACK(J) 
C energy
        TrkP4(NT,4) = SQRT(TrkP4(NT,6)**2+TMASS(NT)**2) 
        DO I=1,4
          P4(I) = P4(I) + TrkP4(Nt,I)
        END DO
      END DO
      NDIM = 0
      Nvx  = 0
      DO Nv=1,MaxVtx
        NVRTX(Nv) = 0
      END DO
      DMASSQ = 0
C Loop over tracks in the fit, find dM
      DO 10 J=1,NTRKQ                   
        Nt = MTRACK(J) 
        NtF = TOFF(Nt)
C find the vertex for this track
        Nv = 1                          
        DMASSQ = DMASSQ + TMASS(Nt)
        DO WHILE (.NOT.TrkVtx(Nt,Nv))
          Nv = Nv+1
        END DO
C count number of vertices involved
        IF (NVRTX(Nv).EQ.0) THEN        
          Nvx = Nvx+1
          NVRTX(Nv) = Nv
        END IF
        NvF = VOFF(Nv)
C maximum matrix dimension in VMAT
        NDIM = MAX0(NDIM,NtF+3)         
C first,dM**2/dXsv and  dM**2/dYsv, treat the displacement to this vertex
        DM(NvF+1) = DM(NvF+1)-2.0*(P4(1)*DDA(Nt,3)+P4(2)*DDA(Nt,7)) 
        DM(NvF+2) = DM(NvF+2)-2.0*(P4(1)*DDA(Nt,4)+P4(2)*DDA(Nt,8))
C then, dM**2/dc, dM**2/dphi, dM**2/dcot(theta), helix parameters for this track
        SUM = -P4(4)*TrkP4(Nt,6)**2/TrkP4(Nt,4) + P4(3)*TrkP4(Nt,3)
        SUM = SUM/PAR(1,Nt) - (P4(1)*DDA(Nt,1) + P4(2)*DDA(Nt,5))
        DM(NtF+1) = 2.0*SUM
        DM(NtF+2) =-2.0*(P4(1)*DDA(Nt,2) + P4(2)*DDA(Nt,6))
        SUM = (P4(4)*TrkP4(Nt,3)/TrkP4(Nt,4) - P4(3)) * TrkP4(Nt,5)
        DM(NtF+3) = 2.0*SUM
   10 CONTINUE
C               mass, mass**2
      QMASS = DSQRT(P4(1)**2+P4(2)**2+P4(3)**2)
      IF (P4(4).LE.QMASS) THEN
        IF (Nvx.EQ.1 .AND. Cvtx(Nv).GT.0) THEN
          QMASS = DMASSQ
        ELSE
C ? crunch
          QMASS = DSQRT((P4(4)+QMASS)*(P4(4)-QMASS)) 
        END IF
      ELSE
        QMASS = DSQRT((P4(4)+QMASS)*(P4(4)-QMASS))
      END IF
      FMASS = QMASS
C Set error in mass squared = 0 for summing
      DMASSQ = 0.0              
      DO I=1,NDIM
      DO J=1,NDIM
        DMASSQ = DMASSQ + DM(I)*VMAT(J,I)*DM(J)
      END DO
      END DO
C check for, and flag, negative error**2
      IF (DMASSQ.GE.0.) THEN    
         DMASS = SQRT(DMASSQ)
      ELSE
        DMASS = -SQRT(-DMASSQ)
      ENDIF
      IF (QMASS.EQ.0.0) THEN
        FMASS = -0.001
      ELSE
        DMASS = DMASS/(2.*QMASS)
      END IF
  100 CONTINUE
      RETURN
      END