VertexFit/interface/MultiVertexFitterC.h

//--------------------------------------------------------------------------------------------------
// $Id: MultiVertexFitterC.h,v 1.2 2009/03/20 13:33:03 loizides Exp $
//
// 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 CLHEP::HepVector &pars, const CLHEP::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     (CLHEP::Hep3Vector pv);
      bool setPrimaryVertexError(const CLHEP::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, CLHEP::HepSymMatrix berr,
                              float xzbslope,float yzbslope);
      bool beamlineConstraint(CLHEP::Hep3Vector pv, CLHEP::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
      FourVector       getTrackP4       (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 CLHEP::Hep3Vector& dir,
                                         float& dlerr) const;
      float            getDecayLength   (vertexNumber nv, vertexNumber mv, const ThreeVector& dir,
                                         float& dlerr) const;                                      
      float            getZDecayLength  (vertexNumber nv, vertexNumber mv,
                                         const CLHEP::Hep3Vector& dir, float& dlerr) const;  
      float            getZDecayLength  (vertexNumber nv, vertexNumber mv,
                                         const ThreeVector& dir, float& dlerr) const;                
      float            getImpactPar     (vertexNumber prdV, vertexNumber dcyV,
                                         const CLHEP::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
      CLHEP::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;
      CLHEP::HepSymMatrix getErrorMatrixHep(vertexNumber nv) const;
      CLHEP::HepSymMatrix getErrorMatrixHep(const int trkId) const;
      void                getPosMomErr     (CLHEP::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:
      double            _bField;                       // B field in Tesla
      int               _currentAllocatedVertexNumber; // index to enum vertexNumber
      ThreeVector       _referencePoint;               // reference point of track
      CLHEP::Hep3Vector _primaryVertex;                // primary vertex relative to _referencePoint
      CLHEP::Hep3Vector _cdfPrimaryVertex;             // primary vertex in CDF coordinate system
      bool              _extrapolateTrackErrors;       // extrapolate track errors _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.3
Committed:	Mon Oct 12 21:39:21 2009 UTC (15 years, 6 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_032, Mit_031, Mit_025c_branch2, Mit_025c_branch1, Mit_030, Mit_029c, Mit_030_pre1, Mit_029a, Mit_029, Mit_029_pre1, Mit_028a, Mit_025c_branch0, Mit_028, Mit_027a, Mit_027, Mit_026, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, V07-05-00, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b, ConvRejection-10-06-09, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013, Mit_013pre1, Mit_012i, Mit_012g, Mit_012f, Mit_012e, Mit_012d, Mit_012c, Mit_012b, Mit_012a, Mit_012, Mit_011a, HEAD
Branch point for:	Mit_025c_branch
Changes since 1.2:	+29 -46 lines
Log Message:	Use CLHEP namespace.
#	User	Rev	Content
1	paus	1.1	//--------------------------------------------------------------------------------------------------
2	loizides	1.3	// $Id: MultiVertexFitterC.h,v 1.2 2009/03/20 13:33:03 loizides Exp $
3	paus	1.1	//
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	loizides	1.3	bool addTrack (const CLHEP::HepVector &pars, const CLHEP::HepSymMatrix &cov,
89			int trackid, float mass, vertexNumber jv);
90	paus	1.1
91	loizides	1.3	bool addTrack (const TVectorD &pars, const TMatrixDSym &cov,
92			int trackid, double mass, vertexNumber jv);
93	paus	1.1
94			bool vertexPoint_2d (vertexNumber jv1, vertexNumber jv2);
95			bool vertexPoint_3d (vertexNumber jv1, vertexNumber jv2);
96			bool vertexPoint_1track (vertexNumber jv1, vertexNumber jv2);
97			bool vertexPoint_0track (vertexNumber jv1, vertexNumber jv2);
98			bool conversion_2d (vertexNumber jv);
99			bool conversion_3d (vertexNumber jv);
100			bool massConstrain (int ntrk, const int trkIds[], float mass);
101
102			void setPrimaryVertex (float xv, float yv, float zv);
103	loizides	1.3	void setPrimaryVertex (CLHEP::Hep3Vector pv);
104			bool setPrimaryVertexError(const CLHEP::HepSymMatrix &xverr);
105	paus	1.1	void setPrimaryVertexError(const float xverr[3][3]);
106
107			bool fit ();
108
109			void print (std::ostream& os) const;
110			void print () const;
111			void printErr (std::ostream& os) const;
112			void printErr () const;
113
114			void restoreFromCommons ();
115
116			// -------------------------------------------------------------------------------------------
117			// Turn ON this option only when fitting the primary with no additional constraints.
118			// The routineis protected and will not function if those options are selected as well as the
119			// beamline constraint. In case this option is chosen the fit will calculate a result also
120			// with just one track in input.
121			//
122			// In order to enable this option the user must do the following: Fit the primary vertex as
123			// VERTEX_1 (without any other vertices). Enable beamlineConstraint by setting the pointing
124			// constraint variable vtxpnt[0][0] to -100 (no pointing constraints when <0) Note that index
125			// 0,0 corresponds to index 1,1 in fortran Provide the beamline parameters:
126			//
127			// --> assume xv = xb + xzbslope*z,
128			// yv = yb + yzbslope*z, where (xv, yv) is the transverse
129			// beam position at z, (xb, yb) is the transverse beam position at z = 0
130			// and (xzbslope, yzbslope) are the slopes of the beamline
131			// --> set primary vertex at z=0 (xb, yb, 0)
132			// --> set the diagonal elements of the primary vertex error matrix to the
133			// sigma**2 of beam spread in x, y and z:
134			// xverr[1][1] = (~30 - 50 um)**2
135			// xverr[2][2] = (~30 - 50 um)**2
136			// xverr[3][3] = (~30 cm) **2
137			// All other elements should be 0
138			// For help email: Joao Guimaraes guima@huhepl.harvard.edu, Franco Bedeschi bed@fnal.gov
139			// See more information in the ctvmft.f
140			// -------------------------------------------------------------------------------------------
141	loizides	1.3	bool beamlineConstraint(float xb, float yb, CLHEP::HepSymMatrix berr,
142			float xzbslope,float yzbslope);
143			bool beamlineConstraint(CLHEP::Hep3Vector pv, CLHEP::HepSymMatrix berr,
144			float xzbslope, float yzbslope);
145	paus	1.1
146			//--------------------------------------------------------------------------------------------
147			// Accessors
148			//--------------------------------------------------------------------------------------------
149
150			void setBField (double bField) { _bField = bField; }
151			double bField () const { return _bField; }
152			void setExcuse ();
153			void setNoExcuse();
154
155			std::string expert () const; // name/email of MultiVertexFitter expert
156			int status () const; // return status of fit
157			// overall fit quality paramters
158			int ndof () const; // number of degrees of freedom
159			float prob () const; // return probability of chi-square
160			float chisq () const; // return chi-square of fit
161			float chisq (const int trkId) const;
162			float chisq_rphi () const;
163			float chisq_rphi (const int trkId) const;
164			float chisq_z () const;
165			float chisq_z (const int trkId) const;
166			float chisq_rphiz() const;
167			float chisq_rphiz(const int trkId) const;
168
169			// return fit track four momentum
170			FourVector getTrackP4 (const int trkId) const;
171
172			// return fit mass and get error
173			float getMass (int ntrk, const int trkIds[], float& dmass) const;
174
175			// return decay length
176	loizides	1.3	float getDecayLength (vertexNumber nv, vertexNumber mv,
177			const CLHEP::Hep3Vector& dir,
178	paus	1.1	float& dlerr) const;
179			float getDecayLength (vertexNumber nv, vertexNumber mv, const ThreeVector& dir,
180			float& dlerr) const;
181			float getZDecayLength (vertexNumber nv, vertexNumber mv,
182	loizides	1.3	const CLHEP::Hep3Vector& dir, float& dlerr) const;
183	paus	1.1	float getZDecayLength (vertexNumber nv, vertexNumber mv,
184	loizides	1.2	const ThreeVector& dir, float& dlerr) const;
185	paus	1.1	float getImpactPar (vertexNumber prdV, vertexNumber dcyV,
186	loizides	1.3	const CLHEP::Hep3Vector &v, float &dxyerr) const;
187	paus	1.1	float getImpactPar (vertexNumber prdV, vertexNumber dcyV,
188	loizides	1.2	const ThreeVector &v, float &dxyerr) const;
189	paus	1.1	float get_dr (vertexNumber nv, vertexNumber mv, float& drerr) const;
190			float get_dz (vertexNumber nv, vertexNumber mv, float& dzerr) const;
191
192			// return location of vertex
193	loizides	1.3	CLHEP::Hep3Vector getVertexHep (vertexNumber nv) const;
194			ThreeVector getVertex (vertexNumber nv) const;
195	paus	1.1
196			// return error matrix element.
197	loizides	1.3	ThreeSymMatrix getErrorMatrix (vertexNumber nv) const;
198			double getErrorMatrixHep(int j, int k) const;
199			CLHEP::HepSymMatrix getErrorMatrixHep(vertexNumber nv) const;
200			CLHEP::HepSymMatrix getErrorMatrixHep(const int trkId) const;
201			void getPosMomErr (CLHEP::HepMatrix& errors) const;
202	paus	1.1	int vOff (vertexNumber jv) const;
203			int tOff (const int trkId) const;
204			int pOff (int lp) const;
205			int cOff (int lc) const;
206			int mOff () const;
207			double VMat (int i, int j) const;
208			float getPtError (const int trkId) const;
209			MultiVertexFitterC::vertexNumber
210	loizides	1.3	allocateVertexNumber();
211	paus	1.1	void resetAllocatedVertexNumber();
212
213	loizides	1.2	// Accessors for getting information relative to ijk errors.
214			// Get the error code from the three ijk indexes into the argument variables
215	paus	1.1	void getIJKErr(int& err0, int& err1, int& err2) const;
216			// Return each error code from the three ijk indexes
217			int getIJKErr0() const;
218			int getIJKErr1() const;
219			int getIJKErr2() const;
220
221	loizides	1.2	// Get the track-id of the track causing a fatal error as indicated
222			// by the corresponding ijk error
223	paus	1.1	int getErrTrackId() const;
224
225			// Set new track reference point
226			void setTrackReferencePoint(const ThreeVector &ref);
227
228			//--------------------------------------------------------------------------------------------
229			// Overload operators
230			//--------------------------------------------------------------------------------------------
231			friend std::ostream& operator << (std::ostream& os, const MultiVertexFitterC& vfit);
232
233			protected:
234			std::string _expert; // string: name and email of expert
235			int _stat; // status returned from fit
236
237			static const int _maxvtx = CCTVMFT_MAXVTX; // Maximum number of vertices
238			static const int _maxmcn = CCTVMFT_MAXMCN; // Maximum number of mass constraints
239			static const int _maxtrk = CCTVMFT_MAXTRK; // Maximum number of tracks
240			static const int _maxitr = CCTVMFT_MAXITR; // Maximum number of iteration steps
241			static const int _maxdim = 5(CCTVMFT_MAXVTX+1)+3CCTVMFT_MAXTRK+CCTVMFT_MAXMCN;
242
243			// FIDDLE must have access to protected data like _maxvtx, etc. It must therefore be a
244			// friend
245			struct FIDDLE;
246			friend struct FIDDLE;
247			struct FIDDLE {
248			int excuse;
249			float chisqmax;
250			};
251
252	loizides	1.2	// CTVMQ must have access to MultiVertexFitterC's (other) protected data like _maxvtx, etc.
253			// It must therefore be a friend.
254	paus	1.1	struct CTVMQ;
255			friend struct CTVMQ;
256			struct CTVMQ {
257			int runnum;
258			int trgnum;
259			int iter;
260			int ntscut;
261			int nvertx;
262			int nmassc;
263			int ntrack;
264			int trkvtx[_maxvtx][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
265			int trkmcn[_maxmcn][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
266			int vtxpnt[2][_maxvtx];
267			float cmass [_maxmcn];
268			int cvtx [_maxvtx];
269			int vtxvtx[_maxvtx][_maxvtx]; // Logical in FORTRAN, but integer here to get the size
270			char tkbank[_maxtrk][4];
271			int list [_maxtrk];
272			float tmass [_maxtrk];
273			int matdim;
274			int tkerr [_maxtrk];
275			int ndof;
276			float chisqr[_maxitr+1];
277			float chit [_maxtrk];
278			float chiv [_maxvtx+1];
279			float chim [_maxmcn];
280			float xyzpv0[3];
281			float exyzpv[3][3];
282			float xzslope;
283			float yzslope;
284			float xyzvrt[_maxvtx+1][3];
285			float dxyzpv[3];
286			float par [_maxtrk][5];
287			float g [_maxtrk][5][5];
288			float trkp4 [6][_maxtrk];
289			float vtxp4 [_maxvtx][4];
290			float mcnp4 [_maxmcn][4];
291			float dda [8][_maxtrk];
292			int voff [_maxvtx];
293			int toff [_maxtrk];
294			int poff [_maxvtx];
295			int coff [_maxvtx];
296			int moff;
297			float par0 [_maxtrk][5];
298			float pardif[_maxtrk][5];
299			float fmcdif[_maxmcn];
300			float pcon [2][_maxvtx];
301			float sang [2][_maxvtx];
302			float drmax;
303			float rvmax;
304			float dzmax;
305			float trnmax;
306			float dsmin;
307			int ijkerr[3];
308			float pscale;
309			};
310
311			// CTVMFR must have access to MultiVertexFitterC's (other) protected data like _maxdim, etc.
312			// It must therefore be a friend.
313			struct CTVMFR;
314			friend struct CTVMFR;
315			struct CTVMFR {
316			double vmat[_maxdim+1][_maxdim];
317			};
318
319			// TRKPRM must have access to MultiVertexFitterC's (other) protected data like _maxtrk, etc.
320			// It must therefore be a friend.
321			struct TRKPRM;
322			friend struct TRKPRM;
323			struct TRKPRM {
324			float trhelix[_maxtrk][5];
325			float trem [_maxtrk][5][5];
326			};
327
328			CTVMQ _ctvmq;
329			CTVMQ* _ctvmq_com;
330			CTVMFR _ctvmfr;
331			CTVMFR* _ctvmfr_com;
332			FIDDLE _fiddle;
333			FIDDLE* _fiddle_com;
334			TRKPRM _trkprm;
335			TRKPRM* _trkprm_com;
336
337			private:
338	loizides	1.3	double _bField; // B field in Tesla
339			int _currentAllocatedVertexNumber; // index to enum vertexNumber
340			ThreeVector _referencePoint; // reference point of track
341			CLHEP::Hep3Vector _primaryVertex; // primary vertex relative to _referencePoint
342			CLHEP::Hep3Vector _cdfPrimaryVertex; // primary vertex in CDF coordinate system
343			bool _extrapolateTrackErrors; // extrapolate track errors _referencePoint
344	paus	1.1	};
345			}
346
347			//--------------------------------------------------------------------------------------------------
348			// Inline functions
349			//--------------------------------------------------------------------------------------------------
350			inline
351			std::ostream& operator << (std::ostream& os, const mithep::MultiVertexFitterC& vfit)
352			{
353			vfit.print(os);
354			return (os);
355			}
356	loizides	1.2
357			//--------------------------------------------------------------------------------------------------
358	paus	1.1	inline
359			void mithep::MultiVertexFitterC::setExcuse()
360			{
361			_fiddle.excuse = 1; // the default
362			}
363	loizides	1.2
364			//--------------------------------------------------------------------------------------------------
365	paus	1.1	inline
366			void mithep::MultiVertexFitterC::setNoExcuse()
367			{
368			_fiddle.excuse = 0; // crash on input error
369			}
370	loizides	1.2
371			//--------------------------------------------------------------------------------------------------
372	paus	1.1	inline
373			void mithep::MultiVertexFitterC::setChisqMax(const float chisqmx)
374			{
375			_fiddle.chisqmax = chisqmx;
376			}
377			#endif