Fork me on GitHub

source: git/examples/Example6.C

Last change on this file was dc883b4, checked in by Pavel Demin <pavel.demin@…>, 4 years ago

add info to TreeWriter
  • Property mode set to 100644
File size: 10.3 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 Z_{0}", 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)", 100, -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 magnetifc field
79 Double_t Bz = treeReader->GetInfo("Bz");
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 // Get associated generated particle
110 GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
111 //
112 // Position of origin in meters
113 Double_t x = 1.0e-3 * gp->X;
114 Double_t y = 1.0e-3 * gp->Y;
115 Double_t z = 1.0e-3 * gp->Z;
116 TVector3 xv(x, y, z);
117 //
118 // Momentum at origin
119 Double_t px = gp->Px;
120 Double_t py = gp->Py;
121 Double_t pz = gp->Pz;
122 TVector3 pv(px, py, pz);
123 //
124 // Unsmeared original parameters
125 Double_t Q = gp->Charge;
126 TVectorD gPar = TrkUtil::XPtoPar(xv, pv, Q, Bz); // Get parameters from position and momentum
127 TVectorD genPar = TrkUtil::ParToMm(gPar); // Convert to mm
128 // Unpack parameters
129 Double_t genD0 = genPar(0);
130 Double_t genPhi = genPar(1);
131 Double_t genC = genPar(2);
132 Double_t genZ0 = genPar(3);
133 Double_t genCtg = genPar(4);
134 Double_t genPt = pv.Pt();
135 //
136 // Calculate and fill pulls
137 Double_t pullD0 = (obsD0 - genD0) / trk->ErrorD0;
138 Double_t pullPhi = (obsPhi - genPhi) / trk->ErrorPhi;
139 Double_t pullC = (obsC - genC) / trk->ErrorC;
140 Double_t pullZ0 = (obsZ0 - genZ0) / trk->ErrorDZ;
141 Double_t pullCtg = (obsCtg - genCtg) / trk->ErrorCtgTheta;
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 histPtobs->Fill(obsPt);
152 if (TMath::Abs(genCtg) < maxCtAcc)histPtobsC->Fill(obsPt); // pt for central tracks
153 if (obsPt > minPtAcc)histCtobsH->Fill(obsCtg); // cot(theta) for high pt tracks
154 } // End loop on tracks
155
156 // If event contains at least 1 generated charged
157 //
158 if (branchGenPart->GetEntries() > 0)
159 {
160 // Loop on generated particles
161 for (Int_t it = 0; it < branchGenPart->GetEntries(); it++) {
162 GenParticle* gpart = (GenParticle*)branchGenPart->At(it);
163 //
164 // Plot charged particle parameters
165 // Only final state particles (Status = 1)
166 if (gpart->Status == 1 && TMath::Abs(gpart->Charge) > 0) {
167 //
168 // Position of origin in meters
169 Double_t x = 1.0e-3 * gpart->X;
170 Double_t y = 1.0e-3 * gpart->Y;
171 Double_t z = 1.0e-3 * gpart->Z;
172 TVector3 xv(x, y, z);
173 //
174 // Momentum at origin
175 Double_t px = gpart->Px;
176 Double_t py = gpart->Py;
177 Double_t pz = gpart->Pz;
178 TVector3 pv(px, py, pz);
179 //
180 // Original parameters
181 Double_t Q = gpart->Charge;
182 TVectorD gPar = TrkUtil::XPtoPar(xv, pv, Q, Bz); // Get parameters from position and momentum
183 TVectorD genPar = TrkUtil::ParToMm(gPar); // Convert to mm
184 // Unpack parameters
185 Double_t genD0 = genPar(0);
186 Double_t genPhi = genPar(1);
187 Double_t genC = genPar(2);
188 Double_t genZ0 = genPar(3);
189 Double_t genCtg = genPar(4);
190 Double_t genPt = pv.Pt();
191 // Fill histograms
192 histDgen ->Fill(genD0);
193 histP0gen->Fill(genPhi);
194 histCgen ->Fill(genC);
195 histZ0gen->Fill(genZ0);
196 histCtgen->Fill(genCtg);
197 //
198 // Acceptance plots
199 histPtgen->Fill(genPt);
200 if (TMath::Abs(genCtg) < maxCtAcc)histPtgenC->Fill(genPt);
201 if (genPt > minPtAcc)histCtgenH->Fill(genCtg);
202 }
203 }
204 }
205 }
206 }
207 //
208 // Show resulting histograms
209 //
210 TCanvas* Cnv = new TCanvas("Cnv", "Delphes generated track plots", 50, 50, 900, 500);
211 Cnv->Divide(3, 2);
212 Cnv->cd(1); histDgen->Draw();
213 Cnv->cd(2); histP0gen->Draw();
214 Cnv->cd(3); histCgen->Draw();
215 Cnv->cd(4); histZ0gen->Draw();
216 Cnv->cd(5); histCtgen->Draw();
217 //
218 TCanvas* Cnv1 = new TCanvas("Cnv1", "Delphes observed track plots", 100, 100, 900, 500);
219 Cnv1->Divide(3, 2);
220 Cnv1->cd(1); histDobs->Draw();
221 Cnv1->cd(2); histP0obs->Draw();
222 Cnv1->cd(3); histCobs->Draw();
223 Cnv1->cd(4); histZ0obs->Draw();
224 Cnv1->cd(5); histCtobs->Draw();
225 //
226 TCanvas* Cnv2 = new TCanvas("Cnv2", "Delphes observed track pulls", 150, 150, 900, 500);
227 Cnv2->Divide(3, 2);
228 Cnv2->cd(1); gPad->SetLogy(1);
229 histDpull->Draw();
230 Cnv2->cd(2); gPad->SetLogy(1);
231 histP0pull->Draw();
232 Cnv2->cd(3); gPad->SetLogy(1);
233 histCpull->Draw();
234 Cnv2->cd(4); gPad->SetLogy(1);
235 histZ0pull->Draw();
236 Cnv2->cd(5); gPad->SetLogy(1);
237 histCtpull->Draw();
238 //
239 // Acceptance plots
240 //
241 TCanvas* Cnv4 = new TCanvas("Cnv4", "Delphes acceptance", 200, 200, 500, 500);
242 Cnv4->Divide(2, 2);
243 Cnv4->cd(1);
244 histCtgen->Draw();
245 histCtobs->SetLineColor(kRed);
246 histCtobs->Draw("SAME");
247 Cnv4->cd(2);
248 Int_t NbCt = histCtgen->GetNbinsX();
249 for (Int_t i = 1; i < NbCt + 1; i++)
250 {
251 Float_t cgen = histCtgen->GetBinContent(i);
252 Float_t cobs = histCtobs->GetBinContent(i);
253 Float_t AccCtg = 0.0;
254 if (cgen > 0.0) AccCtg = cobs / cgen;
255 histAccCtg->SetBinContent(i, AccCtg);
256 }
257 histAccCtg->Draw();
258 Cnv4->cd(3); gPad->SetLogy(1);
259 histPtgen->Draw();
260 histPtobs->SetLineColor(kRed);
261 histPtobs->Draw("SAME");
262 Cnv4->cd(4);
263 Int_t NbPt = histPtgen->GetNbinsX();
264 for (Int_t i = 1; i < NbPt + 1; i++)
265 {
266 Float_t pgen = histPtgen->GetBinContent(i);
267 Float_t pobs = histPtobs->GetBinContent(i);
268 Float_t AccPt = 0.0;
269 if (pgen > 0.0) AccPt = pobs / pgen;
270 histAccPt->SetBinContent(i, AccPt);
271 }
272 histAccPt->Draw();
273 //
274 TCanvas* Cnv5 = new TCanvas("Cnv5", "Delphes acceptance (constrained)", 250, 250, 500, 500);
275 Cnv5->Divide(2, 2);
276 Cnv5->cd(1);
277 histCtgenH->Draw();
278 histCtobsH->SetLineColor(kRed);
279 histCtobsH->Draw("SAME");
280 Cnv5->cd(2);
281 NbCt = histCtgenH->GetNbinsX();
282 for (Int_t i = 1; i < NbCt + 1; i++)
283 {
284 Float_t cgen = histCtgenH->GetBinContent(i);
285 Float_t cobs = histCtobsH->GetBinContent(i);
286 Float_t AccCtg = 0.0;
287 if (cgen > 0.0) AccCtg = cobs / cgen;
288 histAccCtgH->SetBinContent(i, AccCtg);
289 }
290 histAccCtgH->Draw();
291 Cnv5->cd(3); gPad->SetLogy(1);
292 histPtgenC->Draw();
293 histPtobsC->SetLineColor(kRed);
294 histPtobsC->Draw("SAME");
295 Cnv5->cd(4);
296 NbPt = histPtgenC->GetNbinsX();
297 for (Int_t i = 1; i < NbPt + 1; i++)
298 {
299 Float_t pgen = histPtgenC->GetBinContent(i);
300 Float_t pobs = histPtobsC->GetBinContent(i);
301 Float_t AccPt = 0.0;
302 if (pgen > 0.0) AccPt = pobs / pgen;
303 histAccPtC->SetBinContent(i, AccPt);
304 }
305 histAccPtC->Draw();
306}
Note: See TracBrowser for help on using the repository browser.