GenericBenchmark.h

Go to the documentation of this file.

#ifndef RecoParticleFlow_Benchmark_GenericBenchmark_h
#define RecoParticleFlow_Benchmark_GenericBenchmark_h

//COLIN: necessary?
#include "RecoParticleFlow/Benchmark/interface/PFBenchmarkAlgo.h"

#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/METReco/interface/MET.h"
#include "DQMServices/Core/interface/DQMStore.h"

#include <string>

#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>

class BenchmarkTree;

//COLIN: this class REALLY needs to be cleaned up and rationalized, on the model of PFCandidateBenchmark.
class GenericBenchmark {
public:
  typedef dqm::legacy::DQMStore DQMStore;
  typedef dqm::legacy::MonitorElement MonitorElement;

  GenericBenchmark();
  virtual ~GenericBenchmark() noexcept(false);

  void setup(DQMStore *DQM = nullptr,
             bool PlotAgainstReco_ = true,
             float minDeltaEt = -100.,
             float maxDeltaEt = 50.,
             float minDeltaPhi = -0.5,
             float maxDeltaPhi = 0.5,
             bool doMetPlots = false);

  template <typename C>
  void fill(const C *RecoCollection,
            const C *GenCollection,
            bool startFromGen = false,
            bool PlotAgainstReco = true,
            bool onlyTwoJets = false,
            double recPt_cut = -1.,
            double minEta_cut = -1.,
            double maxEta_cut = -1.,
            double deltaR_cut = -1.);

  void write(std::string Filename);

  void setfile(TFile *file);

private:
  bool accepted(const reco::Candidate *particle, double ptCut, double minEtaCut, double maxEtaCut) const;

  void fillHistos(const reco::Candidate *genParticle,
                  const reco::Candidate *recParticle,
                  double deltaR_cut,
                  bool plotAgainstReco);

  TFile *file_;

  TH1F *hDeltaEt;
  TH1F *hDeltaEx;
  TH1F *hDeltaEy;
  TH2F *hDeltaEtvsEt;
  TH2F *hDeltaEtOverEtvsEt;
  TH2F *hDeltaEtvsEta;
  TH2F *hDeltaEtOverEtvsEta;
  TH2F *hDeltaEtvsPhi;
  TH2F *hDeltaEtOverEtvsPhi;
  TH2F *hDeltaEtvsDeltaR;
  TH2F *hDeltaEtOverEtvsDeltaR;

  TH2F *hEtRecvsEt;
  TH2F *hEtRecOverTrueEtvsTrueEt;

  TH1F *hDeltaEta;
  TH2F *hDeltaEtavsEt;
  TH2F *hDeltaEtavsEta;

  TH1F *hDeltaPhi;
  TH2F *hDeltaPhivsEt;
  TH2F *hDeltaPhivsEta;

  TH1F *hDeltaR;
  TH2F *hDeltaRvsEt;
  TH2F *hDeltaRvsEta;

  TH1F *hNRec;

  TH1F *hEtGen;
  TH1F *hEtaGen;
  TH2F *hEtvsEtaGen;
  TH2F *hEtvsPhiGen;
  TH1F *hPhiGen;

  TH1F *hNGen;

  TH1F *hEtSeen;
  TH1F *hEtaSeen;
  TH1F *hPhiSeen;
  TH2F *hEtvsEtaSeen;
  TH2F *hEtvsPhiSeen;

  TH1F *hEtRec;
  TH1F *hExRec;
  TH1F *hEyRec;
  TH1F *hPhiRec;

  TH1F *hSumEt;
  TH1F *hTrueSumEt;
  TH2F *hDeltaSetvsSet;
  TH2F *hDeltaMexvsSet;
  TH2F *hDeltaSetOverSetvsSet;
  TH2F *hRecSetvsTrueSet;
  TH2F *hRecSetOverTrueSetvsTrueSet;
  TH2F *hTrueMexvsTrueSet;

  BenchmarkTree *tree_;

  bool doMetPlots_;

protected:
  DQMStore *dbe_;
  PFBenchmarkAlgo *algo_;
};

template <typename C>
void GenericBenchmark::fill(const C *RecoCollection,
                            const C *GenCollection,
                            bool startFromGen,
                            bool PlotAgainstReco,
                            bool onlyTwoJets,
                            double recPt_cut,
                            double minEta_cut,
                            double maxEta_cut,
                            double deltaR_cut) {
  //if (doMetPlots_)
  //{
  //  const reco::MET* met1=static_cast<const reco::MET*>(&((*RecoCollection)[0]));
  //  if (met1!=NULL) std::cout << "FL: met1.sumEt() = " << (*met1).sumEt() << std::endl;
  //}

  // loop over reco particles

  if (!startFromGen) {
    int nRec = 0;
    for (unsigned int i = 0; i < RecoCollection->size(); i++) {
      // generate histograms comparing the reco and truth candidate (truth = closest in delta-R)
      const reco::Candidate *particle = &(*RecoCollection)[i];

      assert(particle != nullptr);
      if (!accepted(particle, recPt_cut, minEta_cut, maxEta_cut))
        continue;

      // Count the number of jets with a larger energy
      if (onlyTwoJets) {
        unsigned highJets = 0;
        for (unsigned j = 0; j < RecoCollection->size(); j++) {
          const reco::Candidate *otherParticle = &(*RecoCollection)[j];
          if (j != i && otherParticle->pt() > particle->pt())
            highJets++;
        }
        if (highJets > 1)
          continue;
      }
      nRec++;

      const reco::Candidate *gen_particle = algo_->matchByDeltaR(particle, GenCollection);
      if (gen_particle == nullptr)
        continue;

      // fill histograms
      fillHistos(gen_particle, particle, deltaR_cut, PlotAgainstReco);
    }

    hNRec->Fill(nRec);
  }

  // loop over gen particles

  //   std::cout<<"Reco size = "<<RecoCollection->size()<<", ";
  //   std::cout<<"Gen size = "<<GenCollection->size()<<std::endl;

  int nGen = 0;
  for (unsigned int i = 0; i < GenCollection->size(); i++) {
    const reco::Candidate *gen_particle = &(*GenCollection)[i];

    if (!accepted(gen_particle, recPt_cut, minEta_cut, maxEta_cut)) {
      continue;
    }

    hEtaGen->Fill(gen_particle->eta());
    hPhiGen->Fill(gen_particle->phi());
    hEtGen->Fill(gen_particle->et());
    hEtvsEtaGen->Fill(gen_particle->eta(), gen_particle->et());
    hEtvsPhiGen->Fill(gen_particle->phi(), gen_particle->et());

    const reco::Candidate *rec_particle = algo_->matchByDeltaR(gen_particle, RecoCollection);
    nGen++;
    if (!rec_particle)
      continue;  // no match
    // must make a cut on delta R - so let's do the cut

    double deltaR = algo_->deltaR(rec_particle, gen_particle);
    if (deltaR > deltaR_cut && deltaR_cut != -1.)
      continue;

    hEtaSeen->Fill(gen_particle->eta());
    hPhiSeen->Fill(gen_particle->phi());
    hEtSeen->Fill(gen_particle->et());
    hEtvsEtaSeen->Fill(gen_particle->eta(), gen_particle->et());
    hEtvsPhiSeen->Fill(gen_particle->phi(), gen_particle->et());

    hPhiRec->Fill(rec_particle->phi());
    hEtRec->Fill(rec_particle->et());
    hExRec->Fill(rec_particle->px());
    hEyRec->Fill(rec_particle->py());

    hEtRecvsEt->Fill(gen_particle->et(), rec_particle->et());
    if (gen_particle->et() != 0.0)
      hEtRecOverTrueEtvsTrueEt->Fill(gen_particle->et(), rec_particle->et() / gen_particle->et());

    if (startFromGen)
      fillHistos(gen_particle, rec_particle, deltaR_cut, PlotAgainstReco);
  }
  hNGen->Fill(nGen);
}

#endif  // RecoParticleFlow_Benchmark_GenericBenchmark_h