Fork me on GitHub

Context Navigation

source: git/external/TrackCovariance/VertexFit.cc@ cb80e6f

Visit:

Last change on this file since cb80e6f was 91ef0b8, checked in by Franco BEDESCHI <bed@…>, 4 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
11	VertexFit::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
22	VertexFit::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
50	VertexFit::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
84	VertexFit::~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	//
115	Double_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	}
157	Double_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
205	TMatrixDSym 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	//
270	TMatrixD 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	//
313	TVectorD 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	//
338	TVectorD 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	//
360	TVectorD 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	//
383	void 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 = (ACA')^-1; W^-1 = ACA'
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 = ACA'
490	Winvi[i] = new TMatrixDSym(Winv); // Store W^-1 matrix
491	//std::cout << "VertexFinder: Winvi stored" << std::endl;
492	TMatrixDSym W = RegInv3(Winv); // W = (ACA')^-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	//
609	TVectorD VertexFit::GetVtx()
610	{
611	//std::cout << "GetVtx: flag set to " << fVtxDone << std::endl;
612	if (!fVtxDone)VertexFinder();
613	return fXv;
614	}
615
616	TMatrixDSym VertexFit::GetVtxCov()
617	{
618	if (!fVtxDone)VertexFinder();
619	return fcovXv;
620	}
621
622	Double_t VertexFit::GetVtxChi2()
623	{
624	if (!fVtxDone)VertexFinder();
625	return fChi2;
626	}
627
628	TVectorD VertexFit::GetVtxChi2List()
629	{
630	if (!fVtxDone)VertexFinder();
631	return fChi2List;
632	}
633	//
634	// Handle tracks/constraints
635	void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov) // Add gaussian vertex constraint
636	{
637	std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
638	}
639	//
640	void VertexFit::AddTrk(TVectorD par, TMatrixDSym Cov) // Add track to input list
641	{
642	std::cout << "VertexFit::AddTrk: Not implemented yet" << std::endl;
643	}
644	void 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.

Download in other formats: