- Timestamp:
- Jan 19, 2022, 12:08:01 PM (3 years ago)
- Branches:
- master
- Children:
- 2c8865f, 56fb0be, 5c03893
- Parents:
- 0c0c9af (diff), 78ce8d1 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the(diff)
links above to see all the changes relative to each parent. - git-author:
- Michele Selvaggi <michele.selvaggi@…> (01/19/22 12:08:01)
- git-committer:
- GitHub <noreply@…> (01/19/22 12:08:01)
- Location:
- examples
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
examples/ExamplePVtxFind.C
r0c0c9af rdd263e4 1 1 /* 2 2 Example of using vertex fitter class to fit primary vertex 3 assumed to be generated in (0,0,0)3 assumed to be generated with Pythia8/ee_zh_smr-shf.cmn 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); 46 47 // 47 48 // Loop over all events … … 87 88 pr[it] = new TVectorD(obsPar); 88 89 cv[it] = new TMatrixDSym(trk->CovarianceMatrix()); 90 // Check covariance 91 /* 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->ErrorD0 102 <<" - "<<TMath::Sqrt(Ctmp(0,0))<<std::endl; 103 } 104 */ 89 105 //std::cout << "Track loaded it= " << it << std::endl; 90 106 // … … 120 136 } // End loop on tracks 121 137 } 138 std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl; 122 139 // 123 140 // Find primary vertex … … 132 149 covpvc(1, 1) = 2.55e-05 * 2.55e-05; 133 150 covpvc(2, 2) = 0.64 * 0.64; 151 if(Nprim == 0){ 152 xpv = xpvc(0); 153 ypv = xpvc(1); 154 zpv = xpvc(2); 155 } 134 156 // 135 157 // … … 137 159 Int_t nSkim = 0; 138 160 Int_t* nSkimmed = new Int_t[NtrG]; 139 VertexFit* Vskim = new VertexFit();140 Vskim->AddVtxConstraint(xpvc, covpvc);141 Double_t MaxChi2 = 2.0+3* sqrt(2.);161 TVectorD** PrSk = new TVectorD * [1]; 162 TMatrixDSym** CvSk = new TMatrixDSym * [1]; 163 Double_t MaxChi2 = 2.0+3*2.; 142 164 for (Int_t n = 0; n < NtrG; n++) { 143 Vskim->AddTrk(pr[n], cv[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); 144 169 Double_t Chi2One = Vskim->GetVtxChi2(); 145 170 //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl; … … 149 174 nSkim++; 150 175 } 151 Vskim->RemoveTrk(0); 152 } 176 // Cleanup 177 delete Vskim; 178 } 179 delete PrSk[0]; 180 delete CvSk[0]; 181 delete[] PrSk; 182 delete[] CvSk; 183 // 153 184 // Load tracks for primary fit 154 185 Int_t MinTrk = 1; // Minumum # tracks for vertex fit 155 TVectorD** PrFit = new TVectorD * [nSkim];156 TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];157 186 if (nSkim >= MinTrk) { 187 TVectorD** PrFit = new TVectorD * [nSkim]; 188 TMatrixDSym** CvFit = new TMatrixDSym * [nSkim]; 158 189 for (Int_t n = 0; n < nSkim; n++) { 159 190 PrFit[n] = new TVectorD(*pr[nSkimmed[n]]); 160 191 CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]); 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) { 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; 170 197 VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit); 171 198 //std::cout << "Vertex fit created " << std::endl; … … 173 200 // 174 201 // Remove tracks with large chi2 175 Double_t MaxChi2Fit = 9.;202 Double_t MaxChi2Fit = 4.5; 176 203 Bool_t Done = kFALSE; 177 204 while (!Done) { … … 180 207 TVectorD Chi2List = Vtx->GetVtxChi2List(); // Get contributions to Chi2 181 208 //std::cout << "After Chi2List. " << std::endl; Chi2List.Print(); 182 Double_t* Chi2L = new Double_t[Nfound];209 //Double_t* Chi2L = new Double_t[Nfound]; 183 210 Chi2L = Chi2List.GetMatrixArray(); 184 211 Int_t iMax = TMath::LocMax(Nfound, Chi2L); … … 186 213 Double_t Chi2Mx = Chi2L[iMax]; 187 214 //std::cout << "Chi2Mx "<<Chi2Mx << std::endl; 188 if (Chi2Mx > MaxChi2Fit && Nfound > 2) {215 if (Chi2Mx > MaxChi2Fit && Nfound > 1) { 189 216 //std::cout << "Before remove. Nfound = "<<Nfound << std::endl; 190 217 Vtx->RemoveTrk(iMax); … … 195 222 Done = kTRUE; 196 223 } 197 //std::cout << "Done = " << Done << std::endl;198 //delete[] Chi2L;199 //std::cout << "Array Chi2L removed " << std::endl;200 224 } 201 225 // … … 203 227 // 204 228 // Require minimum number of tracks in vertex 205 Int_t Nmin = 2;229 Int_t Nmin = 1; 206 230 if (Nfound >= Nmin) { 207 231 TVectorD xvtx = Vtx->GetVtx(); 208 //std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;232 std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl; 209 233 TMatrixDSym covX = Vtx->GetVtxCov(); 210 234 Double_t Chi2 = Vtx->GetVtxChi2(); … … 222 246 hTrPrim->Fill(Nprim); 223 247 hTrFound->Fill((Double_t)Nfound); 248 hTrDiff->Fill((Double_t)Nfound-Nprim); 224 249 //std::cout << "Histograms filled " << std::endl; 225 250 } 226 251 // 252 // Clean 227 253 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; 228 258 } 229 259 … … 233 263 for (Int_t i = 0; i < NtrG; i++) delete pr[i]; 234 264 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];237 265 delete[] pr; 238 266 delete[] cv; 239 delete[] PrFit;240 delete[] CvFit;241 267 } 242 268 // … … 254 280 // 255 281 TCanvas* CnvN = new TCanvas("CnvN", "Primary tracks found", 100, 100, 900, 500); 256 CnvN->Divide(2, 1);282 CnvN->Divide(2, 2); 257 283 CnvN->cd(1); 258 284 hTrPrim->Draw(); … … 260 286 hTrFound->SetLineColor(kRed); 261 287 hTrFound->Draw(); 288 CnvN->cd(3); 289 hTrDiff->Draw(); 262 290 } -
examples/ExamplePVtxFitEC.C
r0c0c9af rdd263e4 1 1 /* 2 2 Example of using vertex fitter class to fit primary vertex 3 assumed to be generated in (0,0,0)3 assumed to be generated with Pythia8/ee_zh_smr-shf.cmn 4 4 */ 5 5 … … 21 21 void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5) 22 22 { 23 // 24 // Beam constraint 25 // (consistent with generation with ee_zh_smr-shf.cmnd) 26 // 27 // Mean beam position 28 TVectorD xpvc(3); 29 xpvc(0) = 1.0; 30 xpvc(1) = -2.0; 31 xpvc(2) = 10.0; 32 // Interaction region covariance 33 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 // 23 39 // Create chain of root trees 24 40 TChain chain("Delphes"); … … 35 51 // Book histograms 36 52 Int_t Nbin = 100; 53 // Vertex fit pulls 37 54 TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.); 38 55 TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.); 39 56 TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.); 40 57 TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.); 58 // Generation check 59 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 41 63 // 42 64 // Loop over all events … … 65 87 // Get associated generated particle 66 88 GenParticle* gp = (GenParticle*)trk->Particle.GetObject(); 67 TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z); 89 TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z); // mm -> meters 68 90 TVector3 pg(gp->Px,gp->Py,gp->Pz); 69 91 Double_t Q = (Double_t)gp->Charge; 70 92 Double_t Bz = 2.0; 71 93 TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz); 72 TVectorD genPar = TrkUtil::ParToMm(genParM); 94 TVectorD genPar = TrkUtil::ParToMm(genParM); // -> back to mm 95 96 // 97 // Check if original track is from primary vertex 73 98 // 74 99 // Position of origin in mm … … 79 104 Int_t mp = gp->M1; // Mother 80 105 while(mp>0){ 81 GenParticle* gm = 82 (GenParticle*)branchGenPart->At(mp); 106 GenParticle* gm = (GenParticle*)branchGenPart->At(mp); 83 107 Double_t xm = gm->X; 84 108 Double_t ym = gm->Y; 85 109 Double_t zm = gm->Z; 86 110 if(x!=xm || y!=ym || z!=zm){ 87 prim = kFALSE; 111 prim = kFALSE; // Not primary 88 112 break; 89 113 }else mp = gm->M1; 90 114 } 115 91 116 // 92 117 // group tracks originating from the primary vertex 93 118 if (prim) 94 119 { 95 // std::cout<<"Event: "<<entry<<", track: "<<it;96 // std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl;120 // 121 // Store true primary vertex for this event 97 122 xpv = x; 98 123 ypv = y; … … 104 129 Double_t obsC = trk->C; 105 130 Double_t obsZ0 = trk->DZ; 106 //std::cout<<"Z0 track = "<< obsZ0107 //<<", gen Z0 = "<<genPar(3)108 //<<", gen cot = "<<genPar(4)<<std::endl;109 131 Double_t obsCtg = trk->CtgTheta; 110 132 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg }; … … 112 134 // 113 135 pr[Ntr] = new TVectorD(obsPar); 114 //std::cout<<"Cov Matrix:"<<std::endl;115 //trk->CovarianceMatrix().Print();116 136 cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix()); 117 137 Ntr++; … … 122 142 // Fit primary vertex with beam constraint 123 143 // 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 144 Int_t MinTrk = 3; // Minumum # tracks for vertex fit 134 145 if (Ntr >= MinTrk) { 135 146 VertexFit* Vtx = new VertexFit(Ntr, pr, cv); … … 139 150 Double_t Chi2 = Vtx->GetVtxChi2(); 140 151 Double_t Ndof = 2 * (Double_t)Ntr; 152 delete Vtx; 153 // 141 154 Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0)); 142 155 Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1)); 143 156 Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2)); 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; 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 // 148 162 // 149 163 // Fill histograms … … 152 166 hZpull->Fill(PullZ); 153 167 hChi2->Fill(Chi2 / Ndof); 168 // 169 hXvpull->Fill(PullXv); 170 hYvpull->Fill(PullYv); 171 hZvpull->Fill(PullZv); 154 172 } 155 173 156 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;157 174 // 158 175 // Cleanup … … 162 179 delete[] cv; 163 180 } 181 164 182 // 165 183 // Show resulting histograms … … 174 192 hZpull->Fit("gaus"); hZpull->Draw(); 175 193 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(); 176 203 }
Note:
See TracChangeset
for help on using the changeset viewer.