subroutine kvtx_known(ntrk, list_trk, update_tracks, luse_vtx,
     *                      vtx, Vvtx, vtx0, Vvtx0, chisq, error)
*
* begin_doc
*
*  Fit several tracks to a vertex which has a prior position and
*  covariance matrix.
*
*  Input Parameters:
*  ntrk            integer variable
*                  # of tracks to fit (length of list_trk)
*
*  list_trk(*)     integer  array
*                  List of tracks in track list to fit (length = ntrk)
*
*  update_tracks   integer variable
*                  0 ==> don't update input track parameters
*                  1 ==> update input track parameters
*                  2 ==> update input track parameters & covariance matrices
*                        and replace the original unconstrainted parameters
*                        by the new ones. ***** This is a dangerous option.
*
*  luse_vtx        logical variable
*                  TRUE  ==> use the value of vtx(1-3) as the starting value
*                  FALSE ==> use the value of vtx0(1-3) as the starting value
*
*  vtx(3)          double precision array
*                  Initial values of (x,y,z)
*
*  vtx0(3)         double precision array
*                  Positions of known vertex
*
*  Vvtx0(3)        double precision array
*                  3x3 covariance matrix of vtx0
*
*  Output Parameters:
*  vtx(3)          double precision array
*                  Updated values of (x,y,z)
*
*  Vvtx(3,3)       double precision array
*                  3x3 covariance matrix of vtx(1-3)
*
*  chisq           double precision variable
*                  Chisquare of fit
*
*  error           integer variable
*                   0  if all OK
*                  >0  if error
*
*  Other routines:
*
*  Notes:
*
*  Author:   Paul R Avery      Created:  Mon Jul 28 07:05:06 EDT 1997
*
*  Major revisions:
*
*  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*
*     The constraint equations for a solenoidal field are:
*
*       0 = dely*Px - delx*Py - 0.5*a*(delx**2 + dely**2)
*       0 = delz - (Pz/a) * arcsin[a*(delx*Px + dely*Py)/Pt**2]
*
*     where a = -0.299792458*B*Q
*           B = magnetic field in Tesla
*           Q = charge of particle
*           delx = x - x0
*           dely = y - y0
*           delz = z - z0
*
*     We expand to first order to evaluate (see Fitting Theory 6)
*
* 1. 0 = (dy-dy0)*Px + dely*dPx - (dx-dx0)*Py - delx*dPy  -
*        a*[delx*(dx-dx0) + dely*(dy-dy0)] + d1
*
* 2. 0 =  dz - dz0 - dPz/a * arcsin[a*S] - Pz/sqrt(1-a**2*S**2)*FF + d2
*
*  S  = (delx*Px + dely*Py)/Pt**2
*  FF = -2*(Px*dPx + Py*dPy)(S/Pt**2 +
*       [(dx-dx0)*Px+delx*dPx+(dy-dy0)*Py+dely*dPy]/Pt**2
*
*     where d1 and d2 are the values of the constraint equations before
*     the fit.
*
*     The fitting algorithm is as follows. Let the track parameters be
*     called alpha and the vertex parameters be called z. We assume that
*     the original alpha parameters and their covariance matrix are
*     available from an earlier fit. The following chisquare equation
*     must be minimized:
*
*        chisq = (alpha-alpha0)(t) * Valpha0(inv) * (alpha-alpha0)
*              + (z - z0)(t) * Vz0(inv) * (z - z0)
*              + 2*lambda(t) * (D*alpha + E*z + d1),
*
*     which has contributions from the tracks, the vertices and the
*     constraints. The variables are defined as follows:
*
*     lambda = Lagrange multipliers
*     alpha0 = initial track parameters
*     alpha  = final track parameters
*    Valpha0 = initial track covariance matrix.
*         z0 = initial vertex parameters
*          z = final vertex parameters
*        Vz0 = initial vertex covariance matrix.
*
*      The solution that minimizes this chisquare is
*
*             lambda0 = VD * (D*alpha0 + d),
*             z = z0 - Vz0*E(t)*lambda
*               = Vz*Vz0(inv)*z0 -  Vz*E(t)*lambda0
*             lambda = VD' * (D*alpha0 + E*z0 + d)
*             alpha = alpha0 - Valpha0*D(t)*lambda,
*                VD = [D * Valpha0 * D(t)](inv)
*                VD'= VD - VD*E*Vz*E(t)*VD
*
*     The chisquare can be written
*
*           chisq = lambda(t) * VD'(inv) * lambda
*
*     The covariance matrices are given by
*
*             Vz = [Vz0(inv) + E(t)*VD*E](inv)
*             Vlambda = VD'
*             Valpha = Valpha0 - Valpha0*D(t)*Vlambda*D*Valpha0.
*                    = Valpha0 - Valpha0*D(t)*VD*D*Valpha0
*                              + Valpha0*D(t)*VD*E*Vz*E(t)*VD*D*Valpha0
*
*     The above solution is completely general for any set of parameters
*     appearing in a constraint equation. However, in vertex fitting the
*     elements of the matrices D, E and d do not "mix" tracks so we only
*     have to invert small 2x2 matrices local to each track and a final
*     3 x 3 matrix.
*  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
* end_doc