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