ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitCommon/VertexFit/interface/MultiVertexFitter.h
Revision: 1.6
Committed: Mon Jul 20 03:12:22 2009 UTC (15 years, 9 months ago) by loizides
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_011, Mit_010a, Mit_010, Mit_008
Changes since 1.5: +3 -1 lines
Log Message:
Changes for docu.

File Contents

# User Rev Content
1 loizides 1.1 //--------------------------------------------------------------------------------------------------
2 loizides 1.6 // $Id: MultiVertexFitter.h,v 1.5 2009/03/20 13:33:03 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    
100     public:
101 paus 1.4
102 loizides 1.1 //--------------------------------------------------------------------------------------------
103     // Enumerations
104     //--------------------------------------------------------------------------------------------
105 paus 1.3 enum vertexNumber { PRIMARY_VERTEX,VERTEX_1 };
106 loizides 1.1 enum vertexIndex { X_INDEX=0, Y_INDEX, Z_INDEX, P1_INDEX, P2_INDEX };
107     enum trackIndex { CURVATURE_INDEX=0, PHI_INDEX, COTTH_INDEX };
108    
109     //--------------------------------------------------------------------------------------------
110     // *structors
111     //--------------------------------------------------------------------------------------------
112     MultiVertexFitter();
113     ~MultiVertexFitter() {}
114    
115     //--------------------------------------------------------------------------------------------
116     // Fundamental funtions
117     //--------------------------------------------------------------------------------------------
118     void init (double bfield = 3.8);
119     void setChisqMax (const float chisqmx);
120    
121     //--------------------------------------------------------------------------------------------
122     // CMS parameter ordering for the vector/matrix, which is assumed here:
123     //
124     // qoverp, cotTheta, phi0, d0, z0;
125     // mapping to CDF is therefore { 1, 0*, 4, 3, 2 }
126     //
127     // Note that the radius of curvature (in cm) is then:
128     //
129     // Rc = cos(theta) / [ 0.0029979.... * (q/p) * B ],
130     //
131     // where B is the magnetic field in Tesla and tht is the angle between the field and the
132     // direction. With p * cos(theta) = pT it follows:
133     //
134     // Rc = pT / [ 0.0029979.... * q * B ],
135     // fullCurvature = 1 / Rc = 0.0029979 * q * B / pT = - 0.0029979 * B / pT.
136     //
137     //--------------------------------------------------------------------------------------------
138    
139     //--------------------------------------------------------------------------------------------
140     // Parameter ordering for the vector/matrix, assumed here:
141     //
142     // cotTheta, curvature, z0, d0, phi0;
143     //
144     // And also: curvature = -0.5 * CurvConst * BFieldT / pT (strictly speaking: half curvature)
145     // CurvConst = 0.0029979, BField in Tesla
146     //
147     // Incidentally this is the ordering used in the CDF experiment.
148     //--------------------------------------------------------------------------------------------
149     bool addTrack (const HepVector &pars, const HepSymMatrix &cov, int trackid,
150     float mass, vertexNumber jv);
151    
152     bool addTrack (const TVectorD &pars, const TMatrixDSym &cov, int trackid,
153     double mass, vertexNumber jv);
154    
155    
156     bool vertexPoint_2d (vertexNumber jv1, vertexNumber jv2);
157     bool vertexPoint_3d (vertexNumber jv1, vertexNumber jv2);
158     bool vertexPoint_1track (vertexNumber jv1, vertexNumber jv2);
159     bool vertexPoint_0track (vertexNumber jv1, vertexNumber jv2);
160     bool conversion_2d (vertexNumber jv);
161     bool conversion_3d (vertexNumber jv);
162     bool massConstrain (int ntrk, const int trkIds[], float mass);
163    
164     void setPrimaryVertex (float xv, float yv, float zv);
165     void setPrimaryVertex (Hep3Vector pv);
166     bool setPrimaryVertexError(const HepSymMatrix &xverr);
167     void setPrimaryVertexError(const float xverr[3][3]);
168    
169     bool fit ();
170    
171     void print (std::ostream& os) const;
172     void print () const;
173     void printErr (std::ostream& os) const;
174     void printErr () const;
175    
176     void restoreFromCommons ();
177    
178     // -------------------------------------------------------------------------------------------
179     // Turn ON this option only when fitting the primary with no additional constraints.
180     // The routineis protected and will not function if those options are selected as well as the
181     // beamline constraint. In case this option is chosen the fit will calculate a result also
182     // with just one track in input.
183     //
184     // In order to enable this option the user must do the following: Fit the primary vertex as
185     // VERTEX_1 (without any other vertices). Enable beamlineConstraint by setting the pointing
186     // constraint variable vtxpnt[0][0] to -100 (no pointing constraints when <0) Note that index
187     // 0,0 corresponds to index 1,1 in fortran Provide the beamline parameters:
188     //
189     // --> assume xv = xb + xzbslope*z,
190     // yv = yb + yzbslope*z, where (xv, yv) is the transverse
191     // beam position at z, (xb, yb) is the transverse beam position at z = 0
192     // and (xzbslope, yzbslope) are the slopes of the beamline
193     // --> set primary vertex at z=0 (xb, yb, 0)
194     // --> set the diagonal elements of the primary vertex error matrix to the
195     // sigma**2 of beam spread in x, y and z:
196     // xverr[1][1] = (~30 - 50 um)**2
197     // xverr[2][2] = (~30 - 50 um)**2
198     // xverr[3][3] = (~30 cm) **2
199     // All other elements should be 0
200     // For help email: Joao Guimaraes guima@huhepl.harvard.edu, Franco Bedeschi bed@fnal.gov
201     // See more information in the ctvmft.f
202     // -------------------------------------------------------------------------------------------
203     bool beamlineConstraint(float xb, float yb, HepSymMatrix berr,float xzbslope,float yzbslope);
204     bool beamlineConstraint(Hep3Vector pv, HepSymMatrix berr, float xzbslope, float yzbslope);
205    
206     //--------------------------------------------------------------------------------------------
207     // Accessors
208     //--------------------------------------------------------------------------------------------
209    
210     void setBField (double bField) { _bField = bField; }
211     double bField () const { return _bField; }
212     void setExcuse ();
213     void setNoExcuse();
214    
215     std::string expert () const; // name/email of MultiVertexFitter expert
216     int status () const; // return status of fit
217     // overall fit quality paramters
218     int ndof () const; // number of degrees of freedom
219     float prob () const; // return probability of chi-square
220     float chisq () const; // return chi-square of fit
221     float chisq (const int trkId) const;
222     float chisq_rphi () const;
223     float chisq_rphi (const int trkId) const;
224     float chisq_z () const;
225     float chisq_z (const int trkId) const;
226     float chisq_rphiz() const;
227     float chisq_rphiz(const int trkId) const;
228    
229     // return fit track four momentum
230     //HepLorentzVector getTrackP4 (const int trkId) const;
231     FourVector getTrackP4 (const int trkId) const;
232    
233     //// return fit track parameters
234     //Helix getHelix (const int trkId) const;
235     // return fit mass and get error
236     float getMass (int ntrk, const int trkIds[], float& dmass) const;
237    
238     // return decay length
239     float getDecayLength(vertexNumber nv, vertexNumber mv, const Hep3Vector& dir,
240     float& dlerr) const;
241 bendavid 1.2 float getDecayLength(vertexNumber nv, vertexNumber mv, const ThreeVector& dir,
242     float& dlerr) const;
243     float getZDecayLength(vertexNumber nv, vertexNumber mv,
244     const Hep3Vector& dir, float& dlerr) const;
245     float getZDecayLength(vertexNumber nv, vertexNumber mv,
246 loizides 1.5 const ThreeVector& dir, float& dlerr) const;
247 bendavid 1.2 float getImpactPar(vertexNumber prdV, vertexNumber dcyV,
248     const Hep3Vector &v, float &dxyerr) const;
249     float getImpactPar(vertexNumber prdV, vertexNumber dcyV,
250 loizides 1.5 const ThreeVector &v, float &dxyerr) const;
251 loizides 1.1 float get_dr(vertexNumber nv, vertexNumber mv, float& drerr) const;
252     float get_dz(vertexNumber nv, vertexNumber mv, float& dzerr) const;
253     // return location of vertex
254     Hep3Vector getVertexHep(vertexNumber nv) const;
255     ThreeVector getVertex(vertexNumber nv) const;
256     // return error matrix element.
257     ThreeSymMatrix getErrorMatrix(vertexNumber nv) const;
258     double getErrorMatrixHep(int j, int k) const;
259     HepSymMatrix getErrorMatrixHep(vertexNumber nv) const;
260     HepSymMatrix getErrorMatrixHep(const int trkId) const;
261     void getPosMomErr (HepMatrix& errors) const;
262     int vOff (vertexNumber jv) const;
263     int tOff (const int trkId) const;
264     int pOff (int lp) const;
265     int cOff (int lc) const;
266     int mOff () const;
267     double VMat (int i, int j) const;
268     float getPtError (const int trkId) const;
269     MultiVertexFitter::vertexNumber
270     allocateVertexNumber();
271     void resetAllocatedVertexNumber();
272    
273 paus 1.4 // Accessors for getting information relative to ijk errors get the error code from the three
274     // ijk indexes into the argument variables
275 loizides 1.1 void getIJKErr(int& err0, int& err1, int& err2) const;
276     // Return each error code from the three ijk indexes
277     int getIJKErr0() const;
278     int getIJKErr1() const;
279     int getIJKErr2() const;
280    
281 paus 1.4 // Get track-id of the track causing a fatal error as indicated by the corresponding ijk error
282 loizides 1.1 int getErrTrackId() const;
283    
284     // Set new track reference point
285     void setTrackReferencePoint(const ThreeVector &ref);
286    
287     //--------------------------------------------------------------------------------------------
288     // Overload operators
289     //--------------------------------------------------------------------------------------------
290     friend std::ostream& operator << (std::ostream& os, const MultiVertexFitter& vfit);
291    
292     protected:
293 paus 1.4 std::string _expert; // string: name, email of MultiVertexFitter expert
294 loizides 1.1 int _stat; // status returned from fit.
295    
296 paus 1.4 static const int _maxvtx = CTVMFT_MAXVTX; // maximum number of vertices
297     static const int _maxmcn = CTVMFT_MAXMCN; // maximum number of mass constraints
298     static const int _maxtrk = CTVMFT_MAXTRK; // maximum number of tracks
299     static const int _maxitr = CTVMFT_MAXITR; // maximum number of iteration steps
300 loizides 1.1 static const int _maxdim = 5*(CTVMFT_MAXVTX+1)+3*CTVMFT_MAXTRK+CTVMFT_MAXMCN;
301    
302     // FIDDLE must have access to MultiVertexFitter's (other) protected data like _maxvtx, etc. It
303     // must therefore be a friend.
304     struct FIDDLE;
305     friend struct FIDDLE;
306     struct FIDDLE {
307     int excuse;
308     float chisqmax;
309     };
310    
311     // CTVMQ must have access to MultiVertexFitter's (other) protected data like _maxvtx, etc. It
312     // must therefore be a friend.
313     struct CTVMQ;
314     friend struct CTVMQ;
315     struct CTVMQ {
316     int runnum;
317     int trgnum;
318     int iter;
319     int ntscut;
320     int nvertx;
321     int nmassc;
322     int ntrack;
323     int trkvtx[_maxvtx][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
324     int trkmcn[_maxmcn][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
325     int vtxpnt[2][_maxvtx];
326     float cmass [_maxmcn];
327     int cvtx [_maxvtx];
328     int vtxvtx[_maxvtx][_maxvtx]; // Logical in FORTRAN, but integer here to get the size
329     char tkbank[_maxtrk][4];
330     int list [_maxtrk];
331     float tmass [_maxtrk];
332     int matdim;
333     int tkerr [_maxtrk];
334     int ndof;
335     float chisqr[_maxitr+1];
336     float chit [_maxtrk];
337     float chiv [_maxvtx+1];
338     float chim [_maxmcn];
339     float xyzpv0[3];
340     float exyzpv[3][3];
341     float xzslope;
342     float yzslope;
343     float xyzvrt[_maxvtx+1][3];
344     float dxyzpv[3];
345     float par [_maxtrk][5];
346     float g [_maxtrk][5][5];
347     float trkp4 [6][_maxtrk];
348     float vtxp4 [_maxvtx][4];
349     float mcnp4 [_maxmcn][4];
350     float dda [8][_maxtrk];
351     int voff [_maxvtx];
352     int toff [_maxtrk];
353     int poff [_maxvtx];
354     int coff [_maxvtx];
355     int moff;
356     float par0 [_maxtrk][5];
357     float pardif[_maxtrk][5];
358     float fmcdif[_maxmcn];
359     float pcon [2][_maxvtx];
360     float sang [2][_maxvtx];
361     float drmax;
362     float rvmax;
363     float dzmax;
364     float trnmax;
365     float dsmin;
366     int ijkerr[3];
367     float pscale;
368     };
369    
370     // CTVMFR must have access to MultiVertexFitter's (other) protected data like _maxdim, etc.
371     // It must therefore be a friend.
372     struct CTVMFR;
373     friend struct CTVMFR;
374     struct CTVMFR {
375     double vmat[_maxdim+1][_maxdim];
376     };
377    
378     // TRKPRM must have access to MultiVertexFitter's (other) protected data like _maxtrk, etc.
379     // It must therefore be a friend.
380     struct TRKPRM;
381     friend struct TRKPRM;
382     struct TRKPRM {
383     float trhelix[_maxtrk][5];
384     float trem [_maxtrk][5][5];
385     };
386    
387     CTVMQ _ctvmq;
388     CTVMQ* _ctvmq_com;
389     CTVMFR _ctvmfr;
390     CTVMFR* _ctvmfr_com;
391     FIDDLE _fiddle;
392     FIDDLE* _fiddle_com;
393     TRKPRM _trkprm;
394     TRKPRM* _trkprm_com;
395    
396     //--------------------------------------------------------------------------------------------
397     // Private functions used by class
398     //--------------------------------------------------------------------------------------------
399    
400     private:
401    
402     // Moves reference point of track parameters and errors to _referencePoint
403     //void moveReferencePoint(HepVector &v, HepSymMatrix &m);
404     // Moves reference point of track parameters and errors to _referencePoint. Additionally
405     //checks if track has already been moved by examining it's "derived" link
406     //void moveReferencePoint(const CdfTrack *trk, HepVector &v, HepSymMatrix &m);
407    
408     //--------------------------------------------------------------------------------------------
409     // Data members of class
410     //--------------------------------------------------------------------------------------------
411 paus 1.3 double _bField; // B field in Tesla
412 loizides 1.1
413 paus 1.3 int _currentAllocatedVertexNumber; // index to enum vertexNumber
414     ThreeVector _referencePoint; // reference point of track
415     Hep3Vector _primaryVertex; // primary vertex relative to _referencePoint
416     Hep3Vector _cdfPrimaryVertex; // primary vertex in CDF coordinate system
417     bool _extrapolateTrackErrors; // extrapolate track errors to _referencePoint
418 loizides 1.1 };
419     }
420    
421     //--------------------------------------------------------------------------------------------------
422     // Inline functions
423     //--------------------------------------------------------------------------------------------------
424     inline
425     std::ostream& operator << (std::ostream& os, const mithep::MultiVertexFitter& vfit)
426     {
427     vfit.print(os);
428     return (os);
429     }
430 loizides 1.5
431     //--------------------------------------------------------------------------------------------------
432 loizides 1.1 inline
433     void mithep::MultiVertexFitter::setExcuse()
434     {
435     _fiddle.excuse = 1; // the default
436     }
437 loizides 1.5
438     //--------------------------------------------------------------------------------------------------
439 loizides 1.1 inline
440     void mithep::MultiVertexFitter::setNoExcuse()
441     {
442     _fiddle.excuse = 0; // crash on input error
443     }
444 loizides 1.5
445     //--------------------------------------------------------------------------------------------------
446 loizides 1.1 inline
447     void mithep::MultiVertexFitter::setChisqMax(const float chisqmx)
448     {
449     _fiddle.chisqmax = chisqmx;
450     }
451     #endif