Fork me on GitHub

source: git/external/TrackCovariance/VertexFit.cc@ 91ef0b8

Last change on this file since 91ef0b8 was 91ef0b8, checked in by Franco BEDESCHI <bed@…>, 3 years ago

Fixed conflicts and Reginv3
  • Property mode set to 100644
File size: 17.0 KB
Line 
1#include <TMath.h>
2#include <TVectorD.h>
3#include <TVector3.h>
4#include <TMatrixD.h>
5#include <TMatrixDSym.h>
6#include "VertexFit.h"
7//
8// Constructors
9//
10// Empty
11VertexFit::VertexFit()
12{
13 fNtr = 0;
14 fVtxDone = kFALSE;
15 fVtxCst = kFALSE;
16 fxCst.ResizeTo(3);
17 fCovCst.ResizeTo(3, 3);
18 fXv.ResizeTo(3);
19 fcovXv.ResizeTo(3, 3);
20}
21// Parameters and covariances
22VertexFit::VertexFit(Int_t Ntr, TVectorD** trkPar, TMatrixDSym** trkCov)
23{
24 fNtr = Ntr;
25 fVtxDone = kFALSE;
26 fVtxCst = kFALSE;
27 fxCst.ResizeTo(3);
28 fCovCst.ResizeTo(3, 3);
29 fXv.ResizeTo(3);
30 fcovXv.ResizeTo(3, 3);
31 //
32 fPar = trkPar;
33 fCov = trkCov;
34 fChi2List.ResizeTo(Ntr);
35 //
36 ffi = new Double_t[Ntr]; // Fit phases
37 fx0i = new TVectorD * [Ntr]; // Track expansion points
38 for (Int_t i = 0; i < Ntr; i++) fx0i[i] = new TVectorD(3);
39 fai = new TVectorD * [Ntr]; // dx/dphi
40 for (Int_t i = 0; i < Ntr; i++) fai[i] = new TVectorD(3);
41 fa2i = new Double_t[Ntr]; // a'Wa
42 fDi = new TMatrixDSym * [Ntr]; // W-WBW
43 for (Int_t i = 0; i < Ntr; i++) fDi[i] = new TMatrixDSym(3);
44 fWi = new TMatrixDSym * [Ntr]; // (ACA')^-1
45 for (Int_t i = 0; i < Ntr; i++) fWi[i] = new TMatrixDSym(3);
46 fWinvi = new TMatrixDSym * [Ntr]; // ACA'
47 for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3);
48}
49// ObsTrk list
50VertexFit::VertexFit(Int_t Ntr, ObsTrk** track)
51{
52 fNtr = Ntr;
53 fVtxDone = kFALSE;
54 fVtxCst = kFALSE;
55 fxCst.ResizeTo(3);
56 fCovCst.ResizeTo(3, 3);
57 fXv.ResizeTo(3);
58 fcovXv.ResizeTo(3, 3);
59 //
60 fPar = new TVectorD * [Ntr];
61 fCov = new TMatrixDSym * [Ntr];
62 fChi2List.ResizeTo(Ntr);
63 for (Int_t i = 0; i < Ntr; i++)
64 {
65 fPar[i] = new TVectorD(track[i]->GetObsPar());
66 fCov[i] = new TMatrixDSym(track[i]->GetCov());
67 }
68 //
69 ffi = new Double_t[Ntr]; // Fit phases
70 fx0i = new TVectorD * [Ntr]; // Track expansion points
71 for (Int_t i = 0; i < Ntr; i++) fx0i[i] = new TVectorD(3);
72 fai = new TVectorD * [Ntr]; // dx/dphi
73 for (Int_t i = 0; i < Ntr; i++) fai[i] = new TVectorD(3);
74 fa2i = new Double_t[Ntr]; // a'Wa
75 fDi = new TMatrixDSym * [Ntr]; // W-WBW
76 for (Int_t i = 0; i < Ntr; i++) fDi[i] = new TMatrixDSym(3);
77 fWi = new TMatrixDSym * [Ntr]; // (ACA')^-1
78 for (Int_t i = 0; i < Ntr; i++) fWi[i] = new TMatrixDSym(3);
79 fWinvi = new TMatrixDSym * [Ntr]; // ACA'
80 for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3);
81}
82//
83// Destructor
84VertexFit::~VertexFit()
85{
86 fxCst.Clear();
87 fCovCst.Clear();
88 fXv.Clear();
89 fcovXv.Clear();
90 fChi2List.Clear();
91 //
92 for (Int_t i = 0; i < fNtr; i++)
93 {
94 fPar[i]->Clear();
95 fCov[i]->Clear();
96 //
97 fx0i[i]->Clear(); delete fx0i[i];
98 fai[i]->Clear(); delete fai[i];
99 fDi[i]->Clear(); delete fDi[i];
100 fWi[i]->Clear(); delete fWi[i];
101 fWinvi[i]->Clear(); delete fWinvi[i];
102 }
103 fNtr = 0;
104 delete[] fPar;
105 delete[] fCov;
106 delete[] ffi;
107 delete[] fa2i;
108 delete[] fx0i;
109 delete[] fai;
110 delete[] fDi;
111 delete[] fWi;
112 delete[] fWinvi;
113}
114//
115Double_t VertexFit::FastRv1(TVectorD p1, TVectorD p2)
116{
117 //
118 // Find radius of intersection between two tracks in the transverse plane
119 //
120 // p = (D,phi, C, z0, ct)
121 //
122 // Define arrays
123 //
124 Double_t r1 = 1.0 / p1(2);
125 Double_t r2 = 1.0 / p2(2);
126 TVectorD x0 = Fill_x0(p1);
127 TVectorD y0 = Fill_x0(p2);
128 TVectorD n = Fill_a(p1, 0.0);
129 n *= r1;
130 TVectorD k = Fill_a(p2, 0.0);
131 k *= r2;
132 //
133 // Setup and solve linear system
134 //
135 Double_t nn = 0; for (Int_t i = 0; i < 3; i++)nn += n(i) * n(i);
136 Double_t nk = 0; for (Int_t i = 0; i < 3; i++)nk += n(i) * k(i);
137 Double_t kk = 0; for (Int_t i = 0; i < 3; i++)kk += k(i) * k(i);
138 Double_t discr = nn * kk - nk * nk;
139 TMatrixDSym H(2);
140 H(0, 0) = kk;
141 H(0, 1) = nk;
142 H(1, 0) = H(0, 1);
143 H(1, 1) = nn;
144 TVectorD c(2);
145 c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) += n(i) * (y0(i) - x0(i));
146 c(1) = 0; for (Int_t i = 0; i < 3; i++)c(1) += -k(i) * (y0(i) - x0(i));
147 TVectorD smin = (H * c);
148 smin *= 1.0 / discr;
149 //
150 TVectorD X = x0 + smin(0) * n;
151 TVectorD Y = y0 + smin(1) * k;
152 Double_t R1 = TMath::Sqrt(X(0) * X(0) + X(1) * X(1));
153 Double_t R2 = TMath::Sqrt(Y(0) * Y(0) + Y(1) * Y(1));
154 //
155 return 0.5 * (R1 + R2);
156}
157Double_t VertexFit::FastRv(TVectorD p1, TVectorD p2)
158{
159 //
160 // Find radius of minimum distance
161 //
162 // p = (D,phi, C)
163 //
164 // Solving matrix
165 TMatrixDSym H(2);
166 H(0, 0) = -TMath::Cos(p2(1));
167 H(0, 1) = TMath::Cos(p1(1));
168 H(1, 0) = -TMath::Sin(p2(1));
169 H(1, 1) = TMath::Sin(p1(1));
170 Double_t Det = TMath::Sin(p2(1) - p1(1));
171 H *= 1.0 / Det;
172 //
173 // Convergence parameters
174 Int_t Ntry = 0;
175 Int_t NtryMax = 100;
176 Double_t eps = 1000.;
177 Double_t epsMin = 1.0e-6;
178 //
179 // Vertex finding loop
180 //
181 TVectorD cterm(2);
182 cterm(0) = p1(0);
183 cterm(1) = p2(0);
184 TVectorD xv(2);
185 Double_t R = 1000.;
186 while (eps > epsMin)
187 {
188 xv = H * cterm;
189 Ntry++;
190 if (Ntry > NtryMax)
191 {
192 std::cout << "FastRv: maximum number of iteration reached" << std::endl;
193 break;
194 }
195 Double_t Rnew = TMath::Sqrt(xv(0) * xv(0) + xv(1) * xv(1));
196 eps = Rnew - R;
197 R = Rnew;
198 cterm(0) = p1(2) * R * R;
199 cterm(1) = p2(2) * R * R;
200 }
201 //
202 return R;
203}
204
205TMatrixDSym VertexFit::RegInv3(TMatrixDSym& Smat0)
206{
207 //
208 // Regularized inversion of symmetric 3x3 matrix with positive diagonal elements
209 //
210 TMatrixDSym Smat = Smat0;
211 Int_t N = Smat.GetNrows();
212 if (N != 3)
213 {
214 std::cout << "RegInv3 called with matrix size != 3. Abort & return standard inversion." << std::endl;
215 return Smat.Invert();
216 }
217 TMatrixDSym D(N); D.Zero();
218 Bool_t dZero = kTRUE; // No elements less or equal 0 on the diagonal
219 for (Int_t i = 0; i < N; i++) if (Smat(i, i) <= 0.0)dZero = kFALSE;
220 if (dZero)
221 {
222 for (Int_t i = 0; i < N; i++) D(i, i) = 1.0 / TMath::Sqrt(Smat(i, i));
223 TMatrixDSym RegMat = Smat.Similarity(D);
224 TMatrixDSym Q(2);
225 for (Int_t i = 0; i < 2; i++)
226 {
227 for (Int_t j = 0; j < 2; j++)Q(i, j) = RegMat(i, j);
228 }
229 Double_t Det = 1 - Q(0, 1) * Q(1, 0);
230 TMatrixDSym H(2);
231 H = Q;
232 H(0, 1) = -Q(0, 1);
233 H(1, 0) = -Q(1, 0);
234 TVectorD p(2);
235 p(0) = RegMat(0, 2);
236 p(1) = RegMat(1, 2);
237 Double_t pHp = H.Similarity(p);
238 Double_t h = pHp - Det;
239 //
240 TMatrixDSym pp(2); pp.Rank1Update(p);
241 TMatrixDSym F = (h * H) - pp.Similarity(H);
242 F *= 1.0 / Det;
243 TVectorD b = H * p;
244 TMatrixDSym InvReg(3);
245 for (Int_t i = 0; i < 2; i++)
246 {
247 InvReg(i, 2) = b(i);
248 InvReg(2, i) = b(i);
249 for (Int_t j = 0; j < 2; j++) InvReg(i, j) = F(i, j);
250 }
251 InvReg(2, 2) = -Det;
252 //
253 InvReg *= 1.0 / h;
254 //
255 //
256 return InvReg.Similarity(D);
257 }
258 else
259 {
260 D.Zero();
261 for (Int_t i = 0; i < N; i++) D(i, i) = 1.0 / TMath::Sqrt(TMath::Abs(Smat(i, i)));
262 TMatrixDSym RegMat = Smat.Similarity(D);
263 RegMat.Invert();
264 return RegMat.Similarity(D);
265 }
266}
267//
268//
269//
270TMatrixD VertexFit::Fill_A(TVectorD par, Double_t phi)
271{
272 //
273 // Derivative of track 3D position vector with respect to track parameters at constant phase
274 //
275 // par = vector of track parameters
276 // phi = phase
277 //
278 TMatrixD A(3, 5);
279 //
280 // Decode input arrays
281 //
282 Double_t D = par(0);
283 Double_t p0 = par(1);
284 Double_t C = par(2);
285 Double_t z0 = par(3);
286 Double_t ct = par(4);
287 //
288 // Fill derivative matrix dx/d alpha
289 // D
290 A(0, 0) = -TMath::Sin(p0);
291 A(1, 0) = TMath::Cos(p0);
292 A(2, 0) = 0.0;
293 // phi0
294 A(0, 1) = -D * TMath::Cos(p0) + (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);
295 A(1, 1) = -D * TMath::Sin(p0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
296 A(2, 1) = 0.0;
297 // C
298 A(0, 2) = -(TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C * C);
299 A(1, 2) = (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C * C);
300 A(2, 2) = -ct * phi / (2 * C * C);
301 // z0
302 A(0, 3) = 0.0;
303 A(1, 3) = 0.0;
304 A(2, 3) = 1.0;
305 // ct = lambda
306 A(0, 4) = 0.0;
307 A(1, 4) = 0.0;
308 A(2, 4) = phi / (2 * C);
309 //
310 return A;
311}
312//
313TVectorD VertexFit::Fill_a(TVectorD par, Double_t phi)
314{
315 //
316 // Derivative of track 3D position vector with respect to phase at constant track parameters
317 //
318 // par = vector of track parameters
319 // phi = phase
320 //
321 TVectorD a(3);
322 //
323 // Decode input arrays
324 //
325 Double_t D = par(0);
326 Double_t p0 = par(1);
327 Double_t C = par(2);
328 Double_t z0 = par(3);
329 Double_t ct = par(4);
330 //
331 a(0) = TMath::Cos(phi + p0) / (2 * C);
332 a(1) = TMath::Sin(phi + p0) / (2 * C);
333 a(2) = ct / (2 * C);
334 //
335 return a;
336}
337//
338TVectorD VertexFit::Fill_x0(TVectorD par)
339{
340 //
341 // Calculate track 3D position at R = |D| (minimum approach to z-axis)
342 //
343 TVectorD x0(3);
344 //
345 // Decode input arrays
346 //
347 Double_t D = par(0);
348 Double_t p0 = par(1);
349 Double_t C = par(2);
350 Double_t z0 = par(3);
351 Double_t ct = par(4);
352 //
353 x0(0) = -D * TMath::Sin(p0);
354 x0(1) = D * TMath::Cos(p0);
355 x0(2) = z0;
356 //
357 return x0;
358}
359//
360TVectorD VertexFit::Fill_x(TVectorD par, Double_t phi)
361{
362 //
363 // Calculate track 3D position for a given phase, phi
364 //
365 TVectorD x(3);
366 //
367 // Decode input arrays
368 //
369 Double_t D = par(0);
370 Double_t p0 = par(1);
371 Double_t C = par(2);
372 Double_t z0 = par(3);
373 Double_t ct = par(4);
374 //
375 TVectorD x0 = Fill_x0(par);
376 x(0) = x0(0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
377 x(1) = x0(1) - (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);
378 x(2) = x0(2) + ct * phi / (2 * C);
379 //
380 return x;
381}
382//
383void VertexFit::VertexFinder()
384{
385 //
386 // Vertex fit (units are meters)
387 //
388 // Initial variable definitions
389 TVectorD x(3);
390 TMatrixDSym covX(3);
391 TVectorD x0(3); for (Int_t v = 0; v < 3; v++)x0(v) = 100.; // set to large value
392 Double_t Chi2 = 0;
393 //
394 // Stored quantities
395 Double_t* fi = new Double_t[fNtr]; // Phases
396 TVectorD** x0i = new TVectorD * [fNtr]; // Track expansion point
397 TVectorD** ai = new TVectorD * [fNtr]; // dx/dphi
398 Double_t* a2i = new Double_t[fNtr]; // a'Wa
399 TMatrixDSym** Di = new TMatrixDSym * [fNtr]; // W-WBW
400 TMatrixDSym** Wi = new TMatrixDSym * [fNtr]; // (ACA')^-1
401 TMatrixDSym** Winvi = new TMatrixDSym * [fNtr]; // ACA'
402 //
403 // vertex radius approximation
404 // Maximum impact parameter
405 Double_t Rd = 0;
406 for (Int_t i = 0; i < fNtr; i++)
407 {
408 //ObsTrk* t = tracks[i];
409 TVectorD par = *fPar[i];
410 Double_t Dabs = TMath::Abs(par(0));
411 if (Dabs > Rd)Rd = Dabs;
412 }
413 //
414 // Find track pair with largest phi difference
415 Int_t isel = 0; Int_t jsel = 0; // selected track indices
416 Double_t dphiMax = 0.0; // Max phi difference
417
418 for (Int_t i = 0; i < fNtr - 1; i++)
419 {
420 //ObsTrk* ti = tracks[i];
421 TVectorD pari = *fPar[i];
422 Double_t phi1 = pari(1);
423
424 for (Int_t j = i + 1; j < fNtr; j++)
425 {
426 //ObsTrk* tj = tracks[j];
427 TVectorD parj = *fPar[j];
428 Double_t phi2 = parj(1);
429 Double_t dphi = TMath::Abs(phi2 - phi1);
430 if (dphi > TMath::Pi())dphi = TMath::TwoPi() - dphi;
431 if (dphi > dphiMax)
432 {
433 isel = i; jsel = j;
434 dphiMax = dphi;
435 }
436 }
437 }
438 //
439 TVectorD p1 = *fPar[isel];
440 TVectorD p2 = *fPar[jsel];
441 Double_t R = FastRv1(p1, p2);
442 if (R > 1000.0) R = Rd;
443 R = 0.9 * R + 0.1 * Rd;
444 //
445 // Iteration properties
446 //
447 Int_t Ntry = 0;
448 Int_t TryMax = 100;
449 Double_t eps = 1.0e-9; // vertex stability
450 Double_t epsi = 1000.;
451 //
452 while (epsi > eps && Ntry < TryMax) // Iterate until found vertex is stable
453 {
454 x.Zero();
455 TVectorD cterm(3); TMatrixDSym H(3); TMatrixDSym DW1D(3);
456 covX.Zero(); // Reset vertex covariance
457 cterm.Zero(); // Reset constant term
458 H.Zero(); // Reset H matrix
459 DW1D.Zero();
460 //
461 //std::cout << "VertexFinder: start loop on tracks" << std::endl;
462 for (Int_t i = 0; i < fNtr; i++)
463 {
464 // Get track helix parameters and their covariance matrix
465 TVectorD par = *fPar[i];
466 TMatrixDSym Cov = *fCov[i];
467
468 Double_t fs;
469 if (Ntry <= 0) // Initialize all phases on first pass
470 {
471 Double_t D = par(0);
472 Double_t C = par(2);
473 Double_t arg = TMath::Max(1.0e-6, (R * R - D * D) / (1 + 2 * C * D));
474 fs = 2 * TMath::ASin(C * TMath::Sqrt(arg));
475 fi[i] = fs;
476 }
477 //
478 // Starting values
479 //
480 fs = fi[i]; // Get phase
481 //std::cout << "VertexFinder: phase fs set" << std::endl;
482 TVectorD xs = Fill_x(par, fs);
483 //std::cout << "VertexFinder: position xs set" << std::endl;
484 x0i[i] = new TVectorD(xs); // Start helix position
485 //std::cout << "VertexFinder: position x0i stored" << std::endl;
486 // W matrix = (A*C*A')^-1; W^-1 = A*C*A'
487 TMatrixD A = Fill_A(par, fs); // A = dx/da = derivatives wrt track parameters
488 //std::cout << "VertexFinder: derivatives A set" << std::endl;
489 TMatrixDSym Winv = Cov.Similarity(A); // W^-1 = A*C*A'
490 Winvi[i] = new TMatrixDSym(Winv); // Store W^-1 matrix
491 //std::cout << "VertexFinder: Winvi stored" << std::endl;
492 TMatrixDSym W = RegInv3(Winv); // W = (A*C*A')^-1
493 Wi[i] = new TMatrixDSym(W); // Store W matrix
494 //std::cout << "VertexFinder: Wi stored" << std::endl;
495 TVectorD a = Fill_a(par, fs); // a = dx/ds = derivatives wrt phase
496 //std::cout << "VertexFinder: derivatives a set" << std::endl;
497 ai[i] = new TVectorD(a); // Store a
498 //std::cout << "VertexFinder: derivatives a stored" << std::endl;
499 Double_t a2 = W.Similarity(a);
500 a2i[i] = a2; // Store a2
501 // Build D matrix
502 TMatrixDSym B(3);
503
504 B.Rank1Update(a, 1.0);
505 B *= -1. / a2;
506 B.Similarity(W);
507 TMatrixDSym Ds = W + B; // D matrix
508 Di[i] = new TMatrixDSym(Ds); // Store D matrix
509 //std::cout << "VertexFinder: matrix Di stored" << std::endl;
510 TMatrixDSym DsW1Ds = Winv.Similarity(Ds); // Service matrix to calculate covX
511 DW1D += DsW1Ds;
512 // Update hessian
513 H += Ds;
514 // update constant term
515 cterm += Ds * xs;
516 } // End loop on tracks
517 //
518 // update vertex position
519 TMatrixDSym H1 = RegInv3(H);
520 x = H1 * cterm;
521 //std::cout << "VertexFinder: x vertex set" << std::endl;
522 // Update vertex covariance
523 covX = DW1D.Similarity(H1);
524 //std::cout << "VertexFinder: cov vertex set" << std::endl;
525 // Update phases and chi^2
526 Chi2 = 0.0;
527 for (Int_t i = 0; i < fNtr; i++)
528 {
529 TVectorD lambda = (*Di[i]) * (*x0i[i] - x);
530 TMatrixDSym Wm1 = *Winvi[i];
531 fChi2List(i) = Wm1.Similarity(lambda);
532 Chi2 += fChi2List(i);
533 TVectorD a = *ai[i];
534 TVectorD b = (*Wi[i]) * (x - *x0i[i]);
535 for (Int_t j = 0; j < 3; j++)fi[i] += a(j) * b(j) / a2i[i];
536 }
537
538 //
539 TVectorD dx = x - x0;
540 x0 = x;
541 // update vertex stability
542 TMatrixDSym Hess = RegInv3(covX);
543 epsi = Hess.Similarity(dx);
544 Ntry++;
545 //
546 // Store result
547 //
548 fXv = x; // Vertex position
549 fcovXv = covX; // Vertex covariance
550 fChi2 = Chi2; // Vertex fit Chi2
551 //
552 // Store intermediate data
553 //
554
555 //std::cout << "VertexFinder: before store intermediate data" << std::endl;
556 for (Int_t i = 0; i < fNtr; i++)
557 {
558 //std::cout << "VertexFinder: inside store intermediate data" << std::endl;
559 //std::cout << "i = " << i << ", fi[i] = " << fi[i] << std::endl;
560 //std::cout << "i = " << i << ", ffi[i] = " << ffi[i] << std::endl;
561 ffi[i] = fi[i]; // Fit phases
562 //std::cout << "VertexFinder: fi stored" << std::endl;
563 fx0i[i] = x0i[i]; // Track expansion points
564 //std::cout << "VertexFinder: x0i stored" << std::endl;
565 fai[i] = ai[i]; // dx/dphi
566 //std::cout << "VertexFinder: ai stored" << std::endl;
567 fa2i[i] = a2i[i]; // a'Wa
568 //std::cout << "VertexFinder: a2i stored" << std::endl;
569 fDi[i] = Di[i]; // W-WBW
570 //std::cout << "VertexFinder: Di stored" << std::endl;
571 fWi[i] = Wi[i]; // (ACA')^-1
572 //std::cout << "VertexFinder: Wi stored" << std::endl;
573 fWinvi[i] = Winvi[i]; // ACA'
574 //std::cout << "VertexFinder: Winvi stored" << std::endl;
575 }
576 //std::cout << "Iteration " << Ntry << " completed - Before cleanup" << std::endl;
577 //
578 // Cleanup
579 //
580 for (Int_t i = 0; i < fNtr; i++)
581 {
582 x0i[i]->Clear();
583 Winvi[i]->Clear();
584 Wi[i]->Clear();
585 ai[i]->Clear();
586 Di[i]->Clear();
587
588 delete x0i[i];
589 delete Winvi[i];
590 delete Wi[i];
591 delete ai[i];
592 delete Di[i];
593 }
594
595 //std::cout << "Iteration " << Ntry << " completed - After cleanup" << std::endl;
596 }
597 //
598 fVtxDone = kTRUE; // Set fitting completion flag
599 //
600 delete[] fi; // Phases
601 delete[] x0i; // Track expansion point
602 delete[] ai; // dx/dphi
603 delete[] a2i; // a'Wa
604 delete[] Di; // W-WBW
605 delete[] Wi; // (ACA')^-1
606 delete[] Winvi; // ACA'
607}
608//
609TVectorD VertexFit::GetVtx()
610{
611 //std::cout << "GetVtx: flag set to " << fVtxDone << std::endl;
612 if (!fVtxDone)VertexFinder();
613 return fXv;
614}
615
616TMatrixDSym VertexFit::GetVtxCov()
617{
618 if (!fVtxDone)VertexFinder();
619 return fcovXv;
620}
621
622Double_t VertexFit::GetVtxChi2()
623{
624 if (!fVtxDone)VertexFinder();
625 return fChi2;
626}
627
628TVectorD VertexFit::GetVtxChi2List()
629{
630 if (!fVtxDone)VertexFinder();
631 return fChi2List;
632}
633//
634// Handle tracks/constraints
635void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov) // Add gaussian vertex constraint
636{
637 std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
638}
639//
640void VertexFit::AddTrk(TVectorD par, TMatrixDSym Cov) // Add track to input list
641{
642 std::cout << "VertexFit::AddTrk: Not implemented yet" << std::endl;
643}
644void VertexFit::RemoveTrk(Int_t iTrk) // Remove iTrk track
645{
646 std::cout << "VertexFit::RemoveTrk: Not implemented yet" << std::endl;
647}
Note: See TracBrowser for help on using the repository browser.