subroutine kfit_back_to_back_cm(ktrk1, ktrk2, update_track,
     *                                velocity, chisq, error)
*
* begin_doc
*
*  Fit two tracks so that they are back to back, i.e., part of a single
*  helix, taking into account the fact that they have been boosted from
*  the center of mass frame. This routine is needed when beams do not
*  collide head on.
*
*  The boost is assumed to be small and is evaluated to first order in
*  the velocity.
*
*  Note: If the track is not updated here, you can update it later by
*  calling kfit_update_tracks.
*
*  Input Parameters:
*  ktrk1         integer variable
*                Track 1 to fit
*
*  ktrk2         integer variable
*                Track 2 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
*
*  Output Parameters:
*  velocity(3)   double precision array
*                Velocity of boost
*
*  chisq         double precision variable
*                Chisquare of fit
*
*  error         integer variable
*                 0 ==> All OK
*                >0 ==> Error
*
*  Other routines:
*
*  Notes:
*
*  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*   The momemtum components in the center of mass frame are:
*
*       Px* = Px - gam * vx [E - gam/(gam+1) * v o P]
*       Py* = Py - gam * vy [E - gam/(gam+1) * v o P]
*       Pz* = Pz - gam * vz [E - gam/(gam+1) * v o P]
*        E* = gam * [E - v o P]
*
*   where v o P is the dot product of the velocity boost and the
*   momentum. To make the constraint equations tractable, I
*   assume here that the boost velocity is very small and work
*   only to first order in the velocity.
*
*   The back to back constraint is imposed by the method of Lagrange
*   multipliers. There are 5 constraints, of which three are best
*   expressed in terms of the track parameters at the point of
*   closest approach to the origin:
*
*      1. pt1*   =  pt2*
*      2. px1'*py2' = px2'*py1'   (opposite direction)
*      3. d0_1   = -d0_2
*      4. pz1*   = -pz2*
*      5. z0_1   =  z0_2
*
*   chisq = (alpha-alpha0)(t) * Valpha0(inv) * (alpha-alpha0)
*             + 2*lambda(t)*(D*dalpha + d)
*
*   where
*    alpha0 are the unconstrained parameters, (2*7 of them)
*    alpha are the final parameters (2*7) and
*    D*dalpha + d = 0 represents the linearized constraint equations (5).
*
*  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 covariance matrix of the resulting tracks can be computed as
*
*       Valpha = Valpha0 - Valpha0*D(t)*VD*D*Valpha0
*
*   The constraints are somehat complicated to show here in differential
*   form.
*  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*
*  Author:   Paul Avery      Created:  Fri Dec  4 22:12:52 EDT 1997
*
*  Major revisions:
*
*
* end_doc