VertexFit/interface/MultiVertexFitter.h

//--------------------------------------------------------------------------------------------------
// $Id: MultiVertexFitter.h,v 1.6 2009/07/20 03:12:22 loizides Exp $
//
// MultiVertexFitter class header file
//
// Author: Craig Blocker, Brandeis CDF Group
//         Christoph Paus, MIT
//           + porting to CMS
//           + major cleanup
//           + removed track dependency
//           - removed CLHEP dependence, now using root for vector/matrix operations
//           - re-enbale the track extrapolation
//
// Description:  Does generic wrapping of the multiple vertex fitter (wrapper of CTVMFT).
//
// Revision History:
//   May 27, 1999  Craig Blocker     Creation
//   Jun 22, 2001  Craig Blocker     Fix a few bugs in checking on vertices (such as pointing to
//                                   primary vertex wasn't allowed)
//   May 07 2002   Mat Martin        Added support for "zero track" vertices: Pointing multi track
//                 Jonas Rademacker  vertices at another vertex which has no tracks directly      
//                                   associated with it.  This is acheived through the routine:   
//                                   vertexPoint_0track. See the implementation in                
//                                   MultiVertexFitter.cc                                         
//   Sep 07, 2002 Joao Guimaraes     Added accessors for ijk errors and track-id of tracks that
//                                   produce some of those errors. These methods are:
//                                        getIJKErr(int&,int&,int&) const
//                                        getIJKErr0() const
//                                        getIJKErr1() const
//                                        getIJKErr2() const
//                                        getErrTrackId() const
//   Sep 15, 2002 Joao Guimaraes     Increased maximum number of tracks to 50 (_maxtrk).
//   Nov 11, 2002 Craig Blocker      Put protection in Fortran for potential bombs [divide by zero,
//                                   sqrt(neg. number]. Added particle specific addTrack routines
//                                   (addElectron, addPion, etc.) Added methods to directly access
//                                   the ctvmft error matrix and pointers.
//
//   Mar 18, 2003 Dmitry Litvintsev  added ability to restore the state of MultiVertexFitter
//                                   object from FORTRAN common blocks function is
//                                   'restoreFromCommons()' use it with care, fortran common blocks
//                                   have to be properly initialised. Contact me if unsure how to
//                                   use it. This routine is not needed for standard use of
//                                   MultiVertexFitter.
//   Oct 15, 2003 Sal Rappoccio      Added accessor to add tracks to ctvmft via a track id and the
//                                   track parameters and covariance
//   Dec 02, 2003 Sal Rappoccio      Added an interface to use extrapolated track errors in a
//                                   "transparent" way. The new interface is used as:
//
//                                     MultiVertexFitter fit;
//                                     fit.init();
//                                     fit.setPrimaryVertex(p1);    // Note: the order matters
//                                     fit.setReferencePoint(p2);
//                                     fit.addTrack(trk1, M_PION, 1);
//                                     fit.addTrack(trk2, M_PION, 1);
//                                     fit.fit();
//                                     Hep3Vector &pos = fit.getVertex(1);
//
//                                   "pos" will now hold the new vertex constructed from tracks
//                                   with errors referenced at the new reference point "p2". This
//                                   track extrapolation code is OFF by default. To turn it ON use:
//
//                                     fit.setReferencePoint(p2);
//
//                                   AddTrack was modified to implement this feature.
//   Jul 18, 2008  Christoph Paus    Works for CMS and any other type of Helix fitting.
//
//   (This class was taken from CDF and does therefore not follow our coding conventions.)
//--------------------------------------------------------------------------------------------------

#ifndef MITCOMMON_VERTEXFIT_MULTIVERTEXFITTER_H
#define MITCOMMON_VERTEXFIT_MULTIVERTEXFITTER_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/dimensions.hh"

//-------------------------------------------------------------------------------------------------
// Fortran routines to get address of the start of the ctvmq and ctvmfr common blocks
//-------------------------------------------------------------------------------------------------
extern "C" {
   int ctvmq_address_ (void);
   int ctvmfr_address_(void);
   int fiddle_address_(void);
   int trkprm_address_(void);
}

namespace mithep {
  class MultiVertexFitter {
    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
      //--------------------------------------------------------------------------------------------
      MultiVertexFitter();
      ~MultiVertexFitter() {}

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

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

      static const int _maxvtx = CTVMFT_MAXVTX; // maximum number of vertices
      static const int _maxmcn = CTVMFT_MAXMCN; // maximum number of mass constraints
      static const int _maxtrk = CTVMFT_MAXTRK; // maximum number of tracks
      static const int _maxitr = CTVMFT_MAXITR; // maximum number of iteration steps
      static const int _maxdim = 5*(CTVMFT_MAXVTX+1)+3*CTVMFT_MAXTRK+CTVMFT_MAXMCN;

      // FIDDLE must have access to MultiVertexFitter's (other) 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 MultiVertexFitter'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 MultiVertexFitter'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 MultiVertexFitter'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 rel. 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::MultiVertexFitter& vfit)
{
  vfit.print(os);
  return    (os);
}

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

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

//--------------------------------------------------------------------------------------------------
inline
void mithep::MultiVertexFitter::setChisqMax(const float chisqmx)
{
  _fiddle.chisqmax = chisqmx;
}
#endif
Revision:	1.7
Committed:	Mon Oct 12 21:39:20 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.6:	+71 -87 lines
Log Message:	Use CLHEP namespace.
#	User	Rev	Content
1	loizides	1.1	//--------------------------------------------------------------------------------------------------
2	loizides	1.7	// $Id: MultiVertexFitter.h,v 1.6 2009/07/20 03:12:22 loizides Exp $
3	loizides	1.1	//
4			// MultiVertexFitter class header file
5			//
6			// Author: Craig Blocker, Brandeis CDF Group
7			// Christoph Paus, MIT
8			// + porting to CMS
9			// + major cleanup
10			// + removed track dependency
11			// - removed CLHEP dependence, now using root for vector/matrix operations
12			// - re-enbale the track extrapolation
13			//
14			// Description: Does generic wrapping of the multiple vertex fitter (wrapper of CTVMFT).
15			//
16			// Revision History:
17			// May 27, 1999 Craig Blocker Creation
18			// Jun 22, 2001 Craig Blocker Fix a few bugs in checking on vertices (such as pointing to
19			// primary vertex wasn't allowed)
20			// May 07 2002 Mat Martin Added support for "zero track" vertices: Pointing multi track
21			// Jonas Rademacker vertices at another vertex which has no tracks directly
22			// associated with it. This is acheived through the routine:
23			// vertexPoint_0track. See the implementation in
24			// MultiVertexFitter.cc
25			// Sep 07, 2002 Joao Guimaraes Added accessors for ijk errors and track-id of tracks that
26			// produce some of those errors. These methods are:
27			// getIJKErr(int&,int&,int&) const
28			// getIJKErr0() const
29			// getIJKErr1() const
30			// getIJKErr2() const
31			// getErrTrackId() const
32			// Sep 15, 2002 Joao Guimaraes Increased maximum number of tracks to 50 (_maxtrk).
33			// Nov 11, 2002 Craig Blocker Put protection in Fortran for potential bombs [divide by zero,
34			// sqrt(neg. number]. Added particle specific addTrack routines
35			// (addElectron, addPion, etc.) Added methods to directly access
36			// the ctvmft error matrix and pointers.
37			//
38			// Mar 18, 2003 Dmitry Litvintsev added ability to restore the state of MultiVertexFitter
39			// object from FORTRAN common blocks function is
40			// 'restoreFromCommons()' use it with care, fortran common blocks
41			// have to be properly initialised. Contact me if unsure how to
42			// use it. This routine is not needed for standard use of
43			// MultiVertexFitter.
44			// Oct 15, 2003 Sal Rappoccio Added accessor to add tracks to ctvmft via a track id and the
45			// track parameters and covariance
46			// Dec 02, 2003 Sal Rappoccio Added an interface to use extrapolated track errors in a
47			// "transparent" way. The new interface is used as:
48			//
49			// MultiVertexFitter fit;
50			// fit.init();
51			// fit.setPrimaryVertex(p1); // Note: the order matters
52			// fit.setReferencePoint(p2);
53			// fit.addTrack(trk1, M_PION, 1);
54			// fit.addTrack(trk2, M_PION, 1);
55			// fit.fit();
56			// Hep3Vector &pos = fit.getVertex(1);
57			//
58			// "pos" will now hold the new vertex constructed from tracks
59			// with errors referenced at the new reference point "p2". This
60			// track extrapolation code is OFF by default. To turn it ON use:
61			//
62			// fit.setReferencePoint(p2);
63			//
64			// AddTrack was modified to implement this feature.
65			// Jul 18, 2008 Christoph Paus Works for CMS and any other type of Helix fitting.
66	loizides	1.6	//
67			// (This class was taken from CDF and does therefore not follow our coding conventions.)
68	loizides	1.1	//--------------------------------------------------------------------------------------------------
69
70			#ifndef MITCOMMON_VERTEXFIT_MULTIVERTEXFITTER_H
71			#define MITCOMMON_VERTEXFIT_MULTIVERTEXFITTER_H
72
73			#include <string>
74			#include <iostream>
75
76			#include <Rtypes.h>
77			#include <TVectorD.h>
78			#include <TMatrixDSym.h>
79			#include <CLHEP/Vector/LorentzVector.h>
80			#include <CLHEP/Vector/ThreeVector.h>
81			#include <CLHEP/Matrix/SymMatrix.h>
82			#include <CLHEP/Matrix/Vector.h>
83
84			#include "MitCommon/DataFormats/interface/Types.h"
85			#include "MitCommon/Ctvmft/interface/dimensions.hh"
86
87			//-------------------------------------------------------------------------------------------------
88			// Fortran routines to get address of the start of the ctvmq and ctvmfr common blocks
89			//-------------------------------------------------------------------------------------------------
90			extern "C" {
91			int ctvmq_address_ (void);
92			int ctvmfr_address_(void);
93			int fiddle_address_(void);
94			int trkprm_address_(void);
95			}
96
97			namespace mithep {
98			class MultiVertexFitter {
99			public:
100	paus	1.4
101	loizides	1.1	//--------------------------------------------------------------------------------------------
102			// Enumerations
103			//--------------------------------------------------------------------------------------------
104	paus	1.3	enum vertexNumber { PRIMARY_VERTEX,VERTEX_1 };
105	loizides	1.1	enum vertexIndex { X_INDEX=0, Y_INDEX, Z_INDEX, P1_INDEX, P2_INDEX };
106			enum trackIndex { CURVATURE_INDEX=0, PHI_INDEX, COTTH_INDEX };
107
108			//--------------------------------------------------------------------------------------------
109			// *structors
110			//--------------------------------------------------------------------------------------------
111			MultiVertexFitter();
112			~MultiVertexFitter() {}
113
114			//--------------------------------------------------------------------------------------------
115			// Fundamental funtions
116			//--------------------------------------------------------------------------------------------
117			void init (double bfield = 3.8);
118			void setChisqMax (const float chisqmx);
119
120			//--------------------------------------------------------------------------------------------
121			// CMS parameter ordering for the vector/matrix, which is assumed here:
122			//
123			// qoverp, cotTheta, phi0, d0, z0;
124			// mapping to CDF is therefore { 1, 0*, 4, 3, 2 }
125			//
126			// Note that the radius of curvature (in cm) is then:
127			//
128			// Rc = cos(theta) / [ 0.0029979.... * (q/p) * B ],
129			//
130			// where B is the magnetic field in Tesla and tht is the angle between the field and the
131			// direction. With p * cos(theta) = pT it follows:
132			//
133			// Rc = pT / [ 0.0029979.... * q * B ],
134			// fullCurvature = 1 / Rc = 0.0029979 * q * B / pT = - 0.0029979 * B / pT.
135			//
136			//--------------------------------------------------------------------------------------------
137
138			//--------------------------------------------------------------------------------------------
139			// Parameter ordering for the vector/matrix, assumed here:
140			//
141			// cotTheta, curvature, z0, d0, phi0;
142			//
143			// And also: curvature = -0.5 * CurvConst * BFieldT / pT (strictly speaking: half curvature)
144			// CurvConst = 0.0029979, BField in Tesla
145			//
146			// Incidentally this is the ordering used in the CDF experiment.
147			//--------------------------------------------------------------------------------------------
148	loizides	1.7	bool addTrack (const CLHEP::HepVector &pars, const CLHEP::HepSymMatrix &cov,
149			int trackid, float mass, vertexNumber jv);
150	loizides	1.1
151			bool addTrack (const TVectorD &pars, const TMatrixDSym &cov, int trackid,
152			double mass, vertexNumber jv);
153
154
155			bool vertexPoint_2d (vertexNumber jv1, vertexNumber jv2);
156			bool vertexPoint_3d (vertexNumber jv1, vertexNumber jv2);
157			bool vertexPoint_1track (vertexNumber jv1, vertexNumber jv2);
158			bool vertexPoint_0track (vertexNumber jv1, vertexNumber jv2);
159			bool conversion_2d (vertexNumber jv);
160			bool conversion_3d (vertexNumber jv);
161			bool massConstrain (int ntrk, const int trkIds[], float mass);
162
163			void setPrimaryVertex (float xv, float yv, float zv);
164	loizides	1.7	void setPrimaryVertex (CLHEP::Hep3Vector pv);
165			bool setPrimaryVertexError(const CLHEP::HepSymMatrix &xverr);
166	loizides	1.1	void setPrimaryVertexError(const float xverr[3][3]);
167
168			bool fit ();
169
170			void print (std::ostream& os) const;
171			void print () const;
172			void printErr (std::ostream& os) const;
173			void printErr () const;
174
175			void restoreFromCommons ();
176
177			// -------------------------------------------------------------------------------------------
178			// Turn ON this option only when fitting the primary with no additional constraints.
179			// The routineis protected and will not function if those options are selected as well as the
180			// beamline constraint. In case this option is chosen the fit will calculate a result also
181			// with just one track in input.
182			//
183			// In order to enable this option the user must do the following: Fit the primary vertex as
184			// VERTEX_1 (without any other vertices). Enable beamlineConstraint by setting the pointing
185			// constraint variable vtxpnt[0][0] to -100 (no pointing constraints when <0) Note that index
186			// 0,0 corresponds to index 1,1 in fortran Provide the beamline parameters:
187			//
188			// --> assume xv = xb + xzbslope*z,
189			// yv = yb + yzbslope*z, where (xv, yv) is the transverse
190			// beam position at z, (xb, yb) is the transverse beam position at z = 0
191			// and (xzbslope, yzbslope) are the slopes of the beamline
192			// --> set primary vertex at z=0 (xb, yb, 0)
193			// --> set the diagonal elements of the primary vertex error matrix to the
194			// sigma**2 of beam spread in x, y and z:
195			// xverr[1][1] = (~30 - 50 um)**2
196			// xverr[2][2] = (~30 - 50 um)**2
197			// xverr[3][3] = (~30 cm) **2
198			// All other elements should be 0
199			// For help email: Joao Guimaraes guima@huhepl.harvard.edu, Franco Bedeschi bed@fnal.gov
200			// See more information in the ctvmft.f
201			// -------------------------------------------------------------------------------------------
202	loizides	1.7	bool beamlineConstraint(float xb, float yb, CLHEP::HepSymMatrix berr,
203			float xzbslope,float yzbslope);
204			bool beamlineConstraint(CLHEP::Hep3Vector pv, CLHEP::HepSymMatrix berr,
205			float xzbslope, float yzbslope);
206	loizides	1.1
207			//--------------------------------------------------------------------------------------------
208			// Accessors
209			//--------------------------------------------------------------------------------------------
210
211	loizides	1.7	void setBField (double bField) { _bField = bField; }
212			double bField () const { return _bField; }
213			void setExcuse ();
214			void setNoExcuse();
215	loizides	1.1
216	loizides	1.7	std::string expert () const; // name/email of MultiVertexFitter expert
217			int status () const; // return status of fit
218	loizides	1.1	// overall fit quality paramters
219	loizides	1.7	int ndof () const; // number of degrees of freedom
220			float prob () const; // return probability of chi-square
221			float chisq () const; // return chi-square of fit
222			float chisq (const int trkId) const;
223			float chisq_rphi () const;
224			float chisq_rphi (const int trkId) const;
225			float chisq_z () const;
226			float chisq_z (const int trkId) const;
227			float chisq_rphiz() const;
228			float chisq_rphiz(const int trkId) const;
229	loizides	1.1
230			// return fit track four momentum
231	loizides	1.7	FourVector getTrackP4 (const int trkId) const;
232	loizides	1.1
233			// return fit mass and get error
234	loizides	1.7	float getMass (int ntrk, const int trkIds[], float& dmass) const;
235	loizides	1.1
236			// return decay length
237	loizides	1.7	float getDecayLength(vertexNumber nv, vertexNumber mv,
238			const CLHEP::Hep3Vector& dir, float& dlerr) const;
239			float getDecayLength(vertexNumber nv, vertexNumber mv, const ThreeVector& dir,
240			float& dlerr) const;
241			float getZDecayLength(vertexNumber nv, vertexNumber mv,
242			const CLHEP::Hep3Vector& dir, float& dlerr) const;
243			float getZDecayLength(vertexNumber nv, vertexNumber mv,
244			const ThreeVector& dir, float& dlerr) const;
245			float getImpactPar(vertexNumber prdV, vertexNumber dcyV,
246			const CLHEP::Hep3Vector &v, float &dxyerr) const;
247			float getImpactPar(vertexNumber prdV, vertexNumber dcyV,
248			const ThreeVector &v, float &dxyerr) const;
249			float get_dr(vertexNumber nv, vertexNumber mv, float& drerr) const;
250			float get_dz(vertexNumber nv, vertexNumber mv, float& dzerr) const;
251	loizides	1.1	// return location of vertex
252	loizides	1.7	CLHEP::Hep3Vector getVertexHep(vertexNumber nv) const;
253			ThreeVector getVertex(vertexNumber nv) const;
254			// return error matrix element
255			ThreeSymMatrix getErrorMatrix(vertexNumber nv) const;
256			double getErrorMatrixHep(int j, int k) const;
257			CLHEP::HepSymMatrix getErrorMatrixHep(vertexNumber nv) const;
258			CLHEP::HepSymMatrix getErrorMatrixHep(const int trkId) const;
259			void getPosMomErr (CLHEP::HepMatrix& errors) const;
260			int vOff (vertexNumber jv) const;
261			int tOff (const int trkId) const;
262			int pOff (int lp) const;
263			int cOff (int lc) const;
264			int mOff () const;
265			double VMat (int i, int j) const;
266			float getPtError (const int trkId) const;
267			MultiVertexFitter::vertexNumber
268			allocateVertexNumber();
269			void resetAllocatedVertexNumber();
270	loizides	1.1
271	paus	1.4	// Accessors for getting information relative to ijk errors get the error code from the three
272			// ijk indexes into the argument variables
273	loizides	1.7	void getIJKErr(int& err0, int& err1, int& err2) const;
274	loizides	1.1	// Return each error code from the three ijk indexes
275	loizides	1.7	int getIJKErr0() const;
276			int getIJKErr1() const;
277			int getIJKErr2() const;
278	loizides	1.1
279	paus	1.4	// Get track-id of the track causing a fatal error as indicated by the corresponding ijk error
280	loizides	1.7	int getErrTrackId() const;
281	loizides	1.1
282			// Set new track reference point
283	loizides	1.7	void setTrackReferencePoint(const ThreeVector &ref);
284	loizides	1.1
285			//--------------------------------------------------------------------------------------------
286			// Overload operators
287			//--------------------------------------------------------------------------------------------
288			friend std::ostream& operator << (std::ostream& os, const MultiVertexFitter& vfit);
289
290			protected:
291	paus	1.4	std::string _expert; // string: name, email of MultiVertexFitter expert
292	loizides	1.1	int _stat; // status returned from fit.
293
294	paus	1.4	static const int _maxvtx = CTVMFT_MAXVTX; // maximum number of vertices
295			static const int _maxmcn = CTVMFT_MAXMCN; // maximum number of mass constraints
296			static const int _maxtrk = CTVMFT_MAXTRK; // maximum number of tracks
297			static const int _maxitr = CTVMFT_MAXITR; // maximum number of iteration steps
298	loizides	1.1	static const int _maxdim = 5(CTVMFT_MAXVTX+1)+3CTVMFT_MAXTRK+CTVMFT_MAXMCN;
299
300			// FIDDLE must have access to MultiVertexFitter's (other) protected data like _maxvtx, etc. It
301			// must therefore be a friend.
302			struct FIDDLE;
303			friend struct FIDDLE;
304			struct FIDDLE {
305			int excuse;
306			float chisqmax;
307			};
308
309			// CTVMQ must have access to MultiVertexFitter's (other) protected data like _maxvtx, etc. It
310			// must therefore be a friend.
311			struct CTVMQ;
312			friend struct CTVMQ;
313			struct CTVMQ {
314			int runnum;
315			int trgnum;
316			int iter;
317			int ntscut;
318			int nvertx;
319			int nmassc;
320			int ntrack;
321			int trkvtx[_maxvtx][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
322			int trkmcn[_maxmcn][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
323			int vtxpnt[2][_maxvtx];
324			float cmass [_maxmcn];
325			int cvtx [_maxvtx];
326			int vtxvtx[_maxvtx][_maxvtx]; // Logical in FORTRAN, but integer here to get the size
327			char tkbank[_maxtrk][4];
328			int list [_maxtrk];
329			float tmass [_maxtrk];
330			int matdim;
331			int tkerr [_maxtrk];
332			int ndof;
333			float chisqr[_maxitr+1];
334			float chit [_maxtrk];
335			float chiv [_maxvtx+1];
336			float chim [_maxmcn];
337			float xyzpv0[3];
338			float exyzpv[3][3];
339			float xzslope;
340			float yzslope;
341			float xyzvrt[_maxvtx+1][3];
342			float dxyzpv[3];
343			float par [_maxtrk][5];
344			float g [_maxtrk][5][5];
345			float trkp4 [6][_maxtrk];
346			float vtxp4 [_maxvtx][4];
347			float mcnp4 [_maxmcn][4];
348			float dda [8][_maxtrk];
349			int voff [_maxvtx];
350			int toff [_maxtrk];
351			int poff [_maxvtx];
352			int coff [_maxvtx];
353			int moff;
354			float par0 [_maxtrk][5];
355			float pardif[_maxtrk][5];
356			float fmcdif[_maxmcn];
357			float pcon [2][_maxvtx];
358			float sang [2][_maxvtx];
359			float drmax;
360			float rvmax;
361			float dzmax;
362			float trnmax;
363			float dsmin;
364			int ijkerr[3];
365			float pscale;
366			};
367
368			// CTVMFR must have access to MultiVertexFitter's (other) protected data like _maxdim, etc.
369			// It must therefore be a friend.
370			struct CTVMFR;
371			friend struct CTVMFR;
372			struct CTVMFR {
373			double vmat[_maxdim+1][_maxdim];
374			};
375
376			// TRKPRM must have access to MultiVertexFitter's (other) protected data like _maxtrk, etc.
377			// It must therefore be a friend.
378			struct TRKPRM;
379			friend struct TRKPRM;
380			struct TRKPRM {
381			float trhelix[_maxtrk][5];
382			float trem [_maxtrk][5][5];
383			};
384
385			CTVMQ _ctvmq;
386			CTVMQ* _ctvmq_com;
387			CTVMFR _ctvmfr;
388			CTVMFR* _ctvmfr_com;
389			FIDDLE _fiddle;
390			FIDDLE* _fiddle_com;
391			TRKPRM _trkprm;
392			TRKPRM* _trkprm_com;
393
394			private:
395	loizides	1.7	double _bField; // B field in Tesla
396	loizides	1.1
397	loizides	1.7	int _currentAllocatedVertexNumber; // index to enum vertexNumber
398			ThreeVector _referencePoint; // reference point of track
399			CLHEP::Hep3Vector _primaryVertex; // primary vertex rel. to _referencePoint
400			CLHEP::Hep3Vector _cdfPrimaryVertex; // primary vertex in CDF coordinate system
401			bool _extrapolateTrackErrors; // extrapolate track errors _referencePoint
402	loizides	1.1	};
403			}
404
405			//--------------------------------------------------------------------------------------------------
406			// Inline functions
407			//--------------------------------------------------------------------------------------------------
408			inline
409			std::ostream& operator << (std::ostream& os, const mithep::MultiVertexFitter& vfit)
410			{
411			vfit.print(os);
412			return (os);
413			}
414	loizides	1.5
415			//--------------------------------------------------------------------------------------------------
416	loizides	1.1	inline
417			void mithep::MultiVertexFitter::setExcuse()
418			{
419			_fiddle.excuse = 1; // the default
420			}
421	loizides	1.5
422			//--------------------------------------------------------------------------------------------------
423	loizides	1.1	inline
424			void mithep::MultiVertexFitter::setNoExcuse()
425			{
426			_fiddle.excuse = 0; // crash on input error
427			}
428	loizides	1.5
429			//--------------------------------------------------------------------------------------------------
430	loizides	1.1	inline
431			void mithep::MultiVertexFitter::setChisqMax(const float chisqmx)
432			{
433			_fiddle.chisqmax = chisqmx;
434			}
435			#endif