CDF KITS Project:/ VertexAlg/ src/ VertexFit.cc		[ source navigation ] [ diff markup ] [ identifier search ] [ freetext search ] [ file search ]
	Architecture: [ i386 ] Version: [ 4.10.4 ] [ 4.10.5 ] [ 4.8.4 ] [ 4.8.4l3s ] [ 4.8.5 ] [ 4.9.0 ] [ 4.9.1 ] [ 4.9.1.hpt3 ] [ 4.9.1hpt3 ] [ 4.9.1top1 ] [ 5.0.0 ] [ 5.1.0 ] [ 5.1.0beamonly ] [ 5.1.1 ] [ 5.2.0 ] [ 5.3.0 ] [ 5.3.1 ] [ 5.3.1dsp ] [ 5.3.3 ] [ 5.3.3_nt ] [ 5.3.4 ] [ 6.1.1 ] [ 6.1.1b ] [ 6.1.2 ] [ 6.1.3 ] [ 6.1.4 ] [ 6.1.4int3 ] [ 6.1.4mc ] [ 6.1.4mc_a ] [ 6.1.6 ] [ development ]

  /************************************************************************

  * VertexFit class implementation file                                   *

  *                                                                       *

  *                                                                       *

  * Author: Craig Blocker, CDF Group                                      *

  *         Phone 781-736-2920                                            *

  *                                                                       *

  * Description: Routines to implement VertexFit class                    *

  *                                                                       *

  * See VertexFit.hh in the VertexAlg/ area and                           *

  * the file VertexFit.txt in the /doc area of VertexAlg for              *

  * documentation.                                                        *

  *                                                                       *

  * Revision History:                                                     *

  *    May 27, 1999   Craig Blocker   Initial creation                    * 

  *    June 22, 2001  Craig Blocker   Fix a few bugs in checking on       *

  *                                     vertices (such as pointing to     *

  *                                     primary vertex wasn't allowed)    *

  *                                                                       *

  *    Sept 07, 2002 Joao Guimaraes   Added accessors for ijk errors and  *

  *                                   track-id of tracks that produce     *

  *                                   some of those errors                *

  *                                   These methods are:                  *

  *                                      getIJKErr(int&,int&,int&) const  *

  *                                      getIJKErr0() const               *

  *                                      getIJKErr1() const               *

  *                                      getIJKErr2() const               *

  *                                      getErrTrackId() const            *

  *                                                                       *

  *    Feb 10, 2002 Joao Guimaraes    Added two methods to allow for a    *

  *                                   beam constrained fit with CTVMFT.   *

  *                                   This only makes sense when fitting  *

  *                                   a primary vertex and should not be  *

  *                                   used for fitting secondary vertices *

  *    Oct 15 2003 Sal Rappoccio (rappocc@fnal.gov)                       *

  *                  Added accessor to add tracks to ctvmft via a track   *

  *                  id and the track parameters and covariance           *

  *                                                                       *

  *  December 2 2003 Sal Rappoccio (rappocc@fnal.gov)                     *

  *                  Added an interface to use extrapolated track errors  *

  *                  in a "transparent" way. The new interface is used as:*

  *                        VertexFit fit;                                 *

  *                        fit.init();                                    *

  *                        fit.setPrimaryVertex ( p1 ); // Note the order *

  *                                                     // matters        *

  *                        fit.setReferencePoint( p2 );                   *

  *                        fit.addTrack( trk1, M_PION, 1 );               *

  *                        fit.addTrack( trk2, M_PION, 1 );               *

  *                        fit.fit();                                     *

  *                        Hep3Vector & pos = fit.getVertex( 1 );         *

  *                  "pos" will now hold the new vertex in CDF coordinates*

  *                  constructed from tracks with errors referenced at the*

  *                  new reference point "p2".                            *

  *                  This track extrapolation code is OFF by default. It  *

  *                  is turned ON by:   fit.setReferencePoint( p2 );      * 

  *                  AddTrack was modified to use this feature.           *

  *  Dec 2 2003 Joao Guimaraes                                            *

  *                  Updated VertexFit::beamlineConstraint to use new     *

  *                  reference point (see entry above).                   *

  *                                                                       *

  *************************************************************************/
 
 
 #ifndef VERTEXFIT_CC
 #define VERTEXFIT_CC
 
 //-----------------------

 // This Class's Header --

 //-----------------------

 #include "VertexAlg/VertexFit.hh"
 
 extern "C" void ctvmft_(int&,int&,int&);
 extern "C" float prob_(float&,int&);
 extern "C" bool mcalc_(int&,int*,float&,float&,double*);
 extern "C" void dcalc_(int&,int&,float*,float&,float&,float*);
 
 // C++ Headers --

 //---------------

 #if defined(DEFECT_OLD_IOSTREAM_HEADERS)
 #include <iostream>
 #else // Standard C++

 #include <iostream>
 #endif
 
 #if defined(DEFECT_NO_STDLIB_NAMESPACES)
 // Do nothing

 #else // Standard C++

 using std::cerr ;
 using std::cout ;
 using std::endl ;
 #endif
 
 #include <algorithm>
 #include <math.h>
 
 //-------------------------------

 // Collaborating Class Headers --

 //-------------------------------

 #include "TrackingObjects/Tracks/CdfTrack.hh"
 #include "CLHEP/Matrix/Matrix.h"
 #include "TrackingUserHL/Utility/TrackExtrapolator.hh"
 
 //for floating point exception handling Joe and Farrukh

 #include <csignal>
 #include <csetjmp>
 
 /******************************************************************************

  * Default constructor                                                        *

  *                                                                            *

  *                                                                            *

  *****************************************************************************/
 
 //the following will have scope throughout.

 
  jmp_buf env;
  extern "C" void VertexFitSetStatus(int i) { 
    std::cout << "Warning, you are handling a severe error in VertexFit" << std::endl;
    longjmp(env,-66);
  }
 
 VertexFit::VertexFit() :
   _currentAllocatedVertexNumber(0),   // facilitates CandNode recursion

   _referencePoint( 0,0,0 ), // Set Reference Point to (0,0,0) initially

   _primaryVertex( 0,0,0 ), // Set primary vertex to (0,0,0) initially

   _cdfPrimaryVertex( 0,0,0 ) // Set p.v. in CDF coords to (0,0,0) initially

 {
 // Set name and email of VertexFit expert.

    _expert="Craig Blocker (blocker@fnal.gov)";
 // First get pointers to various FORTAN common blocks.

    _ctvmq_com = (CTVMQ*) ctvmq_address_();
    _ctvmfr_com = (CTVMFR*) ctvmfr_address_();
    _fiddle_com = (FIDDLE*) fiddle_address_();
    _trkprm_com = (TRKPRM*) trkprm_address_();
 // Initialize various arrays

    init();
 // Don't bomb program on error.

    _fiddle.excuse=1;
 // Don't extrapolate track errors by default

    _extrapolateTrackErrors = false;
 }
 
 /******************************************************************************

  * Destructor                                                                 *

  *        Use default, i.e. do nothing                                        *

  *****************************************************************************/
 
 VertexFit::~VertexFit(){ }
 
 
 
 // Particle masses

 const float VertexFit::EL_MASS = 0.000510998902;
 const float VertexFit::MU_MASS = 0.105658357;
 const float VertexFit::PI_MASS = 0.13957018;
 const float VertexFit::K_MASS = 0.493677;
 const float VertexFit::PROTON_MASS = 0.938272;
 
 /*=======================================================================

                          VertexFit::init

   =======================================================================*/
 
 void VertexFit::init()
 {
   double bmag = 1.4116; // Magnetic field in Tesla

    init(bmag);
 }
 
 void VertexFit::init(double bmag)
 {
 // Now initialize CTVMQ common.  Run and trigger numbers are

 // dummies - they are not used.

    _ctvmq.runnum=1;
    _ctvmq.trgnum=100;
 // Eventually, we have to get the magnetic field from the right place.

 // Origignal with      bmag = 14.116:

 // _ctvmq.pscale=0.000149896*bmag;

 // New bfield in Tesla bmag = 1.4116 [T]:

    _ctvmq.pscale=0.00149896*bmag;
 // Set the default maximum chi-square per dof.

    _fiddle.chisqmax = 225.;
 // We also need to get the primary vertex from the right place,

 // but for now we put in (0,0,0).

    setPrimaryVertex(0.,0.,0.);
    float xverr[3][3];
    for(int j = 0; j<3; j++) {
       for(int k = 0; k<3; k++) {
          xverr[j][k] = 0.;
       }
    }
    xverr[0][0]=.005;
    xverr[1][1]=.005;
    xverr[2][2]=1.;
    setPrimaryVertexError(xverr);
 // Zero number of tracks, vertices and mass constraints.

    _ctvmq.ntrack=0;
    _ctvmq.nvertx=0;
    _ctvmq.nmassc=0;
 // Zero track list and arrays containing the vertex and

 // mass constraint configuration.

    for(int j=0; j<_maxtrk; ++j) {
       _ctvmq.list[j]=0;
       for(int jv=0; jv<_maxvtx; ++jv) {
          _ctvmq.trkvtx[jv][j]=false;
       }
       for(int jmc=0; jmc<_maxmcn; ++jmc) {
          _ctvmq.trkmcn[jmc][j]=false;
       }
    }
 // Initialize the conversion and vertex pointing arrays.

    for(int jv=0; jv<_maxvtx; ++jv) {
       _ctvmq.cvtx[jv]=0;
       _ctvmq.vtxpnt[0][jv]=-1;
       _ctvmq.vtxpnt[1][jv]=0;
    }
    _ctvmq.drmax=2.;
    _ctvmq.dzmax=20.;
    _ctvmq.rvmax=70.;
    _ctvmq.trnmax=0.5;
    _ctvmq.dsmin=-2.;   
 // Set stat to -999  and chisqq to -1 to symbolize that no fit has yet been done.

    stat=-999;
    _ctvmq.chisqr[0]=-1.;
 
    _primaryVertex = Hep3Vector(0,0,0);
    _cdfPrimaryVertex = Hep3Vector(0,0,0);
    _referencePoint = Hep3Vector(0,0,0);
 }
 
 /*=======================================================================

                          VertexFit::addTrack

   =======================================================================*/
 
 bool VertexFit::addTrack(const CdfTrack* trk, float mass,
                          vertexNumber jv)
 {
 // Check that this vertex number is within the allowed range.

    if(jv<VERTEX_1 || jv>_maxvtx) return false;
    _ctvmq.nvertx = jv>_ctvmq.nvertx ? jv : _ctvmq.nvertx;
 
 // Check that we have not exceeded the maximum number of tracks.

    if(_ctvmq.ntrack>=_maxtrk) return false;
 
 
 // Sal Rappoccio: Add track extrapolation, if the user set it

 
    HepVector w = trk->par();
    HepSymMatrix m = trk->getHelixFit().getErrorMatrix();
    if ( _extrapolateTrackErrors ) {
      // This will move the reference point

      moveReferencePoint( w, m );
    }
 
 // Add this track.

    _ctvmq.list[_ctvmq.ntrack]=trk->id().value();
    _ctvmq.tkbank[_ctvmq.ntrack][0]='Q';
    _ctvmq.tkbank[_ctvmq.ntrack][1]='T';
    _ctvmq.tkbank[_ctvmq.ntrack][2]='R';
    _ctvmq.tkbank[_ctvmq.ntrack][3]='K';
    _ctvmq.tmass[_ctvmq.ntrack]=mass;
    _ctvmq.trkvtx[jv-1][_ctvmq.ntrack]=true;
 
 // Put this track's helix parameters and error matrix into

 // a fortran common block so that they can  be accessed

 // by gettrk.  This is a dummy for now.

    _trkprm.trhelix[_ctvmq.ntrack][0]=w[0];
    _trkprm.trhelix[_ctvmq.ntrack][1]=w[1];
    _trkprm.trhelix[_ctvmq.ntrack][2]=w[2];
    _trkprm.trhelix[_ctvmq.ntrack][3]=w[3];
    _trkprm.trhelix[_ctvmq.ntrack][4]=w[4];
 
    for(int j=0; j<5; ++j) {
       for(int k=0; k<5; ++k) {
          _trkprm.trem[_ctvmq.ntrack][j][k]=m[j][k];
       }
    }
    _ctvmq.ntrack++;
 
    return true;
 }
 
 // Sal Rappoccio: 11/04/03: Added this method to add track using

 // their track id and the explicit helix parameters.

 // This is useful when for adding tracks with helical parameters relative 

 // to a point other than (0,0)

 
 /*=======================================================================

                          VertexFit::addTrack

   =======================================================================*/
 
 bool VertexFit::addTrack( const HepVector & v,
                           const HepSymMatrix & cov,
                           int trackId,
                           float mass,
                           vertexNumber jv)
 {
 // Check that this vertex number is within the allowed range.

    if(jv<VERTEX_1 || jv>_maxvtx) return false;
    _ctvmq.nvertx = jv>_ctvmq.nvertx ? jv : _ctvmq.nvertx;
 
 // Sal Rappoccio: Add track extrapolation, if the user set it

 
    HepVector w = v;
    HepSymMatrix m = cov;
    if ( _extrapolateTrackErrors ) {
      // This will move the reference point

      moveReferencePoint( w, m );
    }
 
 // Check that we have not exceeded the maximum number of tracks.

    if(_ctvmq.ntrack>=_maxtrk) return false;
 
 // Add this track.

    _ctvmq.list[_ctvmq.ntrack]=trackId;
    _ctvmq.tkbank[_ctvmq.ntrack][0]='Q';
    _ctvmq.tkbank[_ctvmq.ntrack][1]='T';
    _ctvmq.tkbank[_ctvmq.ntrack][2]='R';
    _ctvmq.tkbank[_ctvmq.ntrack][3]='K';
    _ctvmq.tmass[_ctvmq.ntrack]=mass;
    _ctvmq.trkvtx[jv-1][_ctvmq.ntrack]=true;
 
 // Put this track's helix parameters and error matrix into

 // a fortran common block so that they can  be accessed

 // by gettrk.  This is a dummy for now.

    _trkprm.trhelix[_ctvmq.ntrack][0]=w[0];
    _trkprm.trhelix[_ctvmq.ntrack][1]=w[1];
    _trkprm.trhelix[_ctvmq.ntrack][2]=w[2];
    _trkprm.trhelix[_ctvmq.ntrack][3]=w[3];
    _trkprm.trhelix[_ctvmq.ntrack][4]=w[4];
 
    for(int j=0; j<5; ++j) {
       for(int k=0; k<5; ++k) {
          _trkprm.trem[_ctvmq.ntrack][j][k]=m[j][k];
 
       }
    }
    _ctvmq.ntrack++;
 
    return true;
 }
 
 /*=======================================================================

                          VertexFit::vertexPoint_2d

   =======================================================================*/
   
 bool VertexFit::vertexPoint_2d(vertexNumber jv1, vertexNumber jv2)
 {
 // Check that these vertex numbers are within allowed range and

 // that the vertices are unique.

    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
    if(jv1 <= jv2) return false;
 
 // Setup vertex pointing.

    _ctvmq.vtxpnt[0][jv1-1]=jv2;
    _ctvmq.vtxpnt[1][jv1-1]=1;           // 2d pointing.

 
    return true;
 }
 
 /*=======================================================================

                          VertexFit::vertexPoint_3d

   =======================================================================*/
 
 bool VertexFit::vertexPoint_3d(vertexNumber jv1, vertexNumber jv2)
 {
 // Check that these vertex numbers are within allowed range and

 // that the vertices are distinct.

    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
    if(jv1 <= jv2) return false;
 // Setup vertex pointing.

    _ctvmq.vtxpnt[0][jv1-1]=jv2;
    _ctvmq.vtxpnt[1][jv1-1]=2;           // 3d pointing.

 
    return true;
 }
 
 /*=======================================================================

                          VertexFit::vertexPoint_1track

   =======================================================================*/
 
 bool VertexFit::vertexPoint_1track(vertexNumber jv1, vertexNumber jv2)
 {
 // Check that these vertex numbers are within allowed range and are distinct.

    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
    if(jv1 <= jv2) return false;
 
 // Setup vertex pointing.

    _ctvmq.vtxpnt[0][jv1-1]=jv2;
    _ctvmq.vtxpnt[1][jv1-1]=3;           // Point to 1 track vertex.

 
    return true;
 }
 /*=======================================================================

                          VertexFit::vertexPoint_0track

   =======================================================================*/
 
 bool VertexFit::vertexPoint_0track(vertexNumber jv1, vertexNumber jv2)
 {
 
   // jv2 is the zero track vertex. jv1 is the multi track vertex which points to jv2

 
   // Note: You must call this routine at least twice in order for ctvmft_zerotrackvtx to work

   // ie there must be at least 2 vertices pointed at a zero track vertex. ctvmft

   // for this, but the error message may not make it to your log file (look at the local

   // variables in the stack frame especially IJKERR(2). The significance of this

   // error code is documented at the top of whatever routine chucked you out (ctvmf00 in

   // this case)

 
   // see ctvmft.f source file for discussion. See especially comments at the top

   // of subroutines: ctvmft and ctvmfa

 
 // Check that these vertex numbers are within allowed range and are distinct.

    if(jv1>_maxvtx || jv1<VERTEX_1) return false;
    if(jv2>_maxvtx || jv2<PRIMARY_VERTEX) return false;
    if(jv1 <= jv2) return false;
 
 // Setup vertex pointing.

    _ctvmq.vtxpnt[0][jv1-1]=jv2;
    _ctvmq.vtxpnt[1][jv1-1]=4;           // Point to 0 track vertex.

 
    return true;
 }//vertexPoint_0track

 
 /*=======================================================================

                          VertexFit::conversion_2d

   =======================================================================*/
 
 bool VertexFit::conversion_2d(vertexNumber jv)
 {
    if(jv<VERTEX_1 || jv>_ctvmq.nvertx) {
       return false;
    }
    _ctvmq.cvtx[jv-1]=1;
    return true;
 }
 
 /*=======================================================================

                          VertexFit::conversion_3d

   =======================================================================*/
 
 bool VertexFit::conversion_3d(vertexNumber jv)
 {
    if(jv<VERTEX_1 || jv>_ctvmq.nvertx) {
       return false;
    }
    _ctvmq.cvtx[jv-1]=2;
    return true;
 }
 
 /*=======================================================================

                          VertexFit::massConstrain

   =======================================================================*/
 
 bool VertexFit::massConstrain(const CdfTrack* trk1,
                               const CdfTrack* trk2, float mass)
 {
    int ntrk=2;
    const CdfTrack* tracks[2];
    tracks[0]=trk1;
    tracks[1]=trk2;
    return massConstrain(ntrk,tracks,mass);
 }
 
 bool VertexFit::massConstrain(const CdfTrack* trk1, const CdfTrack* trk2,
                                 const CdfTrack* trk3, float mass)
 {
    int ntrk=3;
    const CdfTrack* tracks[3];
    tracks[0]=trk1;
    tracks[1]=trk2;
    tracks[2]=trk3;
    return massConstrain(ntrk,tracks,mass);
 }
 
 bool VertexFit::massConstrain(const CdfTrack* trk1, const CdfTrack* trk2,
                     const CdfTrack* trk3, const CdfTrack* trk4,float mass)
 {
    int ntrk=4;
    const CdfTrack* tracks[4];
    tracks[0]=trk1;
    tracks[1]=trk2;
    tracks[2]=trk3;
    tracks[3]=trk4;
    return massConstrain(ntrk,tracks,mass);
 }
 
 bool VertexFit::massConstrain(int ntrk, const CdfTrack* tracks[], float mass)
 {
 // Check that we have not exceeded the allowed number of mass constraints.

    if(_ctvmq.nmassc>=_maxmcn) return false;
 
 // Set constraint mass

    _ctvmq.cmass[_ctvmq.nmassc]=mass;
      
 // For each track in contraint, set trkmcn true.  Since the number in

 // tracks[] is the track number, we have to find each track in the

 // list of tracks.

    for(int jt=0; jt<ntrk; ++jt) {
       bool found=false;
       for(int kt=0; kt<_ctvmq.ntrack; ++kt) {
          if(tracks[jt]->id().value()==_ctvmq.list[kt]) {
             _ctvmq.trkmcn[_ctvmq.nmassc][kt]=true;
             found=true;
          }
       }
       if(!found) return false;
    }
 // Increment number of mass constraints.   

    _ctvmq.nmassc++;
 
    return true;
 }
 
 
 /*=======================================================================

                      VertexFit::beamlineConstraint

   =======================================================================*/
 //  Turn ON this option only when fitting the primary with no additional 

 // constraints. The routine is protected and will not function

 // if those options are selected as well as the beamline constraint.

 // In case this option is chosen the fit will calculate a result also with

 // just one track in input.

 // In order to enable this option the user must do the following:

 // Fit the primary vertex as VERTEX_1 (without any other vertices).

 // Enable beamlineConstraint by setting the pointing constraint variable

 // vtxpnt[0][0] to -100 (no pointing constraints when <0)

 // Note that index 0,0 corresponds to index 1,1 in fortran

 // Provide the beamline parameters:

 //  --> assume xv = xb + xzbslope*z, 

 //             yv = yb + yzbslope*z, where (xv, yv) is the transverse

 //      beam position at z, (xb, yb) is the transverse beam position at z = 0 

 //      and (xzbslope, yzbslope) are the slopes of the beamline

 //  --> set primary vertex at z=0 (xb, yb, 0)

 //  --> set the diagonal elements of the primary vertex error matrix to the

 //      sigma**2 of beam spread in x, y and z:

 //      xverr[1][1] =  (~30 - 50 um)**2

 //      xverr[2][2] =  (~30 - 50 um)**2

 //      xverr[3][3] =  (~30 cm) **2

 //      All other elements should be 0

 // For help email: Joao Guimaraes  guima@huhepl.harvard.edu

 //                 Franco Bedeschi bed@fnal.gov 

 // See more information in the ctvmft.f

 bool VertexFit::beamlineConstraint(float xb, float yb, HepSymMatrix berr,
                                    float xzbslope, float yzbslope)
 {
   // Set beam position at z=0

   setPrimaryVertex(xb,yb,0);
   if (_extrapolateTrackErrors) {
     float newXb = xb - _referencePoint.x() + 
       _referencePoint.z() * xzbslope;
     float newYb = yb - _referencePoint.y() + 
       _referencePoint.z() * yzbslope;
     setPrimaryVertex( newXb, newYb, 0);
   }
 
   bool success = setPrimaryVertexError(berr);
   //Set the beamline slope values

   _ctvmq.xzslope= xzbslope;
   _ctvmq.yzslope= yzbslope;
   // Turn ON beamline constraint

   _ctvmq.vtxpnt[0][0]=-100; 
 
   return success;
 }
 
 bool VertexFit::beamlineConstraint(Hep3Vector pv, HepSymMatrix berr,
                                    float xzbslope, float yzbslope)
 {
   // Check if input beam position coordinates are at z=0

   if (pv.z() != 0) return false;
 
   return beamlineConstraint(pv.x(),pv.y(),berr,xzbslope,yzbslope);
 }
 
 /*=======================================================================

                          VertexFit::setPrimaryVertex

   =======================================================================*/
 
 void VertexFit::setPrimaryVertex(float xv, float yv, float zv)
 {
 // Set x,y,z position of the primary vertex.

    _ctvmq.xyzpv0[0]=xv;
    _ctvmq.xyzpv0[1]=yv;
    _ctvmq.xyzpv0[2]=zv;
 
    _primaryVertex = Hep3Vector( xv, yv, zv );
    return;
 }
 
 void VertexFit::setPrimaryVertex(Hep3Vector pv)
 {
 // Set x,y,z position of the primary vertex.

    _ctvmq.xyzpv0[0]=pv.x();
    _ctvmq.xyzpv0[1]=pv.y();
    _ctvmq.xyzpv0[2]=pv.z();
 
    _primaryVertex = pv;
    return;
 }
 
 /*=======================================================================

                          VertexFit::setPrimaryVertexError

   =======================================================================*/
 
 void VertexFit::setPrimaryVertexError(const float xverr[3][3])
 {
 // Set the error matrix for the primary vertex.

    for(int j=0; j<3; ++j) {
       for(int k=0; k<3; ++k) {
          _ctvmq.exyzpv[j][k]=xverr[j][k];
       }
    }
    return;
 }
 
 bool VertexFit::setPrimaryVertexError(const HepSymMatrix &xverr)
 {
 // Set the error matrix for the primary vertex using a HepSymMatrix.

 // First check that the matrix is the correct size.

    if(xverr.num_row() != 3) return false;
    for(int j=0; j<3; j++) {
      for(int k=0; k<3; k++) {
        _ctvmq.exyzpv[j][k]=xverr[j][k];
      }
    }
    return true;
 }
 
 /*=======================================================================

                          VertexFit::fit

   =======================================================================*/
 
 bool VertexFit::fit()
 {
 // Check that the diagonal elements of all the track error matrices

 //are positive.

   bool mstat = true;
   for(int trk=0; trk<_ctvmq.ntrack; ++trk) {
     for(int j=0; j<5; ++j) {
 // Check diagonal elements of error matrix.

       if(_trkprm.trem[trk][j][j] < 0.) {
         mstat = false;
 // Record offending track

         _ctvmq.ijkerr[2] = trk + 1;
       }
     }
   }
   if(!mstat) {
 // Tell user that the covariance matrix could not be inverted.

     _ctvmq.ijkerr[0] = 3;
     _ctvmq.ijkerr[1] = 2;
     stat = 1;
     return false;
   }
 
 // First copy information into CTVMFT common blocks

    *_ctvmq_com=_ctvmq;
    *_ctvmfr_com=_ctvmfr;
    *_fiddle_com=_fiddle;
    *_trkprm_com=_trkprm;
 // Do the vertex fit.

    int print=0;
    int level=0;
    ctvmft_(print,level,stat);
 // Now copy information from CTVMFT common blocks to local storage

    _ctvmq=*_ctvmq_com;
    _ctvmfr=*_ctvmfr_com;
    _fiddle=*_fiddle_com;
    _trkprm=*_trkprm_com;
    
    return stat==0;
 }
 
 /*=======================================================================

                          VertexFit::print

   =======================================================================*/
 
 void VertexFit::print() const
 {
    print(std::cout);
 }
 
 void VertexFit::print(std::ostream& os) const
 {
    os << "****************************** VertexFit "
       << "******************************" << std::endl;
    os << "Number of tracks: " << _ctvmq.ntrack << std::endl;
    os << "   Tracks: ";
    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
       if(jt != 0) os << ",  ";
       CdfTrack::Id trkId(_ctvmq.list[jt]);
       os << trkId;
    }
    os << std::endl;
    os << "Number of vertices: " << _ctvmq.nvertx << std::endl;
    for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
       os << "   Vertex " << jv+1 << " tracks: ";
       for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
          if(_ctvmq.trkvtx[jv][jt]) {
             os << " " << _ctvmq.list[jt];
          }
       }
       os << std::endl;
    }
    for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
       if(_ctvmq.vtxpnt[0][jv]==0) {
          os << "   Vertex " << jv+1 << " points to the primary vertex ";
       }
       else if(_ctvmq.vtxpnt[0][jv]>0) {
          os << "   Vertex " << jv+1 << " points to vertex "
               << _ctvmq.vtxpnt[0][jv];
       }
       if(_ctvmq.vtxpnt[1][jv]==1) {
          os << " in 2 dimensions" << std::endl;
       }
       else if(_ctvmq.vtxpnt[1][jv]==2) {
          os << " in 3 dimensions" << std::endl;
       }
       else if(_ctvmq.vtxpnt[1][jv]==3) {
          os << ", a single track vertex" << std::endl;
       }
       if(_ctvmq.cvtx[jv]>0) {
          os << "   Vertex " << jv+1 << " is a conversion" << std::endl;
       }
    }
    os << "Number of mass constraints: " << _ctvmq.nmassc << std::endl;
    for(int jmc=0; jmc<_ctvmq.nmassc; ++jmc) {
       os << "   Tracks ";
       for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
          if(_ctvmq.trkmcn[jmc][jt]) {
             os << " " << _ctvmq.list[jt];
          }
       }
       os << " constrained to mass " << _ctvmq.cmass[jmc]
          << " Gev/c^2" << std::endl;
    }
    if(stat==-999) {
       os << "No fit has been done." << std::endl;
    }
    else {
       os << "***** Results of Fit *****" << std::endl;
       printErr(os);
       os << "   Status = " << stat << std::endl;
       os << "   Chi-square = " << _ctvmq.chisqr[0]
          << " for " << _ctvmq.ndof << " degrees of freedom." << std::endl;
       os << "   => probability = " << prob() << std::endl;
       for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
          os << "Vertex " << jv+1 << " position: "
             << _ctvmq.xyzvrt[jv+1][0] << " "
             << _ctvmq.xyzvrt[jv+1][1] << " "
             << _ctvmq.xyzvrt[jv+1][2] << std::endl;
       }
       for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
          os << "Track " << _ctvmq.list[jt] << " P4: " 
             << _ctvmq.trkp4[0][jt] << " "
             << _ctvmq.trkp4[1][jt] << " "
             << _ctvmq.trkp4[2][jt] << " "
             << _ctvmq.trkp4[3][jt] << " " << std::endl;
       }
    }
    os << "****************************************"
       << "**************************" << std::endl;
 
    return;
 }
 
 /*=======================================================================

                          VertexFit::printErr

   =======================================================================*/
 
 void VertexFit::printErr() const
 {
    printErr(std::cout);
 }
 
 void VertexFit::printErr(std::ostream& os) const
 {
    os << "VertexFit: IJKERR = " << _ctvmq.ijkerr[0] << ", "
                         << _ctvmq.ijkerr[1] << ", "
                         << _ctvmq.ijkerr[2] << std::endl;
    if(status()==0 && _ctvmq.ijkerr[0]==0) return;
    if(_ctvmq.ijkerr[0] == -1) {
       os << "   Problem with GETTRK:  track requested is not in list."
          << std::endl
          << "   This should not happen - Contact VertexFit expert "
          << _expert << "." <<std::endl;
    }
    else if(_ctvmq.ijkerr[0]==1) {
       os << "   Problem in CTVM00:" << std::endl;
       if(_ctvmq.ijkerr[1]==1) {
          os << "      Number of tracks is " << _ctvmq.ntrack
             << "." << std::endl;
          if(_ctvmq.ntrack < 2) {
             os << ", which is too few (must be at least 2)." << std::endl;
          }
          else if(_ctvmq.ntrack > _maxtrk) {
             os << ", which is too many (maximum is " << _maxtrk
                << ")." << std::endl;
          }
          else {
             os << "      Problem with number of tracks"
                << " for unknown reasons." << std::endl;
          }
       }
       else if(_ctvmq.ijkerr[1]==2) {
          os << "      Number of vertices is " << _ctvmq.nvertx
             << "." << std::endl;
          if(_ctvmq.nvertx < 1) {
             os << ", which is too few (must be at least 1)." << std::endl;
          }
          else if(_ctvmq.nvertx > _maxvtx) {
             os << ", which is too many (maximum is " << _maxvtx
                << ")." << std::endl;
          }
          else {
             os << std::endl << "      Problem with number of vertices"
                << " for unknown reasons." << std::endl;
          }
       }
       else if(_ctvmq.ijkerr[1]==3) {
          os << "      Number of mass constraints is " << _ctvmq.nmassc
             << "." << std::endl;
          if(_ctvmq.nmassc < 0) {
             os << ", which is negative." << std::endl;
          }
          else if(_ctvmq.nmassc > _maxmcn) {
             os << ", which is too many (maximum is " << _maxmcn
                << ")." << std::endl;
          }
          else {
             os << std::endl << "      Problem with number of mass"
                << " constraints for unknown reasons." << std::endl;
          }
       }
       else if(_ctvmq.ijkerr[1]==11) {
          os << "      Vertex " << _ctvmq.ijkerr[2]
             << " has less than one track." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==12) {
          os << "      Vertex " << _ctvmq.ijkerr[2]
             << " is a conversion vertex with a number of tracks"
             << " different than two." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==13) {
          os << "      Vertex " << _ctvmq.ijkerr[2]
             << " is a one track vertex that has no multi-track"
             << " descendents." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==14) {
          os << "      Vertex " << _ctvmq.ijkerr[2]
             << " does not point at a vertex with a lower number."
             << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==15) {
          os << "      Vertex " << _ctvmq.ijkerr[2]
             << " has a parent vertex that is a conversion." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==16) {
          os << "      Vertex " << _ctvmq.ijkerr[2]
             << " does 1 track pointing to a vertex with"
             << " more than 1 track." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==17) {
          os << "      Vertex " << _ctvmq.ijkerr[2]
             << " does 0 track pointing to a vertex with"
             << " more than 0 track (?)." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==19) {
          os << "      Primary vertex error matrix is singular." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==21) {
          os << "      Track with Id " << _ctvmq.ijkerr[2]
             << "is not in any vertex." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==22) {
          os << "      Track with Id " << _ctvmq.ijkerr[2]
             << "is in multiple vertices." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==23) {
          os << "      Track with Id " << _ctvmq.ijkerr[2]
             << "occurs more than once." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==31) {
          os << "      A mass constraint has less than 2 tracks." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==32) {
          os << "      The sum masses of the tracks in a mass constraint"
             << " exceeds the constraint mass." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==33) {
          os << "      Beamline constraint. Beam covariance not set properly."
             << " Negative diagonal elements." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==34) {
          os << "      Beamline constraint. Beam covariance not set properly."
             << " Off-diagonal elements not zero." << std::endl;
       }
       else if(_ctvmq.ijkerr[1]==36) {
          os << "      Beamline constraint. Number of vertices = " 
             << _ctvmq.nvertx << " Should be 1." << std::endl;
       }
    }
    else if(_ctvmq.ijkerr[0] == 2) {
       if(_ctvmq.ijkerr[1] == 20) {
          os << "   Problem in CTVM00: " << std::endl;
          os << "      Track has negative Id = "
             << _ctvmq.list[_ctvmq.ijkerr[2]-1] << "." << std::endl;
       }
       else {
          os << "   Problem in CTVMFA with vertex "
             << _ctvmq.ijkerr[2] << ": " << std::endl;
          os << "      Failure in vertex first approximation." << std::endl;
          if(_ctvmq.ijkerr[1] == 1) {
             os << "      Tracks are concentric circles." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 2) {
             os << "      Conversion vertex has widely separated"
                << " exterior circles at midpoint." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 3) {
             os << "      Conversion vertex has widely separated"
                << " interior circles at midpoint." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 4) {
             os << "      Vertex has widely separated"
                << " exterior circles at approximate vertex." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 5) {
             os << "      Vertex has widely separated"
                << " interior circles at approximate vertex." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 6) {
             os << "      Rv is too large at the chosen"
                << " intersection point." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 7) {
             os << "      Delta z is too large at the chosen"
                << " intersection point." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 8) {
             os << "      A track's turning to the chosen vertex"
                << " is too large." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 9) {
             os << "      There is no solution with an adequately"
                << " positive arc length." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 21) {
             os << "      zero-track vertexing: either/both vertex "
                << " momenta are too small (<0.01 MeV)." << std::endl;
          }
          if(_ctvmq.ijkerr[1] == 22) {
             os << "      zero-track vertexing: Two lines (tracks) are "
                << " parallel/antiparallel." << std::endl;
          }
 
       }
    }
    else if(_ctvmq.ijkerr[0] == 3) {
       os << "   Problem in CTVM01 with track with Id = "
          << _ctvmq.list[_ctvmq.ijkerr[2]-1] << ": " << std::endl;
       if(_ctvmq.ijkerr[1] == 1) {
          os << "      GETTRK cannot find Id in list." << std::endl;
       }
       if(_ctvmq.ijkerr[1] == 2) {
          os << "      Covariance matrix could not be inverted."  << endl 
             << "      Offending track number (in order addded) is "
             << _ctvmq.ijkerr[2] << "." << std::endl;
       }
       if(_ctvmq.ijkerr[1] == 3) {
          os << "      Track turns through too large an angle"
             << " to the vertex." << std::endl;
       }
       if(_ctvmq.ijkerr[1] == 4) {
          os << "      Track moves too far backward to vertex." << std::endl;
       }
    }
    else if(status() == 9) {
       os << "   General fit problem: " << std::endl;
       if(_ctvmq.ijkerr[1] == 1) {
          os << "      Singular solution matrix." << std::endl;
       }
       if(_ctvmq.ijkerr[1] == 2 || _ctvmq.ijkerr[1] == 3) {
          os << "      Too many iterations ( "
             << _ctvmq.ijkerr[2] << "(." << std::endl;
       }
       if(_ctvmq.ijkerr[1] == 4) {
          os << "      Convergence failure." << std::endl;
       }
       if(_ctvmq.ijkerr[1] == 5) {
          os << "      Bad convergence." << std::endl;
       }
       if(_ctvmq.ijkerr[1] == 9) {
          os << "      Ill-formed  covariance matrix." << std::endl;
       }
    }
    else {
       os << "   The error codes above are not recognized." << std::endl
          << "   Contact VertexFit expert " << _expert << "." << std::endl;
    }
    return;
 }
 
 
 /*=======================================================================

                          VertexFit::getIJKErr

   =======================================================================*/
 
 void VertexFit::getIJKErr(int& err0, int& err1, int& err2) const
 {
   err0 = _ctvmq.ijkerr[0];
   err1 = _ctvmq.ijkerr[1];
   err2 = _ctvmq.ijkerr[2];
   return;
 }
 
 int VertexFit::getIJKErr0() const
 {
   return _ctvmq.ijkerr[0];
 }
 int VertexFit::getIJKErr1() const
 {
   return _ctvmq.ijkerr[1];
 }
 int VertexFit::getIJKErr2() const
 {
   return _ctvmq.ijkerr[2];
 }
 
 /*=======================================================================

                        VertexFit::getErrTrackId()

   =======================================================================*/
 
 int VertexFit::getErrTrackId() const
  {
    if (status() == 0) return 0;
    int trkId = 0;
    // Problems with track in CTVM01 or track has negative id in CTVM00

    // See PrintErr() for a more detailed list of error codes.

    if ((_ctvmq.ijkerr[0] == 2 && _ctvmq.ijkerr[1] == 20) ||
        _ctvmq.ijkerr[0] == 3) {
      trkId = _ctvmq.list[_ctvmq.ijkerr[2]-1];
    }
    
    return trkId;
  }
 
 /*=======================================================================

                          VertexFit::expert

   =======================================================================*/
 
 string VertexFit::expert() const
 {
    return _expert;
 }
 
 /*=======================================================================

                          VertexFit::status

   =======================================================================*/
 
 int VertexFit::status() const
 {
 
    return stat;
 }
 
 /*=======================================================================

                          VertexFit::chisq

   =======================================================================*/
 
 float VertexFit::chisq() const
 {
 // Chi-square of fit

       return _ctvmq.chisqr[0];
 }
 
 /*=======================================================================

                          VertexFit::ndof

   =======================================================================*/
 
 int VertexFit::ndof() const
 {
 // Number of degrees of freedom of fit.

    if(_ctvmq.chisqr[0] >= 0) {
       return _ctvmq.ndof; }
    else {
       return 0;}
 }
 
 /*=======================================================================

                          VertexFit::prob

   =======================================================================*/
 
 float VertexFit::prob() const
 {
 // Probability of chi-square of fit

    if(_ctvmq.chisqr[0]>=0.) {
       float chisq=_ctvmq.chisqr[0];
       int nd=_ctvmq.ndof;
       return prob_(chisq,nd);
    }
    else {
       return -999.;
    }
 }
 
 /*=======================================================================

                          VertexFit::chisq(CdfTrack*)

   =======================================================================*/
 
 float VertexFit::chisq(const CdfTrack* trk) const
 {
 // This method returns the chisquare contribution for one track 

 // If fit not successful or not done yet, return -1.

    if(_ctvmq.chisqr[0] < 0) return -1.;
 // Look for this track in the list of tracks.

    for(int jt = 0; jt < _ctvmq.ntrack; ++jt) {
      if(trk->id().value() == _ctvmq.list[jt]) {
 // Found the track, so return its chisquare contribution.

        return _ctvmq.chit[jt];
      }
    }
 // If track is not in list, return -1.

    return -1.;
 }
 
 /*=======================================================================

                          VertexFit::chisq_rphi()

   =======================================================================*/
 
 float VertexFit::chisq_rphi() const
 {
 // This method returns the chisquare contribution in the r-phi plane.

    int index[3] = {0,1,3};
 // If fit not successful or not done yet, return -1.

    if(_ctvmq.chisqr[0] < 0) return -1.;
 // Loop over the tracks in the event.

    float chisq = 0.;
    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
 // Double loop over the relevant parameter indices.

      for(int k1=0; k1<3; ++k1) {
        for(int k2=0; k2<3; ++k2) {
 // Add contribution to chisquare.

          chisq += _ctvmq.pardif[jt][index[k1]] *
                   _ctvmq.g[jt][index[k1]][index[k2]] *
                   _ctvmq.pardif[jt][index[k2]];
        }
      }
    }
 // Return the chisquare.

    return chisq;
 }
 
 
 /*=======================================================================

                          VertexFit::chisq_z()

   =======================================================================*/
 
 float VertexFit::chisq_z() const
 {
 // This method returns the chisquare contribution in the z direction.

    int index[2] = {2,4};
 // If fit not successful or not done yet, return -1.

    if(_ctvmq.chisqr[0] < 0) return -1.;
 // Loop over the tracks in the event.

    float chisq = 0.;
    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
 // Double loop over the relevant parameter indices.

      for(int k1=0; k1<2; ++k1) {
        for(int k2=0; k2<2; ++k2) {
 // Add contribution to chisquare.

          chisq += _ctvmq.pardif[jt][index[k1]] *
                   _ctvmq.g[jt][index[k1]][index[k2]] *
                   _ctvmq.pardif[jt][index[k2]];
        }
      }
    }
 // Return the chisquare.

    return chisq;
 }
 
 
 /*=======================================================================

                          VertexFit::chisq_rphiz()

   =======================================================================*/
 
 float VertexFit::chisq_rphiz() const
 {
 // This method returns the chisquare contribution of the cross

 // terms in the r-phi and z directions.

    int index1[2] = {2,4};
    int index2[3] = {0,1,3};
 // If fit not successful or not done yet, return -1.

    if(_ctvmq.chisqr[0] < 0) return -1.;
 // Loop over the tracks in the event.

    float chisq = 0.;
    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
 // Double loop over the relevant parameter indices.

      for(int k1=0; k1<2; ++k1) {
        for(int k2=0; k2<3; ++k2) {
 // Add contribution to chisquare.

          chisq += _ctvmq.pardif[jt][index1[k1]] *
                   _ctvmq.g[jt][index1[k1]][index2[k2]] *
                   _ctvmq.pardif[jt][index2[k2]];
        }
      }
    }
 // Return the chisquare.

    return 2.*chisq;
 }
 
 /*=======================================================================

                          VertexFit::chisq_rphi(CdfTrack*)

   =======================================================================*/
 
 float VertexFit::chisq_rphi(const CdfTrack* trk) const
 {
 // This method returns the chisquare contribution in the r-phi plane.

    int index[3] = {0,1,3};
 // If fit not successful or not done yet, return -1.

    if(_ctvmq.chisqr[0] < 0) return -1.;
 // Loop over the tracks in the event, looking for the one we want

    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
      if(trk->id().value() == _ctvmq.list[jt]) {
 // Found the track, so calculate its chisquare contribution.

        float chisq = 0.;
 // Double loop over the relevant parameter indices.

        for(int k1=0; k1<3; ++k1) {
          for(int k2=0; k2<3; ++k2) {
 // Add contribution to chisquare.

            chisq += _ctvmq.pardif[jt][index[k1]] *
                     _ctvmq.g[jt][index[k1]][index[k2]] *
                     _ctvmq.pardif[jt][index[k2]];
          }
        }
        return chisq;
      }
    }
 // Track not found, return -1.

    return -1.;
 }
 
 
 /*=======================================================================

                          VertexFit::chisq_z(CdfTrack*)

   =======================================================================*/
 
 float VertexFit::chisq_z(const CdfTrack* trk) const
 {
 // This method returns the chisquare contribution in the z direction.

    int index[2] = {2,4};
 // If fit not successful or not done yet, return -1.

    if(_ctvmq.chisqr[0] < 0) return -1.;
 // Loop over the tracks in the event, looking for the one we want.

    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
      if(trk->id().value() == _ctvmq.list[jt]) {
 // Found the track, so calculate its chisquare contribution.

        float chisq = 0.;
 // Double loop over the relevant parameter indices.

        for(int k1=0; k1<2; ++k1) {
          for(int k2=0; k2<2; ++k2) {
 // Add contribution to chisquare.

            chisq += _ctvmq.pardif[jt][index[k1]] *
                     _ctvmq.g[jt][index[k1]][index[k2]] *
                     _ctvmq.pardif[jt][index[k2]];
          }
        }
        return chisq;
      }
    }
 // Track not found - return -1.

    return -1.;
 }
 
 
 /*=======================================================================

                          VertexFit::chisq_rphiz(CdfTrack*)

   =======================================================================*/
 
 float VertexFit::chisq_rphiz(const CdfTrack* trk) const
 {
 // This method returns the chisquare contribution of the cross

 // terms in the r-phi and z directions.

    int index1[2] = {2,4};
    int index2[3] = {0,1,3};
 // If fit not successful or not done yet, return -1.

    if(_ctvmq.chisqr[0] < 0) return -1.;
 // Loop over the tracks in the event.

    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
      if(trk->id().value() == _ctvmq.list[jt]) {
 // Found the track, so calculate its chisquare contribution.

        float chisq = 0.;
 // Double loop over the relevant parameter indices.

        for(int k1=0; k1<2; ++k1) {
          for(int k2=0; k2<3; ++k2) {
 // Add contribution to chisquare.

            chisq += _ctvmq.pardif[jt][index1[k1]] *
                     _ctvmq.g[jt][index1[k1]][index2[k2]] *
                     _ctvmq.pardif[jt][index2[k2]];
          }
        }
        return 2.*chisq;
      }
    }
 // Track not found so return -1.

    return -1.;
 }
 
 /*=======================================================================

                          VertexFit::getTrackP4

   =======================================================================*/
 
 HepLorentzVector VertexFit::getTrackP4(const CdfTrack* trk) const
 {
    if(stat != 0) return HepLorentzVector(0,0,0,0);
 // return four momentum of fit track 

    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
 //   Find which track matches this Id.

       if(trk->id().value()==_ctvmq.list[jt]) {
          HepLorentzVector p4((double)_ctvmq.trkp4[0][jt],
                              (double)_ctvmq.trkp4[1][jt],
                              (double)_ctvmq.trkp4[2][jt],
                              (double)_ctvmq.trkp4[3][jt]);
          return p4;
       }
    }
    return HepLorentzVector(0,0,0,0);
 }
 
 /*=======================================================================

                          VertexFit::getHelix

   =======================================================================*/
 
 Helix VertexFit::getHelix(const CdfTrack* trk) const
 {
    if(stat != 0) return Helix(0,0,0,0,0);
 // return helix parameters of fit track 

    for(int jt=0; jt<_ctvmq.ntrack; ++jt) {
 //   Find which track matches this Id.

       if(trk->id().value()==_ctvmq.list[jt]) {
          Helix trhelix((double)_ctvmq.par[jt][2],
                        (double)_ctvmq.par[jt][0],
                        (double)_ctvmq.par[jt][4],
                        (double)_ctvmq.par[jt][3],
                        (double)_ctvmq.par[jt][1]);
          return trhelix;
       }
    }
    return Helix(0,0,0,0,0);
 }
 
 /*=======================================================================

                          VertexFit::getMass

   =======================================================================*/
 
 float VertexFit::getMass(const CdfTrack* trk1,
                          const CdfTrack* trk2, float& dmass) const
 {
 
 //

 // Vertex Fitting is prone to floating point errors.  Here, we

 // turn those off during the duration of this routine: Joe and Farrukh

 //

 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
   sigaction(SIGFPE, &myaction, &oldaction);
   if (setjmp(env)!=0) {
     sigaction(SIGFPE, &oldaction,0);
     return -999;
   }
 #endif
 
    dmass=-999.;
    if(stat!=0) return -999.;
 // Get the fit invariant mass of tracks trk1 and trk2.

 // dmass is the error on the mass.

    int ntrk=2;
    const CdfTrack* trks[2];
    trks[0]=trk1;
    trks[1]=trk2;
 
   double result = getMass(ntrk,trks,dmass);
    //return something for FPE Protection-> Joe and Farrukh

 
 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
   sigaction(SIGFPE, &oldaction,0);
 #endif
 
   return result;
 }
 
 float VertexFit::getMass(const CdfTrack* trk1, const CdfTrack* trk2,
                                 const CdfTrack* trk3,float& dmass) const
 {
 
 //

 // Vertex Fitting is prone to floating point errors.  Here, we

 // turn those off during the duration of this routine: Joe and Farrukh

 //

 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
   sigaction(SIGFPE, &myaction, &oldaction);
   if (setjmp(env)!=0) {
     sigaction(SIGFPE, &oldaction,0);
     return -999;
   }
 #endif
 
 
 
    dmass=-999.;
    if(stat!=0) return -999.;
 // Get the fit invariant mass of tracks trk1, trk2, and trk3.

 // dmass is the error on the mass.

    int ntrk=3;
    const CdfTrack* trks[3];
    trks[0]=trk1;
    trks[1]=trk2;
    trks[2]=trk3;
 
 
    double result = getMass(ntrk,trks,dmass);
 
 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
   sigaction(SIGFPE, &oldaction,0);
 #endif
 
   return result;
 }
 
 float VertexFit::getMass(const CdfTrack* trk1, const CdfTrack* trk2,
      const CdfTrack* trk3, const CdfTrack* trk4, float& dmass) const
 {
    dmass=-999.;
    if(stat!=0) return -999.;
 // Get the fit invariant mass of tracks trk1, trk2, trk3, and trk4.

 // dmass is the error on the mass.

    int ntrk=4;
    const CdfTrack* trks[4];
    trks[0]=trk1;
    trks[1]=trk2;
    trks[2]=trk3;
    trks[3]=trk4;
 //

 // Vertex Fitting is prone to floating point errors.  Here, we

 // turn those off during the duration of this routine: Joe and Farrukh

 //

 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
   sigaction(SIGFPE, &myaction, &oldaction);
   if (setjmp(env)!=0) {
     sigaction(SIGFPE, &oldaction,0);
     return -999;
   }
 #endif
 
   double result = getMass(ntrk,trks,dmass);
 
    //return something for FPE Protection-> Joe and Farrukh

 
 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
   sigaction(SIGFPE, &oldaction,0);
 #endif
 
   return result;
 }
 
 float VertexFit::getMass(int ntrk, const CdfTrack* tracks[],
                                                 float& dmass) const
 {
 
 // Vertex Fitting is prone to floating point errors.  Here, we

 // turn those off during the duration of this routine: Joe and Farrukh

 //

 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
   struct sigaction myaction = {VertexFitSetStatus, 0, 0, 0}, oldaction;
   sigaction(SIGFPE, &myaction, &oldaction);
   if (setjmp(env)!=0) {
     sigaction(SIGFPE, &oldaction,0);
     return -999;
   }
 #endif
 
    dmass=-999.;
    if(stat!=0) return -999.;
 // Get fit invariant mass of ntrk tracks listed in array tracks.

 // dmass is the error on the mass.

    dmass=-999.;
    if(ntrk <= 0) return 0;
    int jtrks[_maxtrk];
    for(int jt=0; jt<ntrk; ++jt) {
       bool found=false;
       for(int kt=0; kt<_ctvmq.ntrack; ++kt) {
          if(tracks[jt]->id().value()==_ctvmq.list[kt]) {
             found=true;
             jtrks[jt]=kt+1;
          }
       }
       if(!found) return 0;
    }
 // Copy information into CTVMFT common blocks

    *_ctvmq_com=_ctvmq;
    *_ctvmfr_com=_ctvmfr;
    int ntr=ntrk;
    float mass;
    double p4[4];
    mcalc_(ntr, jtrks, mass, dmass, p4);
 
    //return something for FPE-> Joe and Farrukh

 #if ( defined(LINUX) && defined(__USE_BSD) ) || defined(OSF1)
   sigaction(SIGFPE, &oldaction,0);
 #endif
 
    return mass;
 }
 
 /*=======================================================================

                          VertexFit::getDecayLength

   =======================================================================*/
 
 float VertexFit::getDecayLength(vertexNumber nv, vertexNumber mv, 
                                 const Hep3Vector& dir, float& dlerr) const
 {
    dlerr=-999.;
    if(stat!=0) return -999.;
 // Get the signed decay length from vertex nv to vertex mv

 // along the x-y direction defined by the 3 vector dir.

 // dlerr is the error on the decay length including the

 // full error matrix.

 // Check that vertices are within range.

    if(nv<0 || nv>=_ctvmq.nvertx) return -999.;
    if(mv<1 || mv>_ctvmq.nvertx) return -999.;
    if(nv>=mv) return -999.;
    float dir_t=sqrt(dir.x()*dir.x()+dir.y()*dir.y());
    if(dir_t <= 0.) return -999.;
    Hep3Vector dv=getVertex(mv)-getVertex(nv);
    float dl=(dv.x()*dir.x()+dv.y()*dir.y())/dir_t;
 // Set up the column matrix of derivatives

    HepMatrix A(4,1);
    A(1,1)=dir.x()/dir_t;
    A(2,1)=dir.y()/dir_t;
    A(3,1)=-dir.x()/dir_t;
    A(4,1)=-dir.y()/dir_t;
 // Check if first vertex (nv) is primary vertex.  If it is,

 // check if it was used in the primary vertex.  If not,

 // all of the corresponding error matrix elements are those

 // supplied for the primary vertex.

    int nvf=0;
    if(nv==0) {
       nvf=-1;
       for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
          if(_ctvmq.vtxpnt[0][jv]==0) nvf=0;
       }
    }
 // Get the relevant submatrix of the full error matrix.

    HepMatrix V(4,4,0);
    if(nvf < 0) {
       V(1,1)=getErrorMatrix(_ctvmq.voff[mv-1]+1,_ctvmq.voff[mv-1]+1);
       V(1,2)=getErrorMatrix(_ctvmq.voff[mv-1]+1,_ctvmq.voff[mv-1]+2);
       V(2,1)=getErrorMatrix(_ctvmq.voff[mv-1]+2,_ctvmq.voff[mv-1]+1);
       V(2,2)=getErrorMatrix(_ctvmq.voff[mv-1]+2,_ctvmq.voff[mv-1]+2);
       V(3,3)=_ctvmq.exyzpv[0][0];
       V(3,4)=_ctvmq.exyzpv[0][1];
       V(4,3)=_ctvmq.exyzpv[1][0];
       V(4,4)=_ctvmq.exyzpv[1][1];
    }
    else {
 //    Get the indices into the error matrix vmat

       int index[4]={_ctvmq.voff[mv-1]+1,_ctvmq.voff[mv-1]+2,0,0};
       if(nv==0) {
          index[2]=1;
          index[3]=2;
       }
       else {
         index[2] = _ctvmq.voff[nv-1]+1;
         index[3] = _ctvmq.voff[nv-1]+2;
       }
       for(int j=0; j<4; ++j) {
          for(int k=0; k<4; ++k) {
             V[j][k]=getErrorMatrix(index[j],index[k]);
          }
       }
    }
 // Calculate square of dlerr

    dlerr=(A.T()*V*A)(1,1);
    if(dlerr >= 0.) {
       dlerr=sqrt(dlerr);
    }
    else {
       dlerr=-sqrt(-dlerr);
    }
 
    return dl;
 }
 
 /*=======================================================================

                          VertexFit::get_dr

   =======================================================================*/
 
 float VertexFit::get_dr(vertexNumber nv, vertexNumber mv,
                                            float& drerr) const
 {
    drerr=-999.;
    if(stat!=0) return -999.;
 // Get the transvese distance between vertices nv and mv

 // and return it as the function value.  drerr is the

 // uncertainty on the transverse distance, calculated from the

 // full error matrix including correlations.

    float dxyz[3];
    float dr;
    float dz;
    float dl[3];
 
    int mvert=mv;
    int nvert=nv;
 // Copy information into CTVMFT common blocks

    *_ctvmq_com=_ctvmq;
    *_ctvmfr_com=_ctvmfr;
 // Do calculation

    dcalc_(mvert,nvert,dxyz,dr,dz,dl);
    drerr=dl[0];
    return dr;
 }
 
 /*=======================================================================

                          VertexFit::get_dz

   =======================================================================*/
 
 float VertexFit::get_dz(vertexNumber nv, vertexNumber mv,
                                              float& dzerr) const
 {
    dzerr=-999.;
    if(stat!=0) return -999.;
 // Get the longitudinal distance between vertices nv and mv

 // and return it as the function value.  dzerr is the

 // uncertainty on the longitudinal distance, calculated from the

 // full error matrix including correlations.

    float dxyz[3];
    float dr;
    float dz;
    float dl[3];
 
    int mvert=mv;
    int nvert=nv;
 // Copy information into CTVMFT common blocks

    *_ctvmq_com=_ctvmq;
    *_ctvmfr_com=_ctvmfr;
 // Do calculation

    dcalc_(mvert,nvert,dxyz,dr,dz,dl);
    dzerr=dl[1];
    return dz;
 }
 
 /*=======================================================================

                          VertexFit::getVertex

   =======================================================================*/
 
 Hep3Vector VertexFit::getVertex(vertexNumber nv) const
 {
    if(stat!=0) return -999.;
 // Return x,y,z position of vertex nv.

    if(nv<0 || nv>_ctvmq.nvertx) return Hep3Vector(0,0,0);
 // Check if first vertex (nv) is primary vertex.  If it is,

 // check if it was used in the primary vertex.  If not,

 // all of the corresponding error matrix elements are those

 // supplied for the primary vertex.

    int nvf=0;
    if(nv==0) {
       nvf=-1;
       for(int jv=0; jv<_ctvmq.nvertx; ++jv) {
          if(_ctvmq.vtxpnt[0][jv]==0) nvf=0;
       }
    }
    Hep3Vector vertex;
 // If primary vertex was not part of fit, take vertex

 // location as supplied.

    if(nvf < 0) {
       vertex.setX(_ctvmq.xyzpv0[0]);
       vertex.setY(_ctvmq.xyzpv0[1]);
       vertex.setZ(_ctvmq.xyzpv0[2]);
    }
    else{
       vertex.setX(_ctvmq.xyzvrt[nv][0]);
       vertex.setY(_ctvmq.xyzvrt[nv][1]);
       vertex.setZ(_ctvmq.xyzvrt[nv][2]);
    }
 // If we have a different reference point, need to add it back in

    vertex += _referencePoint;
    return vertex;
 }
 
 /*=======================================================================

                          VertexFit::getErrorMatrix

   =======================================================================*/
 
 double VertexFit::getErrorMatrix(int j, int k) const
 {
    if(stat!=0) return -999.;
 // Return the j,k element of the full error matrix VMAT.

 // Since the CTVMFT documentation assumes the indices

 // start from 1 (ala Fortran), we will also assume this

 // and convert the C++ indices.  Note also the that order of

 // Fortran and C++ indices is different.  We assume that

 // j and k are in the Fortran order.

    if(j<1 || k<1 || j>_maxdim+1 || k>_maxdim) {
       return -999.;
    }
    return _ctvmfr.vmat[k-1][j-1];
 }
 
 HepSymMatrix VertexFit::getErrorMatrix(VertexFit::vertexNumber nv) const 
 {  
 // Return errors for vertex nv.

 
   HepSymMatrix cov(3,0);
   
 // If this is the primary vertex, return the error matrix the

 // user supplied.

   if(nv==PRIMARY_VERTEX) {
     for(int j=0; j<3; j++)
       for(int k=0; k<3; k++) 
         cov[j][k] = _ctvmq.exyzpv[j][k];
     return cov;
   }
   
   if(stat!=0) return cov;
 // Return x,y,z position of vertex nv.

   if(nv<VERTEX_1 || nv>_ctvmq.nvertx) return cov;
 // get offset for vertex nv

   int voff = _ctvmq.voff[nv-1];
 // fill matrix

   for(int i = 0 ; i < 3 ; ++i)
     for(int j = i ; j < 3 ; ++j)
       cov[i][j] = _ctvmfr.vmat[voff+i][voff+j];
   return cov;
 }
 
 HepSymMatrix VertexFit::getErrorMatrix(const CdfTrack* trk) const 
 {
   HepSymMatrix cov(3,0);
   if (stat != 0) return cov;
 
   for (int nt=0; nt<_ctvmq.ntrack; ++nt) {
 
     // Find which track matches this Id

     if (trk->id().value() == _ctvmq.list[nt]) {
 
       // Position of track nt get offset for track nt

       int toff = _ctvmq.toff[nt];
 
       // Fill matrix -- Crv,Phi,Ctg

       for(int i = 0 ; i < 3 ; ++i)
         for(int j = i ; j < 3 ; ++j)
           cov[i][j] = _ctvmfr.vmat[toff+i][toff+j];
     }
   }
   return cov;
 }
 
 float VertexFit::getPtError(const CdfTrack* trk) const
 {  
   if (stat != 0) return 0;
   
   int   toff;
   float pt,curv,curvErr,ptErr;
 
   for (int nt=0; nt<_ctvmq.ntrack; ++nt) {
 
     // Find which track matches this Id

     if (trk->id().value() == _ctvmq.list[nt]) {
 
       // Position of track nt get offset for track nt

       toff = _ctvmq.toff[nt];
 
       // Curvature error

       pt      = sqrt(_ctvmq.trkp4[0][nt]*_ctvmq.trkp4[0][nt] +
                      _ctvmq.trkp4[1][nt]*_ctvmq.trkp4[1][nt]);
       curv    = _ctvmq.pscale/pt;
       curvErr = sqrt(_ctvmfr.vmat[toff+0][toff+0]);
       ptErr   = _ctvmq.pscale/curv/curv*curvErr;
       return ptErr;
     }
   }
 
   return 0;
 }
 
 /*=======================================================================

                          VertexFit::getPosMomErr

   =======================================================================*/
 void VertexFit::getPosMomErr(HepMatrix& errors) const
 {
   // A c++ rewrite of the FORTRAN MASSAG function

   // The result of this function is an error matrix in position-momentum basis.

   // A 7x7 matrix of errors where the rows/columns are x, y, z, px, py, pz, e.

 
   double cosph[_maxtrk];
   double sinph[_maxtrk];
   double cosdph[_maxtrk];
   double sindph[_maxtrk];
   Hep3Vector mom3[_maxtrk];
   HepLorentzVector pmom[_maxtrk];
   HepLorentzVector total_mom;
 
   for(int lvtx = 0; lvtx < _ctvmq.nvertx; lvtx++){
     for(int ltrk = 0; ltrk < _ctvmq.ntrack; ltrk++){
       if(!_ctvmq.trkvtx[lvtx][ltrk]) continue;
       cosph[ltrk] = cos(_ctvmq.par[ltrk][1]);
       sinph[ltrk] = sin(_ctvmq.par[ltrk][1]);
       double dphi = 0;
       sindph[ltrk] = 2 * _ctvmq.par[ltrk][0] * 
        (_ctvmq.xyzvrt[lvtx + 1][0] * cosph[ltrk] + 
         _ctvmq.xyzvrt[lvtx + 1][1] * sinph[ltrk]);
       cosdph[ltrk] = sqrt(1.0 - sindph[ltrk] * sindph[ltrk]);
       if(fabs(sindph[ltrk]) <= 1.0){
         dphi = asin(sindph[ltrk]);
       }
       double pt = _ctvmq.pscale * fabs(1./_ctvmq.par[ltrk][0]);
       mom3[ltrk].setX(pt * cos(_ctvmq.par[ltrk][1] + dphi));
       mom3[ltrk].setY(pt * sin(_ctvmq.par[ltrk][1] + dphi));
       mom3[ltrk].setZ(pt * _ctvmq.par[ltrk][2]);
       double e  = sqrt(_ctvmq.tmass[ltrk] * _ctvmq.tmass[ltrk] 
                        + mom3[ltrk].mag2());
       pmom[ltrk].setVect(mom3[ltrk]);
       pmom[ltrk].setE(e);
 
       total_mom += pmom[ltrk];
     }
   }
 
 // Easy so far, but now it gets ugly.

 // Fill dp_dpar with the derivatives of the position and momentum with respect

 // to the parameters.

 
   int ctvmft_dim = 3 * (_ctvmq.nvertx + _ctvmq.ntrack);
   HepMatrix deriv(ctvmft_dim, 7, 0);
   HepMatrix dp_dpar[_maxvtx] = {deriv, deriv, deriv};
 
 // Fill the x, y, z rows: 

 
   for(int nvtx = 0; nvtx < _ctvmq.nvertx; nvtx++){
     for(int lcomp = 0; lcomp < 3; lcomp++){
       dp_dpar[nvtx][(3 * nvtx) + lcomp][lcomp] = 1.0;
     }
 
 // Fill the px, py, pz, e rows:

 
     for(int lvtx = 0; lvtx < _ctvmq.nvertx; lvtx++){
       for(int ltrk = 0; ltrk < _ctvmq.ntrack; ltrk++){
         if(!_ctvmq.trkvtx[lvtx][ltrk]) continue;
 
 // Find the derivatives of dphi with respect to x, y, curvature, and phi0:

 
         double dphi_dx = 2.0 * _ctvmq.par[ltrk][0] * cosph[ltrk]/cosdph[ltrk];
         double dphi_dy = 2.0 * _ctvmq.par[ltrk][0] * sinph[ltrk]/cosdph[ltrk];
         double dphi_dc = 2.0 * 
           (_ctvmq.xyzvrt[lvtx + 1][0] * cosph[ltrk] + 
            _ctvmq.xyzvrt[lvtx + 1][1] * sinph[ltrk])/cosdph[ltrk];
         double dphi_dp = 2.0 * _ctvmq.par[ltrk][0] *
           (-_ctvmq.xyzvrt[lvtx + 1][0] * sinph[ltrk] + 
             _ctvmq.xyzvrt[lvtx + 1][1] * cosph[ltrk])/cosdph[ltrk];
 
 // Now load the derivative matrix

 
         int lvele = 3 * lvtx;
 // dPx/dx:

         dp_dpar[nvtx][lvele][3] += -pmom[ltrk].y() * dphi_dx;
 // dPy/dx:

         dp_dpar[nvtx][lvele][4] +=  pmom[ltrk].x() * dphi_dx;
 // dPz/dx:

         dp_dpar[nvtx][lvele][5]  = 0.; 
 // dPy/dx:

         dp_dpar[nvtx][lvele][6]  = 0.; 
 
         lvele++;
 // dPx/dy:

         dp_dpar[nvtx][lvele][3] += -pmom[ltrk].y() * dphi_dy;
 // dPy/dy:

         dp_dpar[nvtx][lvele][4] +=  pmom[ltrk].x() * dphi_dy;
 // dPz/dy:

         dp_dpar[nvtx][lvele][5]  = 0.; 
 // dE/dy:

         dp_dpar[nvtx][lvele][6]  = 0.; 
 
         lvele++;
 // dPx/dz:

         dp_dpar[nvtx][lvele][3]  = 0.;
 // dPy/dz:

         dp_dpar[nvtx][lvele][4]  = 0.;
 // dPz/dz:

         dp_dpar[nvtx][lvele][5]  = 0.; 
 // dE/dz:

         dp_dpar[nvtx][lvele][6]  = 0.; 
 
         lvele = 3 * (ltrk + _ctvmq.nvertx);
 // dPx/dcurv[ltrk]:

         dp_dpar[nvtx][lvele][3]  = -(pmom[ltrk].x()/_ctvmq.par[ltrk][0])
                                    - pmom[ltrk].y() * dphi_dc;
 // dPy/dcurv[ltrk]:

         dp_dpar[nvtx][lvele][4]  = -(pmom[ltrk].y()/_ctvmq.par[ltrk][0])
                                    + pmom[ltrk].x() * dphi_dc;
 // dPz/dcurv[ltrk]:

         dp_dpar[nvtx][lvele][5]  = -(pmom[ltrk].z()/_ctvmq.par[ltrk][0]); 
 // dE/dcurv[ltrk]:

         dp_dpar[nvtx][lvele][6]  = 
           -mom3[ltrk].mag2()/(_ctvmq.par[ltrk][0] * pmom[ltrk].e()); 
 
         lvele++;
 // dPx/dphi[ltrk]:

         dp_dpar[nvtx][lvele][3]  = -pmom[ltrk].y() * (1.0 + dphi_dp);
 // dPy/dphi[ltrk]:

         dp_dpar[nvtx][lvele][4]  =  pmom[ltrk].x() * (1.0 + dphi_dp);
 // dPz/dphi[ltrk]:

         dp_dpar[nvtx][lvele][5]  = 0.;
 // dE/dphi[ltrk]:

         dp_dpar[nvtx][lvele][6]  = 0.;
 
         lvele++;
 // dPx/dcot[ltrk]:

         dp_dpar[nvtx][lvele][3]  = 0;
 // dPy/dcot[ltrk]:

         dp_dpar[nvtx][lvele][4]  = 0;
 // dPz/dcot[ltrk]:

         dp_dpar[nvtx][lvele][5]  = pmom[ltrk].perp();
 // dE/dcot[ltrk]:

         dp_dpar[nvtx][lvele][6]  = 
          pmom[ltrk].perp2() * _ctvmq.par[ltrk][2] / pmom[ltrk].e();
       }
     }
   }
 // Now calculate the new error matrix

 
   // Extract the interesting bits from VMAT.

 
   HepMatrix vmat(ctvmft_dim,ctvmft_dim,0);
   for(int lpar = 0; lpar < ctvmft_dim; lpar++){
     int l = lpar%3;
     int lvele = 0;
     if(lpar < 3  * _ctvmq.nvertx){
       int lvtx = lpar/3;
       lvele = _ctvmq.voff[lvtx] + l;
     }
     else{
       int ltrk = (lpar - 3 * _ctvmq.nvertx)/3;
       lvele = _ctvmq.toff[ltrk] + l;
     }
     for(int kpar = 0; kpar < ctvmft_dim; kpar++){ 
       int k = kpar%3;
       int kvele = 0;
       if(kpar < 3  * _ctvmq.nvertx){
         int kvtx = kpar/3;
         kvele = _ctvmq.voff[kvtx] + k;
       }
       else{
         int ktrk = (kpar - 3 * _ctvmq.nvertx)/3;
         kvele = _ctvmq.toff[ktrk] + k;
       }
       vmat[kpar][lpar] = _ctvmfr.vmat[kvele][lvele];
     }
   }
   //  Just about done

   HepMatrix ans(7,7,0);
   HepMatrix answer[_maxvtx] = {ans, ans, ans};
   for(int nvtx = 0; nvtx < _ctvmq.nvertx; nvtx++){
     answer[nvtx] = (dp_dpar[nvtx].T() * vmat) * dp_dpar[nvtx];
   }
   errors = answer[0];
 }
 
 /*=======================================================================

                          VertexFit::vOff

   =======================================================================*/
 int VertexFit::vOff(vertexNumber jv) const
 { if(jv < VERTEX_1 || jv > _maxvtx) {return -999;}
   else {return _ctvmq.voff[jv-1];}
 }
 
 /*=======================================================================

                          VertexFit::tOff

   =======================================================================*/
 int VertexFit::tOff(const CdfTrack* trk) const
 { for(int kt=0; kt<_ctvmq.ntrack; ++kt) {
     if(trk->id().value()==_ctvmq.list[kt]) return _ctvmq.toff[kt];
   }
   return -999;
 }
 
 /*=======================================================================

                          VertexFit::pOff

   =======================================================================*/
 int VertexFit::pOff(int lp) const
 { if(lp < 1) {
     return -999;
   }
   else return _ctvmq.poff[lp-1];
 }
 
 /*=======================================================================

                          VertexFit::cOff

   =======================================================================*/
 int VertexFit::cOff(int lc) const
 { if(lc < 1) {
     return -999;
   }
   else return _ctvmq.coff[lc-1];
 }
 
 /*=======================================================================

                          VertexFit::mOff

   =======================================================================*/
 int VertexFit::mOff() const
 { return _ctvmq.moff; }
 
 /*=======================================================================

                          VertexFit::VMat

   =======================================================================*/
 double VertexFit::VMat(int i, int j) const
 { if(i <0 || j < 0) {return -999.;}
   else return _ctvmfr.vmat[i][j];
 }
 
 /*=======================================================================

                          VertexFit::allocateVertexNumber

   =======================================================================*/
 
 // Facilitates CandNode recursion.  CandNodes have no way of deciding

 // which vertex they are, and these trivial functions help them do that.

 VertexFit::vertexNumber VertexFit::allocateVertexNumber()
 {
   if ((_currentAllocatedVertexNumber < PRIMARY_VERTEX) ||
       (_currentAllocatedVertexNumber > _maxvtx)) {
     std::cout << "VertexFit::allocateVertexNumber: out of range!" << std::endl;
     return PRIMARY_VERTEX;
   }
   return vertexNumber(++_currentAllocatedVertexNumber);
 }
 
 void VertexFit::resetAllocatedVertexNumber()
 {
   _currentAllocatedVertexNumber = 0;
 }
 
 /*=======================================================================

                          VertexFit::restoreFromCommons()

   =======================================================================*/
 
 void VertexFit::restoreFromCommons() {
         stat=0;
         _ctvmq_com  = (CTVMQ*) ctvmq_address_();
         _ctvmfr_com = (CTVMFR*) ctvmfr_address_();
         _fiddle_com = (FIDDLE*) fiddle_address_();
         _trkprm_com = (TRKPRM*) trkprm_address_();
         _ctvmq  = *_ctvmq_com;
         _ctvmfr = *_ctvmfr_com;
         _fiddle = *_fiddle_com;
         _trkprm = *_trkprm_com;
 }
 
 
 /*=======================================================================

                          VertexFit::setTrackReferencePoint

   =======================================================================*/
 void VertexFit::setTrackReferencePoint( const Hep3Vector & ref )
 {
   // Set extrapolation flag

   _extrapolateTrackErrors = true; 
   // Keep new reference point

   _referencePoint = ref;
   // Keep track of old primary vertex in CDF coordinates

   _cdfPrimaryVertex = _primaryVertex;
   // Set new primary vertex relative to new reference point

   setPrimaryVertex( _cdfPrimaryVertex - _referencePoint);
 }
 
 /*=======================================================================

                          VertexFit::moveReferencePoint()

   =======================================================================*/
 
 void VertexFit::moveReferencePoint( HepVector & v,
                                     HepSymMatrix & m ) {
 
   // Move the reference point of the track to _referencePoint

   if ( !_extrapolateTrackErrors ) 
     return;
   
   TrackExtrapolator e;
 
   // Modifies v and m

   e.moveReferencePoint( v, m, _referencePoint );
 }
 
 
 #endif  /* VERTEXFIT_CC */
 
 

This page was automatically generated by the LXR engine. The LXR team