Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:51:13

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // JetResolution_t
0004 // ---------------
0005 //
0006 //            11/05/2010 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
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   // SIMPLE TEST
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   // GENERATE SPECTRA AND SMEAR
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   // RUN ROOT APPLICATION AND DRAW BOTH DISTRIBUTIONS
0155   argc = 1;
0156   TApplication* app = new TApplication("JetResolution_t", &argc, argv);
0157   gROOT->SetStyle("Plain");
0158 
0159   // PLOT RESOLUTION FOR DIFFERENT ETA BINS
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   // PLOT GEN VS SMEARED DISTRIBUTIONS
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 }