subroutine kfit_mass(ktrk, update_track, mass, chisq, error)
*
* begin_doc
*
*  Fit a track by constraining it to a fixed mass.
*
*  Note: If the track is not updated here, you can update it later by
*  calling kfit_update_tracks.
*
*  Input Parameters:
*  ktrk          integer variable
*                Track number to fit
*
*  update_track  integer variable
*                 0 ==> Calculate chisquare only (do not update track)
*                 1 ==> Calculate chisquare, update track
*                 2 ==> Calculate chisquare, update track & cov matrix
*
*  mass          double precision variable
*                Mass to constrain to
*
*  Output Parameters:
*  chisq         double precision variable
*                Chisquare of fit
*
*  error         integer variable
*                 0 ==> All OK
*                >0 ==> Error
*
*  Other routines:
*
*  Notes:
*
*  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*   The mass constraint is imposed by the method of Lagrange multipliers
*
*   chisq = (alpha-alpha0)(t) * Valpha0(inv) * (alpha-alpha0)
*             + 2*lambda(t)*(D*dalpha + d)
*
*   where
*    alpha0 are the unconstrained parameters,
*    alpha are the final parameters and
*    D*dalpha + d = 0 is the linearized constraint equation.
*
*  The solution is
*
*       alpha = alpha0 - Valpha0*D(t)*lambda
*      lambda = VD * (D*dalpha0 + d)
*          VD = [D * Valpha0 * D(t)](inv)
*      Valpha = Valpha0 - Valpha0*D(t)*VD*D*Valpha0
*
*   The mass constraint has the form
*
*         E**2 - Px**2 - Py**2 - Pz**2 - Mc**2 = 0
*
*   This equation can be linearized by expanding around an approximate
*   solution x1 = (px1, py1, pz1, E1):
*
*      E1*dE1 - Px1*dpx1 - Py1*dpy1 - Pz1*dpz1 + (M1**2 - Mc**2)/2 = 0
*
*   where M1 is the mass obtained with the vector (px1, py1, pz1, E1).
*   This means that D can be written
*
*       D = (-Px1, -Py1, -Pz1, E1)
*       d = (M1**2 - Mc**2) / 2
*
*   The covariance matrix of the resulting track can be computed as
*
*       Valpha = Valpha0 - Valpha0*D(t)*VD*D*Valpha0
*  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*
*  Author:   Paul Avery      Created:  Fri Aug 29 16:12:52 EDT 1997
*
*  Major revisions:
*
*
* end_doc