ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitCommon/VertexFit/interface/MultiVertexFitter.h
Revision: 1.3
Committed: Thu Nov 13 16:34:28 2008 UTC (16 years, 5 months ago) by paus
Content type: text/plain
Branch: MAIN
Changes since 1.2: +8 -8 lines
Log Message:
Submitting the ugly three versions. But no elegant solution found.

File Contents

# User Rev Content
1 loizides 1.1 //--------------------------------------------------------------------------------------------------
2 paus 1.3 // $Id: MultiVertexFitter.h,v 1.2 2008/09/30 12:48:12 bendavid 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     //--------------------------------------------------------------------------------------------
100     // Enumerations
101     //--------------------------------------------------------------------------------------------
102 paus 1.3 enum vertexNumber { PRIMARY_VERTEX,VERTEX_1 };
103 loizides 1.1 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 bendavid 1.2 float getDecayLength(vertexNumber nv, vertexNumber mv, const ThreeVector& dir,
239     float& dlerr) const;
240     float getZDecayLength(vertexNumber nv, vertexNumber mv,
241     const Hep3Vector& dir, float& dlerr) const;
242     float getZDecayLength(vertexNumber nv, vertexNumber mv,
243     const ThreeVector& dir, float& dlerr) const;
244     float getImpactPar(vertexNumber prdV, vertexNumber dcyV,
245     const Hep3Vector &v, float &dxyerr) const;
246     float getImpactPar(vertexNumber prdV, vertexNumber dcyV,
247     const ThreeVector &v, float &dxyerr) const;
248 loizides 1.1 float get_dr(vertexNumber nv, vertexNumber mv, float& drerr) const;
249     float get_dz(vertexNumber nv, vertexNumber mv, float& dzerr) const;
250     // return location of vertex
251     Hep3Vector getVertexHep(vertexNumber nv) const;
252     ThreeVector getVertex(vertexNumber nv) const;
253     // return error matrix element.
254     ThreeSymMatrix getErrorMatrix(vertexNumber nv) const;
255     double getErrorMatrixHep(int j, int k) const;
256     HepSymMatrix getErrorMatrixHep(vertexNumber nv) const;
257     HepSymMatrix getErrorMatrixHep(const int trkId) const;
258     void getPosMomErr (HepMatrix& errors) const;
259     int vOff (vertexNumber jv) const;
260     int tOff (const int trkId) const;
261     int pOff (int lp) const;
262     int cOff (int lc) const;
263     int mOff () const;
264     double VMat (int i, int j) const;
265     float getPtError (const int trkId) const;
266     MultiVertexFitter::vertexNumber
267     allocateVertexNumber();
268     void resetAllocatedVertexNumber();
269    
270     // Accessors for getting information relative to ijk errors Get the error code from the three ijk
271     // indexes into the argument variables
272     void getIJKErr(int& err0, int& err1, int& err2) const;
273     // Return each error code from the three ijk indexes
274     int getIJKErr0() const;
275     int getIJKErr1() const;
276     int getIJKErr2() const;
277    
278     // Get the track-id of the track causing a fatal error as indicated by the corresponding ijk error
279     int getErrTrackId() const;
280    
281     // Set new track reference point
282     void setTrackReferencePoint(const ThreeVector &ref);
283    
284     //--------------------------------------------------------------------------------------------
285     // Overload operators
286     //--------------------------------------------------------------------------------------------
287     friend std::ostream& operator << (std::ostream& os, const MultiVertexFitter& vfit);
288    
289     protected:
290     std::string _expert; // string: name and email of MultiVertexFitter expert
291     int _stat; // status returned from fit.
292    
293     static const int _maxvtx = CTVMFT_MAXVTX; // 6; // Maximum number of vertices
294     static const int _maxmcn = CTVMFT_MAXMCN; // 4; // Maximum number of mass constraints
295     static const int _maxtrk = CTVMFT_MAXTRK; // 50; // Maximum number of tracks
296     static const int _maxitr = CTVMFT_MAXITR; // 10; // Maximum number of iteration steps
297     static const int _maxdim = 5*(CTVMFT_MAXVTX+1)+3*CTVMFT_MAXTRK+CTVMFT_MAXMCN;
298    
299     // FIDDLE must have access to MultiVertexFitter's (other) protected data like _maxvtx, etc. It
300     // must therefore be a friend.
301     struct FIDDLE;
302     friend struct FIDDLE;
303     struct FIDDLE {
304     int excuse;
305     float chisqmax;
306     };
307    
308     // CTVMQ must have access to MultiVertexFitter's (other) protected data like _maxvtx, etc. It
309     // must therefore be a friend.
310     struct CTVMQ;
311     friend struct CTVMQ;
312     struct CTVMQ {
313     int runnum;
314     int trgnum;
315     int iter;
316     int ntscut;
317     int nvertx;
318     int nmassc;
319     int ntrack;
320     int trkvtx[_maxvtx][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
321     int trkmcn[_maxmcn][_maxtrk]; // Logical in FORTRAN, but integer here to get the size
322     int vtxpnt[2][_maxvtx];
323     float cmass [_maxmcn];
324     int cvtx [_maxvtx];
325     int vtxvtx[_maxvtx][_maxvtx]; // Logical in FORTRAN, but integer here to get the size
326     char tkbank[_maxtrk][4];
327     int list [_maxtrk];
328     float tmass [_maxtrk];
329     int matdim;
330     int tkerr [_maxtrk];
331     int ndof;
332     float chisqr[_maxitr+1];
333     float chit [_maxtrk];
334     float chiv [_maxvtx+1];
335     float chim [_maxmcn];
336     float xyzpv0[3];
337     float exyzpv[3][3];
338     float xzslope;
339     float yzslope;
340     float xyzvrt[_maxvtx+1][3];
341     float dxyzpv[3];
342     float par [_maxtrk][5];
343     float g [_maxtrk][5][5];
344     float trkp4 [6][_maxtrk];
345     float vtxp4 [_maxvtx][4];
346     float mcnp4 [_maxmcn][4];
347     float dda [8][_maxtrk];
348     int voff [_maxvtx];
349     int toff [_maxtrk];
350     int poff [_maxvtx];
351     int coff [_maxvtx];
352     int moff;
353     float par0 [_maxtrk][5];
354     float pardif[_maxtrk][5];
355     float fmcdif[_maxmcn];
356     float pcon [2][_maxvtx];
357     float sang [2][_maxvtx];
358     float drmax;
359     float rvmax;
360     float dzmax;
361     float trnmax;
362     float dsmin;
363     int ijkerr[3];
364     float pscale;
365     };
366    
367     // CTVMFR must have access to MultiVertexFitter's (other) protected data like _maxdim, etc.
368     // It must therefore be a friend.
369     struct CTVMFR;
370     friend struct CTVMFR;
371     struct CTVMFR {
372     double vmat[_maxdim+1][_maxdim];
373     };
374    
375     // TRKPRM must have access to MultiVertexFitter's (other) protected data like _maxtrk, etc.
376     // It must therefore be a friend.
377     struct TRKPRM;
378     friend struct TRKPRM;
379     struct TRKPRM {
380     float trhelix[_maxtrk][5];
381     float trem [_maxtrk][5][5];
382     };
383    
384     CTVMQ _ctvmq;
385     CTVMQ* _ctvmq_com;
386     CTVMFR _ctvmfr;
387     CTVMFR* _ctvmfr_com;
388     FIDDLE _fiddle;
389     FIDDLE* _fiddle_com;
390     TRKPRM _trkprm;
391     TRKPRM* _trkprm_com;
392    
393     //--------------------------------------------------------------------------------------------
394     // Private functions used by class
395     //--------------------------------------------------------------------------------------------
396    
397     private:
398    
399     // Moves reference point of track parameters and errors to _referencePoint
400     //void moveReferencePoint(HepVector &v, HepSymMatrix &m);
401     // Moves reference point of track parameters and errors to _referencePoint. Additionally
402     //checks if track has already been moved by examining it's "derived" link
403     //void moveReferencePoint(const CdfTrack *trk, HepVector &v, HepSymMatrix &m);
404    
405     //--------------------------------------------------------------------------------------------
406     // Data members of class
407     //--------------------------------------------------------------------------------------------
408 paus 1.3 double _bField; // B field in Tesla
409 loizides 1.1
410 paus 1.3 int _currentAllocatedVertexNumber; // index to enum vertexNumber
411     ThreeVector _referencePoint; // reference point of track
412     Hep3Vector _primaryVertex; // primary vertex relative to _referencePoint
413     Hep3Vector _cdfPrimaryVertex; // primary vertex in CDF coordinate system
414     bool _extrapolateTrackErrors; // extrapolate track errors to _referencePoint
415 loizides 1.1 };
416     }
417    
418     //--------------------------------------------------------------------------------------------------
419     // Inline functions
420     //--------------------------------------------------------------------------------------------------
421     inline
422     std::ostream& operator << (std::ostream& os, const mithep::MultiVertexFitter& vfit)
423     {
424     vfit.print(os);
425     return (os);
426     }
427     inline
428     void mithep::MultiVertexFitter::setExcuse()
429     {
430     _fiddle.excuse = 1; // the default
431     }
432     inline
433     void mithep::MultiVertexFitter::setNoExcuse()
434     {
435     _fiddle.excuse = 0; // crash on input error
436     }
437     inline
438     void mithep::MultiVertexFitter::setChisqMax(const float chisqmx)
439     {
440     _fiddle.chisqmax = chisqmx;
441     }
442     #endif