ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitCommon/VertexFit/interface/MultiVertexFitter.h
Revision: 1.1
Committed: Wed Sep 17 04:01:51 2008 UTC (16 years, 7 months ago) by loizides
Content type: text/plain
Branch: MAIN
Log Message:
Moved MitVertex contents to MitCommon. MitVertex therefore is obsolute and should not be touched anymore.

File Contents

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