Fork me on GitHub

source: git/examples/Example6.C@ 825beb7

Last change on this file since 825beb7 was 40e890c, checked in by michele <michele.selvaggi@…>, 4 years ago

added Example6 for validating track covariance module (F. Bedeschi)
  • Property mode set to 100644
File size: 10.8 KB
Line 
1/*
2Simple macro showing how to access branches from the delphes output root file,
3loop over events, and plot track pulls and geometrical acceptance.
4
5root -l examples/Example6.C'("delphes_output.root")'
6*/
7
8#ifdef __CLING__
9R__LOAD_LIBRARY(libDelphes)
10#include "classes/DelphesClasses.h"
11#include "external/ExRootAnalysis/ExRootTreeReader.h"
12#include "modules/TrackCovariance.h"
13#include "external/TrackCovariance/TrkUtil.h"
14#endif
15
16
17//------------------------------------------------------------------------------
18
19void Example6(const char* inputFile)
20{
21 gSystem->Load("libDelphes");
22
23 // Create chain of root trees
24 TChain chain("Delphes");
25 chain.Add(inputFile);
26
27 // Create object of class ExRootTreeReader
28 ExRootTreeReader* treeReader = new ExRootTreeReader(&chain);
29 Long64_t numberOfEntries = treeReader->GetEntries();
30
31 // Get pointers to branches used in this analysis
32 TClonesArray* branchGenPart = treeReader->UseBranch("Particle");
33 TClonesArray* branchTrack = treeReader->UseBranch("Track");
34
35 // Book histograms
36 //
37 // Generated track parameters
38 TH1* histDgen = new TH1F("h_Dgen", "Generated impact parameter", 100, -2., 2.); // mm
39 TH1* histP0gen = new TH1F("h_P0gen", "Generated #phi_{0}", 100, -TMath::Pi(), TMath::Pi()); // rad
40 TH1* histCgen = new TH1F("h_Cgen", "Generated half curvature", 200, -.01, .01); // mm
41 TH1* histZ0gen = new TH1F("h_Z0gen", "Generated Z_{0}", 100, -2., 2.); // mm
42 TH1* histCtgen = new TH1F("h_Ctgen", "Generated cotg(#theta)", 100, -10., 10.);
43 //
44 // Reconstructed impact parameters
45 TH1* histDobs = new TH1F("h_Dobs", "Reconstructed impact parameter", 100, -2., 2.); // mm
46 TH1* histP0obs = new TH1F("h_P0obs", "Reconstructed #phi_{0}", 100, -TMath::Pi(), TMath::Pi()); // rad
47 TH1* histCobs = new TH1F("h_Cobs", "Reconstructed half curvature", 200, -.01, .01); // mm
48 TH1* histZ0obs = new TH1F("h_Z0obs", "Reconstructed Z_{0}", 100, -2., 2.); // mm
49 TH1* histCtobs = new TH1F("h_Ctobs", "Reconstructed cotg(#theta)", 100, -10., 10.);
50 //
51 // Track parameter pulls
52 TH1* histDpull = new TH1F("h_Dpull", "Pull impact parameter", 100, -10., 10.);
53 TH1* histP0pull = new TH1F("h_P0pull", "Pull phi_{0}", 100, -10., 10.);
54 TH1* histCpull = new TH1F("h_Cpull", "Pull half curvature", 100, -10., 10.);
55 TH1* histZ0pull = new TH1F("h_Z0pull", "Pull Z0", 100, -10., 10.);
56 TH1* histCtpull = new TH1F("h_Ctpull", "Pull cotg(#theta)", 100, -10., 10.);
57 //
58 // Accptance plots
59 Double_t minPtAcc = 1.0; // 1 GeV
60 Double_t maxCtAcc = 1.0; // 45 degrees
61 // All tracks
62 TH1* histPtgen = new TH1F("h_Ptgen", "Generated Pt", 500, 0., 50.);
63 TH1* histPtobs = new TH1F("h_Ptobs", "Reconstructed Pt", 500, 0., 50.);
64 // CEntral track over min Pt
65 TH1* histPtgenC = new TH1F("h_PtgenC", "Generated Pt - Central", 500, 0., 50.); // pt for central tracks;
66 TH1* histPtobsC = new TH1F("h_PtobsC", "Reconstructed Pt - Central", 500, 0., 50.); // pt for central
67 TH1* histCtgenH = new TH1F("h_CtgenH", "Generateded Cotangent", 100, -10., 10.); // cot(theta) for high pt tracks;
68 TH1* histCtobsH = new TH1F("h_CtobsH", "Reconstructed Cotangent", 100, -10., 10.); // cot(theta) for high pt tracks
69 // All tracks
70 TH1* histAccCtg = new TH1F("h_AccCtg", "Cotangent acceptance", 100, -10., 10.);
71 TH1* histAccPt = new TH1F("h_AccPt", "Pt acceptance", 500, 0., 50.);
72 // Tracks over min pt
73 TH1* histAccCtgH = new TH1F("h_AccCtgH", "Cotangent acceptance (high pt)", 500, -10., 10.);
74 // tracks above min angle
75 TH1* histAccPtC = new TH1F("h_AccPtC", "Pt acceptance (central tracks)", 500, 0., 50.);
76
77 //
78 // Get magnetic field
79 Double_t Bz = 2.0;
80
81 // Loop over all events
82 for (Int_t entry = 0; entry < numberOfEntries; ++entry)
83 {
84 // Load selected branches with data from specified event
85 treeReader->ReadEntry(entry);
86
87 // If event contains at least 1 track
88 //
89 if (branchTrack->GetEntries() > 0)
90 {
91 // Loop on tracks
92 for (Int_t it = 0; it < branchTrack->GetEntries(); it++)
93 {
94 Track* trk = (Track*)branchTrack->At(it);
95 //
96 // Reconstructed track parameters
97 Double_t obsD0 = trk->D0;
98 Double_t obsPhi = trk->Phi;
99 Double_t obsC = trk->C;
100 Double_t obsZ0 = trk->DZ;
101 Double_t obsCtg = trk->CtgTheta;
102 // Fill histograms
103 histDobs ->Fill(obsD0);
104 histP0obs->Fill(obsPhi);
105 histCobs ->Fill(obsC);
106 histZ0obs->Fill(obsZ0);
107 histCtobs->Fill(obsCtg);
108
109 GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
110
111 // Position of origin in meters
112 Double_t x = 1.0e-3 * gp->X;
113 Double_t y = 1.0e-3 * gp->Y;
114 Double_t z = 1.0e-3 * gp->Z;
115 TVector3 xv(x, y, z);
116 //
117 // Momentum at origin
118 Double_t px = gp->Px;
119 Double_t py = gp->Py;
120 Double_t pz = gp->Pz;
121 TVector3 pv(px, py, pz);
122 //
123 // Unsmeared original parameters
124 Double_t Q = gp->Charge;
125 TVectorD gPar = TrkUtil::XPtoPar(xv, pv, Q, Bz); // Get parameters from position and momentum
126 TVectorD genPar = TrkUtil::ParToMm(gPar); // Convert to mm
127 // Unpack parameters
128 Double_t genD0 = genPar(0);
129 Double_t genPhi = genPar(1);
130 Double_t genC = genPar(2);
131 Double_t genZ0 = genPar(3);
132 Double_t genCtg = genPar(4);
133 Double_t genPt = pv.Pt();
134 //
135 // Calculate and fill pulls
136 Double_t pullD0 = (obsD0 - genD0) / trk->ErrorD0;
137 Double_t pullPhi = (obsPhi - genPhi) / trk->ErrorPhi;
138 Double_t pullC = (obsC - genC) / trk->ErrorC;
139 Double_t pullZ0 = (obsZ0 - genZ0) / trk->ErrorDZ;
140 Double_t pullCtg = (obsCtg - genCtg) / trk->ErrorCtgTheta;
141
142 if (fabs(pullZ0) > 5)
143 {
144 // Get associated generated particle
145 Int_t unique = gp->GetUniqueID();
146 //cout<<unique<<endl;
147
148 if (branchGenPart->GetEntries() > 0)
149 {
150
151 // Loop on tracks
152 for (Int_t itp = 0; itp < branchGenPart->GetEntries(); itp++)
153 {
154 GenParticle* part = (GenParticle*)branchGenPart->At(itp);
155 if (part->GetUniqueID() == unique)
156 {
157 cout<<" ------------ found: ------------------------"<<unique<<endl;
158
159 cout<<gp->X<<","<<part->X<<","<<TMath::Power((gp->X - part->X),2)*1e6<<endl;
160 cout<<gp->Y<<","<<part->Y<<","<<TMath::Power((gp->Y - part->Y),2)*1e6<<endl;
161 cout<<gp->Z<<","<<part->Z<<","<<TMath::Power((gp->Z - part->Z),2)*1e6<<endl;
162 cout<<gp->Px<<","<<part->Px<<","<<TMath::Power((gp->Px - part->Px),2)*1e6<<endl;
163 cout<<gp->Py<<","<<part->Py<<","<<TMath::Power((gp->Py - part->Py),2)*1e6<<endl;
164 cout<<gp->Pz<<","<<part->Pz<<","<<TMath::Power((gp->Pz - part->Pz),2)*1e6<<endl;
165
166 }
167 }
168 }
169
170
171 }
172
173
174 //
175 histDpull->Fill(pullD0);
176 histP0pull->Fill(pullPhi);
177 histCpull->Fill(pullC);
178 histZ0pull->Fill(pullZ0);
179 histCtpull->Fill(pullCtg);
180 //
181 // Acceptance plots
182 Double_t obsPt = trk->PT;
183 if (TMath::Abs(genCtg) < maxCtAcc)histPtobsC->Fill(obsPt); // pt for central tracks
184 if (genPt > minPtAcc)histCtobsH->Fill(obsCtg); // cot(theta) for high pt tracks
185 } // End loop on tracks
186
187 // If event contains at least 1 generated charged
188 //
189 if (branchGenPart->GetEntries() > 0)
190 {
191 // Loop on generated particles
192 for (Int_t it = 0; it < branchGenPart->GetEntries(); it++) {
193 GenParticle* gpart = (GenParticle*)branchGenPart->At(it);
194 //
195 // Plot charged particle parameters
196 // Only final state particles (Status = 1)
197 if (gpart->Status == 1 && TMath::Abs(gpart->Charge) > 0) {
198 //
199 // Position of origin
200 Double_t x = gpart->X;
201 Double_t y = gpart->Y;
202 Double_t z = gpart->Z;
203 TVector3 xv(x, y, z);
204 //
205 // Momentum at origin
206 Double_t px = gpart->Px;
207 Double_t py = gpart->Py;
208 Double_t pz = gpart->Pz;
209 TVector3 pv(px, py, pz);
210 //
211 // Original parameters
212 Double_t Q = gpart->Charge;
213 TVectorD genPar = TrkUtil::XPtoPar(xv, pv, Q, Bz); // Get parameters from position and momentum
214 // Unpack parameters
215 Double_t genD0 = genPar(0);
216 Double_t genPhi = genPar(1);
217 Double_t genC = genPar(2);
218 Double_t genZ0 = genPar(3);
219 Double_t genCtg = genPar(4);
220 // Fill histograms
221 histDgen->Fill(genD0);
222 histP0gen->Fill(genPhi);
223 histCgen->Fill(genC);
224 histZ0gen->Fill(genZ0);
225 histCtgen->Fill(genCtg);
226 //
227 // Acceptance plots
228 if (TMath::Abs(genCtg) < maxCtAcc)histPtgenC->Fill(gpart->PT);
229 if (gpart->PT > minPtAcc)histCtgenH->Fill(genCtg);
230 }
231 }
232 }
233 }
234 }
235 //
236 // Show resulting histograms
237 //
238 TCanvas* Cnv = new TCanvas("Cnv", "Delphes generated track plots", 50, 50, 900, 500);
239 Cnv->Divide(3, 2);
240 Cnv->cd(1); histDgen->Draw();
241 Cnv->cd(2); histP0gen->Draw();
242 Cnv->cd(3); histCgen->Draw();
243 Cnv->cd(4); histZ0gen->Draw();
244 Cnv->cd(5); histCtgen->Draw();
245 //
246 TCanvas* Cnv1 = new TCanvas("Cnv1", "Delphes observed track plots", 100, 100, 900, 500);
247 Cnv1->Divide(3, 2);
248 Cnv1->cd(1); histDobs->Draw();
249 Cnv1->cd(2); histP0obs->Draw();
250 Cnv1->cd(3); histCobs->Draw();
251 Cnv1->cd(4); histZ0obs->Draw();
252 Cnv1->cd(5); histCtobs->Draw();
253 //
254 TCanvas* Cnv2 = new TCanvas("Cnv2", "Delphes observed track pulls", 150, 150, 900, 500);
255 Cnv2->Divide(3, 2);
256 Cnv2->cd(1); histDpull->Draw();
257 Cnv2->cd(2); histP0pull->Draw();
258 Cnv2->cd(3); histCpull->Draw();
259 Cnv2->cd(4); histZ0pull->Draw();
260 Cnv2->cd(5); histCtpull->Draw();
261 //
262 // Acceptance plots
263 //
264 TCanvas* Cnv4 = new TCanvas("Cnv4", "Delphes acceptance", 350, 350, 500, 500);
265 Cnv4->Divide(2, 2);
266 Cnv4->cd(1);
267 histCtgen->Draw();
268 histCtobs->Draw("SAME");
269 Cnv4->cd(2);
270 Int_t NbCt = histCtgen->GetNbinsX();
271 for (Int_t i = 1; i < NbCt + 1; i++)
272 {
273 Float_t cgen = histCtgen->GetBinContent(i);
274 Float_t cobs = histCtobs->GetBinContent(i);
275 Float_t AccCtg = 0.0;
276 if (cgen > 0.0) AccCtg = cobs / cgen;
277 histAccCtg->SetBinContent(i, AccCtg);
278 }
279 histAccCtg->Draw();
280 Cnv4->cd(3);
281 histPtgen->Draw();
282 histPtobs->Draw("SAME");
283 Cnv4->cd(4);
284 Int_t NbPt = histPtgen->GetNbinsX();
285 for (Int_t i = 1; i < NbPt + 1; i++)
286 {
287 Float_t pgen = histPtgen->GetBinContent(i);
288 Float_t pobs = histPtobs->GetBinContent(i);
289 Float_t AccPt = 0.0;
290 if (pgen > 0.0) AccPt = pobs / pgen;
291 histAccPt->SetBinContent(i, AccPt);
292 }
293 histAccPt->Draw();
294 //
295 TCanvas* Cnv5 = new TCanvas("Cnv5", "Delphes acceptance (constrained)", 400, 400, 500, 500);
296 Cnv5->Divide(2, 2);
297 Cnv5->cd(1);
298 histCtgenH->Draw();
299 histCtobsH->Draw("SAME");
300 Cnv5->cd(2);
301 NbCt = histCtgenH->GetNbinsX();
302 for (Int_t i = 1; i < NbCt + 1; i++)
303 {
304 Float_t cgen = histCtgenH->GetBinContent(i);
305 Float_t cobs = histCtobsH->GetBinContent(i);
306 Float_t AccCtg = 0.0;
307 if (cgen > 0.0) AccCtg = cobs / cgen;
308 histAccCtgH->SetBinContent(i, AccCtg);
309 }
310 histAccCtgH->Draw();
311 Cnv5->cd(3);
312 histPtgenC->Draw();
313 histPtobsC->Draw("SAME");
314 Cnv5->cd(4);
315 NbPt = histPtgenC->GetNbinsX();
316 for (Int_t i = 1; i < NbPt + 1; i++)
317 {
318 Float_t pgen = histPtgenC->GetBinContent(i);
319 Float_t pobs = histPtobsC->GetBinContent(i);
320 Float_t AccPt = 0.0;
321 if (pgen > 0.0) AccPt = pobs / pgen;
322 histAccPtC->SetBinContent(i, AccPt);
323 }
324 histAccPtC->Draw();
325}
Note: See TracBrowser for help on using the repository browser.