Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-10-16 05:03:42

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 #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   // SIMPLE TEST
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   // GENERATE SPECTRA AND SMEAR
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   // RUN ROOT APPLICATION AND DRAW BOTH DISTRIBUTIONS
0148   argc = 1;
0149   TApplication* app = new TApplication("JetResolution_t", &argc, argv);
0150   gROOT->SetStyle("Plain");
0151 
0152   // PLOT RESOLUTION FOR DIFFERENT ETA BINS
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   // PLOT GEN VS SMEARED DISTRIBUTIONS
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 }