00001
00002
00003
00004
00005
00006
00008
00009
00010 #include "CondFormats/JetMETObjects/interface/JetResolution.h"
00011
00012
00013 #include <TROOT.h>
00014 #include <TApplication.h>
00015 #include <TRandom3.h>
00016 #include <TCanvas.h>
00017 #include <TLegend.h>
00018 #include <TH1F.h>
00019 #include <TMath.h>
00020
00021 #include <cstdlib>
00022 #include <iostream>
00023 #include <sstream>
00024 #include <string>
00025 #include <sys/stat.h>
00026
00027
00028 using namespace std;
00029
00030
00031
00032 int main(int argc,char**argv)
00033 {
00034 if (argc>1&&string(argv[1])=="--help") {
00035 cout<<"USAGE: JetResolution_t --era <era> --alg <alg> --nevts <n> --gaussian"<<endl;
00036 return 0;
00037 }
00038
00039 string era("Spring10");
00040 string alg("AK5Calo");
00041 unsigned nevts(10000);
00042 bool doGaussian(false);
00043
00044 for (int i=1;i<argc;i++) {
00045 if (string(argv[i])=="--era") { era = argv[i+1]; i++; }
00046 else if (string(argv[i])=="--alg") { alg = argv[i+1]; i++; }
00047 else if (string(argv[i])=="--nevts") { stringstream ss; ss<<argv[i+1]; ss>>nevts; i++; }
00048 else if (string(argv[i])=="--gaussian") { doGaussian = true; }
00049 else {
00050 cerr<<"ERROR: unknown option '"<<argv[i]<<"'"<<endl;
00051 return 0;
00052 }
00053 }
00054
00055 cout<<"era: "<<era<<endl;
00056 cout<<"alg: "<<alg<<endl;
00057 cout<<"nevts: "<<nevts<<endl;
00058 cout<<"gaussian: "<<doGaussian<<endl<<endl;
00059
00060
00061 string cmssw_base(getenv("CMSSW_BASE"));
00062 string cmssw_release_base(getenv("CMSSW_RELEASE_BASE"));
00063 string path = cmssw_base + "/src/CondFormats/JetMETObjects/data";
00064 struct stat st;
00065 if (stat(path.c_str(),&st)!=0)
00066 path = cmssw_release_base + "/src/CondFormats/JetMETObjects/data";
00067 if (stat(path.c_str(),&st)!=0) {
00068 cerr<<"ERROR: tried to set path but failed, abort."<<endl;
00069 return 0;
00070 }
00071
00072
00073 string ptFileName = path + "/" + era + "_PtResolution_" +alg+".txt";
00074 string etaFileName = path + "/" + era + "_EtaResolution_"+alg+".txt";
00075 string phiFileName = path + "/" + era + "_PhiResolution_"+alg+".txt";
00076
00077 cout<<ptFileName<<endl;
00078 cout<<etaFileName<<endl;
00079 cout<<phiFileName<<endl;
00080 cout<<endl;
00081
00082 JetResolution ptResol(ptFileName,doGaussian);
00083 JetResolution etaResol(etaFileName,doGaussian);
00084 JetResolution phiResol(phiFileName,doGaussian);
00085
00086
00087
00088 float pt = 200.*gRandom->Exp(0.1);
00089 float eta = gRandom->Uniform(-5.0,5.0);
00090 float phi = gRandom->Uniform(-TMath::Pi(),TMath::Pi());
00091
00092 cout<<"pT="<<pt<<" eta="<<eta<<" phi="<<phi<<endl;
00093
00094 TF1* fPtResol = ptResol.resolutionEtaPt(eta,pt);
00095 TF1* fEtaResol = etaResol.resolutionEtaPt(eta,pt);
00096 TF1* fPhiResol = phiResol.resolutionEtaPt(eta,pt);
00097
00098 fPtResol ->Print(); cout<<endl;
00099 fEtaResol->Print(); cout<<endl;
00100 fPhiResol->Print(); cout<<endl;
00101
00102
00103
00104 TH1F* hRndPt = new TH1F("RndPt",";random number",200,0.0,5.0);
00105 TH1F* hGenPt = new TH1F("GenPt",";p_{T} [GeV]",100,0.,250.);
00106 TH1F* hJetPt = new TH1F("JetPt",";p_{T} [GeV]",100,0.,250.);
00107
00108 TH1F* hRndEta = new TH1F("RndEta",";random number",200,-5.0,5.0);
00109 TH1F* hGenEta = new TH1F("GenEta",";#eta",51,-5.,5.);
00110 TH1F* hJetEta = new TH1F("JetEta",";#eta",51,-5.,5.);
00111
00112 TH1F* hRndPhi = new TH1F("RndPhi",";random number",200,-3.15,3.15);
00113 TH1F* hGenPhi = new TH1F("GenPhi",";#varphi",41,-3.15,3.15);
00114 TH1F* hJetPhi = new TH1F("JetPhi",";#varphi",41,-3.15,3.15);
00115
00116 for (unsigned ievt=0;ievt<nevts;ievt++) {
00117
00118 if ((ievt+1)%1000==0) cout<<ievt+1<<" events processed."<<endl;
00119
00120 float genpt = 200.*gRandom->Exp(0.1); if (genpt<1.0) continue;
00121 float geneta = gRandom->Uniform(-5.0,5.0);
00122 float genphi = gRandom->Uniform(-TMath::Pi(),TMath::Pi());
00123
00124 fPtResol = ptResol.resolutionEtaPt(geneta,genpt);
00125 fEtaResol = etaResol.resolutionEtaPt(geneta,genpt);
00126 fPhiResol = phiResol.resolutionEtaPt(geneta,genpt);
00127
00128 float rndpt = fPtResol ->GetRandom();
00129 float rndeta = fEtaResol->GetRandom();
00130 float rndphi = fPhiResol->GetRandom();
00131
00132 float jetpt = rndpt*genpt;
00133 float jeteta = rndeta+geneta;
00134 float jetphi = rndphi+genphi;
00135
00136 hRndPt->Fill(rndpt);
00137 hGenPt->Fill(genpt);
00138 hJetPt->Fill(jetpt);
00139
00140 hRndEta->Fill(rndeta);
00141 hGenEta->Fill(geneta);
00142 hJetEta->Fill(jeteta);
00143
00144 hRndPhi->Fill(rndphi);
00145 hGenPhi->Fill(genphi);
00146 hJetPhi->Fill(jetphi);
00147 }
00148
00149
00150 argc=1;
00151 TApplication* app = new TApplication("JetResolution_t",&argc,argv);
00152 gROOT->SetStyle("Plain");
00153
00154
00155
00156 TCanvas* cResolution = new TCanvas("Resolution","Resolution",0,0,1000,400);
00157 cResolution->Divide(3,1);
00158
00159 cResolution->cd(1);
00160 gPad->SetLogx();
00161
00162 TF1* fPtEta1 = ptResol.parameterEta("sigma",0.25);
00163 TF1* fPtEta2 = ptResol.parameterEta("sigma",1.75);
00164 TF1* fPtEta3 = ptResol.parameterEta("sigma",2.75);
00165
00166 fPtEta1->SetLineWidth(1); fPtEta2->SetLineWidth(1); fPtEta3->SetLineWidth(1);
00167 fPtEta1->SetNpx(500); fPtEta2->SetNpx(500); fPtEta3->SetNpx(500);
00168 fPtEta1->SetLineColor(kRed); fPtEta2->SetLineColor(kBlue); fPtEta3->SetLineColor(kGreen);
00169 fPtEta1->SetRange(5.0,500.); fPtEta2->SetRange(5.0,500.); fPtEta3->SetRange(5.0,500.);
00170 fPtEta1->Draw(); fPtEta2->Draw("SAME"); fPtEta3->Draw("SAME");
00171
00172 fPtEta1->SetTitle("p_{T} - Resolution");
00173 fPtEta1->GetHistogram()->SetXTitle("p_{T} [GeV]");
00174 fPtEta1->GetHistogram()->SetYTitle("#sigma(p_{T})/p_{T}");
00175 fPtEta1->GetHistogram()->GetYaxis()->CenterTitle();
00176 fPtEta1->GetHistogram()->GetYaxis()->SetTitleOffset(1.33);
00177
00178 TLegend* legPtRes = new TLegend(0.6,0.85,0.85,0.6);
00179 legPtRes->SetFillColor(10);
00180 legPtRes->AddEntry(fPtEta1,"#eta=0.25","l");
00181 legPtRes->AddEntry(fPtEta2,"#eta=1.75","l");
00182 legPtRes->AddEntry(fPtEta3,"#eta=2.75","l");
00183 legPtRes->Draw();
00184
00185
00186 cResolution->cd(2);
00187 gPad->SetLogx();
00188
00189 TF1* fEtaEta1 = etaResol.parameterEta("sigma",0.25);
00190 TF1* fEtaEta2 = etaResol.parameterEta("sigma",1.75);
00191 TF1* fEtaEta3 = etaResol.parameterEta("sigma",2.75);
00192
00193 fEtaEta1->SetLineWidth(1); fEtaEta2->SetLineWidth(1); fEtaEta3->SetLineWidth(1);
00194 fEtaEta1->SetNpx(500); fEtaEta2->SetNpx(500); fEtaEta3->SetNpx(500);
00195 fEtaEta1->SetLineColor(kRed); fEtaEta2->SetLineColor(kBlue); fEtaEta3->SetLineColor(kGreen);
00196 fEtaEta1->SetRange(5.0,500.); fEtaEta2->SetRange(5.0,500.); fEtaEta3->SetRange(5.0,500.);
00197 fEtaEta1->Draw(); fEtaEta2->Draw("SAME"); fEtaEta3->Draw("SAME");
00198
00199 fEtaEta1->SetTitle("#eta - Resolution");
00200 fEtaEta1->GetHistogram()->SetXTitle("p_{T} [GeV]");
00201 fEtaEta1->GetHistogram()->SetYTitle("#sigma(#eta)");
00202 fEtaEta1->GetHistogram()->GetYaxis()->CenterTitle();
00203 fEtaEta1->GetHistogram()->GetYaxis()->SetTitleOffset(1.33);
00204
00205 TLegend* legEtaRes = new TLegend(0.6,0.85,0.85,0.6);
00206 legEtaRes->SetFillColor(10);
00207 legEtaRes->AddEntry(fEtaEta1,"#eta=0.25","l");
00208 legEtaRes->AddEntry(fEtaEta2,"#eta=1.75","l");
00209 legEtaRes->AddEntry(fEtaEta3,"#eta=2.75","l");
00210 legEtaRes->Draw();
00211
00212
00213 cResolution->cd(3);
00214 gPad->SetLogx();
00215
00216 TF1* fPhiEta1 = phiResol.parameterEta("sigma",0.25);
00217 TF1* fPhiEta2 = phiResol.parameterEta("sigma",1.75);
00218 TF1* fPhiEta3 = phiResol.parameterEta("sigma",2.75);
00219
00220 fPhiEta1->SetLineWidth(1); fPhiEta2->SetLineWidth(1); fPhiEta3->SetLineWidth(1);
00221 fPhiEta1->SetNpx(500); fPhiEta2->SetNpx(500); fPhiEta3->SetNpx(500);
00222 fPhiEta1->SetLineColor(kRed); fPhiEta2->SetLineColor(kBlue); fPhiEta3->SetLineColor(kGreen);
00223 fPhiEta1->SetRange(5.0,500.); fPhiEta2->SetRange(5.0,500.); fPhiEta3->SetRange(5.0,500.);
00224 fPhiEta1->Draw(); fPhiEta2->Draw("SAME"); fPhiEta3->Draw("SAME");
00225
00226 fPhiEta1->SetTitle("#varphi - Resolution");
00227 fPhiEta1->GetHistogram()->SetXTitle("p_{T} [GeV]");
00228 fPhiEta1->GetHistogram()->SetYTitle("#sigma(#varphi)");
00229 fPhiEta1->GetHistogram()->GetYaxis()->CenterTitle();
00230 fPhiEta1->GetHistogram()->GetYaxis()->SetTitleOffset(1.33);
00231
00232 TLegend* legPhiRes = new TLegend(0.6,0.85,0.85,0.6);
00233 legPhiRes->SetFillColor(10);
00234 legPhiRes->AddEntry(fPhiEta1,"#eta=0.25","l");
00235 legPhiRes->AddEntry(fPhiEta2,"#eta=1.75","l");
00236 legPhiRes->AddEntry(fPhiEta3,"#eta=2.75","l");
00237 legPhiRes->Draw();
00238
00239
00240
00241
00242 TCanvas* cSmearing = new TCanvas("Smearing","Smearing",100,100,1000,600);
00243 cSmearing->Divide(3,2);
00244
00245 cSmearing->cd(1);
00246 gPad->SetLogy();
00247 hRndPt->Draw();
00248
00249 cSmearing->cd(2);
00250 gPad->SetLogy();
00251 hRndEta->Draw();
00252
00253 cSmearing->cd(3);
00254 gPad->SetLogy();
00255 hRndPhi->Draw();
00256
00257 cSmearing->cd(4);
00258 gPad->SetLogy();
00259 hGenPt->Draw();
00260 hJetPt->Draw("SAME");
00261 hJetPt->SetLineColor(kRed);
00262 if (hGenPt->GetMaximum()<hJetPt->GetMaximum())
00263 hGenPt->SetMaximum(1.1*hJetPt->GetMaximum());
00264 TLegend* legPt = new TLegend(0.6,0.8,0.85,0.65);
00265 legPt->SetFillColor(10);
00266 legPt->AddEntry(hGenPt,"generated","l");
00267 legPt->AddEntry(hJetPt,"smeared","l");
00268 legPt->Draw();
00269
00270 cSmearing->cd(5);
00271 hGenEta->Draw();
00272 hJetEta->Draw("SAME");
00273 hJetEta->SetLineColor(kRed);
00274 if (hGenEta->GetMaximum()<hJetEta->GetMaximum())
00275 hGenEta->SetMaximum(1.1*hJetEta->GetMaximum());
00276 hGenEta->SetMinimum(0.0);
00277 hGenEta->SetMaximum(1.5*hGenEta->GetMaximum());
00278 TLegend* legEta = new TLegend(0.6,0.8,0.85,0.65);
00279 legEta->SetFillColor(10);
00280 legEta->AddEntry(hGenEta,"generated","l");
00281 legEta->AddEntry(hJetEta,"smeared","l");
00282 legEta->Draw();
00283
00284 cSmearing->cd(6);
00285 hGenPhi->Draw();
00286 hJetPhi->Draw("SAME");
00287 hJetPhi->SetLineColor(kRed);
00288 if (hGenPhi->GetMaximum()<hJetPhi->GetMaximum())
00289 hGenPhi->SetMaximum(1.1*hJetPhi->GetMaximum());
00290 hGenPhi->SetMinimum(0.0);
00291 hGenPhi->SetMaximum(1.5*hGenPhi->GetMaximum());
00292 TLegend* legPhi = new TLegend(0.6,0.8,0.85,0.65);
00293 legPhi->SetFillColor(10);
00294 legPhi->AddEntry(hGenPhi,"generated","l");
00295 legPhi->AddEntry(hJetPhi,"smeared","l");
00296 legPhi->Draw();
00297
00298 app->Run();
00299
00300 return 0;
00301 }
00302
1.7.3