ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitCommon/VertexFit/interface/MultiVertexFitter.h
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
Error occurred while calculating annotation data.
Log Message:
Use CLHEP namespace.

File Contents

# Content
1 //--------------------------------------------------------------------------------------------------
2 // $Id: MultiVertexFitter.h,v 1.6 2009/07/20 03:12:22 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 // (This class was taken from CDF and does therefore not follow our coding conventions.)
68 //--------------------------------------------------------------------------------------------------
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
101 //--------------------------------------------------------------------------------------------
102 // Enumerations
103 //--------------------------------------------------------------------------------------------
104 enum vertexNumber { PRIMARY_VERTEX,VERTEX_1 };
105 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 bool addTrack (const CLHEP::HepVector &pars, const CLHEP::HepSymMatrix &cov,
149 int trackid, float mass, vertexNumber jv);
150
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 void setPrimaryVertex (CLHEP::Hep3Vector pv);
165 bool setPrimaryVertexError(const CLHEP::HepSymMatrix &xverr);
166 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 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
207 //--------------------------------------------------------------------------------------------
208 // Accessors
209 //--------------------------------------------------------------------------------------------
210
211 void setBField (double bField) { _bField = bField; }
212 double bField () const { return _bField; }
213 void setExcuse ();
214 void setNoExcuse();
215
216 std::string expert () const; // name/email of MultiVertexFitter expert
217 int status () const; // return status of fit
218 // overall fit quality paramters
219 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
230 // return fit track four momentum
231 FourVector getTrackP4 (const int trkId) const;
232
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,
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 // return location of vertex
252 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
271 // Accessors for getting information relative to ijk errors get the error code from the three
272 // ijk indexes into the argument variables
273 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 // Get track-id of the track causing a fatal error as indicated by the corresponding ijk error
280 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 std::string _expert; // string: name, email of MultiVertexFitter expert
292 int _stat; // status returned from fit.
293
294 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 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 private:
395 double _bField; // B field in Tesla
396
397 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 };
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
415 //--------------------------------------------------------------------------------------------------
416 inline
417 void mithep::MultiVertexFitter::setExcuse()
418 {
419 _fiddle.excuse = 1; // the default
420 }
421
422 //--------------------------------------------------------------------------------------------------
423 inline
424 void mithep::MultiVertexFitter::setNoExcuse()
425 {
426 _fiddle.excuse = 0; // crash on input error
427 }
428
429 //--------------------------------------------------------------------------------------------------
430 inline
431 void mithep::MultiVertexFitter::setChisqMax(const float chisqmx)
432 {
433 _fiddle.chisqmax = chisqmx;
434 }
435 #endif