Fork me on GitHub

source: git/examples/ExamplePVtxFind.C@ c1780a5

Last change on this file since c1780a5 was 46d3442, checked in by Franco BEDESCHI <bed@…>, 4 years ago

Vertex contrained fit and primary vertex finder added
  • Property mode set to 100644
File size: 8.6 KB
RevLine 
[82db145]1/*
2Example of using vertex fitter class to fit primary vertex
3assumed to be generated in (0,0,0)
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);
46 //
47 // Loop over all events
48 Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries);
49 for (Int_t entry = 0; entry < Nev; ++entry)
50 {
51 // Load selected branches with data from specified event
52 treeReader->ReadEntry(entry);
53 Int_t NtrG = branchTrack->GetEntries();
[46d3442]54 if(entry%500 ==0)std::cout << "Event "<<entry<<" opened containing " << NtrG << " tracks" << std::endl;
[82db145]55 TVectorD** pr = new TVectorD * [NtrG];
56 TMatrixDSym** cv = new TMatrixDSym * [NtrG];
57 //
58 // test Particle branch
59 Int_t Ngen = branchGenPart->GetEntries();
60 //std::cout << "Nr. of generated particles: " << Ngen << std::endl;
61 // If event contains at least 1 track
62 //
63 Double_t Nprim = 0.0;
[46d3442]64 Double_t xpv = 0.0; // Init true primary vertex position
65 Double_t ypv = 0.0;
66 Double_t zpv = 0.0;
[82db145]67 if (branchTrack->GetEntries() > 0)
68 {
69 // Loop on tracks
70 for (Int_t it = 0; it < branchTrack->GetEntries(); it++)
71 {
72 Track* trk = (Track*)branchTrack->At(it);
73 //
74 // Start fitting all available tracks
75 //
76 // Reconstructed track parameters
77 Double_t obsD0 = trk->D0;
78 Double_t obsPhi = trk->Phi;
79 Double_t obsC = trk->C;
80 Double_t obsZ0 = trk->DZ;
81 Double_t obsCtg = trk->CtgTheta;
82 //std::cout << "Got track parameters for track " << it << std::endl;
83 //
[46d3442]84 // Load all tracks for vertex fit
85 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
86 TVectorD obsPar(5, oPar); // Fill observed parameters
87 pr[it] = new TVectorD(obsPar);
88 cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
89 //std::cout << "Track loaded it= " << it << std::endl;
90 //
91 // Find true primary vertex
[82db145]92 GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
93 //std::cout << "GenParticle pointer "<<gp << std::endl;
94 //
95 // Position of origin in mm
96 Double_t x = gp->X;
97 Double_t y = gp->Y;
98 Double_t z = gp->Z;
[46d3442]99 Bool_t prim = kTRUE; // Is primary?
100 Int_t mp = gp->M1; // Mother
101 while (mp > 0) {
102 GenParticle* gm =
103 (GenParticle*)branchGenPart->At(mp);
104 Double_t xm = gm->X;
105 Double_t ym = gm->Y;
106 Double_t zm = gm->Z;
107 if (x != xm || y != ym || z != zm) {
108 prim = kFALSE;
109 break;
110 }
111 else mp = gm->M1;
112 }
113 if (prim) { // It's a primary track
114 Nprim++;
115 xpv = x; // Store true primary
116 ypv = y;
117 zpv = z;
118 }
[82db145]119
120 } // End loop on tracks
121 }
122 //
[46d3442]123 // Find primary vertex
[82db145]124 //
[46d3442]125 //Beam constraint
126 TVectorD xpvc(3);
127 xpvc(0) = 1.0;
128 xpvc(1) = -2.0;
129 xpvc(2) = 10.0;
130 TMatrixDSym covpvc(3); covpvc.Zero();
131 covpvc(0, 0) = 0.0097 * 0.0097;
132 covpvc(1, 1) = 2.55e-05 * 2.55e-05;
133 covpvc(2, 2) = 0.64 * 0.64;
134 //
135 //
136 // Skim tracks
137 Int_t nSkim = 0;
138 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.);
142 for (Int_t n = 0; n < NtrG; n++) {
143 Vskim->AddTrk(pr[n], cv[n]);
144 Double_t Chi2One = Vskim->GetVtxChi2();
145 //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
146 if (Chi2One < MaxChi2) {
147 nSkimmed[nSkim] = n;
148 //std::cout << "nSkimmed[" << nSkim << "] = " << n << std::endl;
149 nSkim++;
150 }
151 Vskim->RemoveTrk(0);
152 }
153 // Load tracks for primary fit
154 Int_t MinTrk = 1; // Minumum # tracks for vertex fit
155 TVectorD** PrFit = new TVectorD * [nSkim];
156 TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
157 if (nSkim >= MinTrk) {
158 for (Int_t n = 0; n < nSkim; n++) {
159 PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
160 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) {
170 VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit);
[82db145]171 //std::cout << "Vertex fit created " << std::endl;
[46d3442]172 Vtx->AddVtxConstraint(xpvc, covpvc);
[82db145]173 //
174 // Remove tracks with large chi2
[46d3442]175 Double_t MaxChi2Fit = 9.;
[82db145]176 Bool_t Done = kFALSE;
[46d3442]177 while (!Done) {
[82db145]178 //std::cout << "After while " << std::endl;
179 // Find largest Chi2 contribution
180 TVectorD Chi2List = Vtx->GetVtxChi2List(); // Get contributions to Chi2
181 //std::cout << "After Chi2List. " << std::endl; Chi2List.Print();
182 Double_t* Chi2L = new Double_t[Nfound];
183 Chi2L = Chi2List.GetMatrixArray();
184 Int_t iMax = TMath::LocMax(Nfound, Chi2L);
185 //std::cout << "iMax = "<<iMax << std::endl;
186 Double_t Chi2Mx = Chi2L[iMax];
187 //std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
[46d3442]188 if (Chi2Mx > MaxChi2Fit && Nfound > 2) {
189 //std::cout << "Before remove. Nfound = "<<Nfound << std::endl;
190 Vtx->RemoveTrk(iMax);
191 //std::cout << "After remove." << std::endl;
[82db145]192 Nfound--;
193 }
194 else {
195 Done = kTRUE;
196 }
197 //std::cout << "Done = " << Done << std::endl;
198 //delete[] Chi2L;
199 //std::cout << "Array Chi2L removed " << std::endl;
200 }
201 //
202 //std::cout << "Before getting vertex " << std::endl;
203 //
204 // Require minimum number of tracks in vertex
[46d3442]205 Int_t Nmin = 2;
[82db145]206 if (Nfound >= Nmin) {
207 TVectorD xvtx = Vtx->GetVtx();
208 //std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
209 TMatrixDSym covX = Vtx->GetVtxCov();
210 Double_t Chi2 = Vtx->GetVtxChi2();
[46d3442]211 Double_t Ndof = 2 * (Double_t)Nfound;
212 Double_t PullX = (xvtx(0) - xpv) / TMath::Sqrt(covX(0, 0));
213 Double_t PullY = (xvtx(1) - ypv) / TMath::Sqrt(covX(1, 1));
214 Double_t PullZ = (xvtx(2) - zpv) / TMath::Sqrt(covX(2, 2));
[82db145]215 //
216 // Fill histograms
217 hXpull->Fill(PullX);
218 hYpull->Fill(PullY);
219 hZpull->Fill(PullZ);
220 hChi2->Fill(Chi2 / Ndof);
221 //
222 hTrPrim->Fill(Nprim);
223 hTrFound->Fill((Double_t)Nfound);
224 //std::cout << "Histograms filled " << std::endl;
225 }
226 //
227 delete Vtx;
228 }
229
230 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
231 //
232 // Cleanup
[46d3442]233 for (Int_t i = 0; i < NtrG; i++) delete pr[i];
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];
[82db145]237 delete[] pr;
238 delete[] cv;
[46d3442]239 delete[] PrFit;
240 delete[] CvFit;
[82db145]241 }
242 //
243 // Show resulting histograms
244 //
245 TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
246 Cnv->Divide(2, 2);
247 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
248 hXpull->Fit("gaus"); hXpull->Draw();
249 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
250 hYpull->Fit("gaus"); hYpull->Draw();
251 Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
252 hZpull->Fit("gaus"); hZpull->Draw();
253 Cnv->cd(4); hChi2->Draw();
254 //
255 TCanvas* CnvN = new TCanvas("CnvN", "Primary tracks found", 100, 100, 900, 500);
256 CnvN->Divide(2, 1);
257 CnvN->cd(1);
258 hTrPrim->Draw();
[46d3442]259 CnvN->cd(2);
[82db145]260 hTrFound->SetLineColor(kRed);
261 hTrFound->Draw();
262}
Note: See TracBrowser for help on using the repository browser.