Fork me on GitHub

source: git/examples/ExamplePVtxFind.C@ 3a105e5

Last change on this file since 3a105e5 was 7bca620, checked in by Franco BEDESCHI <bed@…>, 3 years ago

Add feature to force starting radius in vertex fit
  • Property mode set to 100644
File size: 10.2 KB
RevLine 
[82db145]1/*
2Example of using vertex fitter class to fit primary vertex
[b750b0a]3assumed to be generated with Pythia8/ee_zh_smr-shf.cmn
[82db145]4*/
5
6#ifdef __CLING__
7R__LOAD_LIBRARY(libDelphes)
8
9#include "classes/DelphesClasses.h"
10#include "external/ExRootAnalysis/ExRootTreeReader.h"
11#include "modules/TrackCovariance.h"
12#include "external/TrackCovariance/TrkUtil.h"
13#include "external/TrackCovariance/VertexFit.h"
14
15#endif
16
17
18//------------------------------------------------------------------------------
19
20void ExamplePVtxFind(const char* inputFile, Int_t Nevent = 5)
21{
22 // Create chain of root trees
23 TChain chain("Delphes");
24 chain.Add(inputFile);
25
26 // Create object of class ExRootTreeReader
27 ExRootTreeReader* treeReader = new ExRootTreeReader(&chain);
28 Long64_t numberOfEntries = treeReader->GetEntries();
29
30 // Get pointers to branches used in this analysis
31 TClonesArray* branchGenPart = treeReader->UseBranch("Particle");
32 TClonesArray* branchTrack = treeReader->UseBranch("Track");
33
34 // Book histograms
35 // Vertex fit pulls
36 Int_t Nbin = 100;
37 TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.);
38 TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.);
39 TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.);
40 TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.);
41 // Track vertex associations
42 TH1D* hTrPrim = new TH1D("hTrPrim", "Available primary tracks", 41, -0.5, 40.5);
43 TH1D* hTrFound = new TH1D("hTrFound", "Found primary tracks", 41, -0.5, 40.5);
44 TH1D* hTrFoundG = new TH1D("hTrFoundG", "Found GOOD primary tracks", 41, -0.5, 40.5);
45 TH1D* hTrFoundB = new TH1D("hTrFoundB", "Found BAD primary tracks", 41, -0.5, 40.5);
[b750b0a]46 TH1D* hTrDiff = new TH1D("hTrDiff", "Found - available tracks", 21, -10.5, 10.5);
[7bca620]47 // Chi2 for cuts
48 TH1D* hChi2SngP = new TH1D("hChi2SngP", "#chi^2 for single tracks", 200, 0., 50.);
49 TH1D* hChi2MaxP = new TH1D("hChi2MaxP", "#chi^2 max contribution", 200, 0., 50.);
50 TH1D* hChi2SngN = new TH1D("hChi2SngN", "#chi^2 for single tracks", 200, 0., 50.);
51 TH1D* hChi2MaxN = new TH1D("hChi2MaxN", "#chi^2 max contribution", 200, 0., 50.);
52
53
[82db145]54 //
55 // Loop over all events
56 Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries);
57 for (Int_t entry = 0; entry < Nev; ++entry)
58 {
59 // Load selected branches with data from specified event
60 treeReader->ReadEntry(entry);
61 Int_t NtrG = branchTrack->GetEntries();
[7bca620]62 TVectorD** pr = new TVectorD * [NtrG]; // Track Parameters
63 TMatrixDSym** cv = new TMatrixDSym * [NtrG]; // Track covariances
64 Bool_t *fprvx = new Bool_t[NtrG]; // Primary vertex flag
[82db145]65 //
66 // test Particle branch
67 Int_t Ngen = branchGenPart->GetEntries();
68 //std::cout << "Nr. of generated particles: " << Ngen << std::endl;
69 // If event contains at least 1 track
70 //
71 Double_t Nprim = 0.0;
[46d3442]72 Double_t xpv = 0.0; // Init true primary vertex position
73 Double_t ypv = 0.0;
74 Double_t zpv = 0.0;
[82db145]75 if (branchTrack->GetEntries() > 0)
76 {
77 // Loop on tracks
78 for (Int_t it = 0; it < branchTrack->GetEntries(); it++)
79 {
80 Track* trk = (Track*)branchTrack->At(it);
81 //
82 // Start fitting all available tracks
83 //
84 // Reconstructed track parameters
85 Double_t obsD0 = trk->D0;
86 Double_t obsPhi = trk->Phi;
87 Double_t obsC = trk->C;
88 Double_t obsZ0 = trk->DZ;
89 Double_t obsCtg = trk->CtgTheta;
90 //std::cout << "Got track parameters for track " << it << std::endl;
91 //
[46d3442]92 // Load all tracks for vertex fit
93 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
94 TVectorD obsPar(5, oPar); // Fill observed parameters
95 pr[it] = new TVectorD(obsPar);
96 cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
97 //
98 // Find true primary vertex
[82db145]99 GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
100 //std::cout << "GenParticle pointer "<<gp << std::endl;
101 //
102 // Position of origin in mm
103 Double_t x = gp->X;
104 Double_t y = gp->Y;
105 Double_t z = gp->Z;
[46d3442]106 Bool_t prim = kTRUE; // Is primary?
[7bca620]107 fprvx[it] = kFALSE;
[46d3442]108 Int_t mp = gp->M1; // Mother
109 while (mp > 0) {
110 GenParticle* gm =
111 (GenParticle*)branchGenPart->At(mp);
112 Double_t xm = gm->X;
113 Double_t ym = gm->Y;
114 Double_t zm = gm->Z;
115 if (x != xm || y != ym || z != zm) {
116 prim = kFALSE;
117 break;
118 }
119 else mp = gm->M1;
120 }
121 if (prim) { // It's a primary track
122 Nprim++;
[7bca620]123 fprvx[it] = kTRUE;
[46d3442]124 xpv = x; // Store true primary
125 ypv = y;
126 zpv = z;
127 }
[82db145]128
129 } // End loop on tracks
130 }
[7bca620]131 if(entry%500 ==0){
132 std::cout << "Event "<<entry<<" opened containing " << NtrG << " / "<< Nprim
133 << " Total / primary tracks"<< std::endl;
134 }
135 //std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl;
[82db145]136 //
[46d3442]137 // Find primary vertex
[82db145]138 //
[46d3442]139 //Beam constraint
140 TVectorD xpvc(3);
141 xpvc(0) = 1.0;
142 xpvc(1) = -2.0;
143 xpvc(2) = 10.0;
144 TMatrixDSym covpvc(3); covpvc.Zero();
145 covpvc(0, 0) = 0.0097 * 0.0097;
146 covpvc(1, 1) = 2.55e-05 * 2.55e-05;
147 covpvc(2, 2) = 0.64 * 0.64;
[b750b0a]148 if(Nprim == 0){
149 xpv = xpvc(0);
150 ypv = xpvc(1);
151 zpv = xpvc(2);
152 }
[46d3442]153 //
154 //
155 // Skim tracks
156 Int_t nSkim = 0;
157 Int_t* nSkimmed = new Int_t[NtrG];
[b750b0a]158 TVectorD** PrSk = new TVectorD * [1];
159 TMatrixDSym** CvSk = new TMatrixDSym * [1];
[7bca620]160 Double_t MaxChi2 = 9.;
[46d3442]161 for (Int_t n = 0; n < NtrG; n++) {
[b750b0a]162 PrSk[0] = new TVectorD(*pr[n]);
163 CvSk[0] = new TMatrixDSym(*cv[n]);
164 VertexFit* Vskim = new VertexFit(1,PrSk, CvSk);
165 Vskim->AddVtxConstraint(xpvc, covpvc);
[46d3442]166 Double_t Chi2One = Vskim->GetVtxChi2();
167 //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
[7bca620]168 if(fprvx[n])hChi2SngP->Fill(Chi2One);
169 else hChi2SngN->Fill(Chi2One);
170 //
[46d3442]171 if (Chi2One < MaxChi2) {
172 nSkimmed[nSkim] = n;
173 //std::cout << "nSkimmed[" << nSkim << "] = " << n << std::endl;
174 nSkim++;
175 }
[b750b0a]176 // Cleanup
177 delete Vskim;
[46d3442]178 }
[b750b0a]179 delete PrSk[0];
180 delete CvSk[0];
181 delete[] PrSk;
182 delete[] CvSk;
183 //
[46d3442]184 // Load tracks for primary fit
185 Int_t MinTrk = 1; // Minumum # tracks for vertex fit
[7bca620]186 std::vector<Int_t> trnum;
[46d3442]187 if (nSkim >= MinTrk) {
[b750b0a]188 TVectorD** PrFit = new TVectorD * [nSkim];
189 TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
[46d3442]190 for (Int_t n = 0; n < nSkim; n++) {
191 PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
192 CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]);
[7bca620]193 trnum.push_back(nSkimmed[n]);
[46d3442]194 }
[b750b0a]195 delete[] nSkimmed;
196 Int_t Nfound = nSkim;
197 if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
198 const Int_t MaxFound = 100; Double_t Chi2LL[MaxFound]; Double_t *Chi2L = Chi2LL;
[46d3442]199 VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit);
[82db145]200 //std::cout << "Vertex fit created " << std::endl;
[46d3442]201 Vtx->AddVtxConstraint(xpvc, covpvc);
[82db145]202 //
203 // Remove tracks with large chi2
[7bca620]204 Double_t MaxChi2Fit = 8.0;
[82db145]205 Bool_t Done = kFALSE;
[46d3442]206 while (!Done) {
[82db145]207 //std::cout << "After while " << std::endl;
208 // Find largest Chi2 contribution
209 TVectorD Chi2List = Vtx->GetVtxChi2List(); // Get contributions to Chi2
210 //std::cout << "After Chi2List. " << std::endl; Chi2List.Print();
[b750b0a]211 //Double_t* Chi2L = new Double_t[Nfound];
[82db145]212 Chi2L = Chi2List.GetMatrixArray();
213 Int_t iMax = TMath::LocMax(Nfound, Chi2L);
214 //std::cout << "iMax = "<<iMax << std::endl;
215 Double_t Chi2Mx = Chi2L[iMax];
216 //std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
[7bca620]217 if(fprvx[trnum[iMax]])hChi2MaxP->Fill(Chi2Mx);
218 else hChi2MaxN->Fill(Chi2Mx);
[b750b0a]219 if (Chi2Mx > MaxChi2Fit && Nfound > 1) {
[46d3442]220 //std::cout << "Before remove. Nfound = "<<Nfound << std::endl;
221 Vtx->RemoveTrk(iMax);
[7bca620]222 trnum.erase(trnum.begin() + iMax);
[46d3442]223 //std::cout << "After remove." << std::endl;
[82db145]224 Nfound--;
225 }
226 else {
227 Done = kTRUE;
228 }
229 }
230 //
231 //std::cout << "Before getting vertex " << std::endl;
232 //
233 // Require minimum number of tracks in vertex
[b750b0a]234 Int_t Nmin = 1;
[82db145]235 if (Nfound >= Nmin) {
236 TVectorD xvtx = Vtx->GetVtx();
[7bca620]237 //std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
[82db145]238 TMatrixDSym covX = Vtx->GetVtxCov();
239 Double_t Chi2 = Vtx->GetVtxChi2();
[46d3442]240 Double_t Ndof = 2 * (Double_t)Nfound;
241 Double_t PullX = (xvtx(0) - xpv) / TMath::Sqrt(covX(0, 0));
242 Double_t PullY = (xvtx(1) - ypv) / TMath::Sqrt(covX(1, 1));
243 Double_t PullZ = (xvtx(2) - zpv) / TMath::Sqrt(covX(2, 2));
[82db145]244 //
245 // Fill histograms
246 hXpull->Fill(PullX);
247 hYpull->Fill(PullY);
248 hZpull->Fill(PullZ);
249 hChi2->Fill(Chi2 / Ndof);
250 //
251 hTrPrim->Fill(Nprim);
252 hTrFound->Fill((Double_t)Nfound);
[b750b0a]253 hTrDiff->Fill((Double_t)Nfound-Nprim);
[82db145]254 //std::cout << "Histograms filled " << std::endl;
255 }
256 //
[b750b0a]257 // Clean
[82db145]258 delete Vtx;
[b750b0a]259 for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
260 for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
261 delete[] PrFit;
262 delete[] CvFit;
[7bca620]263 delete[] fprvx;
264 trnum.clear();
[82db145]265 }
266
267 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
268 //
269 // Cleanup
[46d3442]270 for (Int_t i = 0; i < NtrG; i++) delete pr[i];
271 for (Int_t i = 0; i < NtrG; i++) delete cv[i];
[82db145]272 delete[] pr;
273 delete[] cv;
274 }
275 //
276 // Show resulting histograms
277 //
278 TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
279 Cnv->Divide(2, 2);
[7bca620]280 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);
[82db145]281 hXpull->Fit("gaus"); hXpull->Draw();
[7bca620]282 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);
[82db145]283 hYpull->Fit("gaus"); hYpull->Draw();
[7bca620]284 Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);
[82db145]285 hZpull->Fit("gaus"); hZpull->Draw();
286 Cnv->cd(4); hChi2->Draw();
287 //
288 TCanvas* CnvN = new TCanvas("CnvN", "Primary tracks found", 100, 100, 900, 500);
[b750b0a]289 CnvN->Divide(2, 2);
[82db145]290 CnvN->cd(1);
291 hTrPrim->Draw();
[46d3442]292 CnvN->cd(2);
[82db145]293 hTrFound->SetLineColor(kRed);
294 hTrFound->Draw();
[b750b0a]295 CnvN->cd(3);
296 hTrDiff->Draw();
[7bca620]297 //
298 TCanvas* CnvCh = new TCanvas("CnvCh", "#chi^2", 200, 200, 900, 500);
299 CnvCh->Divide(2,1);
300 CnvCh->cd(1);
301 hChi2SngP->Draw();
302 hChi2SngN->SetLineColor(kRed);
303 hChi2SngN->Draw("SAME");
304 CnvCh->cd(2);
305 hChi2MaxP->Draw();
306 hChi2MaxN->SetLineColor(kRed);
307 hChi2MaxN->Draw("SAME");
[82db145]308}
Note: See TracBrowser for help on using the repository browser.