Fork me on GitHub

source: git/examples/Example6.C@ bfc6722

Last change on this file since bfc6722 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
RevLine 
[40e890c]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.);
[df408d2]53 TH1* histP0pull = new TH1F("h_P0pull", "Pull #phi_{0}", 100, -10., 10.);
[40e890c]54 TH1* histCpull = new TH1F("h_Cpull", "Pull half curvature", 100, -10., 10.);
[df408d2]55 TH1* histZ0pull = new TH1F("h_Z0pull", "Pull Z_{0}", 100, -10., 10.);
[40e890c]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;
[dc883b4]66 TH1* histPtobsC = new TH1F("h_PtobsC", "Reconstructed Pt - Central", 500, 0., 50.); // pt for central
[40e890c]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
[df408d2]73 TH1* histAccCtgH = new TH1F("h_AccCtgH", "Cotangent acceptance (high pt)", 100, -10., 10.);
[40e890c]74 // tracks above min angle
75 TH1* histAccPtC = new TH1F("h_AccPtC", "Pt acceptance (central tracks)", 500, 0., 50.);
76
77 //
[df408d2]78 // Get magnetifc field
[dc883b4]79 Double_t Bz = treeReader->GetInfo("Bz");
[40e890c]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);
[df408d2]108 //
[dc883b4]109 // Get associated generated particle
[40e890c]110 GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
[df408d2]111 //
[40e890c]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 //
[df408d2]143 histDpull ->Fill(pullD0);
[40e890c]144 histP0pull->Fill(pullPhi);
[df408d2]145 histCpull ->Fill(pullC);
[40e890c]146 histZ0pull->Fill(pullZ0);
147 histCtpull->Fill(pullCtg);
148 //
149 // Acceptance plots
150 Double_t obsPt = trk->PT;
[df408d2]151 histPtobs->Fill(obsPt);
[40e890c]152 if (TMath::Abs(genCtg) < maxCtAcc)histPtobsC->Fill(obsPt); // pt for central tracks
[df408d2]153 if (obsPt > minPtAcc)histCtobsH->Fill(obsCtg); // cot(theta) for high pt tracks
[40e890c]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 //
[dc883b4]164 // Plot charged particle parameters
[40e890c]165 // Only final state particles (Status = 1)
166 if (gpart->Status == 1 && TMath::Abs(gpart->Charge) > 0) {
167 //
[df408d2]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;
[40e890c]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;
[df408d2]182 TVectorD gPar = TrkUtil::XPtoPar(xv, pv, Q, Bz); // Get parameters from position and momentum
183 TVectorD genPar = TrkUtil::ParToMm(gPar); // Convert to mm
[40e890c]184 // Unpack parameters
[df408d2]185 Double_t genD0 = genPar(0);
[40e890c]186 Double_t genPhi = genPar(1);
[df408d2]187 Double_t genC = genPar(2);
188 Double_t genZ0 = genPar(3);
[40e890c]189 Double_t genCtg = genPar(4);
[df408d2]190 Double_t genPt = pv.Pt();
[40e890c]191 // Fill histograms
[df408d2]192 histDgen ->Fill(genD0);
[40e890c]193 histP0gen->Fill(genPhi);
[df408d2]194 histCgen ->Fill(genC);
[40e890c]195 histZ0gen->Fill(genZ0);
196 histCtgen->Fill(genCtg);
197 //
198 // Acceptance plots
[df408d2]199 histPtgen->Fill(genPt);
200 if (TMath::Abs(genCtg) < maxCtAcc)histPtgenC->Fill(genPt);
201 if (genPt > minPtAcc)histCtgenH->Fill(genCtg);
[40e890c]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);
[df408d2]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();
[40e890c]238 //
239 // Acceptance plots
240 //
[df408d2]241 TCanvas* Cnv4 = new TCanvas("Cnv4", "Delphes acceptance", 200, 200, 500, 500);
[40e890c]242 Cnv4->Divide(2, 2);
243 Cnv4->cd(1);
244 histCtgen->Draw();
[df408d2]245 histCtobs->SetLineColor(kRed);
[40e890c]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();
[df408d2]258 Cnv4->cd(3); gPad->SetLogy(1);
[40e890c]259 histPtgen->Draw();
[df408d2]260 histPtobs->SetLineColor(kRed);
[40e890c]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 //
[df408d2]274 TCanvas* Cnv5 = new TCanvas("Cnv5", "Delphes acceptance (constrained)", 250, 250, 500, 500);
[40e890c]275 Cnv5->Divide(2, 2);
276 Cnv5->cd(1);
277 histCtgenH->Draw();
[df408d2]278 histCtobsH->SetLineColor(kRed);
[40e890c]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();
[df408d2]291 Cnv5->cd(3); gPad->SetLogy(1);
[dc883b4]292 histPtgenC->Draw();
[df408d2]293 histPtobsC->SetLineColor(kRed);
[40e890c]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.