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