Changes in / [78ce8d1:0c0c9af] in git
- Files:
-
- 6 edited
-
examples/ExamplePVtxFind.C (modified) (16 diffs)
-
examples/ExamplePVtxFitEC.C (modified) (12 diffs)
-
external/TrackCovariance/TrkUtil.h (modified) (2 diffs)
-
external/TrackCovariance/VertexFit.cc (modified) (16 diffs)
-
external/TrackCovariance/VertexFit.h (modified) (5 diffs)
-
modules/TreeWriter.cc (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
examples/ExamplePVtxFind.C
r78ce8d1 r0c0c9af 1 1 /* 2 2 Example of using vertex fitter class to fit primary vertex 3 assumed to be generated with Pythia8/ee_zh_smr-shf.cmn3 assumed to be generated in (0,0,0) 4 4 */ 5 5 … … 44 44 TH1D* hTrFoundG = new TH1D("hTrFoundG", "Found GOOD primary tracks", 41, -0.5, 40.5); 45 45 TH1D* hTrFoundB = new TH1D("hTrFoundB", "Found BAD primary tracks", 41, -0.5, 40.5); 46 TH1D* hTrDiff = new TH1D("hTrDiff", "Found - available tracks", 21, -10.5, 10.5);47 46 // 48 47 // Loop over all events … … 88 87 pr[it] = new TVectorD(obsPar); 89 88 cv[it] = new TMatrixDSym(trk->CovarianceMatrix()); 90 // Check covariance91 /*92 TMatrixDSymEigen Eign(*cv[it]);93 TVectorD lambda = Eign.GetEigenValues();94 Bool_t good = kTRUE;95 for(Int_t i=0; i<5; i++) if(lambda(i)<0.0) good = kFALSE;96 if(!good){97 std::cout<<"Cov["<<it<<"]"<<" eigenvalues: "98 <<lambda(0)<<", "<<lambda(1)<<", "<<lambda(2)<<", "99 <<lambda(3)<<", "<<lambda(4)<<std::endl;100 TMatrixDSym Ctmp = *cv[it];101 std::cout<<"Error D (dir - cov) "<<trk->ErrorD0102 <<" - "<<TMath::Sqrt(Ctmp(0,0))<<std::endl;103 }104 */105 89 //std::cout << "Track loaded it= " << it << std::endl; 106 90 // … … 136 120 } // End loop on tracks 137 121 } 138 std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl;139 122 // 140 123 // Find primary vertex … … 149 132 covpvc(1, 1) = 2.55e-05 * 2.55e-05; 150 133 covpvc(2, 2) = 0.64 * 0.64; 151 if(Nprim == 0){152 xpv = xpvc(0);153 ypv = xpvc(1);154 zpv = xpvc(2);155 }156 134 // 157 135 // … … 159 137 Int_t nSkim = 0; 160 138 Int_t* nSkimmed = new Int_t[NtrG]; 161 TVectorD** PrSk = new TVectorD * [1];162 TMatrixDSym** CvSk = new TMatrixDSym * [1];163 Double_t MaxChi2 = 2.0+3* 2.;139 VertexFit* Vskim = new VertexFit(); 140 Vskim->AddVtxConstraint(xpvc, covpvc); 141 Double_t MaxChi2 = 2.0+3*sqrt(2.); 164 142 for (Int_t n = 0; n < NtrG; n++) { 165 PrSk[0] = new TVectorD(*pr[n]); 166 CvSk[0] = new TMatrixDSym(*cv[n]); 167 VertexFit* Vskim = new VertexFit(1,PrSk, CvSk); 168 Vskim->AddVtxConstraint(xpvc, covpvc); 143 Vskim->AddTrk(pr[n], cv[n]); 169 144 Double_t Chi2One = Vskim->GetVtxChi2(); 170 145 //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl; … … 174 149 nSkim++; 175 150 } 176 // Cleanup 177 delete Vskim; 178 } 179 delete PrSk[0]; 180 delete CvSk[0]; 181 delete[] PrSk; 182 delete[] CvSk; 183 // 151 Vskim->RemoveTrk(0); 152 } 184 153 // Load tracks for primary fit 185 154 Int_t MinTrk = 1; // Minumum # tracks for vertex fit 155 TVectorD** PrFit = new TVectorD * [nSkim]; 156 TMatrixDSym** CvFit = new TMatrixDSym * [nSkim]; 186 157 if (nSkim >= MinTrk) { 187 TVectorD** PrFit = new TVectorD * [nSkim];188 TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];189 158 for (Int_t n = 0; n < nSkim; n++) { 190 159 PrFit[n] = new TVectorD(*pr[nSkimmed[n]]); 191 160 CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]); 192 } 193 delete[] nSkimmed; 194 Int_t Nfound = nSkim; 195 if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl; 196 const Int_t MaxFound = 100; Double_t Chi2LL[MaxFound]; Double_t *Chi2L = Chi2LL; 161 TVectorD test = *PrFit[n]; 162 //std::cout << "Test *PrFit[" << n << "](0) = " << test(0) << std::endl; 163 } 164 } 165 delete[] nSkimmed; 166 delete Vskim; 167 Int_t Nfound = nSkim; 168 if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl; 169 if (nSkim >= MinTrk) { 197 170 VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit); 198 171 //std::cout << "Vertex fit created " << std::endl; … … 200 173 // 201 174 // Remove tracks with large chi2 202 Double_t MaxChi2Fit = 4.5;175 Double_t MaxChi2Fit = 9.; 203 176 Bool_t Done = kFALSE; 204 177 while (!Done) { … … 207 180 TVectorD Chi2List = Vtx->GetVtxChi2List(); // Get contributions to Chi2 208 181 //std::cout << "After Chi2List. " << std::endl; Chi2List.Print(); 209 //Double_t* Chi2L = new Double_t[Nfound];182 Double_t* Chi2L = new Double_t[Nfound]; 210 183 Chi2L = Chi2List.GetMatrixArray(); 211 184 Int_t iMax = TMath::LocMax(Nfound, Chi2L); … … 213 186 Double_t Chi2Mx = Chi2L[iMax]; 214 187 //std::cout << "Chi2Mx "<<Chi2Mx << std::endl; 215 if (Chi2Mx > MaxChi2Fit && Nfound > 1) {188 if (Chi2Mx > MaxChi2Fit && Nfound > 2) { 216 189 //std::cout << "Before remove. Nfound = "<<Nfound << std::endl; 217 190 Vtx->RemoveTrk(iMax); … … 222 195 Done = kTRUE; 223 196 } 197 //std::cout << "Done = " << Done << std::endl; 198 //delete[] Chi2L; 199 //std::cout << "Array Chi2L removed " << std::endl; 224 200 } 225 201 // … … 227 203 // 228 204 // Require minimum number of tracks in vertex 229 Int_t Nmin = 1;205 Int_t Nmin = 2; 230 206 if (Nfound >= Nmin) { 231 207 TVectorD xvtx = Vtx->GetVtx(); 232 std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;208 //std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl; 233 209 TMatrixDSym covX = Vtx->GetVtxCov(); 234 210 Double_t Chi2 = Vtx->GetVtxChi2(); … … 246 222 hTrPrim->Fill(Nprim); 247 223 hTrFound->Fill((Double_t)Nfound); 248 hTrDiff->Fill((Double_t)Nfound-Nprim);249 224 //std::cout << "Histograms filled " << std::endl; 250 225 } 251 226 // 252 // Clean253 227 delete Vtx; 254 for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];255 for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];256 delete[] PrFit;257 delete[] CvFit;258 228 } 259 229 … … 263 233 for (Int_t i = 0; i < NtrG; i++) delete pr[i]; 264 234 for (Int_t i = 0; i < NtrG; i++) delete cv[i]; 235 for (Int_t i = 0; i < nSkim; i++) delete PrFit[i]; 236 for (Int_t i = 0; i < nSkim; i++) delete CvFit[i]; 265 237 delete[] pr; 266 238 delete[] cv; 239 delete[] PrFit; 240 delete[] CvFit; 267 241 } 268 242 // … … 280 254 // 281 255 TCanvas* CnvN = new TCanvas("CnvN", "Primary tracks found", 100, 100, 900, 500); 282 CnvN->Divide(2, 2);256 CnvN->Divide(2, 1); 283 257 CnvN->cd(1); 284 258 hTrPrim->Draw(); … … 286 260 hTrFound->SetLineColor(kRed); 287 261 hTrFound->Draw(); 288 CnvN->cd(3);289 hTrDiff->Draw();290 262 } -
examples/ExamplePVtxFitEC.C
r78ce8d1 r0c0c9af 1 1 /* 2 2 Example of using vertex fitter class to fit primary vertex 3 assumed to be generated with Pythia8/ee_zh_smr-shf.cmn3 assumed to be generated in (0,0,0) 4 4 */ 5 5 … … 21 21 void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5) 22 22 { 23 //24 // Beam constraint25 // (consistent with generation with ee_zh_smr-shf.cmnd)26 //27 // Mean beam position28 TVectorD xpvc(3);29 xpvc(0) = 1.0;30 xpvc(1) = -2.0;31 xpvc(2) = 10.0;32 // Interaction region covariance33 TMatrixDSym covpvc(3); covpvc.Zero();34 covpvc(0, 0) = 0.0097 * 0.0097;35 covpvc(1, 1) = 2.55e-05 * 2.55e-05;36 covpvc(2, 2) = 0.64 * 0.64;37 38 //39 23 // Create chain of root trees 40 24 TChain chain("Delphes"); … … 51 35 // Book histograms 52 36 Int_t Nbin = 100; 53 // Vertex fit pulls54 37 TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.); 55 38 TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.); 56 39 TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.); 57 40 TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.); 58 // Generation check59 TH1D* hXvpull = new TH1D("hXvpull", "Pull X generated vertex", Nbin, -10., 10.);60 TH1D* hYvpull = new TH1D("hYvpull", "Pull Y generated vertex", Nbin, -10., 10.);61 TH1D* hZvpull = new TH1D("hZvpull", "Pull Z generated vertex", Nbin, -10., 10.);62 63 41 // 64 42 // Loop over all events … … 87 65 // Get associated generated particle 88 66 GenParticle* gp = (GenParticle*)trk->Particle.GetObject(); 89 TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z); // mm -> meters67 TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z); 90 68 TVector3 pg(gp->Px,gp->Py,gp->Pz); 91 69 Double_t Q = (Double_t)gp->Charge; 92 70 Double_t Bz = 2.0; 93 71 TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz); 94 TVectorD genPar = TrkUtil::ParToMm(genParM); // -> back to mm 95 96 // 97 // Check if original track is from primary vertex 72 TVectorD genPar = TrkUtil::ParToMm(genParM); 98 73 // 99 74 // Position of origin in mm … … 104 79 Int_t mp = gp->M1; // Mother 105 80 while(mp>0){ 106 GenParticle* gm = (GenParticle*)branchGenPart->At(mp); 81 GenParticle* gm = 82 (GenParticle*)branchGenPart->At(mp); 107 83 Double_t xm = gm->X; 108 84 Double_t ym = gm->Y; 109 85 Double_t zm = gm->Z; 110 86 if(x!=xm || y!=ym || z!=zm){ 111 prim = kFALSE; // Not primary87 prim = kFALSE; 112 88 break; 113 89 }else mp = gm->M1; 114 90 } 115 116 91 // 117 92 // group tracks originating from the primary vertex 118 93 if (prim) 119 94 { 120 // 121 // Store true primary vertex for this event95 //std::cout<<"Event: "<<entry<<", track: "<<it; 96 //std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl; 122 97 xpv = x; 123 98 ypv = y; … … 129 104 Double_t obsC = trk->C; 130 105 Double_t obsZ0 = trk->DZ; 106 //std::cout<<"Z0 track = "<< obsZ0 107 //<<", gen Z0 = "<<genPar(3) 108 //<<", gen cot = "<<genPar(4)<<std::endl; 131 109 Double_t obsCtg = trk->CtgTheta; 132 110 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg }; … … 134 112 // 135 113 pr[Ntr] = new TVectorD(obsPar); 114 //std::cout<<"Cov Matrix:"<<std::endl; 115 //trk->CovarianceMatrix().Print(); 136 116 cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix()); 137 117 Ntr++; … … 142 122 // Fit primary vertex with beam constraint 143 123 // 144 Int_t MinTrk = 3; // Minumum # tracks for vertex fit 124 //Beam constraint 125 TVectorD xpvc(3); 126 xpvc(0) = 1.0; 127 xpvc(1) = -2.0; 128 xpvc(2) = 10.0; 129 TMatrixDSym covpvc(3); covpvc.Zero(); 130 covpvc(0, 0) = 0.0097 * 0.0097; 131 covpvc(1, 1) = 2.55e-05 * 2.55e-05; 132 covpvc(2, 2) = 0.64 * 0.64; 133 Int_t MinTrk = 2; // Minumum # tracks for vertex fit 145 134 if (Ntr >= MinTrk) { 146 135 VertexFit* Vtx = new VertexFit(Ntr, pr, cv); … … 150 139 Double_t Chi2 = Vtx->GetVtxChi2(); 151 140 Double_t Ndof = 2 * (Double_t)Ntr; 152 delete Vtx;153 //154 141 Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0)); 155 142 Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1)); 156 143 Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2)); 157 // 158 Double_t PullXv = (xpvc(0)-xpv) / TMath::Sqrt(covpvc(0, 0)); 159 Double_t PullYv = (xpvc(1)-ypv) / TMath::Sqrt(covpvc(1, 1)); 160 Double_t PullZv = (xpvc(2)-zpv) / TMath::Sqrt(covpvc(2, 2)); 161 // 144 //std::cout<<"**** True vertex (x, y, z) "<<xpv<<", " 145 //<<ypv<<", "<<zpv<<std::endl; 146 //std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", " 147 //<<xvtx(1)<<", "<<xvtx(2)<<std::endl; 162 148 // 163 149 // Fill histograms … … 166 152 hZpull->Fill(PullZ); 167 153 hChi2->Fill(Chi2 / Ndof); 168 //169 hXvpull->Fill(PullXv);170 hYvpull->Fill(PullYv);171 hZvpull->Fill(PullZv);172 154 } 173 155 156 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl; 174 157 // 175 158 // Cleanup … … 179 162 delete[] cv; 180 163 } 181 182 164 // 183 165 // Show resulting histograms … … 192 174 hZpull->Fit("gaus"); hZpull->Draw(); 193 175 Cnv->cd(4); hChi2->Draw(); 194 //195 TCanvas* Cnv1 = new TCanvas("Cnv1", "Generated primary vertex pulls", 100, 100, 900, 500);196 Cnv1->Divide(3, 1);197 Cnv1->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);198 hXvpull->Fit("gaus"); hXvpull->Draw();199 Cnv1->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);200 hYvpull->Fit("gaus"); hYvpull->Draw();201 Cnv1->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);202 hZvpull->Fit("gaus"); hZvpull->Draw();203 176 } -
external/TrackCovariance/TrkUtil.h
r78ce8d1 r0c0c9af 68 68 // 69 69 static TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q, Double_t Bz); 70 static TVector3 ParToX(TVectorD Par); // position of minimum distance from z axis70 static TVector3 ParToX(TVectorD Par); // position of minimum distance from z axis 71 71 static TVector3 ParToP(TVectorD Par, Double_t Bz); // Get Momentum from track parameters 72 static Double_t ParToQ(TVectorD Par); // Get track charge72 static Double_t ParToQ(TVectorD Par); // Get track charge 73 73 static void LineDistance(TVector3 x0, TVector3 y0, TVector3 dirx, TVector3 diry, Double_t &sx, Double_t &sy, Double_t &distance); 74 74 // 75 75 // Track trajectory 76 76 // 77 static TVector3 Xtrack(TVectorD par, Double_t s); // Parametric track trajectory78 TVectorD derRphi_R(TVectorD par, Double_t R); 77 static TVector3 Xtrack(TVectorD par, Double_t s); // Parametric track trajectory 78 TVectorD derRphi_R(TVectorD par, Double_t R); // Derivatives of R-phi at constant R 79 79 TVectorD derZ_R(TVectorD par, Double_t R); // Derivatives of z at constant R 80 TVectorD derRphi_Z(TVectorD par, Double_t z); 80 TVectorD derRphi_Z(TVectorD par, Double_t z); // Derivatives of R-phi at constant z 81 81 TVectorD derR_Z(TVectorD par, Double_t z); // Derivatives of R at constant z 82 82 // … … 88 88 // 89 89 static TVectorD ParToMm(TVectorD Par); // Parameter conversion 90 static TMatrixDSym CovToMm(TMatrixDSym Cov); 90 static TMatrixDSym CovToMm(TMatrixDSym Cov); // Covariance conversion 91 91 // 92 92 // Inside cylindrical volume -
external/TrackCovariance/VertexFit.cc
r78ce8d1 r0c0c9af 1 /* 2 Vertex fitting code 3 */ 4 #include <TMath.h> 1 #include <TMath.h> 5 2 #include <TVectorD.h> 6 3 #include <TVector3.h> 7 4 #include <TMatrixD.h> 8 5 #include <TMatrixDSym.h> 9 #include <TMatrixDSymEigen.h>10 6 #include "VertexFit.h" 11 7 // … … 42 38 for (Int_t i = 0; i < fNtr; i++) 43 39 { 44 fPar.push_back(new TVectorD(*trkPar[i])); 45 fParNew.push_back(new TVectorD(*trkPar[i])); 46 fCov.push_back(new TMatrixDSym(*trkCov[i])); 47 fCovNew.push_back(new TMatrixDSym(*trkCov[i])); 40 TVectorD pr = *trkPar[i]; 41 fPar.push_back(new TVectorD(pr)); 42 fParNew.push_back(new TVectorD(pr)); 43 TMatrixDSym cv = *trkCov[i]; 44 fCov.push_back(new TMatrixDSym(cv)); 45 fCovNew.push_back(new TMatrixDSym(cv)); 48 46 } 49 47 fChi2List.ResizeTo(fNtr); … … 78 76 void VertexFit::ResetWrkArrays() 79 77 { 80 Int_t N = (Int_t)fdi.size(); 81 if(N > 0){ 82 for (Int_t i = 0; i < N; i++) 83 { 84 if (fx0i[i]) { fx0i[i]->Clear(); delete fx0i[i]; } 85 if (fai[i]) { fai[i]->Clear(); delete fai[i]; } 86 if (fdi[i]) { fdi[i]->Clear(); delete fdi[i]; } 87 if (fAti[i]) { fAti[i]->Clear(); delete fAti[i]; } 88 if (fDi[i]) { fDi[i]->Clear(); delete fDi[i]; } 89 if (fWi[i]) { fWi[i]->Clear(); delete fWi[i]; } 90 if (fWinvi[i]){ fWinvi[i]->Clear(); delete fWinvi[i]; } 91 } 92 fa2i.clear(); 93 fx0i.clear(); 94 fai.clear(); 95 fdi.clear(); 96 fAti.clear(); 97 fDi.clear(); 98 fWi.clear(); 99 fWinvi.clear(); 100 } 78 Int_t N = (Int_t)ffi.size(); 79 for (Int_t i = 0; i < N; i++) 80 { 81 if (fx0i[i]) { fx0i[i]->Clear(); delete fx0i[i]; } 82 if (fai[i]) { fai[i]->Clear(); delete fai[i]; } 83 if (fdi[i]) { fdi[i]->Clear(); delete fdi[i]; } 84 if (fAti[i]) { fAti[i]->Clear(); delete fAti[i]; } 85 if (fDi[i]) { fDi[i]->Clear(); delete fDi[i]; } 86 if (fWi[i]) { fWi[i]->Clear(); delete fWi[i]; } 87 if (fWinvi[i]){ fWinvi[i]->Clear(); delete fWinvi[i]; } 88 } 89 fa2i.clear(); 90 fx0i.clear(); 91 fai.clear(); 92 fdi.clear(); 93 fAti.clear(); 94 fDi.clear(); 95 fWi.clear(); 96 fWinvi.clear(); 101 97 } 102 98 VertexFit::~VertexFit() 103 { 99 { 104 100 fxCst.Clear(); 105 101 fCovCst.Clear(); … … 107 103 fXv.Clear(); 108 104 fcovXv.Clear(); 109 fChi2List.Clear(); 105 fChi2List.Clear(); 110 106 // 111 107 for (Int_t i = 0; i < fNtr; i++) … … 115 111 fCov[i]->Clear(); delete fCov[i]; 116 112 fCovNew[i]->Clear(); delete fCovNew[i]; 117 } 113 } 118 114 fPar.clear(); 119 115 fParNew.clear(); … … 122 118 // 123 119 ResetWrkArrays(); 124 ffi.clear(); 120 ffi.clear(); 125 121 fNtr = 0; 126 122 } 127 123 // 128 // 129 // 130 TVectorD VertexFit::Fill_x0(TVectorD par) 131 { 132 // 133 // Calculate track 3D position at R = |D| (minimum approach to z-axis) 134 // 135 TVectorD x0(3); 124 Double_t VertexFit::FastRv(TVectorD p1, TVectorD p2) 125 { 126 // 127 // Find radius of minimum distance between two tracks 128 // 129 // p = (D,phi, C, z0, ct) 130 // 131 // Define arrays 132 // 133 Double_t C1 = p1(2); 134 Double_t C2 = p2(2); 135 Double_t ph1 = p1(1); 136 Double_t ph2 = p2(1); 137 TVectorD x0 = Fill_x0(p1); 138 TVectorD y0 = Fill_x0(p2); 139 TVectorD n = Fill_a(p1, 0.0); 140 n *= (2*C1); 141 TVectorD k = Fill_a(p2, 0.0); 142 k *= (2*C2); 143 // 144 // Setup and solve linear system 145 // 146 Double_t nn = 0; for (Int_t i = 0; i < 3; i++)nn += n(i) * n(i); 147 Double_t nk = 0; for (Int_t i = 0; i < 3; i++)nk += n(i) * k(i); 148 Double_t kk = 0; for (Int_t i = 0; i < 3; i++)kk += k(i) * k(i); 149 Double_t discr = nn * kk - nk * nk; 150 TMatrixDSym H(2); 151 H(0, 0) = kk; 152 H(0, 1) = nk; 153 H(1, 0) = H(0, 1); 154 H(1, 1) = nn; 155 TVectorD c(2); 156 c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) += n(i) * (y0(i) - x0(i)); 157 c(1) = 0; for (Int_t i = 0; i < 3; i++)c(1) += -k(i) * (y0(i) - x0(i)); 158 TVectorD smin = (H * c); 159 smin *= 1.0 / discr; 160 // 161 TVectorD X = x0 + smin(0) * n; 162 TVectorD Y = y0 + smin(1) * k; 163 // 164 // Higher order corrections 165 X(0) += -C1 * smin(0) * smin(0) * TMath::Sin(ph1); 166 X(1) += C1 * smin(0) * smin(0) * TMath::Cos(ph1); 167 Y(0) += -C2 * smin(1) * smin(1) * TMath::Sin(ph2); 168 Y(1) += C2 * smin(1) * smin(1) * TMath::Cos(ph2); 169 // 170 TVectorD Xavg = 0.5 * (X + Y); 171 // 172 // 173 return TMath::Sqrt(Xavg(0)*Xavg(0)+Xavg(1)*Xavg(1)); 174 } 175 // 176 // Starting radius determination 177 Double_t VertexFit::StartRadius() 178 { 179 // 180 // Maximum impact parameter 181 Double_t Rd = 0; 182 for (Int_t i = 0; i < fNtr; i++) 183 { 184 TVectorD par = *fPar[i]; 185 Double_t Dabs = TMath::Abs(par(0)); 186 if (Dabs > Rd)Rd = Dabs; 187 } 188 //----------------------------- 189 // 190 // Find track pair with phi difference closest to pi/2 191 Int_t isel = 0; Int_t jsel = 0; // selected track indices 192 Double_t dSinMax = 0.0; // Max phi difference 193 for (Int_t i = 0; i < fNtr - 1; i++) 194 { 195 TVectorD pari = *fPar[i]; 196 Double_t phi1 = pari(1); 197 198 for (Int_t j = i + 1; j < fNtr; j++) 199 { 200 TVectorD parj = *fPar[j]; 201 Double_t phi2 = parj(1); 202 Double_t Sindphi = TMath::Abs(TMath::Sin(phi2 - phi1)); 203 if (Sindphi > dSinMax) 204 { 205 isel = i; jsel = j; 206 dSinMax = Sindphi; 207 } 208 } 209 } 210 // 211 //------------------------------------------ 212 // 213 // Find radius of minimum distrance between tracks 214 TVectorD p1 = *fPar[isel]; 215 TVectorD p2 = *fPar[jsel]; 216 Double_t R = FastRv(p1, p2); 217 // 218 R = 0.9 * R + 0.1 * Rd; // Protect for overshoot 219 // 220 return R; 221 } 222 // 223 // Regularized symmetric matrix inversion 224 // 225 TMatrixDSym VertexFit::RegInv(TMatrixDSym& Min) 226 { 227 TMatrixDSym M = Min; // Decouple from input 228 Int_t N = M.GetNrows(); // Matrix size 229 TMatrixDSym D(N); D.Zero(); // Normaliztion matrix 230 TMatrixDSym R(N); // Normarized matrix 231 TMatrixDSym Rinv(N); // Inverse of R 232 TMatrixDSym Minv(N); // Inverse of M 233 // 234 // Check for 0's and normalize 235 for (Int_t i = 0; i < N; i++) 236 { 237 if (M(i, i) != 0.0) D(i, i) = 1. / TMath::Sqrt(TMath::Abs(M(i, i))); 238 else D(i, i) = 1.0; 239 } 240 R = M.Similarity(D); 241 // 242 // Recursive algorithms stops when N = 2 243 // 244 //**************** 245 // case N = 2 *** 246 //**************** 247 if (N == 2) 248 { 249 Double_t det = R(0, 0) * R(1, 1) - R(0, 1) * R(1, 0); 250 if (det == 0) 251 { 252 std::cout << "VertexFit::RegInv: null determinant for N = 2" << std::endl; 253 Rinv.Zero(); // Return null matrix 254 } 255 else 256 { 257 // invert matrix 258 Rinv(0, 0) = R(1, 1); 259 Rinv(0, 1) = -R(0, 1); 260 Rinv(1, 0) = Rinv(0, 1); 261 Rinv(1, 1) = R(0, 0); 262 Rinv *= 1. / det; 263 } 264 } 265 //**************** 266 // case N > 2 *** 267 //**************** 268 else 269 { 270 // Break up matrix 271 TMatrixDSym Q = R.GetSub(0, N - 2, 0, N - 2); // Upper left 272 TVectorD p(N - 1); 273 for (Int_t i = 0; i < N - 1; i++)p(i) = R(N - 1, i); 274 Double_t q = R(N - 1, N - 1); 275 //Invert pieces and re-assemble 276 TMatrixDSym Ainv(N - 1); 277 TMatrixDSym A(N - 1); 278 if (TMath::Abs(q) > 1.0e-15) 279 { 280 // Case |q| > 0 281 Ainv.Rank1Update(p, -1.0 / q); 282 Ainv += Q; 283 A = RegInv(Ainv); // Recursive call 284 TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A; 285 // 286 TVectorD b = (-1.0 / q) * (A * p); 287 for (Int_t i = 0; i < N - 1; i++) 288 { 289 Rinv(N - 1, i) = b(i); 290 Rinv(i, N - 1) = b(i); 291 } 292 // 293 Double_t pdotb = 0.; 294 for (Int_t i = 0; i < N - 1; i++)pdotb += p(i) * b(i); 295 Double_t c = (1.0 - pdotb) / q; 296 Rinv(N - 1, N - 1) = c; 297 } 298 else 299 { 300 // case q = 0 301 TMatrixDSym Qinv = RegInv(Q); // Recursive call 302 Double_t a = Qinv.Similarity(p); 303 Double_t c = -1.0 / a; 304 Rinv(N - 1, N - 1) = c; 305 // 306 TVectorD b = (1.0 / a) * (Qinv * p); 307 for (Int_t i = 0; i < N - 1; i++) 308 { 309 Rinv(N - 1, i) = b(i); 310 Rinv(i, N - 1) = b(i); 311 } 312 // 313 A.Rank1Update(p, -1 / a); 314 A += Q; 315 A.Similarity(Qinv); 316 TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A; 317 } 318 } 319 Minv = Rinv.Similarity(D); 320 return Minv; 321 } 322 // 323 // 324 // 325 TMatrixD VertexFit::Fill_A(TVectorD par, Double_t phi) 326 { 327 // 328 // Derivative of track 3D position vector with respect to track parameters at constant phase 329 // 330 // par = vector of track parameters 331 // phi = phase 332 // 333 TMatrixD A(3, 5); 136 334 // 137 335 // Decode input arrays … … 143 341 Double_t ct = par(4); 144 342 // 145 x0(0) = -D * TMath::Sin(p0); 146 x0(1) = D * TMath::Cos(p0); 147 x0(2) = z0; 148 // 149 return x0; 150 } 151 // 152 TVectorD VertexFit::Fill_x(TVectorD par, Double_t phi) 153 { 154 // 155 // Calculate track 3D position for a given phase, phi 156 // 157 TVectorD x(3); 343 // Fill derivative matrix dx/d alpha 344 // D 345 A(0, 0) = -TMath::Sin(p0); 346 A(1, 0) = TMath::Cos(p0); 347 A(2, 0) = 0.0; 348 // phi0 349 A(0, 1) = -D * TMath::Cos(p0) + (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C); 350 A(1, 1) = -D * TMath::Sin(p0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C); 351 A(2, 1) = 0.0; 352 // C 353 A(0, 2) = -(TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C * C); 354 A(1, 2) = (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C * C); 355 A(2, 2) = -ct * phi / (2 * C * C); 356 // z0 357 A(0, 3) = 0.0; 358 A(1, 3) = 0.0; 359 A(2, 3) = 1.0; 360 // ct = lambda 361 A(0, 4) = 0.0; 362 A(1, 4) = 0.0; 363 A(2, 4) = phi / (2 * C); 364 // 365 return A; 366 } 367 // 368 TVectorD VertexFit::Fill_a(TVectorD par, Double_t phi) 369 { 370 // 371 // Derivative of track 3D position vector with respect to phase at constant track parameters 372 // 373 // par = vector of track parameters 374 // phi = phase 375 // 376 TVectorD a(3); 158 377 // 159 378 // Decode input arrays … … 165 384 Double_t ct = par(4); 166 385 // 386 a(0) = TMath::Cos(phi + p0) / (2 * C); 387 a(1) = TMath::Sin(phi + p0) / (2 * C); 388 a(2) = ct / (2 * C); 389 // 390 return a; 391 } 392 // 393 TVectorD VertexFit::Fill_x0(TVectorD par) 394 { 395 // 396 // Calculate track 3D position at R = |D| (minimum approach to z-axis) 397 // 398 TVectorD x0(3); 399 // 400 // Decode input arrays 401 // 402 Double_t D = par(0); 403 Double_t p0 = par(1); 404 Double_t C = par(2); 405 Double_t z0 = par(3); 406 Double_t ct = par(4); 407 // 408 x0(0) = -D * TMath::Sin(p0); 409 x0(1) = D * TMath::Cos(p0); 410 x0(2) = z0; 411 // 412 return x0; 413 } 414 // 415 TVectorD VertexFit::Fill_x(TVectorD par, Double_t phi) 416 { 417 // 418 // Calculate track 3D position for a given phase, phi 419 // 420 TVectorD x(3); 421 // 422 // Decode input arrays 423 // 424 Double_t D = par(0); 425 Double_t p0 = par(1); 426 Double_t C = par(2); 427 Double_t z0 = par(3); 428 Double_t ct = par(4); 429 // 167 430 TVectorD x0 = Fill_x0(par); 168 431 x(0) = x0(0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C); … … 176 439 { 177 440 // 178 // Get track parameters, covariance and phase 179 Double_t fs = ffi[i]; // Get phase 441 // Get track parameters and their covariance 180 442 TVectorD par = *fParNew[i]; 181 443 TMatrixDSym Cov = *fCov[i]; 182 444 // 183 445 // Fill all track related work arrays arrays 184 TMatrixD A = derXdPar(par, fs); // A = dx/da = derivatives wrt track parameters 185 TMatrixDSym Winv = Cov; 186 Winv.Similarity(A); // W^-1 = A*C*A' 187 188 TMatrixD At(TMatrixD::kTransposed, A); // A transposed 189 fAti.push_back(new TMatrixD(At)); // Store A' 190 fWinvi.push_back(new TMatrixDSym(Winv)); // Store W^-1 matrix 191 // 446 Double_t fs = ffi[i]; // Get phase 192 447 TVectorD xs = Fill_x(par, fs); 193 448 fx0i.push_back(new TVectorD(xs)); // Start helix position 194 // 195 TVectorD di = A * (par - *fPar[i]); // x-shift 196 fdi.push_back(new TVectorD(di)); // Store x-shift 449 // 450 TMatrixD A = Fill_A(par, fs); // A = dx/da = derivatives wrt track parameters 451 TMatrixD At(TMatrixD::kTransposed, A); // A transposed 452 fAti.push_back(new TMatrixD(At)); // Store A' 453 TVectorD di = A * (par - *fPar[i]); // x-shift 454 fdi.push_back(new TVectorD(di)); // Store x-shift 455 TMatrixDSym Winv = Cov.Similarity(A); // W^-1 = A*C*A' 456 fWinvi.push_back(new TMatrixDSym(Winv)); // Store W^-1 matrix 197 457 // 198 458 TMatrixDSym W = RegInv(Winv); // W = (A*C*A')^-1 199 459 fWi.push_back(new TMatrixDSym(W)); // Store W matrix 200 460 // 201 TVectorD a = derXds(par, fs); // a = dx/ds = derivatives wrt phase461 TVectorD a = Fill_a(par, fs); // a = dx/ds = derivatives wrt phase 202 462 fai.push_back(new TVectorD(a)); // Store a 203 463 // … … 213 473 } 214 474 // 215 void VertexFit::VtxFitNoSteer()216 {217 //218 // Initialize219 //220 std::vector<TVectorD*> x0i; // Tracks at ma221 std::vector<TVectorD*> ni; // Track derivative wrt phase222 std::vector<TMatrixDSym*> Ci; // Position error matrix at fixed phase223 std::vector<TVectorD*> wi; // Ci*ni224 //225 // Track loop226 for (Int_t i = 0; i < fNtr; i++)227 {228 Double_t s = 0.;229 TVectorD par = *fPar[i];230 TMatrixDSym Cov = *fCov[i];231 x0i.push_back(new TVectorD(Fill_x0(par)));232 ni.push_back(new TVectorD(derXds(par, s)));233 TMatrixD A = derXdPar(par, s);234 Ci.push_back(new TMatrixDSym(Cov.Similarity(A)));235 TMatrixDSym Cinv = RegInv(*Ci[i]);236 wi.push_back(new TVectorD(Cinv * (*ni[i])));237 }238 //std::cout << "Vtx init completed. fNtr = "<<fNtr << std::endl;239 //240 // Get fit vertex241 //242 TMatrixDSym D(3); D.Zero();243 TVectorD Dx(3); Dx.Zero();244 for (Int_t i = 0; i < fNtr; i++)245 {246 TMatrixDSym Cinv = RegInv(*Ci[i]);247 TMatrixDSym W(3);248 W.Rank1Update(*wi[i], 1. / Ci[i]->Similarity(*wi[i]));249 TMatrixDSym Dd = Cinv - W;250 D += Dd;251 Dx += Dd * (*x0i[i]);252 }253 if(fVtxCst){254 D += fCovCstInv;255 Dx += fCovCstInv*fxCst;256 }257 fXv = RegInv(D) * Dx;258 //std::cout << "Fast vertex (x, y, z) = "<<fXv(0)<<", "<<fXv(1)<<", "<<fXv(2) << std::endl;259 //260 // Get fit phases261 //262 for (Int_t i = 0; i < fNtr; i++){263 Double_t si = Dot(*wi[i], fXv - (*x0i[i])) / Ci[i]->Similarity(*wi[i]);264 ffi.push_back(si);265 //TVectorD xvi = Fill_x(*fPar[i],si);266 //std::cout << "Fast vertex "<<i<<": xvi = "<<xvi(0)<<", "<<xvi(1)<<", "<<xvi(2)267 // <<", si = "<<si<< std::endl;268 }269 //270 // Cleanup271 for (Int_t i = 0; i < fNtr; i++)272 {273 x0i[i]->Clear(); delete x0i[i];274 ni[i]->Clear(); delete ni[i];275 Ci[i]->Clear(); delete Ci[i];276 wi[i]->Clear(); delete wi[i];277 }278 x0i.clear();;279 ni.clear();280 Ci.clear();281 wi.clear();282 }283 //284 475 void VertexFit::VertexFitter() 285 476 { … … 296 487 TMatrixDSym covX(3); 297 488 Double_t Chi2 = 0; 298 // 299 VtxFitNoSteer(); // Fast vertex finder on first pass (set ffi and fXv) 300 TVectorD x0 = fXv; 489 TVectorD x0 = fXv; // If previous fit done 490 if (fRold < 0.0) { 491 // External constraint 492 if (fVtxCst) fRold = TMath::Sqrt(fxCst(0) * fxCst(0) + fxCst(1) * fxCst(1)); 493 // No constraint 494 else for (Int_t i = 0; i < 3; i++)x0(i) = 1000.; // Set to arbitrary large value 495 } 496 // 497 // Starting vertex radius approximation 498 // 499 Double_t R = fRold; // Use previous fit if available 500 if (R < 0.0) R = StartRadius(); // Rough vertex estimate 501 //std::cout << "Start radius: " << R << std::endl; 301 502 // 302 503 // Iteration properties … … 314 515 TVectorD cterm(3); TMatrixDSym H(3); TMatrixDSym DW1D(3); 315 516 covX.Zero(); // Reset vertex covariance 316 cterm.Zero(); 517 cterm.Zero(); // Reset constant term 317 518 H.Zero(); // Reset H matrix 318 519 DW1D.Zero(); … … 329 530 TVectorD par = *fPar[i]; 330 531 TMatrixDSym Cov = *fCov[i]; 532 // 533 // For first iteration only 534 Double_t fs; 535 if (Ntry <= 0) // Initialize all phases on first pass 536 { 537 Double_t D = par(0); 538 Double_t C = par(2); 539 Double_t arg = TMath::Max(1.0e-6, (R * R - D * D) / (1 + 2 * C * D)); 540 fs = 2 * TMath::ASin(C * TMath::Sqrt(arg)); 541 ffi.push_back(fs); 542 } 331 543 // 332 544 // Update track related arrays … … 343 555 // update constant term 344 556 TVectorD xs = *fx0i[i] - *fdi[i]; 345 //TVectorD xx0 = *fx0i[i];346 347 //std::cout << "Iter. " << Ntry << ", trk " << i << ", xs= "348 // << xs(0) << ", " << xs(1) << ", " << xs(2)<<349 //", ph0= "<<par(1)<< std::endl;350 557 TVectorD xx0 = *fx0i[i]; 558 /* 559 std::cout << "Iter. " << Ntry << ", trk " << i << ", x= " 560 << xx0(0) << ", " << xx0(1) << ", " << xx0(2)<< 561 ", ph0= "<<par(1)<< std::endl; 562 */ 351 563 cterm += Ds * xs; 352 564 } // End loop on tracks … … 372 584 TMatrixDSym Wm1 = *fWinvi[i]; 373 585 fChi2List(i) = Wm1.Similarity(lambda); 374 if(fChi2List(i) < 0.0){375 //std::cout<<"# "<<i<<", Chi2= "<<fChi2List(i)<<", Wm1:"<<std::endl; Wm1.Print();376 //std::cout<<"Lambda= "<<std::endl; lambda.Print();377 }378 586 Chi2 += fChi2List(i); 379 587 TVectorD a = *fai[i]; 380 TVectorD b = (*fWi[i]) * (x - *fx0i[i] + *fdi[i]);381 f fi[i] += Dot(a, b) / fa2i[i];588 TVectorD b = (*fWi[i]) * (x - (*fx0i[i])); 589 for (Int_t j = 0; j < 3; j++)ffi[i] += a(j) * b(j) / fa2i[i]; 382 590 TVectorD newPar = *fPar[i] - ((*fCov[i]) * (*fAti[i])) * lambda; 383 591 fParNew[i] = new TVectorD(newPar); … … 402 610 // 403 611 fVtxDone = kTRUE; // Set fit completion flag 404 fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1)); // Store fit 405 //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2) 406 // << ", after "<<Ntry<<" iterations"<<std::endl; 612 fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1)); // Store fit radius 407 613 // 408 614 } -
external/TrackCovariance/VertexFit.h
r78ce8d1 r0c0c9af 6 6 #include <TVectorD.h> 7 7 #include <TMatrixDSym.h> 8 #include "TrkUtil.h"9 8 #include "ObsTrk.h" 10 9 #include <vector> … … 13 12 // Class for vertex fitting 14 13 15 class VertexFit: public TrkUtil 16 { 14 class VertexFit { 17 15 // 18 16 // Vertex fitting with track parameters steering … … 23 21 // 24 22 // Inputs 25 Int_t fNtr; 23 Int_t fNtr; // Number of tracks 26 24 std::vector<TVectorD*> fPar; // Input parameter array 27 25 std::vector<TVectorD*> fParNew; // Updated parameter array … … 29 27 std::vector<TMatrixDSym*> fCovNew; // Updated parameter covariances 30 28 // Constraints 31 Bool_t fVtxCst; 32 TVectorD fxCst; 29 Bool_t fVtxCst; // Vertex constraint flag 30 TVectorD fxCst; // Constraint value 33 31 TMatrixDSym fCovCst; // Constraint 34 32 TMatrixDSym fCovCstInv; // Inverse of constraint covariance 35 33 // 36 34 // Results 37 Bool_t fVtxDone; // Flag vertex fit completed35 Bool_t fVtxDone; // Flag vertex fit completed 38 36 Double_t fRold; // Current value of vertex radius 39 37 TVectorD fXv; // Found vertex … … 54 52 // 55 53 // Service routines 54 //void InitWrkArrays(); // Initializations 56 55 void ResetWrkArrays(); // Clear work arrays 57 //Double_t StartRadius(); // Starting vertex radius determination58 //Double_t FastRv(TVectorD p1, TVectorD p2); // Fast vertex radius determination59 //TMatrixDSym RegInv(TMatrixDSym& Smat0);// Regularized 3D matrix inversion60 //TMatrixD Fill_A(TVectorD par, Double_t phi);// Derivative of track position wrt track parameters61 //TVectorD Fill_a(TVectorD par, Double_t phi);// Derivative of track position wrt track phase56 Double_t StartRadius(); // Starting vertex radius determination 57 Double_t FastRv(TVectorD p1, TVectorD p2); // Fast vertex radius determination 58 TMatrixDSym RegInv(TMatrixDSym& Smat0); // Regularized 3D matrix inversion 59 TMatrixD Fill_A(TVectorD par, Double_t phi); // Derivative of track position wrt track parameters 60 TVectorD Fill_a(TVectorD par, Double_t phi); // Derivative of track position wrt track phase 62 61 TVectorD Fill_x0(TVectorD par); // Track position at dma to z-axis 63 TVectorD Fill_x(TVectorD par, Double_t phi); // Track position at given phase62 TVectorD Fill_x(TVectorD par, Double_t phi); // Track position at given phase 64 63 void UpdateTrkArrays(Int_t i); // Fill track realted arrays 65 void VtxFitNoSteer(); // Vertex fitter routine w/o parameter steering 66 void VertexFitter(); // Vertex fitter routine w/ parameter steering 64 void VertexFitter(); // Vertex finder routine 67 65 public: 68 66 // -
modules/TreeWriter.cc
r78ce8d1 r0c0c9af 369 369 370 370 // add some offdiagonal covariance matrix elements 371 entry->ErrorD0Phi = candidate->TrackCovariance(0,1) *1.e3;371 entry->ErrorD0Phi = candidate->TrackCovariance(0,1); 372 372 entry->ErrorD0C = candidate->TrackCovariance(0,2); 373 entry->ErrorD0DZ = candidate->TrackCovariance(0,3) *1.e6;374 entry->ErrorD0CtgTheta = candidate->TrackCovariance(0,4) *1.e3;375 entry->ErrorPhiC = candidate->TrackCovariance(1,2) *1.e-3;376 entry->ErrorPhiDZ = candidate->TrackCovariance(1,3) *1.e3;373 entry->ErrorD0DZ = candidate->TrackCovariance(0,3); 374 entry->ErrorD0CtgTheta = candidate->TrackCovariance(0,4); 375 entry->ErrorPhiC = candidate->TrackCovariance(1,2); 376 entry->ErrorPhiDZ = candidate->TrackCovariance(1,3); 377 377 entry->ErrorPhiCtgTheta = candidate->TrackCovariance(1,4); 378 378 entry->ErrorCDZ = candidate->TrackCovariance(2,3); 379 entry->ErrorCCtgTheta = candidate->TrackCovariance(2,4) *1.e-3;380 entry->ErrorDZCtgTheta = candidate->TrackCovariance(3,4) *1.e3;379 entry->ErrorCCtgTheta = candidate->TrackCovariance(2,4); 380 entry->ErrorDZCtgTheta = candidate->TrackCovariance(3,4); 381 381 382 382 entry->Xd = candidate->Xd;
Note:
See TracChangeset
for help on using the changeset viewer.
