1 | /* |
---|
2 | root -l examples/Example3.C\(\"delphes_output.root\"\) |
---|
3 | */ |
---|
4 | |
---|
5 | //------------------------------------------------------------------------------ |
---|
6 | |
---|
7 | struct TestPlots |
---|
8 | { |
---|
9 | TH1 *fElectronDeltaPT; |
---|
10 | TH1 *fElectronDeltaEta; |
---|
11 | |
---|
12 | TH1 *fPhotonDeltaPT; |
---|
13 | TH1 *fPhotonDeltaEta; |
---|
14 | TH1 *fPhotonDeltaE; |
---|
15 | |
---|
16 | TH1 *fMuonDeltaPT; |
---|
17 | TH1 *fMuonDeltaEta; |
---|
18 | |
---|
19 | TH1 *fTrackDeltaPT; |
---|
20 | TH1 *fTrackDeltaEta; |
---|
21 | |
---|
22 | TH1 *fTowerDeltaEem; |
---|
23 | TH1 *fTowerDeltaEhad; |
---|
24 | |
---|
25 | TH1 *fJetDeltaPT; |
---|
26 | }; |
---|
27 | |
---|
28 | //------------------------------------------------------------------------------ |
---|
29 | |
---|
30 | class ExRootResult; |
---|
31 | class ExRootTreeReader; |
---|
32 | |
---|
33 | //------------------------------------------------------------------------------ |
---|
34 | |
---|
35 | void BookHistograms(ExRootResult *result, TestPlots *plots) |
---|
36 | { |
---|
37 | TLegend *legend; |
---|
38 | TPaveText *comment; |
---|
39 | |
---|
40 | plots->fElectronDeltaPT = result->AddHist1D( |
---|
41 | "electron delta pt", "(p_{T}^{particle} - p_{T}^{electron})/p_{T}^{particle}", |
---|
42 | "(p_{T}^{particle} - p_{T}^{electron})/p_{T}^{particle}", "number of electrons", |
---|
43 | 100, -0.1, 0.1); |
---|
44 | |
---|
45 | plots->fElectronDeltaEta = result->AddHist1D( |
---|
46 | "electron delta eta", "(#eta^{particle} - #eta^{electron})/#eta^{particle}", |
---|
47 | "(#eta^{particle} - #eta^{electron})/#eta^{particle}", "number of electrons", |
---|
48 | 100, -0.1, 0.1); |
---|
49 | |
---|
50 | plots->fPhotonDeltaPT = result->AddHist1D( |
---|
51 | "photon delta pt", "(p_{T}^{particle} - p_{T}^{photon})/p_{T}^{particle}", |
---|
52 | "(p_{T}^{particle} - p_{T}^{photon})/p_{T}^{particle}", "number of photons", |
---|
53 | 100, -0.1, 0.1); |
---|
54 | |
---|
55 | plots->fPhotonDeltaEta = result->AddHist1D( |
---|
56 | "photon delta eta", "(#eta^{particle} - #eta^{photon})/#eta^{particle}", |
---|
57 | "(#eta^{particle} - #eta^{photon})/#eta^{particle}", "number of photons", |
---|
58 | 100, -0.1, 0.1); |
---|
59 | |
---|
60 | plots->fPhotonDeltaE = result->AddHist1D( |
---|
61 | "photon delta energy", "(E^{particle} - E^{photon})/E^{particle}", |
---|
62 | "(E^{particle} - E^{photon})/E^{particle}", "number of photons", |
---|
63 | 100, -0.1, 0.1); |
---|
64 | |
---|
65 | plots->fMuonDeltaPT = result->AddHist1D( |
---|
66 | "muon delta pt", "(p_{T}^{particle} - p_{T}^{muon})/p_{T}^{particle}", |
---|
67 | "(p_{T}^{particle} - p_{T}^{muon})/p_{T}^{particle}", "number of muons", |
---|
68 | 100, -0.1, 0.1); |
---|
69 | |
---|
70 | plots->fMuonDeltaEta = result->AddHist1D( |
---|
71 | "muon delta eta", "(#eta^{particle} - #eta^{muon})/#eta^{particle}", |
---|
72 | "(#eta^{particle} - #eta^{muon})/#eta^{particle}", "number of muons", |
---|
73 | 100, -0.1, 0.1); |
---|
74 | |
---|
75 | plots->fTrackDeltaPT = result->AddHist1D( |
---|
76 | "track delta pt", "(p_{T}^{particle} - p_{T}^{track})/p_{T}^{particle}", |
---|
77 | "(p_{T}^{particle} - p_{T}^{track})/p_{T}^{particle}", "number of tracks", |
---|
78 | 100, -0.1, 0.1); |
---|
79 | |
---|
80 | plots->fTrackDeltaEta = result->AddHist1D( |
---|
81 | "track delta eta", "(#eta^{particle} - #eta^{track})/#eta^{particle}", |
---|
82 | "(#eta^{particle} - #eta^{track})/#eta^{particle}", "number of tracks", |
---|
83 | 100, -0.1, 0.1); |
---|
84 | |
---|
85 | plots->fTowerDeltaEem = result->AddHist1D( |
---|
86 | "tower delta Eem", "(Eem^{particles} - Eem^{tower})/Eem^{particles}", |
---|
87 | "(Eem^{particles} - Eem^{tower})/Eem^{particles}", "number of tower", |
---|
88 | 100, -0.5, 0.5); |
---|
89 | |
---|
90 | plots->fTowerDeltaEhad = result->AddHist1D( |
---|
91 | "tower delta Ehad", "(Ehad^{particles} - Ehad^{tower})/Ehad^{particles}", |
---|
92 | "(Ehad^{particles} - Ehad^{tower})/Ehad^{particles}", "number of tower", |
---|
93 | 100, -5.0, 5.0); |
---|
94 | |
---|
95 | plots->fJetDeltaPT = result->AddHist1D( |
---|
96 | "jet delta pt", "(p_{T}^{jet} - p_{T}^{constituents})/p_{T}^{jet}", |
---|
97 | "(p_{T}^{jet} - p_{T}^{constituents})/p_{T}^{jet}", "number of jets", |
---|
98 | 100, -1.0e-7, 1.0e-7); |
---|
99 | |
---|
100 | } |
---|
101 | |
---|
102 | //------------------------------------------------------------------------------ |
---|
103 | |
---|
104 | void AnalyseEvents(ExRootTreeReader *treeReader, TestPlots *plots) |
---|
105 | { |
---|
106 | TClonesArray *branchParticle = treeReader->UseBranch("Particle"); |
---|
107 | TClonesArray *branchElectron = treeReader->UseBranch("Electron"); |
---|
108 | TClonesArray *branchPhoton = treeReader->UseBranch("Photon"); |
---|
109 | TClonesArray *branchMuon = treeReader->UseBranch("Muon"); |
---|
110 | |
---|
111 | TClonesArray *branchEFlowTrack = treeReader->UseBranch("EFlowTrack"); |
---|
112 | TClonesArray *branchEFlowTower = treeReader->UseBranch("EFlowTower"); |
---|
113 | TClonesArray *branchEFlowMuon = treeReader->UseBranch("EFlowMuon"); |
---|
114 | TClonesArray *branchJet = treeReader->UseBranch("Jet"); |
---|
115 | |
---|
116 | Long64_t allEntries = treeReader->GetEntries(); |
---|
117 | |
---|
118 | cout << "** Chain contains " << allEntries << " events" << endl; |
---|
119 | |
---|
120 | GenParticle *particle; |
---|
121 | Electron *electron; |
---|
122 | Photon *photon; |
---|
123 | Muon *muon; |
---|
124 | |
---|
125 | Track *track; |
---|
126 | Tower *tower; |
---|
127 | |
---|
128 | Jet *jet; |
---|
129 | TObject *object; |
---|
130 | |
---|
131 | TLorentzVector momentum; |
---|
132 | |
---|
133 | Float_t Eem, Ehad; |
---|
134 | Bool_t skip; |
---|
135 | |
---|
136 | Long64_t entry; |
---|
137 | |
---|
138 | Int_t i, j, pdgCode; |
---|
139 | |
---|
140 | // Loop over all events |
---|
141 | for(entry = 0; entry < allEntries; ++entry) |
---|
142 | { |
---|
143 | // Load selected branches with data from specified event |
---|
144 | treeReader->ReadEntry(entry); |
---|
145 | |
---|
146 | // Loop over all electrons in event |
---|
147 | for(i = 0; i < branchElectron->GetEntriesFast(); ++i) |
---|
148 | { |
---|
149 | electron = (Electron*) branchElectron->At(i); |
---|
150 | particle = (GenParticle*) electron->Particle.GetObject(); |
---|
151 | |
---|
152 | plots->fElectronDeltaPT->Fill((particle->PT - electron->PT)/particle->PT); |
---|
153 | plots->fElectronDeltaEta->Fill((particle->Eta - electron->Eta)/particle->Eta); |
---|
154 | } |
---|
155 | |
---|
156 | // Loop over all photons in event |
---|
157 | for(i = 0; i < branchPhoton->GetEntriesFast(); ++i) |
---|
158 | { |
---|
159 | photon = (Photon*) branchPhoton->At(i); |
---|
160 | |
---|
161 | // skip photons with references to multiple particles |
---|
162 | if(photon->Particles.GetEntriesFast() != 1) continue; |
---|
163 | |
---|
164 | particle = (GenParticle*) photon->Particles.At(0); |
---|
165 | |
---|
166 | plots->fPhotonDeltaPT->Fill((particle->PT - photon->PT)/particle->PT); |
---|
167 | plots->fPhotonDeltaEta->Fill((particle->Eta - photon->Eta)/particle->Eta); |
---|
168 | plots->fPhotonDeltaE->Fill((particle->E - photon->E)/particle->E); |
---|
169 | } |
---|
170 | |
---|
171 | // Loop over all muons in event |
---|
172 | for(i = 0; i < branchMuon->GetEntriesFast(); ++i) |
---|
173 | { |
---|
174 | muon = (Muon*) branchMuon->At(i); |
---|
175 | particle = (GenParticle*) muon->Particle.GetObject(); |
---|
176 | |
---|
177 | plots->fMuonDeltaPT->Fill((particle->PT - muon->PT)/particle->PT); |
---|
178 | plots->fMuonDeltaEta->Fill((particle->Eta - muon->Eta)/particle->Eta); |
---|
179 | } |
---|
180 | |
---|
181 | // Loop over all tracks in event |
---|
182 | for(i = 0; i < branchEFlowTrack->GetEntriesFast(); ++i) |
---|
183 | { |
---|
184 | track = (Track*) branchEFlowTrack->At(i); |
---|
185 | particle = (GenParticle*) track->Particle.GetObject(); |
---|
186 | |
---|
187 | plots->fTrackDeltaPT->Fill((particle->PT - track->PT)/particle->PT); |
---|
188 | plots->fTrackDeltaEta->Fill((particle->Eta - track->Eta)/particle->Eta); |
---|
189 | } |
---|
190 | |
---|
191 | // Loop over all towers in event |
---|
192 | for(i = 0; i < branchEFlowTower->GetEntriesFast(); ++i) |
---|
193 | { |
---|
194 | tower = (Tower*) branchEFlowTower->At(i); |
---|
195 | |
---|
196 | Eem = 0.0; |
---|
197 | Ehad = 0.0; |
---|
198 | skip = kFALSE; |
---|
199 | for(j = 0; j < tower->Particles.GetEntriesFast(); ++j) |
---|
200 | { |
---|
201 | particle = (GenParticle*) tower->Particles.At(j); |
---|
202 | pdgCode = TMath::Abs(particle->PID); |
---|
203 | |
---|
204 | // skip muons and neutrinos |
---|
205 | if(pdgCode == 12 || pdgCode == 13 || pdgCode == 14 || pdgCode == 16) |
---|
206 | { |
---|
207 | continue; |
---|
208 | } |
---|
209 | |
---|
210 | // skip K0short and Lambda |
---|
211 | if(pdgCode == 310 || pdgCode == 3122) |
---|
212 | { |
---|
213 | skip = kTRUE; |
---|
214 | } |
---|
215 | |
---|
216 | if(pdgCode == 11 || pdgCode == 22) |
---|
217 | { |
---|
218 | Eem += particle->E; |
---|
219 | } |
---|
220 | else |
---|
221 | { |
---|
222 | Ehad += particle->E; |
---|
223 | } |
---|
224 | } |
---|
225 | if(skip) continue; |
---|
226 | if(Eem > 0.0 && tower->Eem > 0.0) plots->fTowerDeltaEem->Fill((Eem - tower->Eem)/Eem); |
---|
227 | if(Ehad > 0.0 && tower->Ehad > 0.0) plots->fTowerDeltaEhad->Fill((Ehad - tower->Ehad)/Ehad); |
---|
228 | } |
---|
229 | |
---|
230 | // Loop over all jets in event |
---|
231 | for(i = 0; i < branchJet->GetEntriesFast(); ++i) |
---|
232 | { |
---|
233 | jet = (Jet*) branchJet->At(i); |
---|
234 | |
---|
235 | momentum.SetPxPyPzE(0.0, 0.0, 0.0, 0.0); |
---|
236 | |
---|
237 | // Loop over all jet's constituents |
---|
238 | for(j = 0; j < jet->Constituents.GetEntriesFast(); ++j) |
---|
239 | { |
---|
240 | object = jet->Constituents.At(j); |
---|
241 | |
---|
242 | // Check if the constituent is accessible |
---|
243 | if(object == 0) continue; |
---|
244 | |
---|
245 | if(object->IsA() == GenParticle::Class()) |
---|
246 | { |
---|
247 | momentum += ((GenParticle*) object)->P4(); |
---|
248 | } |
---|
249 | else if(object->IsA() == Track::Class()) |
---|
250 | { |
---|
251 | momentum += ((Track*) object)->P4(); |
---|
252 | } |
---|
253 | else if(object->IsA() == Tower::Class()) |
---|
254 | { |
---|
255 | momentum += ((Tower*) object)->P4(); |
---|
256 | } |
---|
257 | else if(object->IsA() == Muon::Class()) |
---|
258 | { |
---|
259 | momentum += ((Muon*) object)->P4(); |
---|
260 | } |
---|
261 | } |
---|
262 | plots->fJetDeltaPT->Fill((jet->PT - momentum.Pt())/jet->PT ); |
---|
263 | } |
---|
264 | } |
---|
265 | } |
---|
266 | |
---|
267 | //------------------------------------------------------------------------------ |
---|
268 | |
---|
269 | void PrintHistograms(ExRootResult *result, TestPlots *plots) |
---|
270 | { |
---|
271 | result->Print("png"); |
---|
272 | } |
---|
273 | |
---|
274 | //------------------------------------------------------------------------------ |
---|
275 | |
---|
276 | void Example3(const char *inputFile) |
---|
277 | { |
---|
278 | gSystem->Load("libDelphes"); |
---|
279 | |
---|
280 | TChain *chain = new TChain("Delphes"); |
---|
281 | chain->Add(inputFile); |
---|
282 | |
---|
283 | ExRootTreeReader *treeReader = new ExRootTreeReader(chain); |
---|
284 | ExRootResult *result = new ExRootResult(); |
---|
285 | |
---|
286 | TestPlots *plots = new TestPlots; |
---|
287 | |
---|
288 | BookHistograms(result, plots); |
---|
289 | |
---|
290 | AnalyseEvents(treeReader, plots); |
---|
291 | |
---|
292 | PrintHistograms(result, plots); |
---|
293 | |
---|
294 | result->Write("results.root"); |
---|
295 | |
---|
296 | cout << "** Exiting..." << endl; |
---|
297 | |
---|
298 | delete plots; |
---|
299 | delete result; |
---|
300 | delete treeReader; |
---|
301 | delete chain; |
---|
302 | } |
---|
303 | |
---|
304 | //------------------------------------------------------------------------------ |
---|