Fork me on GitHub

source: git/examples/Example3.C@ 282f591

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 282f591 was 27bf162, checked in by mselvaggi <mselvaggi@…>, 11 years ago

separated Photon and Neutral Hadron components in particle algorithm

  • Property mode set to 100644
File size: 8.2 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 *fJetDeltaPT;
23};
24
25//------------------------------------------------------------------------------
26
27class ExRootResult;
28class ExRootTreeReader;
29
30//------------------------------------------------------------------------------
31
32void BookHistograms(ExRootResult *result, TestPlots *plots)
33{
34 TLegend *legend;
35 TPaveText *comment;
36
37 plots->fElectronDeltaPT = result->AddHist1D(
38 "electron delta pt", "(p_{T}^{particle} - p_{T}^{electron})/p_{T}^{particle}",
39 "(p_{T}^{particle} - p_{T}^{electron})/p_{T}^{particle}", "number of electrons",
40 100, -0.1, 0.1);
41
42 plots->fElectronDeltaEta = result->AddHist1D(
43 "electron delta eta", "(#eta^{particle} - #eta^{electron})/#eta^{particle}",
44 "(#eta^{particle} - #eta^{electron})/#eta^{particle}", "number of electrons",
45 100, -0.1, 0.1);
46
47 plots->fPhotonDeltaPT = result->AddHist1D(
48 "photon delta pt", "(p_{T}^{particle} - p_{T}^{photon})/p_{T}^{particle}",
49 "(p_{T}^{particle} - p_{T}^{photon})/p_{T}^{particle}", "number of photons",
50 100, -0.1, 0.1);
51
52 plots->fPhotonDeltaEta = result->AddHist1D(
53 "photon delta eta", "(#eta^{particle} - #eta^{photon})/#eta^{particle}",
54 "(#eta^{particle} - #eta^{photon})/#eta^{particle}", "number of photons",
55 100, -0.1, 0.1);
56
57 plots->fPhotonDeltaE = result->AddHist1D(
58 "photon delta energy", "(E^{particle} - E^{photon})/E^{particle}",
59 "(E^{particle} - E^{photon})/E^{particle}", "number of photons",
60 100, -0.1, 0.1);
61
62 plots->fMuonDeltaPT = result->AddHist1D(
63 "muon delta pt", "(p_{T}^{particle} - p_{T}^{muon})/p_{T}^{particle}",
64 "(p_{T}^{particle} - p_{T}^{muon})/p_{T}^{particle}", "number of muons",
65 100, -0.1, 0.1);
66
67 plots->fMuonDeltaEta = result->AddHist1D(
68 "muon delta eta", "(#eta^{particle} - #eta^{muon})/#eta^{particle}",
69 "(#eta^{particle} - #eta^{muon})/#eta^{particle}", "number of muons",
70 100, -0.1, 0.1);
71
72 plots->fTrackDeltaPT = result->AddHist1D(
73 "track delta pt", "(p_{T}^{particle} - p_{T}^{track})/p_{T}^{particle}",
74 "(p_{T}^{particle} - p_{T}^{track})/p_{T}^{particle}", "number of tracks",
75 100, -0.1, 0.1);
76
77 plots->fTrackDeltaEta = result->AddHist1D(
78 "track delta eta", "(#eta^{particle} - #eta^{track})/#eta^{particle}",
79 "(#eta^{particle} - #eta^{track})/#eta^{particle}", "number of tracks",
80 100, -0.1, 0.1);
81
82 plots->fJetDeltaPT = result->AddHist1D(
83 "jet delta pt", "(p_{T}^{jet} - p_{T}^{constituents})/p_{T}^{jet}",
84 "(p_{T}^{jet} - p_{T}^{constituents})/p_{T}^{jet}", "number of jets",
85 100, -1.0e-1, 1.0e-1);
86
87}
88
89//------------------------------------------------------------------------------
90
91void AnalyseEvents(ExRootTreeReader *treeReader, TestPlots *plots)
92{
93 TClonesArray *branchParticle = treeReader->UseBranch("Particle");
94 TClonesArray *branchElectron = treeReader->UseBranch("Electron");
95 TClonesArray *branchPhoton = treeReader->UseBranch("Photon");
96 TClonesArray *branchMuon = treeReader->UseBranch("Muon");
97
98 TClonesArray *branchEFlowTrack = treeReader->UseBranch("EFlowTrack");
99 TClonesArray *branchEFlowPhoton = treeReader->UseBranch("EFlowPhoton");
100 TClonesArray *branchEFlowNeutralHadron = treeReader->UseBranch("EFlowNeutralHadron");
101 TClonesArray *branchJet = treeReader->UseBranch("Jet");
102
103 Long64_t allEntries = treeReader->GetEntries();
104
105 cout << "** Chain contains " << allEntries << " events" << endl;
106
107 GenParticle *particle;
108 Electron *electron;
109 Photon *photon;
110 Muon *muon;
111
112 Track *track;
113 Tower *tower;
114
115 Jet *jet;
116 TObject *object;
117
118 TLorentzVector momentum;
119
120 Float_t Eem, Ehad;
121 Bool_t skip;
122
123 Long64_t entry;
124
125 Int_t i, j, pdgCode;
126
127 // Loop over all events
128 for(entry = 0; entry < allEntries; ++entry)
129 {
130 // Load selected branches with data from specified event
131 treeReader->ReadEntry(entry);
132
133 // Loop over all electrons in event
134 for(i = 0; i < branchElectron->GetEntriesFast(); ++i)
135 {
136 electron = (Electron*) branchElectron->At(i);
137 particle = (GenParticle*) electron->Particle.GetObject();
138
139 plots->fElectronDeltaPT->Fill((particle->PT - electron->PT)/particle->PT);
140 plots->fElectronDeltaEta->Fill((particle->Eta - electron->Eta)/particle->Eta);
141 }
142
143 // Loop over all photons in event
144 for(i = 0; i < branchPhoton->GetEntriesFast(); ++i)
145 {
146 photon = (Photon*) branchPhoton->At(i);
147
148 // skip photons with references to multiple particles
149 if(photon->Particles.GetEntriesFast() != 1) continue;
150
151 particle = (GenParticle*) photon->Particles.At(0);
152
153 plots->fPhotonDeltaPT->Fill((particle->PT - photon->PT)/particle->PT);
154 plots->fPhotonDeltaEta->Fill((particle->Eta - photon->Eta)/particle->Eta);
155 plots->fPhotonDeltaE->Fill((particle->E - photon->E)/particle->E);
156 }
157
158 // Loop over all muons in event
159 for(i = 0; i < branchMuon->GetEntriesFast(); ++i)
160 {
161 muon = (Muon*) branchMuon->At(i);
162 particle = (GenParticle*) muon->Particle.GetObject();
163
164 plots->fMuonDeltaPT->Fill((particle->PT - muon->PT)/particle->PT);
165 plots->fMuonDeltaEta->Fill((particle->Eta - muon->Eta)/particle->Eta);
166 }
167
168 // Loop over all tracks in event
169 for(i = 0; i < branchEFlowTrack->GetEntriesFast(); ++i)
170 {
171 track = (Track*) branchEFlowTrack->At(i);
172 particle = (GenParticle*) track->Particle.GetObject();
173
174 plots->fTrackDeltaPT->Fill((particle->PT - track->PT)/particle->PT);
175 plots->fTrackDeltaEta->Fill((particle->Eta - track->Eta)/particle->Eta);
176 }
177
178 cout<<"-- New event -- "<<endl;
179
180 // Loop over all jets in event
181 for(i = 0; i < branchJet->GetEntriesFast(); ++i)
182 {
183 jet = (Jet*) branchJet->At(i);
184
185 momentum.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
186
187 cout<<"Looping over jet constituents. Jet pt: "<<jet->PT<<", eta: "<<jet->Eta<<", phi: "<<jet->Phi<<endl;
188
189 // Loop over all jet's constituents
190 for(j = 0; j < jet->Constituents.GetEntriesFast(); ++j)
191 {
192 object = jet->Constituents.At(j);
193
194 // Check if the constituent is accessible
195 if(object == 0) continue;
196
197 if(object->IsA() == GenParticle::Class())
198 {
199 cout<<" GenPart pt: "<<((GenParticle*) object)->PT<<", eta: "<<((GenParticle*) object)->Eta<<", phi: "<<((GenParticle*) object)->Phi<<endl;
200 momentum += ((GenParticle*) object)->P4();
201 }
202 else if(object->IsA() == Track::Class())
203 {
204 cout<<" Track pt: "<<((Track*) object)->PT<<", eta: "<<((Track*) object)->Eta<<", phi: "<<((Track*) object)->Phi<<endl;
205 momentum += ((Track*) object)->P4();
206 }
207 else if(object->IsA() == Tower::Class())
208 {
209 cout<<" Tower pt: "<<((Tower*) object)->ET<<", eta: "<<((Tower*) object)->Eta<<", phi: "<<((Tower*) object)->Phi<<endl;
210 momentum += ((Tower*) object)->P4();
211 }
212 else if(object->IsA() == Muon::Class())
213 {
214 cout<<" Muon pt: "<<((Muon*) object)->PT<<", eta: "<<((Muon*) object)->Eta<<", phi: "<<((Muon*) object)->Phi<<endl;
215 momentum += ((Muon*) object)->P4();
216 }
217 }
218 plots->fJetDeltaPT->Fill((jet->PT - momentum.Pt())/jet->PT);
219 }
220 }
221}
222
223//------------------------------------------------------------------------------
224
225void PrintHistograms(ExRootResult *result, TestPlots *plots)
226{
227 result->Print("png");
228}
229
230//------------------------------------------------------------------------------
231
232void Example3(const char *inputFile)
233{
234 gSystem->Load("libDelphes");
235
236 TChain *chain = new TChain("Delphes");
237 chain->Add(inputFile);
238
239 ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
240 ExRootResult *result = new ExRootResult();
241
242 TestPlots *plots = new TestPlots;
243
244 BookHistograms(result, plots);
245
246 AnalyseEvents(treeReader, plots);
247
248 PrintHistograms(result, plots);
249
250 result->Write("results.root");
251
252 cout << "** Exiting..." << endl;
253
254 delete plots;
255 delete result;
256 delete treeReader;
257 delete chain;
258}
259
260//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.