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