Fork me on GitHub

source: git/examples/Example6.C@ c9e7466

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

fixed trailing printouts
  • Property mode set to 100644
File size: 9.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 //
143 histDpull->Fill(pullD0);
144 histP0pull->Fill(pullPhi);
145 histCpull->Fill(pullC);
146 histZ0pull->Fill(pullZ0);
147 histCtpull->Fill(pullCtg);
148 //
149 // Acceptance plots
150 Double_t obsPt = trk->PT;
151 if (TMath::Abs(genCtg) < maxCtAcc)histPtobsC->Fill(obsPt); // pt for central tracks
152 if (genPt > minPtAcc)histCtobsH->Fill(obsCtg); // cot(theta) for high pt tracks
153 } // End loop on tracks
154
155 // If event contains at least 1 generated charged
156 //
157 if (branchGenPart->GetEntries() > 0)
158 {
159 // Loop on generated particles
160 for (Int_t it = 0; it < branchGenPart->GetEntries(); it++) {
161 GenParticle* gpart = (GenParticle*)branchGenPart->At(it);
162 //
163 // Plot charged particle parameters
164 // Only final state particles (Status = 1)
165 if (gpart->Status == 1 && TMath::Abs(gpart->Charge) > 0) {
166 //
167 // Position of origin
168 Double_t x = gpart->X;
169 Double_t y = gpart->Y;
170 Double_t z = gpart->Z;
171 TVector3 xv(x, y, z);
172 //
173 // Momentum at origin
174 Double_t px = gpart->Px;
175 Double_t py = gpart->Py;
176 Double_t pz = gpart->Pz;
177 TVector3 pv(px, py, pz);
178 //
179 // Original parameters
180 Double_t Q = gpart->Charge;
181 TVectorD genPar = TrkUtil::XPtoPar(xv, pv, Q, Bz); // Get parameters from position and momentum
182 // Unpack parameters
183 Double_t genD0 = genPar(0);
184 Double_t genPhi = genPar(1);
185 Double_t genC = genPar(2);
186 Double_t genZ0 = genPar(3);
187 Double_t genCtg = genPar(4);
188 // Fill histograms
189 histDgen->Fill(genD0);
190 histP0gen->Fill(genPhi);
191 histCgen->Fill(genC);
192 histZ0gen->Fill(genZ0);
193 histCtgen->Fill(genCtg);
194 //
195 // Acceptance plots
196 if (TMath::Abs(genCtg) < maxCtAcc)histPtgenC->Fill(gpart->PT);
197 if (gpart->PT > minPtAcc)histCtgenH->Fill(genCtg);
198 }
199 }
200 }
201 }
202 }
203 //
204 // Show resulting histograms
205 //
206 TCanvas* Cnv = new TCanvas("Cnv", "Delphes generated track plots", 50, 50, 900, 500);
207 Cnv->Divide(3, 2);
208 Cnv->cd(1); histDgen->Draw();
209 Cnv->cd(2); histP0gen->Draw();
210 Cnv->cd(3); histCgen->Draw();
211 Cnv->cd(4); histZ0gen->Draw();
212 Cnv->cd(5); histCtgen->Draw();
213 //
214 TCanvas* Cnv1 = new TCanvas("Cnv1", "Delphes observed track plots", 100, 100, 900, 500);
215 Cnv1->Divide(3, 2);
216 Cnv1->cd(1); histDobs->Draw();
217 Cnv1->cd(2); histP0obs->Draw();
218 Cnv1->cd(3); histCobs->Draw();
219 Cnv1->cd(4); histZ0obs->Draw();
220 Cnv1->cd(5); histCtobs->Draw();
221 //
222 TCanvas* Cnv2 = new TCanvas("Cnv2", "Delphes observed track pulls", 150, 150, 900, 500);
223 Cnv2->Divide(3, 2);
224 Cnv2->cd(1); histDpull->Draw();
225 Cnv2->cd(2); histP0pull->Draw();
226 Cnv2->cd(3); histCpull->Draw();
227 Cnv2->cd(4); histZ0pull->Draw();
228 Cnv2->cd(5); histCtpull->Draw();
229 //
230 // Acceptance plots
231 //
232 TCanvas* Cnv4 = new TCanvas("Cnv4", "Delphes acceptance", 350, 350, 500, 500);
233 Cnv4->Divide(2, 2);
234 Cnv4->cd(1);
235 histCtgen->Draw();
236 histCtobs->Draw("SAME");
237 Cnv4->cd(2);
238 Int_t NbCt = histCtgen->GetNbinsX();
239 for (Int_t i = 1; i < NbCt + 1; i++)
240 {
241 Float_t cgen = histCtgen->GetBinContent(i);
242 Float_t cobs = histCtobs->GetBinContent(i);
243 Float_t AccCtg = 0.0;
244 if (cgen > 0.0) AccCtg = cobs / cgen;
245 histAccCtg->SetBinContent(i, AccCtg);
246 }
247 histAccCtg->Draw();
248 Cnv4->cd(3);
249 histPtgen->Draw();
250 histPtobs->Draw("SAME");
251 Cnv4->cd(4);
252 Int_t NbPt = histPtgen->GetNbinsX();
253 for (Int_t i = 1; i < NbPt + 1; i++)
254 {
255 Float_t pgen = histPtgen->GetBinContent(i);
256 Float_t pobs = histPtobs->GetBinContent(i);
257 Float_t AccPt = 0.0;
258 if (pgen > 0.0) AccPt = pobs / pgen;
259 histAccPt->SetBinContent(i, AccPt);
260 }
261 histAccPt->Draw();
262 //
263 TCanvas* Cnv5 = new TCanvas("Cnv5", "Delphes acceptance (constrained)", 400, 400, 500, 500);
264 Cnv5->Divide(2, 2);
265 Cnv5->cd(1);
266 histCtgenH->Draw();
267 histCtobsH->Draw("SAME");
268 Cnv5->cd(2);
269 NbCt = histCtgenH->GetNbinsX();
270 for (Int_t i = 1; i < NbCt + 1; i++)
271 {
272 Float_t cgen = histCtgenH->GetBinContent(i);
273 Float_t cobs = histCtobsH->GetBinContent(i);
274 Float_t AccCtg = 0.0;
275 if (cgen > 0.0) AccCtg = cobs / cgen;
276 histAccCtgH->SetBinContent(i, AccCtg);
277 }
278 histAccCtgH->Draw();
279 Cnv5->cd(3);
280 histPtgenC->Draw();
281 histPtobsC->Draw("SAME");
282 Cnv5->cd(4);
283 NbPt = histPtgenC->GetNbinsX();
284 for (Int_t i = 1; i < NbPt + 1; i++)
285 {
286 Float_t pgen = histPtgenC->GetBinContent(i);
287 Float_t pobs = histPtobsC->GetBinContent(i);
288 Float_t AccPt = 0.0;
289 if (pgen > 0.0) AccPt = pobs / pgen;
290 histAccPtC->SetBinContent(i, AccPt);
291 }
292 histAccPtC->Draw();
293}
Note: See TracBrowser for help on using the repository browser.