CMS 3D CMS Logo
JetResolution_t.cc
Go to the documentation of this file.
1 //
3 // JetResolution_t
4 // ---------------
5 //
6 // 11/05/2010 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
8 
9 #include "CondFormats/JetMETObjects/interface/JetResolution.h"
10 #include "FWCore/ParameterSet/interface/FileInPath.h"
11 
12 #include <TROOT.h>
13 #include <TApplication.h>
14 #include <TRandom3.h>
15 #include <TCanvas.h>
16 #include <TLegend.h>
17 #include <TH1F.h>
18 #include <TMath.h>
19 
20 #include <cstdlib>
21 #include <iostream>
22 #include <sstream>
23 #include <string>
24 #include <sys/stat.h>
25 
26 using namespace std;
27 
28 //______________________________________________________________________________
29 int main(int argc, char** argv) {
30  if (argc > 1 && string(argv[1]) == "--help") {
31  cout << "USAGE: JetResolution_t --era <era> --alg <alg> --nevts <n> --gaussian" << endl;
32  return 0;
33  }
34 
35  string era("Spring10");
36  string alg("AK5Calo");
37  unsigned nevts(10000);
38  bool doGaussian(false);
39 
40  for (int i = 1; i < argc; i++) {
41  if (string(argv[i]) == "--era") {
42  era = argv[i + 1];
43  i++;
44  } else if (string(argv[i]) == "--alg") {
45  alg = argv[i + 1];
46  i++;
47  } else if (string(argv[i]) == "--nevts") {
48  stringstream ss;
49  ss << argv[i + 1];
50  ss >> nevts;
51  i++;
52  } else if (string(argv[i]) == "--gaussian") {
53  doGaussian = true;
54  } else {
55  cerr << "ERROR: unknown option '" << argv[i] << "'" << endl;
56  return 0;
57  }
58  }
59 
60  cout << "era: " << era << endl;
61  cout << "alg: " << alg << endl;
62  cout << "nevts: " << nevts << endl;
63  cout << "gaussian: " << doGaussian << endl << endl;
64 
65  string ptFileName =
66  edm::FileInPath("CondFormats/JetMETObjects/data/" + era + "_PtResolution_" + alg + ".txt").fullPath();
67  string etaFileName =
68  edm::FileInPath("CondFormats/JetMETObjects/data/" + era + "_EtaResolution_" + alg + ".txt").fullPath();
69  string phiFileName =
70  edm::FileInPath("CondFormats/JetMETObjects/data/" + era + "_PhiResolution_" + alg + ".txt").fullPath();
71 
72  cout << ptFileName << endl;
73  cout << etaFileName << endl;
74  cout << phiFileName << endl;
75  cout << endl;
76 
77  JetResolution ptResol(ptFileName, doGaussian);
78  JetResolution etaResol(etaFileName, doGaussian);
79  JetResolution phiResol(phiFileName, doGaussian);
80 
81  // SIMPLE TEST
82  float pt = 200. * gRandom->Exp(0.1);
83  float eta = gRandom->Uniform(-5.0, 5.0);
84  float phi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
85 
86  cout << "pT=" << pt << " eta=" << eta << " phi=" << phi << endl;
87 
88  TF1* fPtResol = ptResol.resolutionEtaPt(eta, pt);
89  TF1* fEtaResol = etaResol.resolutionEtaPt(eta, pt);
90  TF1* fPhiResol = phiResol.resolutionEtaPt(eta, pt);
91 
92  fPtResol->Print();
93  cout << endl;
94  fEtaResol->Print();
95  cout << endl;
96  fPhiResol->Print();
97  cout << endl;
98 
99  // GENERATE SPECTRA AND SMEAR
100  TH1F* hRndPt = new TH1F("RndPt", ";random number", 200, 0.0, 5.0);
101  TH1F* hGenPt = new TH1F("GenPt", ";p_{T} [GeV]", 100, 0., 250.);
102  TH1F* hJetPt = new TH1F("JetPt", ";p_{T} [GeV]", 100, 0., 250.);
103 
104  TH1F* hRndEta = new TH1F("RndEta", ";random number", 200, -5.0, 5.0);
105  TH1F* hGenEta = new TH1F("GenEta", ";#eta", 51, -5., 5.);
106  TH1F* hJetEta = new TH1F("JetEta", ";#eta", 51, -5., 5.);
107 
108  TH1F* hRndPhi = new TH1F("RndPhi", ";random number", 200, -3.15, 3.15);
109  TH1F* hGenPhi = new TH1F("GenPhi", ";#varphi", 41, -3.15, 3.15);
110  TH1F* hJetPhi = new TH1F("JetPhi", ";#varphi", 41, -3.15, 3.15);
111 
112  for (unsigned ievt = 0; ievt < nevts; ievt++) {
113  if ((ievt + 1) % 1000 == 0)
114  cout << ievt + 1 << " events processed." << endl;
115 
116  float genpt = 200. * gRandom->Exp(0.1);
117  if (genpt < 1.0)
118  continue;
119  float geneta = gRandom->Uniform(-5.0, 5.0);
120  float genphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
121 
122  fPtResol = ptResol.resolutionEtaPt(geneta, genpt);
123  fEtaResol = etaResol.resolutionEtaPt(geneta, genpt);
124  fPhiResol = phiResol.resolutionEtaPt(geneta, genpt);
125 
126  float rndpt = fPtResol->GetRandom();
127  float rndeta = fEtaResol->GetRandom();
128  float rndphi = fPhiResol->GetRandom();
129 
130  float jetpt = rndpt * genpt;
131  float jeteta = rndeta + geneta;
132  float jetphi = rndphi + genphi;
133 
134  hRndPt->Fill(rndpt);
135  hGenPt->Fill(genpt);
136  hJetPt->Fill(jetpt);
137 
138  hRndEta->Fill(rndeta);
139  hGenEta->Fill(geneta);
140  hJetEta->Fill(jeteta);
141 
142  hRndPhi->Fill(rndphi);
143  hGenPhi->Fill(genphi);
144  hJetPhi->Fill(jetphi);
145  }
146 
147  // RUN ROOT APPLICATION AND DRAW BOTH DISTRIBUTIONS
148  argc = 1;
149  TApplication* app = new TApplication("JetResolution_t", &argc, argv);
150  gROOT->SetStyle("Plain");
151 
152  // PLOT RESOLUTION FOR DIFFERENT ETA BINS
153  TCanvas* cResolution = new TCanvas("Resolution", "Resolution", 0, 0, 1000, 400);
154  cResolution->Divide(3, 1);
155 
156  cResolution->cd(1);
157  gPad->SetLogx();
158 
159  TF1* fPtEta1 = ptResol.parameterEta("sigma", 0.25);
160  TF1* fPtEta2 = ptResol.parameterEta("sigma", 1.75);
161  TF1* fPtEta3 = ptResol.parameterEta("sigma", 2.75);
162 
163  fPtEta1->SetLineWidth(1);
164  fPtEta2->SetLineWidth(1);
165  fPtEta3->SetLineWidth(1);
166  fPtEta1->SetNpx(500);
167  fPtEta2->SetNpx(500);
168  fPtEta3->SetNpx(500);
169  fPtEta1->SetLineColor(kRed);
170  fPtEta2->SetLineColor(kBlue);
171  fPtEta3->SetLineColor(kGreen);
172  fPtEta1->SetRange(5.0, 500.);
173  fPtEta2->SetRange(5.0, 500.);
174  fPtEta3->SetRange(5.0, 500.);
175  fPtEta1->Draw();
176  fPtEta2->Draw("SAME");
177  fPtEta3->Draw("SAME");
178 
179  fPtEta1->SetTitle("p_{T} - Resolution");
180  fPtEta1->GetHistogram()->SetXTitle("p_{T} [GeV]");
181  fPtEta1->GetHistogram()->SetYTitle("#sigma(p_{T})/p_{T}");
182  fPtEta1->GetHistogram()->GetYaxis()->CenterTitle();
183  fPtEta1->GetHistogram()->GetYaxis()->SetTitleOffset(1.33);
184 
185  TLegend* legPtRes = new TLegend(0.6, 0.85, 0.85, 0.6);
186  legPtRes->SetFillColor(10);
187  legPtRes->AddEntry(fPtEta1, "#eta=0.25", "l");
188  legPtRes->AddEntry(fPtEta2, "#eta=1.75", "l");
189  legPtRes->AddEntry(fPtEta3, "#eta=2.75", "l");
190  legPtRes->Draw();
191 
192  cResolution->cd(2);
193  gPad->SetLogx();
194 
195  TF1* fEtaEta1 = etaResol.parameterEta("sigma", 0.25);
196  TF1* fEtaEta2 = etaResol.parameterEta("sigma", 1.75);
197  TF1* fEtaEta3 = etaResol.parameterEta("sigma", 2.75);
198 
199  fEtaEta1->SetLineWidth(1);
200  fEtaEta2->SetLineWidth(1);
201  fEtaEta3->SetLineWidth(1);
202  fEtaEta1->SetNpx(500);
203  fEtaEta2->SetNpx(500);
204  fEtaEta3->SetNpx(500);
205  fEtaEta1->SetLineColor(kRed);
206  fEtaEta2->SetLineColor(kBlue);
207  fEtaEta3->SetLineColor(kGreen);
208  fEtaEta1->SetRange(5.0, 500.);
209  fEtaEta2->SetRange(5.0, 500.);
210  fEtaEta3->SetRange(5.0, 500.);
211  fEtaEta1->Draw();
212  fEtaEta2->Draw("SAME");
213  fEtaEta3->Draw("SAME");
214 
215  fEtaEta1->SetTitle("#eta - Resolution");
216  fEtaEta1->GetHistogram()->SetXTitle("p_{T} [GeV]");
217  fEtaEta1->GetHistogram()->SetYTitle("#sigma(#eta)");
218  fEtaEta1->GetHistogram()->GetYaxis()->CenterTitle();
219  fEtaEta1->GetHistogram()->GetYaxis()->SetTitleOffset(1.33);
220 
221  TLegend* legEtaRes = new TLegend(0.6, 0.85, 0.85, 0.6);
222  legEtaRes->SetFillColor(10);
223  legEtaRes->AddEntry(fEtaEta1, "#eta=0.25", "l");
224  legEtaRes->AddEntry(fEtaEta2, "#eta=1.75", "l");
225  legEtaRes->AddEntry(fEtaEta3, "#eta=2.75", "l");
226  legEtaRes->Draw();
227 
228  cResolution->cd(3);
229  gPad->SetLogx();
230 
231  TF1* fPhiEta1 = phiResol.parameterEta("sigma", 0.25);
232  TF1* fPhiEta2 = phiResol.parameterEta("sigma", 1.75);
233  TF1* fPhiEta3 = phiResol.parameterEta("sigma", 2.75);
234 
235  fPhiEta1->SetLineWidth(1);
236  fPhiEta2->SetLineWidth(1);
237  fPhiEta3->SetLineWidth(1);
238  fPhiEta1->SetNpx(500);
239  fPhiEta2->SetNpx(500);
240  fPhiEta3->SetNpx(500);
241  fPhiEta1->SetLineColor(kRed);
242  fPhiEta2->SetLineColor(kBlue);
243  fPhiEta3->SetLineColor(kGreen);
244  fPhiEta1->SetRange(5.0, 500.);
245  fPhiEta2->SetRange(5.0, 500.);
246  fPhiEta3->SetRange(5.0, 500.);
247  fPhiEta1->Draw();
248  fPhiEta2->Draw("SAME");
249  fPhiEta3->Draw("SAME");
250 
251  fPhiEta1->SetTitle("#varphi - Resolution");
252  fPhiEta1->GetHistogram()->SetXTitle("p_{T} [GeV]");
253  fPhiEta1->GetHistogram()->SetYTitle("#sigma(#varphi)");
254  fPhiEta1->GetHistogram()->GetYaxis()->CenterTitle();
255  fPhiEta1->GetHistogram()->GetYaxis()->SetTitleOffset(1.33);
256 
257  TLegend* legPhiRes = new TLegend(0.6, 0.85, 0.85, 0.6);
258  legPhiRes->SetFillColor(10);
259  legPhiRes->AddEntry(fPhiEta1, "#eta=0.25", "l");
260  legPhiRes->AddEntry(fPhiEta2, "#eta=1.75", "l");
261  legPhiRes->AddEntry(fPhiEta3, "#eta=2.75", "l");
262  legPhiRes->Draw();
263 
264  // PLOT GEN VS SMEARED DISTRIBUTIONS
265  TCanvas* cSmearing = new TCanvas("Smearing", "Smearing", 100, 100, 1000, 600);
266  cSmearing->Divide(3, 2);
267 
268  cSmearing->cd(1);
269  gPad->SetLogy();
270  hRndPt->Draw();
271 
272  cSmearing->cd(2);
273  gPad->SetLogy();
274  hRndEta->Draw();
275 
276  cSmearing->cd(3);
277  gPad->SetLogy();
278  hRndPhi->Draw();
279 
280  cSmearing->cd(4);
281  gPad->SetLogy();
282  hGenPt->Draw();
283  hJetPt->Draw("SAME");
284  hJetPt->SetLineColor(kRed);
285  if (hGenPt->GetMaximum() < hJetPt->GetMaximum())
286  hGenPt->SetMaximum(1.1 * hJetPt->GetMaximum());
287  TLegend* legPt = new TLegend(0.6, 0.8, 0.85, 0.65);
288  legPt->SetFillColor(10);
289  legPt->AddEntry(hGenPt, "generated", "l");
290  legPt->AddEntry(hJetPt, "smeared", "l");
291  legPt->Draw();
292 
293  cSmearing->cd(5);
294  hGenEta->Draw();
295  hJetEta->Draw("SAME");
296  hJetEta->SetLineColor(kRed);
297  if (hGenEta->GetMaximum() < hJetEta->GetMaximum())
298  hGenEta->SetMaximum(1.1 * hJetEta->GetMaximum());
299  hGenEta->SetMinimum(0.0);
300  hGenEta->SetMaximum(1.5 * hGenEta->GetMaximum());
301  TLegend* legEta = new TLegend(0.6, 0.8, 0.85, 0.65);
302  legEta->SetFillColor(10);
303  legEta->AddEntry(hGenEta, "generated", "l");
304  legEta->AddEntry(hJetEta, "smeared", "l");
305  legEta->Draw();
306 
307  cSmearing->cd(6);
308  hGenPhi->Draw();
309  hJetPhi->Draw("SAME");
310  hJetPhi->SetLineColor(kRed);
311  if (hGenPhi->GetMaximum() < hJetPhi->GetMaximum())
312  hGenPhi->SetMaximum(1.1 * hJetPhi->GetMaximum());
313  hGenPhi->SetMinimum(0.0);
314  hGenPhi->SetMaximum(1.5 * hGenPhi->GetMaximum());
315  TLegend* legPhi = new TLegend(0.6, 0.8, 0.85, 0.65);
316  legPhi->SetFillColor(10);
317  legPhi->AddEntry(hGenPhi, "generated", "l");
318  legPhi->AddEntry(hJetPhi, "smeared", "l");
319  legPhi->Draw();
320 
321  app->Run();
322 
323  return 0;
324 }
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
mps_fire.i
i
Definition: mps_fire.py:429
main
int main(int argc, char **argv)
Definition: JetResolution_t.cc:29
JetResolution::parameterEta
TF1 * parameterEta(const std::string &parameterName, float eta)
Definition: JetResolution.cc:135
JetResolution::resolutionEtaPt
TF1 * resolutionEtaPt(float eta, float pt) const
Definition: JetResolution.cc:108
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
dumpRecoGeometry_cfg.era
era
Definition: dumpRecoGeometry_cfg.py:210
DMR_cfg.cerr
cerr
Definition: DMR_cfg.py:78
jetmet_cfg.nevts
nevts
Definition: jetmet_cfg.py:3
edm::FileInPath::fullPath
const std::string & fullPath() const
Definition: FileInPath.cc:144
gather_cfg.cout
cout
Definition: gather_cfg.py:144
Generated for CMSSW Reference Manual by doxygen 1.8.14