Benchmark/src/GenericBenchmark.cc

0001 #include "RecoParticleFlow/Benchmark/interface/GenericBenchmark.h"
0002 #include "RecoParticleFlow/Benchmark/interface/BenchmarkTree.h"
0003
0004 #include <TROOT.h>
0005 #include <TFile.h>
0006
0007 //Colin: it seems a bit strange to use the preprocessor for that kind of
0008 //thing. Looks like all these macros could be replaced by plain functions.
0009
0010 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
0011 #define BOOK1D(name, title, nbinsx, lowx, highx) \
0012   h##name =                                      \
0013       DQM ? DQM->book1D(#name, title, nbinsx, lowx, highx)->getTH1F() : new TH1F(#name, title, nbinsx, lowx, highx)
0014
0015 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
0016 #define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)                            \
0017   h##name = DQM ? DQM->book2D(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)->getTH2F() \
0018                 : new TH2F(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
0019
0020 // all versions OK
0021 // preprocesor macro for setting axis titles
0022
0023 #define SETAXES(name, xtitle, ytitle)    \
0024   h##name->GetXaxis()->SetTitle(xtitle); \
0025   h##name->GetYaxis()->SetTitle(ytitle)
0026
0027 #define ET (PlotAgainstReco_) ? "reconstructed E_{T} [GeV]" : "generated E_{T} [GeV]"
0028 #define ETA (PlotAgainstReco_) ? "reconstructed #eta" : "generated #eta"
0029 #define PHI (PlotAgainstReco_) ? "reconstructed #phi (rad)" : "generated #phi (rad)"
0030
0031 using namespace std;
0032
0033 GenericBenchmark::GenericBenchmark() {}
0034
0035 GenericBenchmark::~GenericBenchmark() noexcept(false) {}
0036
0037 void GenericBenchmark::setup(DQMStore* DQM,
0038                              bool PlotAgainstReco_,
0039                              float minDeltaEt,
0040                              float maxDeltaEt,
0041                              float minDeltaPhi,
0042                              float maxDeltaPhi,
0043                              bool doMetPlots) {
0044   //std::cout << "minDeltaPhi = " << minDeltaPhi << std::endl;
0045
0046   // CMSSW_2_X_X
0047   // use bare Root if no DQM (FWLite applications)
0048   //if (!DQM)
0049   //{
0050   //  file_ = new TFile("pfmetBenchmark.root", "recreate");
0051   //  cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. "<<endl;
0052   //}
0053   // Book Histograms
0054
0055   //std::cout << "FL : pwd = ";
0056   //gDirectory->pwd();
0057   //std::cout << std::endl;
0058
0059   int nbinsEt = 1000;
0060   float minEt = 0;
0061   float maxEt = 1000;
0062
0063   //float minDeltaEt = -100;
0064   //float maxDeltaEt = 50;
0065
0066   int nbinsEta = 200;
0067   float minEta = -5;
0068   float maxEta = 5;
0069
0070   int nbinsDeltaEta = 1000;
0071   float minDeltaEta = -0.5;
0072   float maxDeltaEta = 0.5;
0073
0074   int nbinsDeltaPhi = 1000;
0075   //float minDeltaPhi = -0.5;
0076   //float maxDeltaPhi = 0.5;
0077
0078   // delta et quantities
0079   BOOK1D(DeltaEt, "#DeltaE_{T}", nbinsEt, minDeltaEt, maxDeltaEt);
0080   BOOK1D(DeltaEx, "#DeltaE_{X}", nbinsEt, minDeltaEt, maxDeltaEt);
0081   BOOK1D(DeltaEy, "#DeltaE_{Y}", nbinsEt, minDeltaEt, maxDeltaEt);
0082   BOOK2D(DeltaEtvsEt, "#DeltaE_{T} vs E_{T}", nbinsEt, minEt, maxEt, nbinsEt, minDeltaEt, maxDeltaEt);
0083   BOOK2D(DeltaEtOverEtvsEt, "#DeltaE_{T}/E_{T} vsE_{T}", nbinsEt, minEt, maxEt, nbinsEt, -1, 1);
0084   BOOK2D(DeltaEtvsEta, "#DeltaE_{T} vs #eta", nbinsEta, minEta, maxEta, nbinsEt, minDeltaEt, maxDeltaEt);
0085   BOOK2D(DeltaEtOverEtvsEta, "#DeltaE_{T}/E_{T} vs #eta", nbinsEta, minEta, maxEta, 100, -1, 1);
0086   BOOK2D(DeltaEtvsPhi, "#DeltaE_{T} vs #phi", 200, -M_PI, M_PI, nbinsEt, minDeltaEt, maxDeltaEt);
0087   BOOK2D(DeltaEtOverEtvsPhi, "#DeltaE_{T}/E_{T} vs #Phi", 200, -M_PI, M_PI, 100, -1, 1);
0088   BOOK2D(DeltaEtvsDeltaR, "#DeltaE_{T} vs #DeltaR", 100, 0, 1, nbinsEt, minDeltaEt, maxDeltaEt);
0089   BOOK2D(DeltaEtOverEtvsDeltaR, "#DeltaE_{T}/E_{T} vs #DeltaR", 100, 0, 1, 100, -1, 1);
0090
0091   // delta eta quantities
0092   BOOK1D(DeltaEta, "#Delta#eta", nbinsDeltaEta, minDeltaEta, maxDeltaEta);
0093   BOOK2D(DeltaEtavsEt, "#Delta#eta vs E_{T}", nbinsEt, minEt, maxEt, nbinsDeltaEta, minDeltaEta, maxDeltaEta);
0094   BOOK2D(DeltaEtavsEta, "#Delta#eta vs #eta", nbinsEta, minEta, maxEta, nbinsDeltaEta, minDeltaEta, maxDeltaEta);
0095
0096   // delta phi quantities
0097   BOOK1D(DeltaPhi, "#Delta#phi", nbinsDeltaPhi, minDeltaPhi, maxDeltaPhi);
0098   BOOK2D(DeltaPhivsEt, "#Delta#phi vs E_{T}", nbinsEt, minEt, maxEt, nbinsDeltaPhi, minDeltaPhi, maxDeltaPhi);
0099   BOOK2D(DeltaPhivsEta, "#Delta#phi vs #eta", nbinsEta, minEta, maxEta, nbinsDeltaPhi, minDeltaPhi, maxDeltaPhi);
0100
0101   // delta R quantities
0102   BOOK1D(DeltaR, "#DeltaR", 100, 0, 1);
0103   BOOK2D(DeltaRvsEt, "#DeltaR vs E_{T}", nbinsEt, minEt, maxEt, 100, 0, 1);
0104   BOOK2D(DeltaRvsEta, "#DeltaR vs #eta", nbinsEta, minEta, maxEta, 100, 0, 1);
0105
0106   BOOK1D(NRec, "Number of reconstructed objects", 20, 0, 20);
0107
0108   // seen and gen distributions, for efficiency computation
0109   BOOK1D(EtaSeen, "seen #eta", 100, -5, 5);
0110   BOOK1D(PhiSeen, "seen #phi", 100, -3.5, 3.5);
0111   BOOK1D(EtSeen, "seen E_{T}", nbinsEt, minEt, maxEt);
0112   BOOK2D(EtvsEtaSeen, "seen E_{T} vs eta", 100, -5, 5, 200, 0, 200);
0113   BOOK2D(EtvsPhiSeen, "seen E_{T} vs seen #phi", 100, -3.5, 3.5, 200, 0, 200);
0114
0115   BOOK1D(PhiRec, "Rec #phi", 100, -3.5, 3.5);
0116   BOOK1D(EtRec, "Rec E_{T}", nbinsEt, minEt, maxEt);
0117   BOOK1D(ExRec, "Rec E_{X}", nbinsEt, -maxEt, maxEt);
0118   BOOK1D(EyRec, "Rec E_{Y}", nbinsEt, -maxEt, maxEt);
0119
0120   BOOK2D(EtRecvsEt, "Rec E_{T} vs E_{T}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
0121   BOOK2D(EtRecOverTrueEtvsTrueEt, "Rec E_{T} / E_{T} vs E_{T}", nbinsEt, minEt, maxEt, 1000, 0., 100.);
0122
0123   BOOK1D(EtaGen, "generated #eta", 100, -5, 5);
0124   BOOK1D(PhiGen, "generated #phi", 100, -3.5, 3.5);
0125   BOOK1D(EtGen, "generated E_{T}", nbinsEt, minEt, maxEt);
0126   BOOK2D(EtvsEtaGen, "generated E_{T} vs generated #eta", 100, -5, 5, 200, 0, 200);
0127   BOOK2D(EtvsPhiGen, "generated E_{T} vs generated #phi", 100, -3.5, 3.5, 200, 0, 200);
0128
0129   BOOK1D(NGen, "Number of generated objects", 20, 0, 20);
0130
0131   if (doMetPlots) {
0132     BOOK1D(SumEt, "SumEt", 1000, 0., 3000.);
0133     BOOK1D(TrueSumEt, "TrueSumEt", 1000, 0., 3000.);
0134     BOOK2D(DeltaSetvsSet, "#DeltaSEt vs trueSEt", 3000, 0., 3000., 1000, -1000., 1000.);
0135     BOOK2D(DeltaMexvsSet, "#DeltaMEX vs trueSEt", 3000, 0., 3000., 1000, -400., 400.);
0136     BOOK2D(DeltaSetOverSetvsSet, "#DeltaSetOverSet vs trueSet", 3000, 0., 3000., 1000, -1., 1.);
0137     BOOK2D(RecSetvsTrueSet, "Set vs trueSet", 3000, 0., 3000., 3000, 0., 3000.);
0138     BOOK2D(RecSetOverTrueSetvsTrueSet, "Set/trueSet vs trueSet", 3000, 0., 3000., 500, 0., 2.);
0139     BOOK2D(TrueMexvsTrueSet, "trueMex vs trueSet", 3000, 0., 3000., nbinsEt, -maxEt, maxEt);
0140   }
0141   // number of truth particles found within given cone radius of reco
0142   //BOOK2D(NumInConeVsConeSize,NumInConeVsConeSize,100,0,1,25,0,25);
0143
0144   // Set Axis Titles
0145
0146   // delta et quantities
0147   SETAXES(DeltaEt, "#DeltaE_{T} [GeV]", "");
0148   SETAXES(DeltaEx, "#DeltaE_{X} [GeV]", "");
0149   SETAXES(DeltaEy, "#DeltaE_{Y} [GeV]", "");
0150   SETAXES(DeltaEtvsEt, ET, "#DeltaE_{T} [GeV]");
0151   SETAXES(DeltaEtOverEtvsEt, ET, "#DeltaE_{T}/E_{T}");
0152   SETAXES(DeltaEtvsEta, ETA, "#DeltaE_{T} [GeV]");
0153   SETAXES(DeltaEtOverEtvsEta, ETA, "#DeltaE_{T}/E_{T}");
0154   SETAXES(DeltaEtvsPhi, PHI, "#DeltaE_{T} [GeV]");
0155   SETAXES(DeltaEtOverEtvsPhi, PHI, "#DeltaE_{T}/E_{T}");
0156   SETAXES(DeltaEtvsDeltaR, "#DeltaR", "#DeltaE_{T} [GeV]");
0157   SETAXES(DeltaEtOverEtvsDeltaR, "#DeltaR", "#DeltaE_{T}/E_{T}");
0158
0159   // delta eta quantities
0160   SETAXES(DeltaEta, "#Delta#eta", "Events");
0161   SETAXES(DeltaEtavsEt, ET, "#Delta#eta");
0162   SETAXES(DeltaEtavsEta, ETA, "#Delta#eta");
0163
0164   // delta phi quantities
0165   SETAXES(DeltaPhi, "#Delta#phi [rad]", "Events");
0166   SETAXES(DeltaPhivsEt, ET, "#Delta#phi [rad]");
0167   SETAXES(DeltaPhivsEta, ETA, "#Delta#phi [rad]");
0168
0169   // delta R quantities
0170   SETAXES(DeltaR, "#DeltaR", "Events");
0171   SETAXES(DeltaRvsEt, ET, "#DeltaR");
0172   SETAXES(DeltaRvsEta, ETA, "#DeltaR");
0173
0174   SETAXES(NRec, "Number of Rec Objects", "");
0175
0176   SETAXES(EtaSeen, "seen #eta", "");
0177   SETAXES(PhiSeen, "seen #phi [rad]", "");
0178   SETAXES(EtSeen, "seen E_{T} [GeV]", "");
0179   SETAXES(EtvsEtaSeen, "seen #eta", "seen E_{T}");
0180   SETAXES(EtvsPhiSeen, "seen #phi [rad]", "seen E_{T}");
0181
0182   SETAXES(PhiRec, "#phi [rad]", "");
0183   SETAXES(EtRec, "E_{T} [GeV]", "");
0184   SETAXES(ExRec, "E_{X} [GeV]", "");
0185   SETAXES(EyRec, "E_{Y} [GeV]", "");
0186   SETAXES(EtRecvsEt, ET, "Rec E_{T} [GeV]");
0187   SETAXES(EtRecOverTrueEtvsTrueEt, ET, "Rec E_{T} / E_{T} [GeV]");
0188
0189   SETAXES(EtaGen, "generated #eta", "");
0190   SETAXES(PhiGen, "generated #phi [rad]", "");
0191   SETAXES(EtGen, "generated E_{T} [GeV]", "");
0192   SETAXES(EtvsPhiGen, "generated #phi [rad]", "generated E_{T} [GeV]");
0193   SETAXES(EtvsEtaGen, "generated #eta", "generated E_{T} [GeV]");
0194
0195   SETAXES(NGen, "Number of Gen Objects", "");
0196
0197   if (doMetPlots) {
0198     SETAXES(SumEt, "SumEt [GeV]", "");
0199     SETAXES(TrueSumEt, "TrueSumEt [GeV]", "");
0200     SETAXES(DeltaSetvsSet, "TrueSumEt", "#DeltaSumEt [GeV]");
0201     SETAXES(DeltaMexvsSet, "TrueSumEt", "#DeltaMEX [GeV]");
0202     SETAXES(DeltaSetOverSetvsSet, "TrueSumEt", "#DeltaSumEt/trueSumEt");
0203     SETAXES(RecSetvsTrueSet, "TrueSumEt", "SumEt");
0204     SETAXES(RecSetOverTrueSetvsTrueSet, "TrueSumEt", "SumEt/trueSumEt");
0205     SETAXES(TrueMexvsTrueSet, "TrueSumEt", "TrueMEX");
0206   }
0207
0208   TDirectory* oldpwd = gDirectory;
0209
0210   TIter next(gROOT->GetListOfFiles());
0211
0212   const bool debug = false;
0213
0214   while (TFile* file = (TFile*)next()) {
0215     if (debug)
0216       cout << "file " << file->GetName() << endl;
0217   }
0218   if (DQM) {
0219     cout << "DQM subdir" << endl;
0220     cout << DQM->pwd().c_str() << endl;
0221
0222     DQM->cd(DQM->pwd());
0223   }
0224
0225   if (debug) {
0226     cout << "current dir" << endl;
0227     gDirectory->pwd();
0228   }
0229
0230   oldpwd->cd();
0231   //gDirectory->pwd();
0232
0233   doMetPlots_ = doMetPlots;
0234
0235   //   tree_ = new BenchmarkTree("genericBenchmark", "Generic Benchmark TTree");
0236 }
0237
0238 bool GenericBenchmark::accepted(const reco::Candidate* particle,
0239                                 double ptCut,
0240                                 double minEtaCut,
0241                                 double maxEtaCut) const {
0242   //skip reconstructed PFJets with p_t < recPt_cut
0243   if (particle->pt() < ptCut and ptCut != -1.)
0244     return false;
0245
0246   if (fabs(particle->eta()) > maxEtaCut and maxEtaCut > 0)
0247     return false;
0248   if (fabs(particle->eta()) < minEtaCut and minEtaCut > 0)
0249     return false;
0250
0251   //accepted!
0252   return true;
0253 }
0254
0255 void GenericBenchmark::fillHistos(const reco::Candidate* genParticle,
0256                                   const reco::Candidate* recParticle,
0257                                   double deltaR_cut,
0258                                   bool plotAgainstReco) {
0259   // get the quantities to place on the denominator and/or divide by
0260   double et = genParticle->et();
0261   double eta = genParticle->eta();
0262   double phi = genParticle->phi();
0263   //std::cout << "FL : et = " << et << std::endl;
0264   //std::cout << "FL : eta = " << eta << std::endl;
0265   //std::cout << "FL : phi = " << phi << std::endl;
0266   //std::cout << "FL : rec et = " << recParticle->et() << std::endl;
0267   //std::cout << "FL : rec eta = " << recParticle->eta() << std::endl;
0268   //std::cout << "FL : rec phi = " <<recParticle-> phi() << std::endl;
0269
0270   if (plotAgainstReco) {
0271     et = recParticle->et();
0272     eta = recParticle->eta();
0273     phi = recParticle->phi();
0274   }
0275
0276   // get the delta quantities
0277   double deltaEt = algo_->deltaEt(recParticle, genParticle);
0278   double deltaR = algo_->deltaR(recParticle, genParticle);
0279   double deltaEta = algo_->deltaEta(recParticle, genParticle);
0280   double deltaPhi = algo_->deltaPhi(recParticle, genParticle);
0281
0282   //std::cout << "FL :deltaR_cut = " << deltaR_cut << std::endl;
0283   //std::cout << "FL :deltaR = " << deltaR << std::endl;
0284
0285   if (deltaR > deltaR_cut && deltaR_cut != -1.)
0286     return;
0287
0288   hDeltaEt->Fill(deltaEt);
0289   hDeltaEx->Fill(recParticle->px() - genParticle->px());
0290   hDeltaEy->Fill(recParticle->py() - genParticle->py());
0291   hDeltaEtvsEt->Fill(et, deltaEt);
0292   hDeltaEtOverEtvsEt->Fill(et, deltaEt / et);
0293   hDeltaEtvsEta->Fill(eta, deltaEt);
0294   hDeltaEtOverEtvsEta->Fill(eta, deltaEt / et);
0295   hDeltaEtvsPhi->Fill(phi, deltaEt);
0296   hDeltaEtOverEtvsPhi->Fill(phi, deltaEt / et);
0297   hDeltaEtvsDeltaR->Fill(deltaR, deltaEt);
0298   hDeltaEtOverEtvsDeltaR->Fill(deltaR, deltaEt / et);
0299
0300   hDeltaEta->Fill(deltaEta);
0301   hDeltaEtavsEt->Fill(et, deltaEta);
0302   hDeltaEtavsEta->Fill(eta, deltaEta);
0303
0304   hDeltaPhi->Fill(deltaPhi);
0305   hDeltaPhivsEt->Fill(et, deltaPhi);
0306   hDeltaPhivsEta->Fill(eta, deltaPhi);
0307
0308   hDeltaR->Fill(deltaR);
0309   hDeltaRvsEt->Fill(et, deltaR);
0310   hDeltaRvsEta->Fill(eta, deltaR);
0311
0312   BenchmarkTreeEntry entry;
0313   entry.deltaEt = deltaEt;
0314   entry.deltaEta = deltaEta;
0315   entry.et = et;
0316   entry.eta = eta;
0317
0318   if (doMetPlots_) {
0319     const reco::MET* met1 = static_cast<const reco::MET*>(genParticle);
0320     const reco::MET* met2 = static_cast<const reco::MET*>(recParticle);
0321     if (met1 != nullptr && met2 != nullptr) {
0322       //std::cout << "FL: met1.sumEt() = " << (*met1).sumEt() << std::endl;
0323       hTrueSumEt->Fill((*met1).sumEt());
0324       hSumEt->Fill((*met2).sumEt());
0325       hDeltaSetvsSet->Fill((*met1).sumEt(), (*met2).sumEt() - (*met1).sumEt());
0326       hDeltaMexvsSet->Fill((*met1).sumEt(), recParticle->px() - genParticle->px());
0327       hDeltaMexvsSet->Fill((*met1).sumEt(), recParticle->py() - genParticle->py());
0328       if ((*met1).sumEt() > 0.01)
0329         hDeltaSetOverSetvsSet->Fill((*met1).sumEt(), ((*met2).sumEt() - (*met1).sumEt()) / (*met1).sumEt());
0330       hRecSetvsTrueSet->Fill((*met1).sumEt(), (*met2).sumEt());
0331       hRecSetOverTrueSetvsTrueSet->Fill((*met1).sumEt(), (*met2).sumEt() / ((*met1).sumEt()));
0332       hTrueMexvsTrueSet->Fill((*met1).sumEt(), (*met1).px());
0333       hTrueMexvsTrueSet->Fill((*met1).sumEt(), (*met1).py());
0334     } else {
0335       std::cout << "Warning : reco::MET* == NULL" << std::endl;
0336     }
0337   }
0338
0339   //     tree_->Fill(entry);
0340 }
0341
0342 void GenericBenchmark::write(std::string Filename) {
0343   if (!Filename.empty() && file_)
0344     file_->Write(Filename.c_str());
0345 }
0346
0347 void GenericBenchmark::setfile(TFile* file) { file_ = file; }