VertexFit/interface/MultiVertexFitterC.h

//--------------------------------------------------------------------------------------------------
// $Id: $
//
// MultiVertexFitterC class header file
//
// Compact version of MultiVertexFitter
//
//--------------------------------------------------------------------------------------------------

#ifndef MITCOMMON_VERTEXFIT_MULTIVERTEXFITTERC_H
#define MITCOMMON_VERTEXFIT_MULTIVERTEXFITTERC_H

#include <string>
#include <iostream>

#include <Rtypes.h>
#include <TVectorD.h>
#include <TMatrixDSym.h>
#include <CLHEP/Vector/LorentzVector.h>
#include <CLHEP/Vector/ThreeVector.h>
#include <CLHEP/Matrix/SymMatrix.h>
#include <CLHEP/Matrix/Vector.h>

#include "MitCommon/DataFormats/interface/Types.h"
#include "MitCommon/Ctvmft/interface/cdimensions.hh"

//-------------------------------------------------------------------------------------------------
// Fortran routines to get address of the start of the ctvmq and ctvmfr common blocks
//-------------------------------------------------------------------------------------------------
extern "C" {
   int cctvmq_address_ (void);
   int cctvmfr_address_(void);
   int cfiddle_address_(void);
   int ctrkprm_address_(void);
}

namespace mithep {
  class MultiVertexFitterC {

    public:
      //--------------------------------------------------------------------------------------------
      // Enumerations
      //--------------------------------------------------------------------------------------------
      enum vertexNumber { PRIMARY_VERTEX,VERTEX_1,VERTEX_2,VERTEX_3,VERTEX_4,VERTEX_5,VERTEX_6 };
      enum vertexIndex  { X_INDEX=0, Y_INDEX, Z_INDEX, P1_INDEX, P2_INDEX };
      enum trackIndex   { CURVATURE_INDEX=0, PHI_INDEX, COTTH_INDEX };

      //--------------------------------------------------------------------------------------------
      // *structors
      //--------------------------------------------------------------------------------------------
      MultiVertexFitterC();
      ~MultiVertexFitterC() {}

      //--------------------------------------------------------------------------------------------
      // Fundamental funtions
      //--------------------------------------------------------------------------------------------
      void init                 (double bfield = 3.8);
      void setChisqMax          (const float chisqmx);
      
      //--------------------------------------------------------------------------------------------
      // CMS parameter ordering for the vector/matrix, which is assumed here:
      //
      //   qoverp, cotTheta, phi0, d0, z0;
      //   mapping to CDF is therefore { 1, 0*, 4, 3, 2 } 
      //
      // Note that the radius of curvature (in cm) is then:
      //
      //   Rc = cos(theta) / [ 0.0029979.... * (q/p) * B ],
      //
      // where B is the magnetic field in Tesla and tht is the angle between the field and the
      // direction. With p * cos(theta) = pT it follows:
      //
      //   Rc = pT / [ 0.0029979.... * q * B ],
      //   fullCurvature = 1 / Rc = 0.0029979 * q * B / pT = - 0.0029979 * B / pT.
      //
      //--------------------------------------------------------------------------------------------

      //--------------------------------------------------------------------------------------------
      // Parameter ordering for the vector/matrix, assumed here:
      //
      //   cotTheta, curvature, z0, d0, phi0;
      //
      // And also: curvature = -0.5 * CurvConst * BFieldT / pT (strictly speaking: half curvature)
      //                        CurvConst = 0.0029979, BField in Tesla
      //
      // Incidentally this is the ordering used in the CDF experiment.
      //--------------------------------------------------------------------------------------------
      bool addTrack             (const HepVector &pars, const HepSymMatrix &cov, int trackid,
                                 float mass, vertexNumber jv);

      bool addTrack             (const TVectorD &pars, const TMatrixDSym &cov, int trackid,
                                 double mass, vertexNumber jv);


      bool vertexPoint_2d       (vertexNumber jv1, vertexNumber jv2);
      bool vertexPoint_3d       (vertexNumber jv1, vertexNumber jv2);
      bool vertexPoint_1track   (vertexNumber jv1, vertexNumber jv2);
      bool vertexPoint_0track   (vertexNumber jv1, vertexNumber jv2);
      bool conversion_2d        (vertexNumber jv);
      bool conversion_3d        (vertexNumber jv);
      bool massConstrain        (int ntrk, const int trkIds[], float mass);

      void setPrimaryVertex     (float xv, float yv, float zv);
      void setPrimaryVertex     (Hep3Vector pv);
      bool setPrimaryVertexError(const HepSymMatrix &xverr);
      void setPrimaryVertexError(const float xverr[3][3]);

      bool fit                  ();

      void print                (std::ostream& os) const;
      void print                () const;
      void printErr             (std::ostream& os) const;
      void printErr             () const;

      void restoreFromCommons   ();

      // -------------------------------------------------------------------------------------------
      // Turn ON this option only when fitting the primary with no additional constraints. 
      // The routineis 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 beamlineConstraint(float xb, float yb, HepSymMatrix berr,float xzbslope,float yzbslope);
      bool beamlineConstraint(Hep3Vector pv, HepSymMatrix berr, float xzbslope, float yzbslope);

      //--------------------------------------------------------------------------------------------
      // Accessors
      //--------------------------------------------------------------------------------------------

      void             setBField  (double bField) { _bField = bField; }
      double           bField     () const { return _bField; }
      void             setExcuse  ();
      void             setNoExcuse();

      std::string      expert     () const;                // name/email of MultiVertexFitter expert
      int              status     () const;                // return status of fit
      // overall fit quality paramters
      int              ndof       () const;                // number of degrees of freedom
      float            prob       () const;                // return probability of chi-square
      float            chisq      () const;                // return chi-square of fit
      float            chisq      (const int trkId) const;
      float            chisq_rphi () const;
      float            chisq_rphi (const int trkId) const;
      float            chisq_z    () const;
      float            chisq_z    (const int trkId) const;
      float            chisq_rphiz() const;
      float            chisq_rphiz(const int trkId) const;

      // return fit track four momentum
      //HepLorentzVector getTrackP4       (const int trkId) const;
      FourVector       getTrackP4       (const int trkId) const;

      //// return fit track parameters
      //Helix            getHelix         (const int trkId) const;

      // return fit mass and get error
      float            getMass          (int ntrk, const int trkIds[], float& dmass) const;

      // return decay length
      float            getDecayLength   (vertexNumber nv, vertexNumber mv, const Hep3Vector& dir,
                                         float& dlerr) const;
      float            getDecayLength   (vertexNumber nv, vertexNumber mv, const ThreeVector& dir,
                                         float& dlerr) const;                                      
      float            getZDecayLength  (vertexNumber nv, vertexNumber mv,
                                         const Hep3Vector& dir, float& dlerr) const;  
      float            getZDecayLength  (vertexNumber nv, vertexNumber mv,
                                         const ThreeVector& dir, float& dlerr) const;                                        
      float            getImpactPar     (vertexNumber prdV, vertexNumber dcyV,
                                         const Hep3Vector &v, float &dxyerr) const;     
      float            getImpactPar     (vertexNumber prdV, vertexNumber dcyV,
                                         const ThreeVector &v, float &dxyerr) const;                      
      float            get_dr           (vertexNumber nv, vertexNumber mv, float& drerr) const;
      float            get_dz           (vertexNumber nv, vertexNumber mv, float& dzerr) const;

      // return location of vertex
      Hep3Vector       getVertexHep     (vertexNumber nv) const;
      ThreeVector      getVertex        (vertexNumber nv) const;

      // return error matrix element.
      ThreeSymMatrix   getErrorMatrix   (vertexNumber nv) const;
      double           getErrorMatrixHep(int j, int k) const;
      HepSymMatrix     getErrorMatrixHep(vertexNumber nv) const;
      HepSymMatrix     getErrorMatrixHep(const int trkId) const;
      void             getPosMomErr     (HepMatrix& errors) const;
      int              vOff             (vertexNumber jv) const;
      int              tOff             (const int trkId) const;
      int              pOff             (int lp) const;
      int              cOff             (int lc) const;
      int              mOff             () const;
      double           VMat             (int i, int j) const;
      float            getPtError       (const int trkId) const;
      MultiVertexFitterC::vertexNumber
        allocateVertexNumber();
      void             resetAllocatedVertexNumber();

      // Accessors for getting information relative to ijk errors Get the error code from the three ijk
      // indexes into the argument variables
      void             getIJKErr(int& err0, int& err1, int& err2) const;
      // Return each error code from the three ijk indexes
      int              getIJKErr0() const;
      int              getIJKErr1() const;
      int              getIJKErr2() const;

      // Get the track-id of the track causing a fatal error as indicated by the corresponding ijk error
      int              getErrTrackId() const;

      // Set new track reference point
      void             setTrackReferencePoint(const ThreeVector &ref);

      //--------------------------------------------------------------------------------------------
      // Overload operators
      //--------------------------------------------------------------------------------------------
      friend std::ostream& operator << (std::ostream& os, const MultiVertexFitterC& vfit);

    protected:
      std::string      _expert;                  // string: name and email of expert
      int              _stat;                    // status returned from fit

      static const int _maxvtx = CCTVMFT_MAXVTX; // Maximum number of vertices
      static const int _maxmcn = CCTVMFT_MAXMCN; // Maximum number of mass constraints
      static const int _maxtrk = CCTVMFT_MAXTRK; // Maximum number of tracks
      static const int _maxitr = CCTVMFT_MAXITR; // Maximum number of iteration steps
      static const int _maxdim = 5*(CCTVMFT_MAXVTX+1)+3*CCTVMFT_MAXTRK+CCTVMFT_MAXMCN;

      // FIDDLE must have access to protected data like _maxvtx, etc. It must therefore be a
      // friend
      struct        FIDDLE;
      friend struct FIDDLE;
      struct FIDDLE {
          int    excuse;
          float  chisqmax;
      };

      // CTVMQ must have access to MultiVertexFitterC's (other) protected data like _maxvtx, etc. It
      // must therefore be a friend.
      struct        CTVMQ;
      friend struct CTVMQ;
      struct        CTVMQ {
          int    runnum;
          int    trgnum;
          int    iter;
          int    ntscut;
          int    nvertx;
          int    nmassc;
          int    ntrack;
          int    trkvtx[_maxvtx][_maxtrk];  // Logical in FORTRAN, but integer here to get the size
          int    trkmcn[_maxmcn][_maxtrk];  // Logical in FORTRAN, but integer here to get the size
          int    vtxpnt[2][_maxvtx];
          float  cmass [_maxmcn];
          int    cvtx  [_maxvtx];
          int    vtxvtx[_maxvtx][_maxvtx];  // Logical in FORTRAN, but integer here to get the size
          char   tkbank[_maxtrk][4];
          int    list  [_maxtrk];
          float  tmass [_maxtrk];
          int    matdim;
          int    tkerr [_maxtrk];
          int    ndof;
          float  chisqr[_maxitr+1];
          float  chit  [_maxtrk];
          float  chiv  [_maxvtx+1];
          float  chim  [_maxmcn];
          float  xyzpv0[3];
          float  exyzpv[3][3];
          float  xzslope;
          float  yzslope;
          float  xyzvrt[_maxvtx+1][3];
          float  dxyzpv[3];
          float  par   [_maxtrk][5];
          float  g     [_maxtrk][5][5];
          float  trkp4 [6][_maxtrk];
          float  vtxp4 [_maxvtx][4];
          float  mcnp4 [_maxmcn][4];
          float  dda   [8][_maxtrk];
          int    voff  [_maxvtx];
          int    toff  [_maxtrk];
          int    poff  [_maxvtx];
          int    coff  [_maxvtx];
          int    moff;
          float  par0  [_maxtrk][5];
          float  pardif[_maxtrk][5];
          float  fmcdif[_maxmcn];
          float  pcon  [2][_maxvtx];
          float  sang  [2][_maxvtx];
          float  drmax;
          float  rvmax;
          float  dzmax;
          float  trnmax;
          float  dsmin;
          int    ijkerr[3];
          float  pscale;
      };

      // CTVMFR must have access to MultiVertexFitterC's (other) protected data like _maxdim, etc. 
      // It must therefore be a friend.
      struct        CTVMFR;
      friend struct CTVMFR;
      struct        CTVMFR {
          double vmat[_maxdim+1][_maxdim];
      };

      // TRKPRM must have access to MultiVertexFitterC's (other) protected data like _maxtrk, etc.
      // It must therefore be a friend.
      struct        TRKPRM;
      friend struct TRKPRM;
      struct        TRKPRM {
          float  trhelix[_maxtrk][5];
          float  trem   [_maxtrk][5][5];
      };

      CTVMQ    _ctvmq;
      CTVMQ*   _ctvmq_com;
      CTVMFR   _ctvmfr;
      CTVMFR*  _ctvmfr_com;
      FIDDLE   _fiddle;
      FIDDLE*  _fiddle_com;
      TRKPRM   _trkprm;
      TRKPRM*  _trkprm_com;

      //--------------------------------------------------------------------------------------------
      // Private functions used by class
      //--------------------------------------------------------------------------------------------

    private:

      // Moves reference point of track parameters and errors to _referencePoint
      //void moveReferencePoint(HepVector &v, HepSymMatrix &m);
      // Moves reference point of track parameters and errors to _referencePoint.  Additionally 
      //checks if track has already been moved by examining it's "derived" link
      //void moveReferencePoint(const CdfTrack *trk, HepVector &v, HepSymMatrix &m);

      //--------------------------------------------------------------------------------------------
      // Data members of class
      //--------------------------------------------------------------------------------------------
      double      _bField;                        // B field in Tesla

      int         _currentAllocatedVertexNumber;  // index to enum vertexNumber
      ThreeVector _referencePoint;                // reference point of track
      Hep3Vector  _primaryVertex;                 // primary vertex relative to _referencePoint
      Hep3Vector  _cdfPrimaryVertex;              // primary vertex in CDF coordinate system
      bool        _extrapolateTrackErrors;        // extrapolate track errors to _referencePoint
  };
}

//--------------------------------------------------------------------------------------------------
// Inline functions
//--------------------------------------------------------------------------------------------------
inline
std::ostream& operator << (std::ostream& os, const mithep::MultiVertexFitterC& vfit)
{
  vfit.print(os);
  return    (os);
}
inline
void mithep::MultiVertexFitterC::setExcuse()
{
  _fiddle.excuse = 1;   // the default
}
inline
void mithep::MultiVertexFitterC::setNoExcuse()
{
  _fiddle.excuse = 0;   // crash on input error
}
inline
void mithep::MultiVertexFitterC::setChisqMax(const float chisqmx)
{
  _fiddle.chisqmax = chisqmx;
}
#endif
Revision:	1.1
Committed:	Thu Nov 13 16:34:28 2008 UTC (16 years, 5 months ago) by paus
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_008pre2, Mit_008pre1, Mit_006b, Mit_006a, Mit_006
Log Message:	Submitting the ugly three versions. But no elegant solution found.
#	User	Rev	Content
1	paus	1.1	//--------------------------------------------------------------------------------------------------
2			// $Id: $
3			//
4			// MultiVertexFitterC class header file
5			//
6			// Compact version of MultiVertexFitter
7			//
8			//--------------------------------------------------------------------------------------------------
9
10			#ifndef MITCOMMON_VERTEXFIT_MULTIVERTEXFITTERC_H
11			#define MITCOMMON_VERTEXFIT_MULTIVERTEXFITTERC_H
12
13			#include <string>
14			#include <iostream>
15
16			#include <Rtypes.h>
17			#include <TVectorD.h>
18			#include <TMatrixDSym.h>
19			#include <CLHEP/Vector/LorentzVector.h>
20			#include <CLHEP/Vector/ThreeVector.h>
21			#include <CLHEP/Matrix/SymMatrix.h>
22			#include <CLHEP/Matrix/Vector.h>
23
24			#include "MitCommon/DataFormats/interface/Types.h"
25			#include "MitCommon/Ctvmft/interface/cdimensions.hh"
26
27			//-------------------------------------------------------------------------------------------------
28			// Fortran routines to get address of the start of the ctvmq and ctvmfr common blocks
29			//-------------------------------------------------------------------------------------------------
30			extern "C" {
31			int cctvmq_address_ (void);
32			int cctvmfr_address_(void);
33			int cfiddle_address_(void);
34			int ctrkprm_address_(void);
35			}
36
37			namespace mithep {
38			class MultiVertexFitterC {
39
40			public:
41			//--------------------------------------------------------------------------------------------
42			// Enumerations
43			//--------------------------------------------------------------------------------------------
44			enum vertexNumber { PRIMARY_VERTEX,VERTEX_1,VERTEX_2,VERTEX_3,VERTEX_4,VERTEX_5,VERTEX_6 };
45			enum vertexIndex { X_INDEX=0, Y_INDEX, Z_INDEX, P1_INDEX, P2_INDEX };
46			enum trackIndex { CURVATURE_INDEX=0, PHI_INDEX, COTTH_INDEX };
47
48			//--------------------------------------------------------------------------------------------
49			// *structors
50			//--------------------------------------------------------------------------------------------
51			MultiVertexFitterC();
52			~MultiVertexFitterC() {}
53
54			//--------------------------------------------------------------------------------------------
55			// Fundamental funtions
56			//--------------------------------------------------------------------------------------------
57			void init (double bfield = 3.8);
58			void setChisqMax (const float chisqmx);
59
60			//--------------------------------------------------------------------------------------------
61			// CMS parameter ordering for the vector/matrix, which is assumed here:
62			//
63			// qoverp, cotTheta, phi0, d0, z0;
64			// mapping to CDF is therefore { 1, 0*, 4, 3, 2 }
65			//
66			// Note that the radius of curvature (in cm) is then:
67			//
68			// Rc = cos(theta) / [ 0.0029979.... * (q/p) * B ],
69			//
70			// where B is the magnetic field in Tesla and tht is the angle between the field and the
71			// direction. With p * cos(theta) = pT it follows:
72			//
73			// Rc = pT / [ 0.0029979.... * q * B ],
74			// fullCurvature = 1 / Rc = 0.0029979 * q * B / pT = - 0.0029979 * B / pT.
75			//
76			//--------------------------------------------------------------------------------------------
77
78			//--------------------------------------------------------------------------------------------
79			// Parameter ordering for the vector/matrix, assumed here:
80			//
81			// cotTheta, curvature, z0, d0, phi0;
82			//
83			// And also: curvature = -0.5 * CurvConst * BFieldT / pT (strictly speaking: half curvature)
84			// CurvConst = 0.0029979, BField in Tesla
85			//
86			// Incidentally this is the ordering used in the CDF experiment.
87			//--------------------------------------------------------------------------------------------
88			bool addTrack (const HepVector &pars, const HepSymMatrix &cov, int trackid,
89			float mass, vertexNumber jv);
90
91			bool addTrack (const TVectorD &pars, const TMatrixDSym &cov, int trackid,
92			double mass, vertexNumber jv);
93
94
95			bool vertexPoint_2d (vertexNumber jv1, vertexNumber jv2);
96			bool vertexPoint_3d (vertexNumber jv1, vertexNumber jv2);
97			bool vertexPoint_1track (vertexNumber jv1, vertexNumber jv2);
98			bool vertexPoint_0track (vertexNumber jv1, vertexNumber jv2);
99			bool conversion_2d (vertexNumber jv);
100			bool conversion_3d (vertexNumber jv);
101			bool massConstrain (int ntrk, const int trkIds[], float mass);
102
103			void setPrimaryVertex (float xv, float yv, float zv);
104			void setPrimaryVertex (Hep3Vector pv);
105			bool setPrimaryVertexError(const HepSymMatrix &xverr);
106			void setPrimaryVertexError(const float xverr[3][3]);
107
108			bool fit ();
109
110			void print (std::ostream& os) const;
111			void print () const;
112			void printErr (std::ostream& os) const;
113			void printErr () const;
114
115			void restoreFromCommons ();
116
117			// -------------------------------------------------------------------------------------------
118			// Turn ON this option only when fitting the primary with no additional constraints.
119			// The routineis protected and will not function if those options are selected as well as the
120			// beamline constraint. In case this option is chosen the fit will calculate a result also
121			// with just one track in input.
122			//
123			// In order to enable this option the user must do the following: Fit the primary vertex as
124			// VERTEX_1 (without any other vertices). Enable beamlineConstraint by setting the pointing
125			// constraint variable vtxpnt[0][0] to -100 (no pointing constraints when <0) Note that index
126			// 0,0 corresponds to index 1,1 in fortran Provide the beamline parameters:
127			//
128			// --> assume xv = xb + xzbslope*z,
129			// yv = yb + yzbslope*z, where (xv, yv) is the transverse
130			// beam position at z, (xb, yb) is the transverse beam position at z = 0
131			// and (xzbslope, yzbslope) are the slopes of the beamline
132			// --> set primary vertex at z=0 (xb, yb, 0)
133			// --> set the diagonal elements of the primary vertex error matrix to the
134			// sigma**2 of beam spread in x, y and z:
135			// xverr[1][1] = (~30 - 50 um)**2
136			// xverr[2][2] = (~30 - 50 um)**2
137			// xverr[3][3] = (~30 cm) **2
138			// All other elements should be 0
139			// For help email: Joao Guimaraes guima@huhepl.harvard.edu, Franco Bedeschi bed@fnal.gov
140			// See more information in the ctvmft.f
141			// -------------------------------------------------------------------------------------------
142			bool beamlineConstraint(float xb, float yb, HepSymMatrix berr,float xzbslope,float yzbslope);
143			bool beamlineConstraint(Hep3Vector pv, HepSymMatrix berr, float xzbslope, float yzbslope);
144
145			//--------------------------------------------------------------------------------------------
146			// Accessors
147			//--------------------------------------------------------------------------------------------
148
149			void setBField (double bField) { _bField = bField; }
150			double bField () const { return _bField; }
151			void setExcuse ();
152			void setNoExcuse();
153
154			std::string expert () const; // name/email of MultiVertexFitter expert
155			int status () const; // return status of fit
156			// overall fit quality paramters
157			int ndof () const; // number of degrees of freedom
158			float prob () const; // return probability of chi-square
159			float chisq () const; // return chi-square of fit
160			float chisq (const int trkId) const;
161			float chisq_rphi () const;
162			float chisq_rphi (const int trkId) const;
163			float chisq_z () const;
164			float chisq_z (const int trkId) const;
165			float chisq_rphiz() const;
166			float chisq_rphiz(const int trkId) const;
167
168			// return fit track four momentum
169			//HepLorentzVector getTrackP4 (const int trkId) const;
170			FourVector getTrackP4 (const int trkId) const;
171
172			//// return fit track parameters
173			//Helix getHelix (const int trkId) const;
174
175			// return fit mass and get error
176			float getMass (int ntrk, const int trkIds[], float& dmass) const;
177
178			// return decay length
179			float getDecayLength (vertexNumber nv, vertexNumber mv, const Hep3Vector& dir,
180			float& dlerr) const;
181			float getDecayLength (vertexNumber nv, vertexNumber mv, const ThreeVector& dir,
182			float& dlerr) const;
183			float getZDecayLength (vertexNumber nv, vertexNumber mv,
184			const Hep3Vector& dir, float& dlerr) const;
185			float getZDecayLength (vertexNumber nv, vertexNumber mv,
186			const ThreeVector& dir, float& dlerr) const;
187			float getImpactPar (vertexNumber prdV, vertexNumber dcyV,
188			const Hep3Vector &v, float &dxyerr) const;
189			float getImpactPar (vertexNumber prdV, vertexNumber dcyV,
190			const ThreeVector &v, float &dxyerr) const;
191			float get_dr (vertexNumber nv, vertexNumber mv, float& drerr) const;
192			float get_dz (vertexNumber nv, vertexNumber mv, float& dzerr) const;
193
194			// return location of vertex
195			Hep3Vector getVertexHep (vertexNumber nv) const;
196			ThreeVector getVertex (vertexNumber nv) const;
197
198			// return error matrix element.
199			ThreeSymMatrix getErrorMatrix (vertexNumber nv) const;
200			double getErrorMatrixHep(int j, int k) const;
201			HepSymMatrix getErrorMatrixHep(vertexNumber nv) const;
202			HepSymMatrix getErrorMatrixHep(const int trkId) const;
203			void getPosMomErr (HepMatrix& errors) const;
204			int vOff (vertexNumber jv) const;
205			int tOff (const int trkId) const;
206			int pOff (int lp) const;
207			int cOff (int lc) const;
208			int mOff () const;
209			double VMat (int i, int j) const;
210			float getPtError (const int trkId) const;
211			MultiVertexFitterC::vertexNumber
212			allocateVertexNumber();
213			void resetAllocatedVertexNumber();
214
215			// Accessors for getting information relative to ijk errors Get the error code from the three ijk
216			// indexes into the argument variables
217			void getIJKErr(int& err0, int& err1, int& err2) const;
218			// Return each error code from the three ijk indexes
219			int getIJKErr0() const;
220			int getIJKErr1() const;
221			int getIJKErr2() const;
222
223			// Get the track-id of the track causing a fatal error as indicated by the corresponding ijk error
224			int getErrTrackId() const;
225
226			// Set new track reference point
227			void setTrackReferencePoint(const ThreeVector &ref);
228
229			//--------------------------------------------------------------------------------------------
230			// Overload operators
231			//--------------------------------------------------------------------------------------------
232			friend std::ostream& operator << (std::ostream& os, const MultiVertexFitterC& vfit);
233
234			protected:
235			std::string _expert; // string: name and email of expert
236			int _stat; // status returned from fit
237
238			static const int _maxvtx = CCTVMFT_MAXVTX; // Maximum number of vertices
239			static const int _maxmcn = CCTVMFT_MAXMCN; // Maximum number of mass constraints
240			static const int _maxtrk = CCTVMFT_MAXTRK; // Maximum number of tracks
241			static const int _maxitr = CCTVMFT_MAXITR; // Maximum number of iteration steps
242			static const int _maxdim = 5(CCTVMFT_MAXVTX+1)+3CCTVMFT_MAXTRK+CCTVMFT_MAXMCN;
243
244			// FIDDLE must have access to protected data like _maxvtx, etc. It must therefore be a
245			// friend
246			struct FIDDLE;
247			friend struct FIDDLE;
248			struct FIDDLE {
249			int excuse;
250			float chisqmax;
251			};
252
253			// CTVMQ must have access to MultiVertexFitterC's (other) protected data like _maxvtx, etc. It
254			// must therefore be a friend.
255			struct CTVMQ;
256			friend struct CTVMQ;
257			struct CTVMQ {
258			int runnum;
259			int trgnum;
260			int iter;
261			int ntscut;
262			int nvertx;
263			int nmassc;
264			int ntrack;
265			int trkvtx[_maxvtx][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
266			int trkmcn[_maxmcn][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
267			int vtxpnt[2][_maxvtx];
268			float cmass [_maxmcn];
269			int cvtx [_maxvtx];
270			int vtxvtx[_maxvtx][_maxvtx]; // Logical in FORTRAN, but integer here to get the size
271			char tkbank[_maxtrk][4];
272			int list [_maxtrk];
273			float tmass [_maxtrk];
274			int matdim;
275			int tkerr [_maxtrk];
276			int ndof;
277			float chisqr[_maxitr+1];
278			float chit [_maxtrk];
279			float chiv [_maxvtx+1];
280			float chim [_maxmcn];
281			float xyzpv0[3];
282			float exyzpv[3][3];
283			float xzslope;
284			float yzslope;
285			float xyzvrt[_maxvtx+1][3];
286			float dxyzpv[3];
287			float par [_maxtrk][5];
288			float g [_maxtrk][5][5];
289			float trkp4 [6][_maxtrk];
290			float vtxp4 [_maxvtx][4];
291			float mcnp4 [_maxmcn][4];
292			float dda [8][_maxtrk];
293			int voff [_maxvtx];
294			int toff [_maxtrk];
295			int poff [_maxvtx];
296			int coff [_maxvtx];
297			int moff;
298			float par0 [_maxtrk][5];
299			float pardif[_maxtrk][5];
300			float fmcdif[_maxmcn];
301			float pcon [2][_maxvtx];
302			float sang [2][_maxvtx];
303			float drmax;
304			float rvmax;
305			float dzmax;
306			float trnmax;
307			float dsmin;
308			int ijkerr[3];
309			float pscale;
310			};
311
312			// CTVMFR must have access to MultiVertexFitterC's (other) protected data like _maxdim, etc.
313			// It must therefore be a friend.
314			struct CTVMFR;
315			friend struct CTVMFR;
316			struct CTVMFR {
317			double vmat[_maxdim+1][_maxdim];
318			};
319
320			// TRKPRM must have access to MultiVertexFitterC's (other) protected data like _maxtrk, etc.
321			// It must therefore be a friend.
322			struct TRKPRM;
323			friend struct TRKPRM;
324			struct TRKPRM {
325			float trhelix[_maxtrk][5];
326			float trem [_maxtrk][5][5];
327			};
328
329			CTVMQ _ctvmq;
330			CTVMQ* _ctvmq_com;
331			CTVMFR _ctvmfr;
332			CTVMFR* _ctvmfr_com;
333			FIDDLE _fiddle;
334			FIDDLE* _fiddle_com;
335			TRKPRM _trkprm;
336			TRKPRM* _trkprm_com;
337
338			//--------------------------------------------------------------------------------------------
339			// Private functions used by class
340			//--------------------------------------------------------------------------------------------
341
342			private:
343
344			// Moves reference point of track parameters and errors to _referencePoint
345			//void moveReferencePoint(HepVector &v, HepSymMatrix &m);
346			// Moves reference point of track parameters and errors to _referencePoint. Additionally
347			//checks if track has already been moved by examining it's "derived" link
348			//void moveReferencePoint(const CdfTrack *trk, HepVector &v, HepSymMatrix &m);
349
350			//--------------------------------------------------------------------------------------------
351			// Data members of class
352			//--------------------------------------------------------------------------------------------
353			double _bField; // B field in Tesla
354
355			int _currentAllocatedVertexNumber; // index to enum vertexNumber
356			ThreeVector _referencePoint; // reference point of track
357			Hep3Vector _primaryVertex; // primary vertex relative to _referencePoint
358			Hep3Vector _cdfPrimaryVertex; // primary vertex in CDF coordinate system
359			bool _extrapolateTrackErrors; // extrapolate track errors to _referencePoint
360			};
361			}
362
363			//--------------------------------------------------------------------------------------------------
364			// Inline functions
365			//--------------------------------------------------------------------------------------------------
366			inline
367			std::ostream& operator << (std::ostream& os, const mithep::MultiVertexFitterC& vfit)
368			{
369			vfit.print(os);
370			return (os);
371			}
372			inline
373			void mithep::MultiVertexFitterC::setExcuse()
374			{
375			_fiddle.excuse = 1; // the default
376			}
377			inline
378			void mithep::MultiVertexFitterC::setNoExcuse()
379			{
380			_fiddle.excuse = 0; // crash on input error
381			}
382			inline
383			void mithep::MultiVertexFitterC::setChisqMax(const float chisqmx)
384			{
385			_fiddle.chisqmax = chisqmx;
386			}
387			#endif