ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitCommon/VertexFit/interface/MultiVertexFitter.h
Revision: 1.5
Committed: Fri Mar 20 13:33:03 2009 UTC (16 years, 1 month ago) by loizides
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_009c, Mit_009b, Mit_009a, Mit_009
Changes since 1.4: +9 -3 lines
Log Message:
Cleanup

File Contents

# User Rev Content
1 loizides 1.1 //--------------------------------------------------------------------------------------------------
2 loizides 1.5 // $Id: MultiVertexFitter.h,v 1.4 2008/11/13 16:54:07 paus 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     //--------------------------------------------------------------------------------------------------
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 paus 1.4
100 loizides 1.1 //--------------------------------------------------------------------------------------------
101     // Enumerations
102     //--------------------------------------------------------------------------------------------
103 paus 1.3 enum vertexNumber { PRIMARY_VERTEX,VERTEX_1 };
104 loizides 1.1 enum vertexIndex { X_INDEX=0, Y_INDEX, Z_INDEX, P1_INDEX, P2_INDEX };
105     enum trackIndex { CURVATURE_INDEX=0, PHI_INDEX, COTTH_INDEX };
106    
107     //--------------------------------------------------------------------------------------------
108     // *structors
109     //--------------------------------------------------------------------------------------------
110     MultiVertexFitter();
111     ~MultiVertexFitter() {}
112    
113     //--------------------------------------------------------------------------------------------
114     // Fundamental funtions
115     //--------------------------------------------------------------------------------------------
116     void init (double bfield = 3.8);
117     void setChisqMax (const float chisqmx);
118    
119     //--------------------------------------------------------------------------------------------
120     // CMS parameter ordering for the vector/matrix, which is assumed here:
121     //
122     // qoverp, cotTheta, phi0, d0, z0;
123     // mapping to CDF is therefore { 1, 0*, 4, 3, 2 }
124     //
125     // Note that the radius of curvature (in cm) is then:
126     //
127     // Rc = cos(theta) / [ 0.0029979.... * (q/p) * B ],
128     //
129     // where B is the magnetic field in Tesla and tht is the angle between the field and the
130     // direction. With p * cos(theta) = pT it follows:
131     //
132     // Rc = pT / [ 0.0029979.... * q * B ],
133     // fullCurvature = 1 / Rc = 0.0029979 * q * B / pT = - 0.0029979 * B / pT.
134     //
135     //--------------------------------------------------------------------------------------------
136    
137     //--------------------------------------------------------------------------------------------
138     // Parameter ordering for the vector/matrix, assumed here:
139     //
140     // cotTheta, curvature, z0, d0, phi0;
141     //
142     // And also: curvature = -0.5 * CurvConst * BFieldT / pT (strictly speaking: half curvature)
143     // CurvConst = 0.0029979, BField in Tesla
144     //
145     // Incidentally this is the ordering used in the CDF experiment.
146     //--------------------------------------------------------------------------------------------
147     bool addTrack (const HepVector &pars, const HepSymMatrix &cov, int trackid,
148     float mass, vertexNumber jv);
149    
150     bool addTrack (const TVectorD &pars, const TMatrixDSym &cov, int trackid,
151     double mass, vertexNumber jv);
152    
153    
154     bool vertexPoint_2d (vertexNumber jv1, vertexNumber jv2);
155     bool vertexPoint_3d (vertexNumber jv1, vertexNumber jv2);
156     bool vertexPoint_1track (vertexNumber jv1, vertexNumber jv2);
157     bool vertexPoint_0track (vertexNumber jv1, vertexNumber jv2);
158     bool conversion_2d (vertexNumber jv);
159     bool conversion_3d (vertexNumber jv);
160     bool massConstrain (int ntrk, const int trkIds[], float mass);
161    
162     void setPrimaryVertex (float xv, float yv, float zv);
163     void setPrimaryVertex (Hep3Vector pv);
164     bool setPrimaryVertexError(const HepSymMatrix &xverr);
165     void setPrimaryVertexError(const float xverr[3][3]);
166    
167     bool fit ();
168    
169     void print (std::ostream& os) const;
170     void print () const;
171     void printErr (std::ostream& os) const;
172     void printErr () const;
173    
174     void restoreFromCommons ();
175    
176     // -------------------------------------------------------------------------------------------
177     // Turn ON this option only when fitting the primary with no additional constraints.
178     // The routineis protected and will not function if those options are selected as well as the
179     // beamline constraint. In case this option is chosen the fit will calculate a result also
180     // with just one track in input.
181     //
182     // In order to enable this option the user must do the following: Fit the primary vertex as
183     // VERTEX_1 (without any other vertices). Enable beamlineConstraint by setting the pointing
184     // constraint variable vtxpnt[0][0] to -100 (no pointing constraints when <0) Note that index
185     // 0,0 corresponds to index 1,1 in fortran Provide the beamline parameters:
186     //
187     // --> assume xv = xb + xzbslope*z,
188     // yv = yb + yzbslope*z, where (xv, yv) is the transverse
189     // beam position at z, (xb, yb) is the transverse beam position at z = 0
190     // and (xzbslope, yzbslope) are the slopes of the beamline
191     // --> set primary vertex at z=0 (xb, yb, 0)
192     // --> set the diagonal elements of the primary vertex error matrix to the
193     // sigma**2 of beam spread in x, y and z:
194     // xverr[1][1] = (~30 - 50 um)**2
195     // xverr[2][2] = (~30 - 50 um)**2
196     // xverr[3][3] = (~30 cm) **2
197     // All other elements should be 0
198     // For help email: Joao Guimaraes guima@huhepl.harvard.edu, Franco Bedeschi bed@fnal.gov
199     // See more information in the ctvmft.f
200     // -------------------------------------------------------------------------------------------
201     bool beamlineConstraint(float xb, float yb, HepSymMatrix berr,float xzbslope,float yzbslope);
202     bool beamlineConstraint(Hep3Vector pv, HepSymMatrix berr, float xzbslope, float yzbslope);
203    
204     //--------------------------------------------------------------------------------------------
205     // Accessors
206     //--------------------------------------------------------------------------------------------
207    
208     void setBField (double bField) { _bField = bField; }
209     double bField () const { return _bField; }
210     void setExcuse ();
211     void setNoExcuse();
212    
213     std::string expert () const; // name/email of MultiVertexFitter expert
214     int status () const; // return status of fit
215     // overall fit quality paramters
216     int ndof () const; // number of degrees of freedom
217     float prob () const; // return probability of chi-square
218     float chisq () const; // return chi-square of fit
219     float chisq (const int trkId) const;
220     float chisq_rphi () const;
221     float chisq_rphi (const int trkId) const;
222     float chisq_z () const;
223     float chisq_z (const int trkId) const;
224     float chisq_rphiz() const;
225     float chisq_rphiz(const int trkId) const;
226    
227     // return fit track four momentum
228     //HepLorentzVector getTrackP4 (const int trkId) const;
229     FourVector getTrackP4 (const int trkId) const;
230    
231     //// return fit track parameters
232     //Helix getHelix (const int trkId) const;
233     // return fit mass and get error
234     float getMass (int ntrk, const int trkIds[], float& dmass) const;
235    
236     // return decay length
237     float getDecayLength(vertexNumber nv, vertexNumber mv, const Hep3Vector& dir,
238     float& dlerr) const;
239 bendavid 1.2 float getDecayLength(vertexNumber nv, vertexNumber mv, const ThreeVector& dir,
240     float& dlerr) const;
241     float getZDecayLength(vertexNumber nv, vertexNumber mv,
242     const Hep3Vector& dir, float& dlerr) const;
243     float getZDecayLength(vertexNumber nv, vertexNumber mv,
244 loizides 1.5 const ThreeVector& dir, float& dlerr) const;
245 bendavid 1.2 float getImpactPar(vertexNumber prdV, vertexNumber dcyV,
246     const Hep3Vector &v, float &dxyerr) const;
247     float getImpactPar(vertexNumber prdV, vertexNumber dcyV,
248 loizides 1.5 const ThreeVector &v, float &dxyerr) const;
249 loizides 1.1 float get_dr(vertexNumber nv, vertexNumber mv, float& drerr) const;
250     float get_dz(vertexNumber nv, vertexNumber mv, float& dzerr) const;
251     // return location of vertex
252     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     HepSymMatrix getErrorMatrixHep(vertexNumber nv) const;
258     HepSymMatrix getErrorMatrixHep(const int trkId) const;
259     void getPosMomErr (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    
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.1 void getIJKErr(int& err0, int& err1, int& err2) const;
274     // Return each error code from the three ijk indexes
275     int getIJKErr0() const;
276     int getIJKErr1() const;
277     int getIJKErr2() const;
278    
279 paus 1.4 // Get track-id of the track causing a fatal error as indicated by the corresponding ijk error
280 loizides 1.1 int getErrTrackId() const;
281    
282     // Set new track reference point
283     void setTrackReferencePoint(const ThreeVector &ref);
284    
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)+3*CTVMFT_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     //--------------------------------------------------------------------------------------------
395     // Private functions used by class
396     //--------------------------------------------------------------------------------------------
397    
398     private:
399    
400     // Moves reference point of track parameters and errors to _referencePoint
401     //void moveReferencePoint(HepVector &v, HepSymMatrix &m);
402     // Moves reference point of track parameters and errors to _referencePoint. Additionally
403     //checks if track has already been moved by examining it's "derived" link
404     //void moveReferencePoint(const CdfTrack *trk, HepVector &v, HepSymMatrix &m);
405    
406     //--------------------------------------------------------------------------------------------
407     // Data members of class
408     //--------------------------------------------------------------------------------------------
409 paus 1.3 double _bField; // B field in Tesla
410 loizides 1.1
411 paus 1.3 int _currentAllocatedVertexNumber; // index to enum vertexNumber
412     ThreeVector _referencePoint; // reference point of track
413     Hep3Vector _primaryVertex; // primary vertex relative to _referencePoint
414     Hep3Vector _cdfPrimaryVertex; // primary vertex in CDF coordinate system
415     bool _extrapolateTrackErrors; // extrapolate track errors to _referencePoint
416 loizides 1.1 };
417     }
418    
419     //--------------------------------------------------------------------------------------------------
420     // Inline functions
421     //--------------------------------------------------------------------------------------------------
422     inline
423     std::ostream& operator << (std::ostream& os, const mithep::MultiVertexFitter& vfit)
424     {
425     vfit.print(os);
426     return (os);
427     }
428 loizides 1.5
429     //--------------------------------------------------------------------------------------------------
430 loizides 1.1 inline
431     void mithep::MultiVertexFitter::setExcuse()
432     {
433     _fiddle.excuse = 1; // the default
434     }
435 loizides 1.5
436     //--------------------------------------------------------------------------------------------------
437 loizides 1.1 inline
438     void mithep::MultiVertexFitter::setNoExcuse()
439     {
440     _fiddle.excuse = 0; // crash on input error
441     }
442 loizides 1.5
443     //--------------------------------------------------------------------------------------------------
444 loizides 1.1 inline
445     void mithep::MultiVertexFitter::setChisqMax(const float chisqmx)
446     {
447     _fiddle.chisqmax = chisqmx;
448     }
449     #endif