Changeset 82db145 in git for external/TrackCovariance/VertexFit.cc
- Timestamp:
- Apr 12, 2021, 6:33:23 PM (4 years ago)
- Branches:
- master
- Children:
- c5696dd
- Parents:
- 3cfe61d
- File:
-
- 1 edited
-
external/TrackCovariance/VertexFit.cc (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/VertexFit.cc
r3cfe61d r82db145 8 8 // Constructors 9 9 // 10 // Empty 10 // 11 // Empty construction (to be used when adding tracks later with AddTrk() ) 11 12 VertexFit::VertexFit() 12 13 { 13 14 fNtr = 0; 15 fRold = -1.0; 14 16 fVtxDone = kFALSE; 15 17 fVtxCst = kFALSE; … … 19 21 fcovXv.ResizeTo(3, 3); 20 22 } 21 // Parameters and covariances 23 // 24 // Build from list of parameters and covariances 22 25 VertexFit::VertexFit(Int_t Ntr, TVectorD** trkPar, TMatrixDSym** trkCov) 23 26 { 24 27 fNtr = Ntr; 28 fRold = -1.0; 25 29 fVtxDone = kFALSE; 26 30 fVtxCst = kFALSE; … … 30 34 fcovXv.ResizeTo(3, 3); 31 35 // 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 36 for (Int_t i = 0; i < fNtr; i++) 37 { 38 TVectorD pr = *trkPar[i]; 39 fPar.push_back(new TVectorD(pr)); 40 TMatrixDSym cv = *trkCov[i]; 41 fCov.push_back(new TMatrixDSym(cv)); 42 } 43 fChi2List.ResizeTo(fNtr); 44 // 45 } 46 // 47 // Build from ObsTrk list of tracks 50 48 VertexFit::VertexFit(Int_t Ntr, ObsTrk** track) 51 49 { 52 50 fNtr = Ntr; 51 fRold = -1.0; 53 52 fVtxDone = kFALSE; 54 53 fVtxCst = kFALSE; … … 58 57 fcovXv.ResizeTo(3, 3); 59 58 // 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); 59 fChi2List.ResizeTo(fNtr); 60 for (Int_t i = 0; i < fNtr; i++) 61 { 62 fPar.push_back(new TVectorD(track[i]->GetObsPar())); 63 fCov.push_back(new TMatrixDSym(track[i]->GetCov())); 64 } 81 65 } 82 66 // 83 67 // Destructor 68 // 69 void VertexFit::ResetWrkArrays() 70 { 71 Int_t N = (Int_t)ffi.size(); 72 for (Int_t i = 0; i < N; i++) 73 { 74 if (fx0i[i]) { fx0i[i]->Clear(); delete fx0i[i]; } 75 if (fai[i]) { fai[i]->Clear(); delete fai[i]; } 76 if (fDi[i]) { fDi[i]->Clear(); delete fDi[i]; } 77 if (fWi[i]) { fWi[i]->Clear(); delete fWi[i]; } 78 if (fWinvi[i]){ fWinvi[i]->Clear(); delete fWinvi[i]; } 79 } 80 fa2i.clear(); 81 fx0i.clear(); 82 fai.clear(); 83 fDi.clear(); 84 fWi.clear(); 85 fWinvi.clear(); 86 } 84 87 VertexFit::~VertexFit() 85 88 { 86 fxCst.Clear(); 87 fCovCst.Clear(); 88 fXv.Clear(); 89 fcovXv.Clear(); 90 fChi2List.Clear(); 89 fxCst.Clear(); 90 fCovCst.Clear(); 91 fXv.Clear(); 92 fcovXv.Clear(); 93 fChi2List.Clear(); 91 94 // 92 95 for (Int_t i = 0; i < fNtr; i++) 93 96 { 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 } 97 fPar[i]->Clear(); delete fPar[i]; 98 fCov[i]->Clear(); delete fCov[i]; 99 } 100 fPar.clear(); 101 fCov.clear(); 102 // 103 ResetWrkArrays(); 104 ffi.clear(); 103 105 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 106 } 107 // 108 Double_t VertexFit::FastRv(TVectorD p1, TVectorD p2) 109 { 110 // 111 // Find radius of minimum distance between two tracks 119 112 // 120 113 // p = (D,phi, C, z0, ct) … … 122 115 // Define arrays 123 116 // 124 Double_t r1 = 1.0 / p1(2); 125 Double_t r2 = 1.0 / p2(2); 117 Double_t C1 = p1(2); 118 Double_t C2 = p2(2); 119 Double_t ph1 = p1(1); 120 Double_t ph2 = p2(1); 126 121 TVectorD x0 = Fill_x0(p1); 127 122 TVectorD y0 = Fill_x0(p2); 128 123 TVectorD n = Fill_a(p1, 0.0); 129 n *= r1;124 n *= (2*C1); 130 125 TVectorD k = Fill_a(p2, 0.0); 131 k *= r2;126 k *= (2*C2); 132 127 // 133 128 // Setup and solve linear system … … 150 145 TVectorD X = x0 + smin(0) * n; 151 146 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) 147 // 148 // Higher order corrections 149 X(0) += -C1 * smin(0) * smin(0) * TMath::Sin(ph1); 150 X(1) += C1 * smin(0) * smin(0) * TMath::Cos(ph1); 151 Y(0) += -C2 * smin(1) * smin(1) * TMath::Sin(ph2); 152 Y(1) += C2 * smin(1) * smin(1) * TMath::Cos(ph2); 153 // 154 TVectorD Xavg = 0.5 * (X + Y); 155 // 156 // 157 return TMath::Sqrt(Xavg(0)*Xavg(0)+Xavg(1)*Xavg(1)); 158 } 159 // 160 // Starting radius determination 161 Double_t VertexFit::StartRadius() 162 { 163 // 164 // Maximum impact parameter 165 Double_t Rd = 0; 166 for (Int_t i = 0; i < fNtr; i++) 167 { 168 TVectorD par = *fPar[i]; 169 Double_t Dabs = TMath::Abs(par(0)); 170 if (Dabs > Rd)Rd = Dabs; 171 } 172 //----------------------------- 173 // 174 // Find track pair with phi difference closest to pi/2 175 Int_t isel = 0; Int_t jsel = 0; // selected track indices 176 Double_t dSinMax = 0.0; // Max phi difference 177 for (Int_t i = 0; i < fNtr - 1; i++) 178 { 179 TVectorD pari = *fPar[i]; 180 Double_t phi1 = pari(1); 181 182 for (Int_t j = i + 1; j < fNtr; j++) 191 183 { 192 std::cout << "FastRv: maximum number of iteration reached" << std::endl; 193 break; 184 TVectorD parj = *fPar[j]; 185 Double_t phi2 = parj(1); 186 Double_t Sindphi = TMath::Abs(TMath::Sin(phi2 - phi1)); 187 if (Sindphi > dSinMax) 188 { 189 isel = i; jsel = j; 190 dSinMax = Sindphi; 191 } 194 192 } 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 } 193 } 194 // 195 //------------------------------------------ 196 // 197 // Find radius of minimum distrance between tracks 198 TVectorD p1 = *fPar[isel]; 199 TVectorD p2 = *fPar[jsel]; 200 Double_t R = FastRv(p1, p2); 201 // 202 R = 0.9 * R + 0.1 * Rd; // Protect for overshoot 201 203 // 202 204 return R; 203 205 } 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++) 206 // 207 // Regularized symmetric matrix inversion 208 // 209 TMatrixDSym VertexFit::RegInv(TMatrixDSym& Min) 210 { 211 TMatrixDSym M = Min; // Decouple from input 212 Int_t N = M.GetNrows(); // Matrix size 213 TMatrixDSym D(N); D.Zero(); // Normaliztion matrix 214 TMatrixDSym R(N); // Normarized matrix 215 TMatrixDSym Rinv(N); // Inverse of R 216 TMatrixDSym Minv(N); // Inverse of M 217 // 218 // Check for 0's and normalize 219 for (Int_t i = 0; i < N; i++) 220 { 221 if (M(i, i) != 0.0) D(i, i) = 1. / TMath::Sqrt(TMath::Abs(M(i, i))); 222 else D(i, i) = 1.0; 223 } 224 R = M.Similarity(D); 225 // 226 // Recursive algorithms stops when N = 2 227 // 228 //**************** 229 // case N = 2 *** 230 //**************** 231 if (N == 2) 232 { 233 Double_t det = R(0, 0) * R(1, 1) - R(0, 1) * R(1, 0); 234 if (det == 0) 226 235 { 227 for (Int_t j = 0; j < 2; j++)Q(i, j) = RegMat(i, j); 236 std::cout << "VertexFit::RegInv: null determinant for N = 2" << std::endl; 237 Rinv.Zero(); // Return null matrix 228 238 } 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++) 239 else 246 240 { 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); 241 // invert matrix 242 Rinv(0, 0) = R(1, 1); 243 Rinv(0, 1) = -R(0, 1); 244 Rinv(1, 0) = Rinv(0, 1); 245 Rinv(1, 1) = R(0, 0); 246 Rinv *= 1. / det; 250 247 } 251 InvReg(2, 2) = -Det; 252 // 253 InvReg *= 1.0 / h; 254 // 255 // 256 return InvReg.Similarity(D); 257 } 248 } 249 //**************** 250 // case N > 2 *** 251 //**************** 258 252 else 259 253 { 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 } 254 // Break up matrix 255 TMatrixDSym Q = R.GetSub(0, N - 2, 0, N - 2); // Upper left 256 TVectorD p(N - 1); 257 for (Int_t i = 0; i < N - 1; i++)p(i) = R(N - 1, i); 258 Double_t q = R(N - 1, N - 1); 259 //Invert pieces and re-assemble 260 TMatrixDSym Ainv(N - 1); 261 TMatrixDSym A(N - 1); 262 if (TMath::Abs(q) > 1.0e-15) 263 { 264 // Case |q| > 0 265 Ainv.Rank1Update(p, -1.0 / q); 266 Ainv += Q; 267 A = RegInv(Ainv); // Recursive call 268 TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A; 269 // 270 TVectorD b = (-1.0 / q) * (A * p); 271 for (Int_t i = 0; i < N - 1; i++) 272 { 273 Rinv(N - 1, i) = b(i); 274 Rinv(i, N - 1) = b(i); 275 } 276 // 277 Double_t pdotb = 0.; 278 for (Int_t i = 0; i < N - 1; i++)pdotb += p(i) * b(i); 279 Double_t c = (1.0 - pdotb) / q; 280 Rinv(N - 1, N - 1) = c; 281 } 282 else 283 { 284 // case q = 0 285 TMatrixDSym Qinv = RegInv(Q); // Recursive call 286 Double_t a = Qinv.Similarity(p); 287 Double_t c = -1.0 / a; 288 Rinv(N - 1, N - 1) = c; 289 // 290 TVectorD b = (1.0 / a) * (Qinv * p); 291 for (Int_t i = 0; i < N - 1; i++) 292 { 293 Rinv(N - 1, i) = b(i); 294 Rinv(i, N - 1) = b(i); 295 } 296 // 297 A.Rank1Update(p, -1 / a); 298 A += Q; 299 A.Similarity(Qinv); 300 TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A; 301 } 302 } 303 Minv = Rinv.Similarity(D); 304 return Minv; 266 305 } 267 306 // … … 381 420 } 382 421 // 383 void VertexFit::VertexFinder() 384 { 385 // 386 // Vertex fit (units are meters) 422 void VertexFit::UpdateTrkArrays(Int_t i) 423 { 424 // 425 // Get track parameters and their covariance 426 TVectorD par = *fPar[i]; 427 TMatrixDSym Cov = *fCov[i]; 428 // 429 // Fill all track related work arrays arrays 430 Double_t fs = ffi[i]; // Get phase 431 TVectorD xs = Fill_x(par, fs); 432 fx0i.push_back(new TVectorD(xs)); // Start helix position 433 // 434 TMatrixD A = Fill_A(par, fs); // A = dx/da = derivatives wrt track parameters 435 TMatrixDSym Winv = Cov.Similarity(A); // W^-1 = A*C*A' 436 fWinvi.push_back(new TMatrixDSym(Winv)); // Store W^-1 matrix 437 // 438 TMatrixDSym W = RegInv(Winv); // W = (A*C*A')^-1 439 fWi.push_back(new TMatrixDSym(W)); // Store W matrix 440 // 441 TVectorD a = Fill_a(par, fs); // a = dx/ds = derivatives wrt phase 442 fai.push_back(new TVectorD(a)); // Store a 443 // 444 Double_t a2 = W.Similarity(a); 445 fa2i.push_back(a2); // Store a2 446 // 447 // Build D matrix 448 TMatrixDSym B(3); 449 B.Rank1Update(a, -1. / a2); 450 B.Similarity(W); 451 TMatrixDSym Ds = W + B; // D matrix 452 fDi.push_back(new TMatrixDSym(Ds)); // Store D matrix 453 } 454 // 455 void VertexFit::VertexFitter() 456 { 457 //std::cout << "VertexFitter: just in" << std::endl; 458 if (fNtr < 2) 459 { 460 std::cout << "VertexFit::VertexFitter - Method called with less than 2 tracks - Aborting " << std::endl; 461 std::exit(1); 462 } 463 // 464 // Vertex fit 387 465 // 388 466 // Initial variable definitions 389 467 TVectorD x(3); 390 468 TMatrixDSym covX(3); 391 TVectorD x0(3); for (Int_t v = 0; v < 3; v++)x0(v) = 100.; // set to large value392 469 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; 470 TVectorD x0 = fXv; // If previous fit done 471 if (fRold < 0.0)for (Int_t i = 0; i < 3; i++)x0(i) = 1000.; // Set to arbitrary large value if not 472 // 473 // Starting vertex radius approximation 474 // 475 Double_t R = fRold; // Use previous fit if available 476 if (R < 0.0) R = StartRadius(); // Rough vertex estimate 444 477 // 445 478 // Iteration properties … … 450 483 Double_t epsi = 1000.; 451 484 // 485 // Iteration loop 452 486 while (epsi > eps && Ntry < TryMax) // Iterate until found vertex is stable 453 487 { 488 // Initialize arrays 454 489 x.Zero(); 455 490 TVectorD cterm(3); TMatrixDSym H(3); TMatrixDSym DW1D(3); 456 covX.Zero(); // Reset vertex covariance491 covX.Zero(); // Reset vertex covariance 457 492 cterm.Zero(); // Reset constant term 458 493 H.Zero(); // Reset H matrix 459 494 DW1D.Zero(); 460 495 // 461 //std::cout << "VertexFinder: start loop on tracks" << std::endl; 496 // Reset work arrays 497 // 498 ResetWrkArrays(); 499 // 500 // Start loop on tracks 501 // 462 502 for (Int_t i = 0; i < fNtr; i++) 463 503 { … … 465 505 TVectorD par = *fPar[i]; 466 506 TMatrixDSym Cov = *fCov[i]; 467 507 // 508 // For first iteration only 468 509 Double_t fs; 469 510 if (Ntry <= 0) // Initialize all phases on first pass … … 473 514 Double_t arg = TMath::Max(1.0e-6, (R * R - D * D) / (1 + 2 * C * D)); 474 515 fs = 2 * TMath::ASin(C * TMath::Sqrt(arg)); 475 f i[i] = fs;516 ffi.push_back(fs); 476 517 } 477 518 // 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; 519 // Update track related arrays 520 // 521 UpdateTrkArrays(i); 522 TMatrixDSym Ds = *fDi[i]; 523 TMatrixDSym Winv = *fWinvi[i]; 510 524 TMatrixDSym DsW1Ds = Winv.Similarity(Ds); // Service matrix to calculate covX 525 // 526 // Update global arrays 511 527 DW1D += DsW1Ds; 512 528 // Update hessian 513 529 H += Ds; 514 530 // update constant term 531 TVectorD xs = *fx0i[i]; 515 532 cterm += Ds * xs; 516 533 } // End loop on tracks 517 534 // 518 535 // update vertex position 519 TMatrixDSym H1 = RegInv 3(H);536 TMatrixDSym H1 = RegInv(H); 520 537 x = H1 * cterm; 521 // std::cout << "VertexFinder: x vertex set" << std::endl;538 // 522 539 // Update vertex covariance 523 540 covX = DW1D.Similarity(H1); 524 // std::cout << "VertexFinder: cov vertex set" << std::endl;541 // 525 542 // Update phases and chi^2 526 543 Chi2 = 0.0; 527 544 for (Int_t i = 0; i < fNtr; i++) 528 545 { 529 TVectorD lambda = (* Di[i]) * (*x0i[i] - x);530 TMatrixDSym Wm1 = * Winvi[i];546 TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x); 547 TMatrixDSym Wm1 = *fWinvi[i]; 531 548 fChi2List(i) = Wm1.Similarity(lambda); 532 549 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++)f i[i] += a(j) * b(j) /a2i[i];550 TVectorD a = *fai[i]; 551 TVectorD b = (*fWi[i]) * (x - (*fx0i[i])); 552 for (Int_t j = 0; j < 3; j++)ffi[i] += a(j) * b(j) / fa2i[i]; 536 553 } 537 538 554 // 539 555 TVectorD dx = x - x0; 540 556 x0 = x; 541 557 // update vertex stability 542 TMatrixDSym Hess = RegInv 3(covX);558 TMatrixDSym Hess = RegInv(covX); 543 559 epsi = Hess.Similarity(dx); 544 560 Ntry++; … … 546 562 // Store result 547 563 // 548 fXv = x; // Vertex position564 fXv = x; // Vertex position 549 565 fcovXv = covX; // Vertex covariance 550 566 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 // 567 } // end of iteration loop 568 // 569 fVtxDone = kTRUE; // Set fit completion flag 570 fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1)); // Store fit radius 571 // 572 } 573 // 574 // Return fit vertex 609 575 TVectorD VertexFit::GetVtx() 610 576 { 611 //std::cout << "GetVtx: flag set to " << fVtxDone << std::endl; 612 if (!fVtxDone)VertexFinder(); 577 if (!fVtxDone) VertexFitter(); 613 578 return fXv; 614 579 } 615 580 // 581 // Return fit vertex covariance 616 582 TMatrixDSym VertexFit::GetVtxCov() 617 583 { 618 if (!fVtxDone) VertexFinder();584 if (!fVtxDone) VertexFitter(); 619 585 return fcovXv; 620 586 } 621 587 // 588 // Return fit vertex chi2 622 589 Double_t VertexFit::GetVtxChi2() 623 590 { 624 if (!fVtxDone) VertexFinder();591 if (!fVtxDone) VertexFitter(); 625 592 return fChi2; 626 593 } 627 594 // 595 // Return array of chi2 contributions from each track 628 596 TVectorD VertexFit::GetVtxChi2List() 629 597 { 630 if (!fVtxDone) VertexFinder();598 if (!fVtxDone) VertexFitter(); 631 599 return fChi2List; 632 600 } … … 638 606 } 639 607 // 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 } 608 // Adding tracks one by one 609 void VertexFit::AddTrk(TVectorD *par, TMatrixDSym *Cov) // Add track to input list 610 { 611 fNtr++; 612 fChi2List.ResizeTo(fNtr); // Resize chi2 array 613 fPar.push_back(par); // add new track 614 fCov.push_back(Cov); 615 // 616 // Reset previous vertex temp arrays 617 ResetWrkArrays(); 618 ffi.clear(); 619 fVtxDone = kFALSE; // Reset vertex done flag 620 } 621 // 622 // Removing tracks one by one 623 void VertexFit::RemoveTrk(Int_t iTrk) // Remove iTrk track 624 { 625 fNtr--; 626 fChi2List.Clear(); 627 fChi2List.ResizeTo(fNtr); // Resize chi2 array 628 fPar.erase(fPar.begin() + iTrk); // Remove track 629 fCov.erase(fCov.begin() + iTrk); 630 // 631 // Reset previous vertex temp arrays 632 ResetWrkArrays(); 633 ffi.clear(); 634 fVtxDone = kFALSE; // Reset vertex done flag 635 }
Note:
See TracChangeset
for help on using the changeset viewer.
