Changes in external/TrackCovariance/VertexFit.cc [537d24f:127644a] in git
- File:
-
- 1 edited
-
external/TrackCovariance/VertexFit.cc (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/VertexFit.cc
r537d24f r127644a 46 46 fWinvi = new TMatrixDSym * [Ntr]; // ACA' 47 47 for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3); 48 48 49 49 } 50 50 // ObsTrk list … … 65 65 { 66 66 fPar[i] = new TVectorD(track[i]->GetObsPar()); 67 // cout << "fPar = " <<endl; fPar[i]->Print();67 //std::cout << "fPar = " << std::endl; fPar[i]->Print(); 68 68 fCov[i] = new TMatrixDSym(track[i]->GetCov()); 69 // cout << "fCov = " <<endl; fCov[i]->Print();69 //std::cout << "fCov = " << std::endl; fCov[i]->Print(); 70 70 } 71 71 // … … 83 83 for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3); 84 84 85 // cout << "VertexFit constructor executed" <<endl;85 //std::cout << "VertexFit constructor executed" << std::endl; 86 86 } 87 87 // … … 89 89 VertexFit::~VertexFit() 90 90 { 91 fxCst.Clear(); 91 fxCst.Clear(); 92 92 fCovCst.Clear(); 93 fXv.Clear(); 94 fcovXv.Clear(); 93 fXv.Clear(); 94 fcovXv.Clear(); 95 95 fChi2List.Clear(); 96 96 // 97 97 for (Int_t i= 0; i < fNtr; i++) 98 98 { 99 fPar[i]->Clear(); 99 fPar[i]->Clear(); 100 100 fCov[i]->Clear(); 101 101 // … … 195 195 if (Ntry > NtryMax) 196 196 { 197 cout << "FastRv: maximum number of iteration reached" <<endl;197 std::cout << "FastRv: maximum number of iteration reached" << std::endl; 198 198 break; 199 199 } … … 217 217 if (N != 3) 218 218 { 219 cout << "RegInv3 called with matrix size != 3. Abort & return standard inversion." <<endl;219 std::cout << "RegInv3 called with matrix size != 3. Abort & return standard inversion." << std::endl; 220 220 return Smat.Invert(); 221 221 } … … 263 263 else 264 264 { 265 cout << "RegInv3: found negative elements in diagonal. Return standard inversion." <<endl;265 std::cout << "RegInv3: found negative elements in diagonal. Return standard inversion." << std::endl; 266 266 return Smat.Invert(); 267 267 } … … 273 273 { 274 274 // 275 // Derivative of track 3D position vector with respect to track parameters at constant phase 275 // Derivative of track 3D position vector with respect to track parameters at constant phase 276 276 // 277 277 // par = vector of track parameters … … 395 395 // 396 396 // Stored quantities 397 Double_t *fi = new Double_t[fNtr]; // Phases 397 Double_t *fi = new Double_t[fNtr]; // Phases 398 398 TVectorD **x0i = new TVectorD*[fNtr]; // Track expansion point 399 399 TVectorD **ai = new TVectorD*[fNtr]; // dx/dphi … … 416 416 // Find track pair with largest phi difference 417 417 Int_t isel=0; Int_t jsel=0; // selected track indices 418 Double_t dphiMax = 0.0; // Max phi difference 418 Double_t dphiMax = 0.0; // Max phi difference 419 419 for (Int_t i = 0; i < fNtr-1; i++) 420 420 { … … 438 438 } 439 439 // 440 // 440 // 441 441 //ObsTrk* t1 = tracks[isel]; 442 442 TVectorD p1 = *fPar[isel]; … … 447 447 R = 0.9*R + 0.1*Rd; 448 448 // 449 // cout << "VertexFinder: fast R = " << R <<endl;449 //std::cout << "VertexFinder: fast R = " << R << std::endl; 450 450 // 451 451 // Iteration properties … … 464 464 H.Zero(); // Reset H matrix 465 465 DW1D.Zero(); 466 // 467 // cout << "VertexFinder: start loop on tracks" <<endl;466 // 467 //std::cout << "VertexFinder: start loop on tracks" << std::endl; 468 468 for (Int_t i = 0; i < fNtr; i++) 469 469 { 470 // Get track helix parameters and their covariance matrix 470 // Get track helix parameters and their covariance matrix 471 471 //ObsTrk *t = tracks[i]; 472 472 TVectorD par = *fPar[i]; 473 TMatrixDSym Cov = *fCov[i]; 473 TMatrixDSym Cov = *fCov[i]; 474 474 Double_t fs; 475 475 if (Ntry <= 0) // Initialize all phases on first pass … … 485 485 // 486 486 fs = fi[i]; // Get phase 487 // cout << "VertexFinder: phase fs set" <<endl;487 //std::cout << "VertexFinder: phase fs set" << std::endl; 488 488 TVectorD xs = Fill_x(par, fs); 489 // cout << "VertexFinder: position xs set" <<endl;489 //std::cout << "VertexFinder: position xs set" << std::endl; 490 490 x0i[i] = new TVectorD(xs); // Start helix position 491 // cout << "VertexFinder: position x0i stored" <<endl;491 //std::cout << "VertexFinder: position x0i stored" << std::endl; 492 492 // W matrix = (A*C*A')^-1; W^-1 = A*C*A' 493 493 TMatrixD A = Fill_A(par, fs); // A = dx/da = derivatives wrt track parameters 494 // cout << "VertexFinder: derivatives A set" <<endl;494 //std::cout << "VertexFinder: derivatives A set" << std::endl; 495 495 TMatrixDSym Winv = Cov.Similarity(A); // W^-1 = A*C*A' 496 496 Winvi[i] = new TMatrixDSym(Winv); // Store W^-1 matrix 497 // cout << "VertexFinder: Winvi stored" <<endl;497 //std::cout << "VertexFinder: Winvi stored" << std::endl; 498 498 TMatrixDSym W = RegInv3(Winv); // W = (A*C*A')^-1 499 499 Wi[i] = new TMatrixDSym(W); // Store W matrix 500 // cout << "VertexFinder: Wi stored" <<endl;500 //std::cout << "VertexFinder: Wi stored" << std::endl; 501 501 TVectorD a = Fill_a(par, fs); // a = dx/ds = derivatives wrt phase 502 // cout << "VertexFinder: derivatives a set" <<endl;502 //std::cout << "VertexFinder: derivatives a set" << std::endl; 503 503 ai[i] = new TVectorD(a); // Store a 504 // cout << "VertexFinder: derivatives a stored" <<endl;504 //std::cout << "VertexFinder: derivatives a stored" << std::endl; 505 505 Double_t a2 = W.Similarity(a); 506 506 a2i[i] = a2; // Store a2 507 507 // Build D matrix 508 TMatrixDSym B(3); 508 TMatrixDSym B(3); 509 509 B.Rank1Update(a, 1.0); 510 510 B *= -1. / a2; … … 512 512 TMatrixDSym Ds = W+B; // D matrix 513 513 Di[i] = new TMatrixDSym(Ds); // Store D matrix 514 // cout << "VertexFinder: matrix Di stored" <<endl;514 //std::cout << "VertexFinder: matrix Di stored" << std::endl; 515 515 TMatrixDSym DsW1Ds = Winv.Similarity(Ds); // Service matrix to calculate covX 516 DW1D += DsW1Ds; 516 DW1D += DsW1Ds; 517 517 // Update hessian 518 518 H += Ds; … … 524 524 TMatrixDSym H1 = RegInv3(H); 525 525 x = H1*cterm; 526 // cout << "VertexFinder: x vertex set" <<endl;526 //std::cout << "VertexFinder: x vertex set" << std::endl; 527 527 // Update vertex covariance 528 528 covX = DW1D.Similarity(H1); 529 // cout << "VertexFinder: cov vertex set" <<endl;529 //std::cout << "VertexFinder: cov vertex set" << std::endl; 530 530 // Update phases and chi^2 531 531 Chi2 = 0.0; … … 541 541 } 542 542 543 // cout << "VertexFinder: end chi2 calculation" <<endl;543 //std::cout << "VertexFinder: end chi2 calculation" << std::endl; 544 544 // 545 545 TVectorD dx = x - x0; … … 559 559 // 560 560 561 // cout << "VertexFinder: before store intermediate data" <<endl;561 //std::cout << "VertexFinder: before store intermediate data" << std::endl; 562 562 for (Int_t i = 0; i < fNtr; i++) 563 563 { 564 // cout << "VertexFinder: inside store intermediate data" <<endl;565 // cout << "i = " << i << ", fi[i] = " << fi[i] <<endl;566 // cout << "i = " << i << ", ffi[i] = " << ffi[i] <<endl;564 //std::cout << "VertexFinder: inside store intermediate data" << std::endl; 565 //std::cout << "i = " << i << ", fi[i] = " << fi[i] << std::endl; 566 //std::cout << "i = " << i << ", ffi[i] = " << ffi[i] << std::endl; 567 567 ffi[i] = fi[i]; // Fit phases 568 // cout << "VertexFinder: fi stored" <<endl;568 //std::cout << "VertexFinder: fi stored" << std::endl; 569 569 fx0i[i] = x0i[i]; // Track expansion points 570 // cout << "VertexFinder: x0i stored" <<endl;570 //std::cout << "VertexFinder: x0i stored" << std::endl; 571 571 fai[i] = ai[i]; // dx/dphi 572 // cout << "VertexFinder: ai stored" <<endl;572 //std::cout << "VertexFinder: ai stored" << std::endl; 573 573 fa2i[i] = a2i[i]; // a'Wa 574 // cout << "VertexFinder: a2i stored" <<endl;574 //std::cout << "VertexFinder: a2i stored" << std::endl; 575 575 fDi[i] = Di[i]; // W-WBW 576 // cout << "VertexFinder: Di stored" <<endl;576 //std::cout << "VertexFinder: Di stored" << std::endl; 577 577 fWi[i] = Wi[i]; // (ACA')^-1 578 // cout << "VertexFinder: Wi stored" <<endl;578 //std::cout << "VertexFinder: Wi stored" << std::endl; 579 579 fWinvi[i] = Winvi[i]; // ACA' 580 // cout << "VertexFinder: Winvi stored" <<endl;580 //std::cout << "VertexFinder: Winvi stored" << std::endl; 581 581 } 582 // cout << "Iteration " << Ntry << " completed - Before cleanup" <<endl;582 //std::cout << "Iteration " << Ntry << " completed - Before cleanup" << std::endl; 583 583 // 584 584 // Cleanup … … 599 599 } 600 600 601 // cout << "Iteration " << Ntry << " completed - After cleanup" <<endl;601 //std::cout << "Iteration " << Ntry << " completed - After cleanup" << std::endl; 602 602 } 603 603 // 604 604 fVtxDone = kTRUE; // Set fitting completion flag 605 605 // 606 delete[] fi; // Phases 606 delete[] fi; // Phases 607 607 delete[] x0i; // Track expansion point 608 608 delete[] ai; // dx/dphi … … 615 615 TVectorD VertexFit::GetVtx() 616 616 { 617 // cout << "GetVtx: flag set to " << fVtxDone <<endl;617 //std::cout << "GetVtx: flag set to " << fVtxDone << std::endl; 618 618 if (!fVtxDone)VertexFinder(); 619 619 return fXv; … … 641 641 void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov) // Add gaussian vertex constraint 642 642 { 643 cout << "VertexFit::AddVtxConstraint: Not implemented yet" <<endl;643 std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl; 644 644 } 645 645 // 646 646 void VertexFit::AddTrk(TVectorD par, TMatrixDSym Cov) // Add track to input list 647 647 { 648 cout << "VertexFit::AddTrk: Not implemented yet" <<endl;648 std::cout << "VertexFit::AddTrk: Not implemented yet" << std::endl; 649 649 } 650 650 void VertexFit::RemoveTrk(Int_t iTrk) // Remove iTrk track 651 651 { 652 cout << "VertexFit::RemoveTrk: Not implemented yet" <<endl;653 } 652 std::cout << "VertexFit::RemoveTrk: Not implemented yet" << std::endl; 653 }
Note:
See TracChangeset
for help on using the changeset viewer.
