Fork me on GitHub

source: git/examples/Example3.C@ 3ccc8586

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 3ccc8586 was 374bbcc, checked in by pavel <pavel@…>, 11 years ago

remove EFlowMuon from examples

  • Property mode set to 100644
File size: 8.9 KB
Line 
1/*
2root -l examples/Example3.C\(\"delphes_output.root\"\)
3*/
4
5//------------------------------------------------------------------------------
6
7struct 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
30class ExRootResult;
31class ExRootTreeReader;
32
33//------------------------------------------------------------------------------
34
35void 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
104void 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 *branchJet = treeReader->UseBranch("Jet");
114
115 Long64_t allEntries = treeReader->GetEntries();
116
117 cout << "** Chain contains " << allEntries << " events" << endl;
118
119 GenParticle *particle;
120 Electron *electron;
121 Photon *photon;
122 Muon *muon;
123
124 Track *track;
125 Tower *tower;
126
127 Jet *jet;
128 TObject *object;
129
130 TLorentzVector momentum;
131
132 Float_t Eem, Ehad;
133 Bool_t skip;
134
135 Long64_t entry;
136
137 Int_t i, j, pdgCode;
138
139 // Loop over all events
140 for(entry = 0; entry < allEntries; ++entry)
141 {
142 // Load selected branches with data from specified event
143 treeReader->ReadEntry(entry);
144
145 // Loop over all electrons in event
146 for(i = 0; i < branchElectron->GetEntriesFast(); ++i)
147 {
148 electron = (Electron*) branchElectron->At(i);
149 particle = (GenParticle*) electron->Particle.GetObject();
150
151 plots->fElectronDeltaPT->Fill((particle->PT - electron->PT)/particle->PT);
152 plots->fElectronDeltaEta->Fill((particle->Eta - electron->Eta)/particle->Eta);
153 }
154
155 // Loop over all photons in event
156 for(i = 0; i < branchPhoton->GetEntriesFast(); ++i)
157 {
158 photon = (Photon*) branchPhoton->At(i);
159
160 // skip photons with references to multiple particles
161 if(photon->Particles.GetEntriesFast() != 1) continue;
162
163 particle = (GenParticle*) photon->Particles.At(0);
164
165 plots->fPhotonDeltaPT->Fill((particle->PT - photon->PT)/particle->PT);
166 plots->fPhotonDeltaEta->Fill((particle->Eta - photon->Eta)/particle->Eta);
167 plots->fPhotonDeltaE->Fill((particle->E - photon->E)/particle->E);
168 }
169
170 // Loop over all muons in event
171 for(i = 0; i < branchMuon->GetEntriesFast(); ++i)
172 {
173 muon = (Muon*) branchMuon->At(i);
174 particle = (GenParticle*) muon->Particle.GetObject();
175
176 plots->fMuonDeltaPT->Fill((particle->PT - muon->PT)/particle->PT);
177 plots->fMuonDeltaEta->Fill((particle->Eta - muon->Eta)/particle->Eta);
178 }
179
180 // Loop over all tracks in event
181 for(i = 0; i < branchEFlowTrack->GetEntriesFast(); ++i)
182 {
183 track = (Track*) branchEFlowTrack->At(i);
184 particle = (GenParticle*) track->Particle.GetObject();
185
186 plots->fTrackDeltaPT->Fill((particle->PT - track->PT)/particle->PT);
187 plots->fTrackDeltaEta->Fill((particle->Eta - track->Eta)/particle->Eta);
188 }
189
190 // Loop over all towers in event
191 for(i = 0; i < branchEFlowTower->GetEntriesFast(); ++i)
192 {
193 tower = (Tower*) branchEFlowTower->At(i);
194
195 Eem = 0.0;
196 Ehad = 0.0;
197 skip = kFALSE;
198 for(j = 0; j < tower->Particles.GetEntriesFast(); ++j)
199 {
200 particle = (GenParticle*) tower->Particles.At(j);
201 pdgCode = TMath::Abs(particle->PID);
202
203 // skip muons and neutrinos
204 if(pdgCode == 12 || pdgCode == 13 || pdgCode == 14 || pdgCode == 16)
205 {
206 continue;
207 }
208
209 // skip K0short and Lambda
210 if(pdgCode == 310 || pdgCode == 3122)
211 {
212 skip = kTRUE;
213 }
214
215 if(pdgCode == 11 || pdgCode == 22)
216 {
217 Eem += particle->E;
218 }
219 else
220 {
221 Ehad += particle->E;
222 }
223 }
224 if(skip) continue;
225 if(Eem > 0.0 && tower->Eem > 0.0) plots->fTowerDeltaEem->Fill((Eem - tower->Eem)/Eem);
226 if(Ehad > 0.0 && tower->Ehad > 0.0) plots->fTowerDeltaEhad->Fill((Ehad - tower->Ehad)/Ehad);
227 }
228
229 // Loop over all jets in event
230 for(i = 0; i < branchJet->GetEntriesFast(); ++i)
231 {
232 jet = (Jet*) branchJet->At(i);
233
234 momentum.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
235
236 // Loop over all jet's constituents
237 for(j = 0; j < jet->Constituents.GetEntriesFast(); ++j)
238 {
239 object = jet->Constituents.At(j);
240
241 // Check if the constituent is accessible
242 if(object == 0) continue;
243
244 if(object->IsA() == GenParticle::Class())
245 {
246 momentum += ((GenParticle*) object)->P4();
247 }
248 else if(object->IsA() == Track::Class())
249 {
250 momentum += ((Track*) object)->P4();
251 }
252 else if(object->IsA() == Tower::Class())
253 {
254 momentum += ((Tower*) object)->P4();
255 }
256 else if(object->IsA() == Muon::Class())
257 {
258 momentum += ((Muon*) object)->P4();
259 }
260 }
261 plots->fJetDeltaPT->Fill((jet->PT - momentum.Pt())/jet->PT );
262 }
263 }
264}
265
266//------------------------------------------------------------------------------
267
268void PrintHistograms(ExRootResult *result, TestPlots *plots)
269{
270 result->Print("png");
271}
272
273//------------------------------------------------------------------------------
274
275void Example3(const char *inputFile)
276{
277 gSystem->Load("libDelphes");
278
279 TChain *chain = new TChain("Delphes");
280 chain->Add(inputFile);
281
282 ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
283 ExRootResult *result = new ExRootResult();
284
285 TestPlots *plots = new TestPlots;
286
287 BookHistograms(result, plots);
288
289 AnalyseEvents(treeReader, plots);
290
291 PrintHistograms(result, plots);
292
293 result->Write("results.root");
294
295 cout << "** Exiting..." << endl;
296
297 delete plots;
298 delete result;
299 delete treeReader;
300 delete chain;
301}
302
303//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.