VertexFit/interface/MultiVertexFitterD.h

//--------------------------------------------------------------------------------------------------
// $Id: MultiVertexFitterD.h,v 1.2 2009/03/20 13:33:03 loizides Exp $
//
// MultiVertexFitterD class header file
//
// Compact version of MultiVertexFitter
//
//--------------------------------------------------------------------------------------------------

#ifndef MITCOMMON_VERTEXFIT_MULTIVERTEXFITTERD_H
#define MITCOMMON_VERTEXFIT_MULTIVERTEXFITTERD_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/ddimensions.hh"

//-------------------------------------------------------------------------------------------------
// Fortran routines to get address of the start of the ctvmq and ctvmfr common blocks
//-------------------------------------------------------------------------------------------------
extern "C" {
   int dctvmq_address_ (void);
   int dctvmfr_address_(void);
   int dfiddle_address_(void);
   int dtrkprm_address_(void);
}

namespace mithep {
  class MultiVertexFitterD {

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

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

      //--------------------------------------------------------------------------------------------
      // 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;
      MultiVertexFitterD::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 MultiVertexFitterD& vfit);

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

      static const int _maxvtx = DCTVMFT_MAXVTX; // Maximum number of vertices
      static const int _maxmcn = DCTVMFT_MAXMCN; // Maximum number of mass constraints
      static const int _maxtrk = DCTVMFT_MAXTRK; // Maximum number of tracks
      static const int _maxitr = DCTVMFT_MAXITR; // Maximum number of iteration steps
      static const int _maxdim = 5*(DCTVMFT_MAXVTX+1)+3*DCTVMFT_MAXTRK+DCTVMFT_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 MultiVertexFitterD'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 MultiVertexFitterD'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 MultiVertexFitterD'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:
      //--------------------------------------------------------------------------------------------
      // Data members of class
      //--------------------------------------------------------------------------------------------
      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::MultiVertexFitterD& vfit)
{
  vfit.print(os);
  return    (os);
}

//--------------------------------------------------------------------------------------------------
inline
void mithep::MultiVertexFitterD::setExcuse()
{
  _fiddle.excuse = 1;   // the default
}

//--------------------------------------------------------------------------------------------------
inline
void mithep::MultiVertexFitterD::setNoExcuse()
{
  _fiddle.excuse = 0;   // crash on input error
}

//--------------------------------------------------------------------------------------------------
inline
void mithep::MultiVertexFitterD::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:	+33 -46 lines
Log Message:	Use CLHEP namespace.
#	User	Rev	Content
1	paus	1.1	//--------------------------------------------------------------------------------------------------
2	loizides	1.3	// $Id: MultiVertexFitterD.h,v 1.2 2009/03/20 13:33:03 loizides Exp $
3	paus	1.1	//
4			// MultiVertexFitterD class header file
5			//
6			// Compact version of MultiVertexFitter
7			//
8			//--------------------------------------------------------------------------------------------------
9
10			#ifndef MITCOMMON_VERTEXFIT_MULTIVERTEXFITTERD_H
11			#define MITCOMMON_VERTEXFIT_MULTIVERTEXFITTERD_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/ddimensions.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 dctvmq_address_ (void);
32			int dctvmfr_address_(void);
33			int dfiddle_address_(void);
34			int dtrkprm_address_(void);
35			}
36
37			namespace mithep {
38			class MultiVertexFitterD {
39
40			public:
41			//--------------------------------------------------------------------------------------------
42			// Enumerations
43			//--------------------------------------------------------------------------------------------
44			enum vertexNumber { PRIMARY_VERTEX,VERTEX_1 };
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			MultiVertexFitterD();
52			~MultiVertexFitterD() {}
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			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	loizides	1.3	void setPrimaryVertex (CLHEP::Hep3Vector pv);
105			bool setPrimaryVertexError(const CLHEP::HepSymMatrix &xverr);
106	paus	1.1	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	loizides	1.3	bool beamlineConstraint(float xb, float yb, CLHEP::HepSymMatrix berr,
143			float xzbslope,float yzbslope);
144			bool beamlineConstraint(CLHEP::Hep3Vector pv, CLHEP::HepSymMatrix berr,
145			float xzbslope, float yzbslope);
146	paus	1.1
147			//--------------------------------------------------------------------------------------------
148			// Accessors
149			//--------------------------------------------------------------------------------------------
150
151			void setBField (double bField) { _bField = bField; }
152			double bField () const { return _bField; }
153			void setExcuse ();
154			void setNoExcuse();
155
156	loizides	1.3	std::string expert () const; // name/email of MultiVertexFitter expert
157			int status () const; // return status of fit
158	paus	1.1	// overall fit quality paramters
159	loizides	1.3	int ndof () const; // number of degrees of freedom
160			float prob () const; // return probability of chi-square
161			float chisq () const; // return chi-square of fit
162	paus	1.1	float chisq (const int trkId) const;
163			float chisq_rphi () const;
164			float chisq_rphi (const int trkId) const;
165			float chisq_z () const;
166			float chisq_z (const int trkId) const;
167			float chisq_rphiz() const;
168			float chisq_rphiz(const int trkId) const;
169
170			// return fit track four momentum
171			FourVector getTrackP4 (const int trkId) const;
172
173			// return fit mass and get error
174			float getMass (int ntrk, const int trkIds[], float& dmass) const;
175
176			// return decay length
177	loizides	1.3	float getDecayLength (vertexNumber nv, vertexNumber mv,
178			const CLHEP::Hep3Vector& dir,
179	paus	1.1	float& dlerr) const;
180			float getDecayLength (vertexNumber nv, vertexNumber mv, const ThreeVector& dir,
181			float& dlerr) const;
182			float getZDecayLength (vertexNumber nv, vertexNumber mv,
183	loizides	1.3	const CLHEP::Hep3Vector& dir, float& dlerr) const;
184	paus	1.1	float getZDecayLength (vertexNumber nv, vertexNumber mv,
185	loizides	1.2	const ThreeVector& dir, float& dlerr) const;
186	paus	1.1	float getImpactPar (vertexNumber prdV, vertexNumber dcyV,
187	loizides	1.3	const CLHEP::Hep3Vector &v, float &dxyerr) const;
188	paus	1.1	float getImpactPar (vertexNumber prdV, vertexNumber dcyV,
189	loizides	1.3	const ThreeVector &v, float &dxyerr) const;
190	paus	1.1	float get_dr (vertexNumber nv, vertexNumber mv, float& drerr) const;
191			float get_dz (vertexNumber nv, vertexNumber mv, float& dzerr) const;
192
193			// return location of vertex
194	loizides	1.3	CLHEP::Hep3Vector getVertexHep (vertexNumber nv) const;
195			ThreeVector getVertex (vertexNumber nv) const;
196	paus	1.1
197			// return error matrix element.
198	loizides	1.3	ThreeSymMatrix getErrorMatrix (vertexNumber nv) const;
199			double getErrorMatrixHep(int j, int k) const;
200			CLHEP::HepSymMatrix getErrorMatrixHep(vertexNumber nv) const;
201			CLHEP::HepSymMatrix getErrorMatrixHep(const int trkId) const;
202			void getPosMomErr (CLHEP::HepMatrix& errors) const;
203	paus	1.1	int vOff (vertexNumber jv) const;
204			int tOff (const int trkId) const;
205			int pOff (int lp) const;
206			int cOff (int lc) const;
207			int mOff () const;
208			double VMat (int i, int j) const;
209			float getPtError (const int trkId) const;
210			MultiVertexFitterD::vertexNumber
211	loizides	1.3	allocateVertexNumber();
212	paus	1.1	void resetAllocatedVertexNumber();
213
214	loizides	1.2	// Accessors for getting information relative to ijk errors.
215			// Get the error code from the three ijk indexes into the argument variables.
216	paus	1.1	void getIJKErr(int& err0, int& err1, int& err2) const;
217			// Return each error code from the three ijk indexes
218			int getIJKErr0() const;
219			int getIJKErr1() const;
220			int getIJKErr2() const;
221
222	loizides	1.2	// Get the track-id of the track causing a fatal error as indicated
223			// by the corresponding ijk error
224	paus	1.1	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 MultiVertexFitterD& 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 = DCTVMFT_MAXVTX; // Maximum number of vertices
239			static const int _maxmcn = DCTVMFT_MAXMCN; // Maximum number of mass constraints
240			static const int _maxtrk = DCTVMFT_MAXTRK; // Maximum number of tracks
241			static const int _maxitr = DCTVMFT_MAXITR; // Maximum number of iteration steps
242			static const int _maxdim = 5(DCTVMFT_MAXVTX+1)+3DCTVMFT_MAXTRK+DCTVMFT_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	loizides	1.2	// CTVMQ must have access to MultiVertexFitterD's (other) protected data like _maxvtx, etc.
254			// It must therefore be a friend.
255	paus	1.1	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 MultiVertexFitterD'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 MultiVertexFitterD'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			private:
339			//--------------------------------------------------------------------------------------------
340			// Data members of class
341			//--------------------------------------------------------------------------------------------
342	loizides	1.3	double _bField; // B field in Tesla
343			int _currentAllocatedVertexNumber; // index to enum vertexNumber
344			ThreeVector _referencePoint; // reference point of track
345			CLHEP::Hep3Vector _primaryVertex; // primary vertex relative to _referencePoint
346			CLHEP::Hep3Vector _cdfPrimaryVertex; // primary vertex in CDF coordinate system
347			bool _extrapolateTrackErrors; // extrapolate track errors _referencePoint
348	paus	1.1	};
349			}
350
351			//--------------------------------------------------------------------------------------------------
352			// Inline functions
353			//--------------------------------------------------------------------------------------------------
354			inline
355			std::ostream& operator << (std::ostream& os, const mithep::MultiVertexFitterD& vfit)
356			{
357			vfit.print(os);
358			return (os);
359			}
360	loizides	1.2
361			//--------------------------------------------------------------------------------------------------
362	paus	1.1	inline
363			void mithep::MultiVertexFitterD::setExcuse()
364			{
365			_fiddle.excuse = 1; // the default
366			}
367	loizides	1.2
368			//--------------------------------------------------------------------------------------------------
369	paus	1.1	inline
370			void mithep::MultiVertexFitterD::setNoExcuse()
371			{
372			_fiddle.excuse = 0; // crash on input error
373			}
374	loizides	1.2
375			//--------------------------------------------------------------------------------------------------
376	paus	1.1	inline
377			void mithep::MultiVertexFitterD::setChisqMax(const float chisqmx)
378			{
379			_fiddle.chisqmax = chisqmx;
380			}
381			#endif