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.
#	Content
1	//--------------------------------------------------------------------------------------------------
2	// $Id: MultiVertexFitterC.h,v 1.2 2009/03/20 13:33:03 loizides Exp $
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 CLHEP::HepVector &pars, const CLHEP::HepSymMatrix &cov,
89	int trackid, float mass, vertexNumber jv);
90
91	bool addTrack (const TVectorD &pars, const TMatrixDSym &cov,
92	int trackid, double mass, vertexNumber jv);
93
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	void setPrimaryVertex (CLHEP::Hep3Vector pv);
104	bool setPrimaryVertexError(const CLHEP::HepSymMatrix &xverr);
105	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	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
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	float getDecayLength (vertexNumber nv, vertexNumber mv,
177	const CLHEP::Hep3Vector& dir,
178	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	const CLHEP::Hep3Vector& dir, float& dlerr) const;
183	float getZDecayLength (vertexNumber nv, vertexNumber mv,
184	const ThreeVector& dir, float& dlerr) const;
185	float getImpactPar (vertexNumber prdV, vertexNumber dcyV,
186	const CLHEP::Hep3Vector &v, float &dxyerr) const;
187	float getImpactPar (vertexNumber prdV, vertexNumber dcyV,
188	const ThreeVector &v, float &dxyerr) const;
189	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	CLHEP::Hep3Vector getVertexHep (vertexNumber nv) const;
194	ThreeVector getVertex (vertexNumber nv) const;
195
196	// return error matrix element.
197	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	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	allocateVertexNumber();
211	void resetAllocatedVertexNumber();
212
213	// Accessors for getting information relative to ijk errors.
214	// Get the error code from the three ijk indexes into the argument variables
215	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	// Get the track-id of the track causing a fatal error as indicated
222	// by the corresponding ijk error
223	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	// CTVMQ must have access to MultiVertexFitterC's (other) protected data like _maxvtx, etc.
253	// It must therefore be a friend.
254	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	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	};
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
357	//--------------------------------------------------------------------------------------------------
358	inline
359	void mithep::MultiVertexFitterC::setExcuse()
360	{
361	_fiddle.excuse = 1; // the default
362	}
363
364	//--------------------------------------------------------------------------------------------------
365	inline
366	void mithep::MultiVertexFitterC::setNoExcuse()
367	{
368	_fiddle.excuse = 0; // crash on input error
369	}
370
371	//--------------------------------------------------------------------------------------------------
372	inline
373	void mithep::MultiVertexFitterC::setChisqMax(const float chisqmx)
374	{
375	_fiddle.chisqmax = chisqmx;
376	}
377	#endif