Fork me on GitHub

source: git/examples/Example4.C@ 88c03bc

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 88c03bc was c9803f7, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

adapt examples to ROOT 6.04
  • Property mode set to 100644
File size: 9.2 KB
Line 
1/*
2This macro shows how to compute jet energy scale.
3root -l examples/Example4.C'("delphes_output.root", "plots.root")'
4
5The output ROOT file contains the pT(MC)/pT(Reco) distributions for
6various pT(Reco) and |eta| bins. The peak value of such distribution is
7interpreted as the jet energy correction to be applied for that
8given pT(Reco), |eta| bin.
9
10This can be done by modifying the "ScaleFormula" input parameter to
11the JetEnergyScale module in the delphes_card_XXX.tcl
12
13e.g a smooth function:
14
15 set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0) }
16
17or a binned function:
18
19 set ScaleFormula {(abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 20.0 && pt <= 50.0) * (1.10) +
20 (abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 50.0 && pt <= 100.0) * (1.05) +
21 (abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 100.0) * (1.00) +
22 (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 20.0 && pt <= 50.0) * (1.10) +
23 (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 50.0 && pt <= 100.0) * (1.05) +
24 (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 100.0) * (1.00)}
25
26Be aware that a binned jet energy scale can produce "steps" in the corrected
27jet pt distribution ...
28*/
29
30#ifdef __CLING__
31R__LOAD_LIBRARY(libDelphes)
32#include "classes/DelphesClasses.h"
33#include "external/ExRootAnalysis/ExRootTreeReader.h"
34#include "external/ExRootAnalysis/ExRootResult.h"
35#else
36class ExRootTreeReader;
37class ExRootResult;
38#endif
39
40//------------------------------------------------------------------------------
41
42struct TestPlots
43{
44 TH1 *fJetPT;
45
46 TH1 *fJetRes_Pt_20_50_Eta_0_25;
47 TH1 *fJetRes_Pt_20_50_Eta_25_5;
48
49 TH1 *fJetRes_Pt_50_100_Eta_0_25;
50 TH1 *fJetRes_Pt_50_100_Eta_25_5;
51
52 TH1 *fJetRes_Pt_100_200_Eta_0_25;
53 TH1 *fJetRes_Pt_100_200_Eta_25_5;
54
55 TH1 *fJetRes_Pt_200_500_Eta_0_25;
56 TH1 *fJetRes_Pt_200_500_Eta_25_5;
57
58 TH1 *fJetRes_Pt_500_inf_Eta_0_25;
59 TH1 *fJetRes_Pt_500_inf_Eta_25_5;
60
61};
62
63//------------------------------------------------------------------------------
64
65class ExRootResult;
66class ExRootTreeReader;
67
68//------------------------------------------------------------------------------
69
70void BookHistograms(ExRootResult *result, TestPlots *plots)
71{
72 TLegend *legend;
73 TPaveText *comment;
74
75 plots->fJetPT = result->AddHist1D(
76 "jet_pt", "p_{T}^{jet}",
77 "p_{T}^{jet} GeV/c", "number of jets",
78 100, 0.0, 1000.0);
79
80 plots->fJetRes_Pt_20_50_Eta_0_25 = result->AddHist1D(
81 "jet_delta_pt_20_50_cen", "p_{T}^{truth,parton}/p_{T}^{jet} , 20 < p_{T} < 50 , 0 < | #eta | < 2.5 ",
82 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
83 100, 0.0, 2.0);
84
85 plots->fJetRes_Pt_20_50_Eta_0_25->SetStats();
86
87 plots->fJetRes_Pt_20_50_Eta_25_5 = result->AddHist1D(
88 "jet_delta_pt_20_50_fwd", "p_{T}^{truth,parton}/p_{T}^{jet} , 20 < p_{T} < 50 , 2.5 < | #eta | < 5 ",
89 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
90 100, 0.0, 2.0);
91
92 plots->fJetRes_Pt_20_50_Eta_25_5->SetStats();
93
94 plots->fJetRes_Pt_50_100_Eta_0_25 = result->AddHist1D(
95 "jet_delta_pt_50_100_cen", "p_{T}^{truth,parton}/p_{T}^{jet} , 50 < p_{T} < 100 , 0 < | #eta | < 2.5 ",
96 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
97 100, 0.0, 2.0);
98
99 plots->fJetRes_Pt_50_100_Eta_0_25->SetStats();
100
101 plots->fJetRes_Pt_50_100_Eta_25_5 = result->AddHist1D(
102 "jet_delta_pt_50_100_fwd", "p_{T}^{truth,parton}/p_{T}^{jet} , 50 < p_{T} < 100 , 2.5 < | #eta | < 5 ",
103 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
104 100, 0.0, 2.0);
105
106 plots->fJetRes_Pt_50_100_Eta_25_5->SetStats();
107
108
109 plots->fJetRes_Pt_100_200_Eta_0_25 = result->AddHist1D(
110 "jet_delta_pt_100_200_cen", "p_{T}^{truth,parton}/p_{T}^{jet} , 100 < p_{T} < 200 , 0 < | #eta | < 2.5 ",
111 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
112 100, 0.0, 2.0);
113
114 plots->fJetRes_Pt_100_200_Eta_0_25->SetStats();
115
116 plots->fJetRes_Pt_100_200_Eta_25_5 = result->AddHist1D(
117 "jet_delta_pt_100_200_fwd", "p_{T}^{truth,parton}/p_{T}^{jet} , 100 < p_{T} < 200 , 2.5 < | #eta | < 5 ",
118 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
119 100, 0.0, 2.0);
120
121 plots->fJetRes_Pt_100_200_Eta_25_5->SetStats();
122
123 plots->fJetRes_Pt_200_500_Eta_0_25 = result->AddHist1D(
124 "jet_delta_pt_200_500_cen", "p_{T}^{truth,parton}/p_{T}^{jet} , 200 < p_{T} < 500 , 0 < | #eta | < 2.5 ",
125 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
126 100, 0.0, 2.0);
127
128 plots->fJetRes_Pt_200_500_Eta_0_25->SetStats();
129
130 plots->fJetRes_Pt_200_500_Eta_25_5 = result->AddHist1D(
131 "jet_delta_pt_200_500_fwd", "p_{T}^{truth,parton}/p_{T}^{jet} , 200 < p_{T} < 500 , 2.5 < | #eta | < 5 ",
132 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
133 100, 0.0, 2.0);
134
135 plots->fJetRes_Pt_200_500_Eta_25_5->SetStats();
136
137 plots->fJetRes_Pt_500_inf_Eta_0_25 = result->AddHist1D(
138 "jet_delta_pt_500_1000_cen", "p_{T}^{truth,parton}/p_{T}^{jet} , 500 < p_{T} < 1000, 0 < | #eta | < 2.5 ",
139 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
140 100, 0.0, 2.0);
141
142 plots->fJetRes_Pt_500_inf_Eta_0_25->SetStats();
143
144 plots->fJetRes_Pt_500_inf_Eta_25_5 = result->AddHist1D(
145 "jet_delta_pt_500_1000_fwd", "p_{T}^{truth,parton}/p_{T}^{jet} , 500 < p_{T} < 1000, 2.5 < | #eta | < 5 ",
146 "p_{T}^{truth,parton}/p_{T}^{jet}", "number of jets",
147 100, 0.0, 2.0);
148
149 plots->fJetRes_Pt_500_inf_Eta_25_5->SetStats();
150
151
152}
153
154//------------------------------------------------------------------------------
155
156void AnalyseEvents(ExRootTreeReader *treeReader, TestPlots *plots)
157{
158 TClonesArray *branchParticle = treeReader->UseBranch("Particle");
159 TClonesArray *branchGenJet = treeReader->UseBranch("GenJet");
160 TClonesArray *branchJet = treeReader->UseBranch("Jet");
161
162 Long64_t allEntries = treeReader->GetEntries();
163
164 cout << "** Chain contains " << allEntries << " events" << endl;
165
166 Jet *jet, *genjet;
167 GenParticle *particle;
168 TObject *object;
169
170 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
171
172 Float_t deltaR;
173 Float_t pt, eta;
174 Long64_t entry;
175
176 Int_t i, j;
177
178 // Loop over all events
179 for(entry = 0; entry < allEntries; ++entry)
180 {
181 // Load selected branches with data from specified event
182 treeReader->ReadEntry(entry);
183
184 if(entry%500 == 0) cout << "Event number: "<< entry <<endl;
185
186 // Loop over all reconstructed jets in event
187 for(i = 0; i < branchJet->GetEntriesFast(); ++i)
188 {
189
190 jet = (Jet*) branchJet->At(i);
191 jetMomentum = jet->P4();
192
193 plots->fJetPT->Fill(jetMomentum.Pt());
194
195 deltaR = 999;
196
197 // Loop over all hard partons in event
198 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
199 {
200 particle = (GenParticle*) branchParticle->At(j);
201
202 genJetMomentum = particle->P4();
203
204 // this is simply to avoid warnings from initial state particle
205 // having infite rapidity ...
206 if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
207
208 // take the closest parton candidate
209 if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
210 {
211 deltaR = genJetMomentum.DeltaR(jetMomentum);
212 bestGenJetMomentum = genJetMomentum;
213 }
214 }
215
216 if(deltaR < 0.3)
217 {
218 pt = jetMomentum.Pt();
219 eta = TMath::Abs(jetMomentum.Eta());
220
221
222 if(pt > 20.0 && pt < 50.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_20_50_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
223 if(pt > 20.0 && pt < 50.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_20_50_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
224
225 if(pt > 50.0 && pt < 100.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_50_100_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
226 if(pt > 50.0 && pt < 100.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_50_100_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
227
228 if(pt > 100.0 && pt < 200.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_100_200_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
229 if(pt > 100.0 && pt < 200.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_100_200_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
230
231 if(pt > 200.0 && pt < 500.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_200_500_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
232 if(pt > 200.0 && pt < 500.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_200_500_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
233
234 if(pt > 500.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_500_inf_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
235 if(pt > 500.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_500_inf_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
236
237 }
238 }
239 }
240}
241
242//------------------------------------------------------------------------------
243
244void Example4(const char *inputFile, const char *outputFile)
245{
246 gSystem->Load("libDelphes");
247
248 TChain *chain = new TChain("Delphes");
249 chain->Add(inputFile);
250
251 ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
252 ExRootResult *result = new ExRootResult();
253
254 TestPlots *plots = new TestPlots;
255
256 BookHistograms(result, plots);
257
258 AnalyseEvents(treeReader, plots);
259
260 result->Write(outputFile);
261
262 cout << "** Exiting..." << endl;
263
264 delete plots;
265 delete result;
266 delete treeReader;
267 delete chain;
268}
269
270//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.