Fork me on GitHub

source: git/examples/Example4.C@ e7e90df

ImprovedOutputFile Timing dual_readout llp
Last change on this file since e7e90df was 35b9204, checked in by Michele <michele.selvaggi@…>, 10 years ago

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