Benchmark/interface/GenericBenchmark.h

0001 #ifndef RecoParticleFlow_Benchmark_GenericBenchmark_h
0002 #define RecoParticleFlow_Benchmark_GenericBenchmark_h
0003
0004 //COLIN: necessary?
0005 #include "RecoParticleFlow/Benchmark/interface/PFBenchmarkAlgo.h"
0006
0007 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0008 #include "DataFormats/Candidate/interface/Candidate.h"
0009 #include "DataFormats/METReco/interface/MET.h"
0010 #include "DQMServices/Core/interface/DQMStore.h"
0011
0012 #include <string>
0013
0014 #include <TH1F.h>
0015 #include <TH2F.h>
0016 #include <TFile.h>
0017
0018 class BenchmarkTree;
0019
0020 //COLIN: this class REALLY needs to be cleaned up and rationalized, on the model of PFCandidateBenchmark.
0021 class GenericBenchmark {
0022 public:
0023   typedef dqm::legacy::DQMStore DQMStore;
0024   typedef dqm::legacy::MonitorElement MonitorElement;
0025
0026   GenericBenchmark();
0027   virtual ~GenericBenchmark() noexcept(false);
0028
0029   void setup(DQMStore *DQM = nullptr,
0030              bool PlotAgainstReco_ = true,
0031              float minDeltaEt = -100.,
0032              float maxDeltaEt = 50.,
0033              float minDeltaPhi = -0.5,
0034              float maxDeltaPhi = 0.5,
0035              bool doMetPlots = false);
0036
0037   template <typename C>
0038   void fill(const C *RecoCollection,
0039             const C *GenCollection,
0040             bool startFromGen = false,
0041             bool PlotAgainstReco = true,
0042             bool onlyTwoJets = false,
0043             double recPt_cut = -1.,
0044             double minEta_cut = -1.,
0045             double maxEta_cut = -1.,
0046             double deltaR_cut = -1.);
0047
0048   void write(std::string Filename);
0049
0050   void setfile(TFile *file);
0051
0052 private:
0053   bool accepted(const reco::Candidate *particle, double ptCut, double minEtaCut, double maxEtaCut) const;
0054
0055   void fillHistos(const reco::Candidate *genParticle,
0056                   const reco::Candidate *recParticle,
0057                   double deltaR_cut,
0058                   bool plotAgainstReco);
0059
0060   TFile *file_;
0061
0062   TH1F *hDeltaEt;
0063   TH1F *hDeltaEx;
0064   TH1F *hDeltaEy;
0065   TH2F *hDeltaEtvsEt;
0066   TH2F *hDeltaEtOverEtvsEt;
0067   TH2F *hDeltaEtvsEta;
0068   TH2F *hDeltaEtOverEtvsEta;
0069   TH2F *hDeltaEtvsPhi;
0070   TH2F *hDeltaEtOverEtvsPhi;
0071   TH2F *hDeltaEtvsDeltaR;
0072   TH2F *hDeltaEtOverEtvsDeltaR;
0073
0074   TH2F *hEtRecvsEt;
0075   TH2F *hEtRecOverTrueEtvsTrueEt;
0076
0077   TH1F *hDeltaEta;
0078   TH2F *hDeltaEtavsEt;
0079   TH2F *hDeltaEtavsEta;
0080
0081   TH1F *hDeltaPhi;
0082   TH2F *hDeltaPhivsEt;
0083   TH2F *hDeltaPhivsEta;
0084
0085   TH1F *hDeltaR;
0086   TH2F *hDeltaRvsEt;
0087   TH2F *hDeltaRvsEta;
0088
0089   TH1F *hNRec;
0090
0091   TH1F *hEtGen;
0092   TH1F *hEtaGen;
0093   TH2F *hEtvsEtaGen;
0094   TH2F *hEtvsPhiGen;
0095   TH1F *hPhiGen;
0096
0097   TH1F *hNGen;
0098
0099   TH1F *hEtSeen;
0100   TH1F *hEtaSeen;
0101   TH1F *hPhiSeen;
0102   TH2F *hEtvsEtaSeen;
0103   TH2F *hEtvsPhiSeen;
0104
0105   TH1F *hEtRec;
0106   TH1F *hExRec;
0107   TH1F *hEyRec;
0108   TH1F *hPhiRec;
0109
0110   TH1F *hSumEt;
0111   TH1F *hTrueSumEt;
0112   TH2F *hDeltaSetvsSet;
0113   TH2F *hDeltaMexvsSet;
0114   TH2F *hDeltaSetOverSetvsSet;
0115   TH2F *hRecSetvsTrueSet;
0116   TH2F *hRecSetOverTrueSetvsTrueSet;
0117   TH2F *hTrueMexvsTrueSet;
0118
0119   BenchmarkTree *tree_;
0120
0121   bool doMetPlots_;
0122
0123 protected:
0124   DQMStore *dbe_;
0125   PFBenchmarkAlgo *algo_;
0126 };
0127
0128 template <typename C>
0129 void GenericBenchmark::fill(const C *RecoCollection,
0130                             const C *GenCollection,
0131                             bool startFromGen,
0132                             bool PlotAgainstReco,
0133                             bool onlyTwoJets,
0134                             double recPt_cut,
0135                             double minEta_cut,
0136                             double maxEta_cut,
0137                             double deltaR_cut) {
0138   //if (doMetPlots_)
0139   //{
0140   //  const reco::MET* met1=static_cast<const reco::MET*>(&((*RecoCollection)[0]));
0141   //  if (met1!=NULL) std::cout << "FL: met1.sumEt() = " << (*met1).sumEt() << std::endl;
0142   //}
0143
0144   // loop over reco particles
0145
0146   if (!startFromGen) {
0147     int nRec = 0;
0148     for (unsigned int i = 0; i < RecoCollection->size(); i++) {
0149       // generate histograms comparing the reco and truth candidate (truth = closest in delta-R)
0150       const reco::Candidate *particle = &(*RecoCollection)[i];
0151
0152       assert(particle != nullptr);
0153       if (!accepted(particle, recPt_cut, minEta_cut, maxEta_cut))
0154         continue;
0155
0156       // Count the number of jets with a larger energy
0157       if (onlyTwoJets) {
0158         unsigned highJets = 0;
0159         for (unsigned j = 0; j < RecoCollection->size(); j++) {
0160           const reco::Candidate *otherParticle = &(*RecoCollection)[j];
0161           if (j != i && otherParticle->pt() > particle->pt())
0162             highJets++;
0163         }
0164         if (highJets > 1)
0165           continue;
0166       }
0167       nRec++;
0168
0169       const reco::Candidate *gen_particle = algo_->matchByDeltaR(particle, GenCollection);
0170       if (gen_particle == nullptr)
0171         continue;
0172
0173       // fill histograms
0174       fillHistos(gen_particle, particle, deltaR_cut, PlotAgainstReco);
0175     }
0176
0177     hNRec->Fill(nRec);
0178   }
0179
0180   // loop over gen particles
0181
0182   //   std::cout<<"Reco size = "<<RecoCollection->size()<<", ";
0183   //   std::cout<<"Gen size = "<<GenCollection->size()<<std::endl;
0184
0185   int nGen = 0;
0186   for (unsigned int i = 0; i < GenCollection->size(); i++) {
0187     const reco::Candidate *gen_particle = &(*GenCollection)[i];
0188
0189     if (!accepted(gen_particle, recPt_cut, minEta_cut, maxEta_cut)) {
0190       continue;
0191     }
0192
0193     hEtaGen->Fill(gen_particle->eta());
0194     hPhiGen->Fill(gen_particle->phi());
0195     hEtGen->Fill(gen_particle->et());
0196     hEtvsEtaGen->Fill(gen_particle->eta(), gen_particle->et());
0197     hEtvsPhiGen->Fill(gen_particle->phi(), gen_particle->et());
0198
0199     const reco::Candidate *rec_particle = algo_->matchByDeltaR(gen_particle, RecoCollection);
0200     nGen++;
0201     if (!rec_particle)
0202       continue;  // no match
0203     // must make a cut on delta R - so let's do the cut
0204
0205     double deltaR = algo_->deltaR(rec_particle, gen_particle);
0206     if (deltaR > deltaR_cut && deltaR_cut != -1.)
0207       continue;
0208
0209     hEtaSeen->Fill(gen_particle->eta());
0210     hPhiSeen->Fill(gen_particle->phi());
0211     hEtSeen->Fill(gen_particle->et());
0212     hEtvsEtaSeen->Fill(gen_particle->eta(), gen_particle->et());
0213     hEtvsPhiSeen->Fill(gen_particle->phi(), gen_particle->et());
0214
0215     hPhiRec->Fill(rec_particle->phi());
0216     hEtRec->Fill(rec_particle->et());
0217     hExRec->Fill(rec_particle->px());
0218     hEyRec->Fill(rec_particle->py());
0219
0220     hEtRecvsEt->Fill(gen_particle->et(), rec_particle->et());
0221     if (gen_particle->et() != 0.0)
0222       hEtRecOverTrueEtvsTrueEt->Fill(gen_particle->et(), rec_particle->et() / gen_particle->et());
0223
0224     if (startFromGen)
0225       fillHistos(gen_particle, rec_particle, deltaR_cut, PlotAgainstReco);
0226   }
0227   hNGen->Fill(nGen);
0228 }
0229
0230 #endif  // RecoParticleFlow_Benchmark_GenericBenchmark_h