File indexing completed on 2024-10-16 05:03:42
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "CondFormats/JetMETObjects/interface/JetResolution.h"
0010 #include "FWCore/ParameterSet/interface/FileInPath.h"
0011
0012 #include <TROOT.h>
0013 #include <TApplication.h>
0014 #include <TRandom3.h>
0015 #include <TCanvas.h>
0016 #include <TLegend.h>
0017 #include <TH1F.h>
0018 #include <TMath.h>
0019
0020 #include <cstdlib>
0021 #include <iostream>
0022 #include <sstream>
0023 #include <string>
0024 #include <sys/stat.h>
0025
0026 using namespace std;
0027
0028
0029 int main(int argc, char** argv) {
0030 if (argc > 1 && string(argv[1]) == "--help") {
0031 cout << "USAGE: JetResolution_t --era <era> --alg <alg> --nevts <n> --gaussian" << endl;
0032 return 0;
0033 }
0034
0035 string era("Spring10");
0036 string alg("AK5Calo");
0037 unsigned nevts(10000);
0038 bool doGaussian(false);
0039
0040 for (int i = 1; i < argc; i++) {
0041 if (string(argv[i]) == "--era") {
0042 era = argv[i + 1];
0043 i++;
0044 } else if (string(argv[i]) == "--alg") {
0045 alg = argv[i + 1];
0046 i++;
0047 } else if (string(argv[i]) == "--nevts") {
0048 stringstream ss;
0049 ss << argv[i + 1];
0050 ss >> nevts;
0051 i++;
0052 } else if (string(argv[i]) == "--gaussian") {
0053 doGaussian = true;
0054 } else {
0055 cerr << "ERROR: unknown option '" << argv[i] << "'" << endl;
0056 return 0;
0057 }
0058 }
0059
0060 cout << "era: " << era << endl;
0061 cout << "alg: " << alg << endl;
0062 cout << "nevts: " << nevts << endl;
0063 cout << "gaussian: " << doGaussian << endl << endl;
0064
0065 string ptFileName =
0066 edm::FileInPath("CondFormats/JetMETObjects/data/" + era + "_PtResolution_" + alg + ".txt").fullPath();
0067 string etaFileName =
0068 edm::FileInPath("CondFormats/JetMETObjects/data/" + era + "_EtaResolution_" + alg + ".txt").fullPath();
0069 string phiFileName =
0070 edm::FileInPath("CondFormats/JetMETObjects/data/" + era + "_PhiResolution_" + alg + ".txt").fullPath();
0071
0072 cout << ptFileName << endl;
0073 cout << etaFileName << endl;
0074 cout << phiFileName << endl;
0075 cout << endl;
0076
0077 JetResolution ptResol(ptFileName, doGaussian);
0078 JetResolution etaResol(etaFileName, doGaussian);
0079 JetResolution phiResol(phiFileName, doGaussian);
0080
0081
0082 float pt = 200. * gRandom->Exp(0.1);
0083 float eta = gRandom->Uniform(-5.0, 5.0);
0084 float phi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
0085
0086 cout << "pT=" << pt << " eta=" << eta << " phi=" << phi << endl;
0087
0088 TF1* fPtResol = ptResol.resolutionEtaPt(eta, pt);
0089 TF1* fEtaResol = etaResol.resolutionEtaPt(eta, pt);
0090 TF1* fPhiResol = phiResol.resolutionEtaPt(eta, pt);
0091
0092 fPtResol->Print();
0093 cout << endl;
0094 fEtaResol->Print();
0095 cout << endl;
0096 fPhiResol->Print();
0097 cout << endl;
0098
0099
0100 TH1F* hRndPt = new TH1F("RndPt", ";random number", 200, 0.0, 5.0);
0101 TH1F* hGenPt = new TH1F("GenPt", ";p_{T} [GeV]", 100, 0., 250.);
0102 TH1F* hJetPt = new TH1F("JetPt", ";p_{T} [GeV]", 100, 0., 250.);
0103
0104 TH1F* hRndEta = new TH1F("RndEta", ";random number", 200, -5.0, 5.0);
0105 TH1F* hGenEta = new TH1F("GenEta", ";#eta", 51, -5., 5.);
0106 TH1F* hJetEta = new TH1F("JetEta", ";#eta", 51, -5., 5.);
0107
0108 TH1F* hRndPhi = new TH1F("RndPhi", ";random number", 200, -3.15, 3.15);
0109 TH1F* hGenPhi = new TH1F("GenPhi", ";#varphi", 41, -3.15, 3.15);
0110 TH1F* hJetPhi = new TH1F("JetPhi", ";#varphi", 41, -3.15, 3.15);
0111
0112 for (unsigned ievt = 0; ievt < nevts; ievt++) {
0113 if ((ievt + 1) % 1000 == 0)
0114 cout << ievt + 1 << " events processed." << endl;
0115
0116 float genpt = 200. * gRandom->Exp(0.1);
0117 if (genpt < 1.0)
0118 continue;
0119 float geneta = gRandom->Uniform(-5.0, 5.0);
0120 float genphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
0121
0122 fPtResol = ptResol.resolutionEtaPt(geneta, genpt);
0123 fEtaResol = etaResol.resolutionEtaPt(geneta, genpt);
0124 fPhiResol = phiResol.resolutionEtaPt(geneta, genpt);
0125
0126 float rndpt = fPtResol->GetRandom();
0127 float rndeta = fEtaResol->GetRandom();
0128 float rndphi = fPhiResol->GetRandom();
0129
0130 float jetpt = rndpt * genpt;
0131 float jeteta = rndeta + geneta;
0132 float jetphi = rndphi + genphi;
0133
0134 hRndPt->Fill(rndpt);
0135 hGenPt->Fill(genpt);
0136 hJetPt->Fill(jetpt);
0137
0138 hRndEta->Fill(rndeta);
0139 hGenEta->Fill(geneta);
0140 hJetEta->Fill(jeteta);
0141
0142 hRndPhi->Fill(rndphi);
0143 hGenPhi->Fill(genphi);
0144 hJetPhi->Fill(jetphi);
0145 }
0146
0147
0148 argc = 1;
0149 TApplication* app = new TApplication("JetResolution_t", &argc, argv);
0150 gROOT->SetStyle("Plain");
0151
0152
0153 TCanvas* cResolution = new TCanvas("Resolution", "Resolution", 0, 0, 1000, 400);
0154 cResolution->Divide(3, 1);
0155
0156 cResolution->cd(1);
0157 gPad->SetLogx();
0158
0159 TF1* fPtEta1 = ptResol.parameterEta("sigma", 0.25);
0160 TF1* fPtEta2 = ptResol.parameterEta("sigma", 1.75);
0161 TF1* fPtEta3 = ptResol.parameterEta("sigma", 2.75);
0162
0163 fPtEta1->SetLineWidth(1);
0164 fPtEta2->SetLineWidth(1);
0165 fPtEta3->SetLineWidth(1);
0166 fPtEta1->SetNpx(500);
0167 fPtEta2->SetNpx(500);
0168 fPtEta3->SetNpx(500);
0169 fPtEta1->SetLineColor(kRed);
0170 fPtEta2->SetLineColor(kBlue);
0171 fPtEta3->SetLineColor(kGreen);
0172 fPtEta1->SetRange(5.0, 500.);
0173 fPtEta2->SetRange(5.0, 500.);
0174 fPtEta3->SetRange(5.0, 500.);
0175 fPtEta1->Draw();
0176 fPtEta2->Draw("SAME");
0177 fPtEta3->Draw("SAME");
0178
0179 fPtEta1->SetTitle("p_{T} - Resolution");
0180 fPtEta1->GetHistogram()->SetXTitle("p_{T} [GeV]");
0181 fPtEta1->GetHistogram()->SetYTitle("#sigma(p_{T})/p_{T}");
0182 fPtEta1->GetHistogram()->GetYaxis()->CenterTitle();
0183 fPtEta1->GetHistogram()->GetYaxis()->SetTitleOffset(1.33);
0184
0185 TLegend* legPtRes = new TLegend(0.6, 0.85, 0.85, 0.6);
0186 legPtRes->SetFillColor(10);
0187 legPtRes->AddEntry(fPtEta1, "#eta=0.25", "l");
0188 legPtRes->AddEntry(fPtEta2, "#eta=1.75", "l");
0189 legPtRes->AddEntry(fPtEta3, "#eta=2.75", "l");
0190 legPtRes->Draw();
0191
0192 cResolution->cd(2);
0193 gPad->SetLogx();
0194
0195 TF1* fEtaEta1 = etaResol.parameterEta("sigma", 0.25);
0196 TF1* fEtaEta2 = etaResol.parameterEta("sigma", 1.75);
0197 TF1* fEtaEta3 = etaResol.parameterEta("sigma", 2.75);
0198
0199 fEtaEta1->SetLineWidth(1);
0200 fEtaEta2->SetLineWidth(1);
0201 fEtaEta3->SetLineWidth(1);
0202 fEtaEta1->SetNpx(500);
0203 fEtaEta2->SetNpx(500);
0204 fEtaEta3->SetNpx(500);
0205 fEtaEta1->SetLineColor(kRed);
0206 fEtaEta2->SetLineColor(kBlue);
0207 fEtaEta3->SetLineColor(kGreen);
0208 fEtaEta1->SetRange(5.0, 500.);
0209 fEtaEta2->SetRange(5.0, 500.);
0210 fEtaEta3->SetRange(5.0, 500.);
0211 fEtaEta1->Draw();
0212 fEtaEta2->Draw("SAME");
0213 fEtaEta3->Draw("SAME");
0214
0215 fEtaEta1->SetTitle("#eta - Resolution");
0216 fEtaEta1->GetHistogram()->SetXTitle("p_{T} [GeV]");
0217 fEtaEta1->GetHistogram()->SetYTitle("#sigma(#eta)");
0218 fEtaEta1->GetHistogram()->GetYaxis()->CenterTitle();
0219 fEtaEta1->GetHistogram()->GetYaxis()->SetTitleOffset(1.33);
0220
0221 TLegend* legEtaRes = new TLegend(0.6, 0.85, 0.85, 0.6);
0222 legEtaRes->SetFillColor(10);
0223 legEtaRes->AddEntry(fEtaEta1, "#eta=0.25", "l");
0224 legEtaRes->AddEntry(fEtaEta2, "#eta=1.75", "l");
0225 legEtaRes->AddEntry(fEtaEta3, "#eta=2.75", "l");
0226 legEtaRes->Draw();
0227
0228 cResolution->cd(3);
0229 gPad->SetLogx();
0230
0231 TF1* fPhiEta1 = phiResol.parameterEta("sigma", 0.25);
0232 TF1* fPhiEta2 = phiResol.parameterEta("sigma", 1.75);
0233 TF1* fPhiEta3 = phiResol.parameterEta("sigma", 2.75);
0234
0235 fPhiEta1->SetLineWidth(1);
0236 fPhiEta2->SetLineWidth(1);
0237 fPhiEta3->SetLineWidth(1);
0238 fPhiEta1->SetNpx(500);
0239 fPhiEta2->SetNpx(500);
0240 fPhiEta3->SetNpx(500);
0241 fPhiEta1->SetLineColor(kRed);
0242 fPhiEta2->SetLineColor(kBlue);
0243 fPhiEta3->SetLineColor(kGreen);
0244 fPhiEta1->SetRange(5.0, 500.);
0245 fPhiEta2->SetRange(5.0, 500.);
0246 fPhiEta3->SetRange(5.0, 500.);
0247 fPhiEta1->Draw();
0248 fPhiEta2->Draw("SAME");
0249 fPhiEta3->Draw("SAME");
0250
0251 fPhiEta1->SetTitle("#varphi - Resolution");
0252 fPhiEta1->GetHistogram()->SetXTitle("p_{T} [GeV]");
0253 fPhiEta1->GetHistogram()->SetYTitle("#sigma(#varphi)");
0254 fPhiEta1->GetHistogram()->GetYaxis()->CenterTitle();
0255 fPhiEta1->GetHistogram()->GetYaxis()->SetTitleOffset(1.33);
0256
0257 TLegend* legPhiRes = new TLegend(0.6, 0.85, 0.85, 0.6);
0258 legPhiRes->SetFillColor(10);
0259 legPhiRes->AddEntry(fPhiEta1, "#eta=0.25", "l");
0260 legPhiRes->AddEntry(fPhiEta2, "#eta=1.75", "l");
0261 legPhiRes->AddEntry(fPhiEta3, "#eta=2.75", "l");
0262 legPhiRes->Draw();
0263
0264
0265 TCanvas* cSmearing = new TCanvas("Smearing", "Smearing", 100, 100, 1000, 600);
0266 cSmearing->Divide(3, 2);
0267
0268 cSmearing->cd(1);
0269 gPad->SetLogy();
0270 hRndPt->Draw();
0271
0272 cSmearing->cd(2);
0273 gPad->SetLogy();
0274 hRndEta->Draw();
0275
0276 cSmearing->cd(3);
0277 gPad->SetLogy();
0278 hRndPhi->Draw();
0279
0280 cSmearing->cd(4);
0281 gPad->SetLogy();
0282 hGenPt->Draw();
0283 hJetPt->Draw("SAME");
0284 hJetPt->SetLineColor(kRed);
0285 if (hGenPt->GetMaximum() < hJetPt->GetMaximum())
0286 hGenPt->SetMaximum(1.1 * hJetPt->GetMaximum());
0287 TLegend* legPt = new TLegend(0.6, 0.8, 0.85, 0.65);
0288 legPt->SetFillColor(10);
0289 legPt->AddEntry(hGenPt, "generated", "l");
0290 legPt->AddEntry(hJetPt, "smeared", "l");
0291 legPt->Draw();
0292
0293 cSmearing->cd(5);
0294 hGenEta->Draw();
0295 hJetEta->Draw("SAME");
0296 hJetEta->SetLineColor(kRed);
0297 if (hGenEta->GetMaximum() < hJetEta->GetMaximum())
0298 hGenEta->SetMaximum(1.1 * hJetEta->GetMaximum());
0299 hGenEta->SetMinimum(0.0);
0300 hGenEta->SetMaximum(1.5 * hGenEta->GetMaximum());
0301 TLegend* legEta = new TLegend(0.6, 0.8, 0.85, 0.65);
0302 legEta->SetFillColor(10);
0303 legEta->AddEntry(hGenEta, "generated", "l");
0304 legEta->AddEntry(hJetEta, "smeared", "l");
0305 legEta->Draw();
0306
0307 cSmearing->cd(6);
0308 hGenPhi->Draw();
0309 hJetPhi->Draw("SAME");
0310 hJetPhi->SetLineColor(kRed);
0311 if (hGenPhi->GetMaximum() < hJetPhi->GetMaximum())
0312 hGenPhi->SetMaximum(1.1 * hJetPhi->GetMaximum());
0313 hGenPhi->SetMinimum(0.0);
0314 hGenPhi->SetMaximum(1.5 * hGenPhi->GetMaximum());
0315 TLegend* legPhi = new TLegend(0.6, 0.8, 0.85, 0.65);
0316 legPhi->SetFillColor(10);
0317 legPhi->AddEntry(hGenPhi, "generated", "l");
0318 legPhi->AddEntry(hJetPhi, "smeared", "l");
0319 legPhi->Draw();
0320
0321 app->Run();
0322
0323 return 0;
0324 }