Fork me on GitHub

Ticket #1079: Example2.C

File Example2.C, 4.8 KB (added by Amin Abou Ibrahim, 8 years ago)
Line 
1/*
2Simple macro showing how to access branches from the delphes output root file,
3loop over events, store histograms in a root file and print them as image files.
4
5root -l examples/Example2.C'("delphes_output.root")'
6*/
7
8#include "TH1.h"
9#include "TSystem.h"
10
11#ifdef __CLING__
12R__LOAD_LIBRARY(libDelphes)
13#include "classes/DelphesClasses.h"
14//R__ADD_INCLUDE_PATH(/home/amin/Documents/Delphes-3.4.0/external)
15#include "external/ExRootAnalysis/ExRootTreeReader.h"
16#include "external/ExRootAnalysis/ExRootResult.h"
17#endif
18
19//------------------------------------------------------------------------------
20
21struct MyPlots
22{
23 TH1 *fJetPT[2];
24 TH1 *fMissingET;
25 TH1 *fElectronPT;
26};
27
28//------------------------------------------------------------------------------
29
30class ExRootResult;
31class ExRootTreeReader;
32
33//------------------------------------------------------------------------------
34
35void BookHistograms(ExRootResult *result, MyPlots *plots)
36{
37 THStack *stack;
38 TLegend *legend;
39 TPaveText *comment;
40
41 // book 2 histograms for PT of 1st and 2nd leading jets
42
43 plots->fJetPT[0] = result->AddHist1D(
44 "jet_pt_0", "leading jet P_{T}",
45 "jet P_{T}, GeV/c", "number of jets",
46 50, 0.0, 100.0);
47
48 plots->fJetPT[1] = result->AddHist1D(
49 "jet_pt_1", "2nd leading jet P_{T}",
50 "jet P_{T}, GeV/c", "number of jets",
51 50, 0.0, 100.0);
52
53 plots->fJetPT[0]->SetLineColor(kRed);
54 plots->fJetPT[1]->SetLineColor(kBlue);
55
56 // book 1 stack of 2 histograms
57
58 stack = result->AddHistStack("jet_pt_all", "1st and 2nd jets P_{T}");
59 stack->Add(plots->fJetPT[0]);
60 stack->Add(plots->fJetPT[1]);
61
62 // book legend for stack of 2 histograms
63
64 legend = result->AddLegend(0.72, 0.86, 0.98, 0.98);
65 legend->AddEntry(plots->fJetPT[0], "leading jet", "l");
66 legend->AddEntry(plots->fJetPT[1], "second jet", "l");
67
68 // attach legend to stack (legend will be printed over stack in .eps file)
69
70 result->Attach(stack, legend);
71
72 // book more histograms
73
74 plots->fElectronPT = result->AddHist1D(
75 "electron_pt", "electron P_{T}",
76 "electron P_{T}, GeV/c", "number of electrons",
77 50, 0.0, 100.0);
78
79 plots->fMissingET = result->AddHist1D(
80 "missing_et", "Missing E_{T}",
81 "Missing E_{T}, GeV", "number of events",
82 60, 0.0, 30.0);
83
84 // book general comment
85
86 comment = result->AddComment(0.64, 0.86, 0.98, 0.98);
87 comment->AddText("demonstration plot");
88 comment->AddText("produced by Example2.C");
89
90 // attach comment to single histograms
91
92 result->Attach(plots->fJetPT[0], comment);
93 result->Attach(plots->fJetPT[1], comment);
94 result->Attach(plots->fElectronPT, comment);
95
96 // show histogram statisics for MissingET
97 plots->fMissingET->SetStats();
98}
99
100//------------------------------------------------------------------------------
101
102void AnalyseEvents(ExRootTreeReader *treeReader, MyPlots *plots)
103{
104 TClonesArray *branchJet = treeReader->UseBranch("Jet");
105 TClonesArray *branchElectron = treeReader->UseBranch("Electron");
106 TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");
107
108 Long64_t allEntries = treeReader->GetEntries();
109
110 cout << "** Chain contains " << allEntries << " events" << endl;
111
112 Jet *jet[2];
113 MissingET *met;
114 Electron *electron;
115
116 Long64_t entry;
117
118 Int_t i;
119
120 // Loop over all events
121 for(entry = 0; entry < allEntries; ++entry)
122 {
123 // Load selected branches with data from specified event
124 treeReader->ReadEntry(entry);
125
126 // Analyse two leading jets
127 if(branchJet->GetEntriesFast() >= 2)
128 {
129 jet[0] = (Jet*) branchJet->At(0);
130 jet[1] = (Jet*) branchJet->At(1);
131
132 plots->fJetPT[0]->Fill(jet[0]->PT);
133 plots->fJetPT[1]->Fill(jet[1]->PT);
134 }
135
136 // Analyse missing ET
137 if(branchMissingET->GetEntriesFast() > 0)
138 {
139 met = (MissingET*) branchMissingET->At(0);
140 plots->fMissingET->Fill(met->MET);
141 }
142
143 // Loop over all electrons in event
144 for(i = 0; i < branchElectron->GetEntriesFast(); ++i)
145 {
146 electron = (Electron*) branchElectron->At(i);
147 plots->fElectronPT->Fill(electron->PT);
148 }
149 }
150}
151
152//------------------------------------------------------------------------------
153
154void PrintHistograms(ExRootResult *result, MyPlots *plots)
155{
156 result->Print("png");
157}
158
159//------------------------------------------------------------------------------
160
161void Example2(const char *inputFile)
162{
163 gSystem->Load("libDelphes");
164
165 TChain *chain = new TChain("Delphes");
166 chain->Add(inputFile);
167
168 ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
169 ExRootResult *result = new ExRootResult();
170
171 MyPlots *plots = new MyPlots;
172
173 BookHistograms(result, plots);
174
175 AnalyseEvents(treeReader, plots);
176
177 PrintHistograms(result, plots);
178
179 result->Write("results.root");
180
181 cout << "** Exiting..." << endl;
182
183 delete plots;
184 delete result;
185 delete treeReader;
186 delete chain;
187}
188
189//------------------------------------------------------------------------------