PFTauElecRejectionBenchmark.cc

Go to the documentation of this file.

#include "RecoParticleFlow/Benchmark/interface/PFTauElecRejectionBenchmark.h"

// preprocessor macro for booking 1d histos with DQMStore -or- bare Root
#define BOOK1D(name, title, nbinsx, lowx, highx) \
  h##name =                                      \
      db_ ? db_->book1D(#name, title, nbinsx, lowx, highx)->getTH1F() : new TH1F(#name, title, nbinsx, lowx, highx)

// preprocessor macro for booking 2d histos with DQMStore -or- bare Root
#define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)                            \
  h##name = db_ ? db_->book2D(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)->getTH2F() \
                : new TH2F(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)

// all versions OK
// preprocesor macro for setting axis titles
#define SETAXES(name, xtitle, ytitle)    \
  h##name->GetXaxis()->SetTitle(xtitle); \
  h##name->GetYaxis()->SetTitle(ytitle)

using namespace reco;
using namespace std;

PFTauElecRejectionBenchmark::PFTauElecRejectionBenchmark() : file_(nullptr) {}

PFTauElecRejectionBenchmark::~PFTauElecRejectionBenchmark() {
  if (file_)
    file_->Close();
}

void PFTauElecRejectionBenchmark::write() {
  // Store the DAQ Histograms
  if (!outputFile_.empty()) {
    if (db_)
      db_->save(outputFile_);
    // use bare Root if no DQM (FWLite applications)
    else if (file_) {
      file_->Write(outputFile_.c_str());
      cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
      file_->Close();
    }
  } else
    cout << "No output file specified (" << outputFile_ << "). Results will not be saved!" << endl;
}

void PFTauElecRejectionBenchmark::setup(string Filename,
                                        string benchmarkLabel,
                                        double maxDeltaR,
                                        double minRecoPt,
                                        double maxRecoAbsEta,
                                        double minMCPt,
                                        double maxMCAbsEta,
                                        string sGenMatchObjectLabel,
                                        bool applyEcalCrackCut,
                                        DQMStore* db_store) {
  maxDeltaR_ = maxDeltaR;
  benchmarkLabel_ = benchmarkLabel;
  outputFile_ = Filename;
  minRecoPt_ = minRecoPt;
  maxRecoAbsEta_ = maxRecoAbsEta;
  minMCPt_ = minMCPt;
  maxMCAbsEta_ = maxMCAbsEta;
  sGenMatchObjectLabel_ = sGenMatchObjectLabel;
  applyEcalCrackCut_ = applyEcalCrackCut;

  file_ = nullptr;
  db_ = db_store;

  // print parameters
  cout << "PFTauElecRejectionBenchmark Setup parameters ==============================================" << endl;
  cout << "Filename to write histograms " << Filename << endl;
  cout << "Benchmark label name " << benchmarkLabel_ << endl;
  cout << "maxDeltaRMax " << maxDeltaR << endl;
  cout << "minRecoPt " << minRecoPt_ << endl;
  cout << "maxRecoAbsEta " << maxRecoAbsEta_ << endl;
  cout << "minMCPt " << minMCPt_ << endl;
  cout << "maxMCAbsEta " << maxMCAbsEta_ << endl;
  cout << "applyEcalCrackCut " << applyEcalCrackCut_ << endl;
  cout << "sGenMatchObjectLabel " << sGenMatchObjectLabel_ << endl;

  // Book histogram

  // Establish DQM Store
  string path = "PFTask/Benchmarks/" + benchmarkLabel_ + "/";
  path += "Gen";
  if (db_) {
    db_->setCurrentFolder(path);
  } else {
    file_ = new TFile(outputFile_.c_str(), "recreate");
    cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. " << endl;
  }

  // E/p
  BOOK1D(EoverP, "E/p", 100, 0., 4.);
  SETAXES(EoverP, "E/p", "Entries");

  BOOK1D(EoverP_barrel, "E/p barrel", 100, 0., 4.);
  SETAXES(EoverP_barrel, "E/p barrel", "Entries");

  BOOK1D(EoverP_endcap, "E/p endcap", 100, 0., 4.);
  SETAXES(EoverP_endcap, "E/p endcap", "Entries");

  BOOK1D(EoverP_preid0, "E/p (preid=0)", 100, 0., 4.);
  SETAXES(EoverP_preid0, "E/p", "Entries");

  BOOK1D(EoverP_preid1, "E/p (preid=1)", 100, 0., 4.);
  SETAXES(EoverP_preid1, "E/p", "Entries");

  // H/p
  BOOK1D(HoverP, "H/p", 100, 0., 2.);
  SETAXES(HoverP, "H/p", "Entries");

  BOOK1D(HoverP_barrel, "H/p barrel", 100, 0., 2.);
  SETAXES(HoverP_barrel, "H/p barrel", "Entries");

  BOOK1D(HoverP_endcap, "H/p endcap", 100, 0., 2.);
  SETAXES(HoverP_endcap, "H/p endcap", "Entries");

  BOOK1D(HoverP_preid0, "H/p (preid=0)", 100, 0., 2.);
  SETAXES(HoverP_preid0, "H/p", "Entries");

  BOOK1D(HoverP_preid1, "H/p (preid=1)", 100, 0., 2.);
  SETAXES(HoverP_preid1, "H/p", "Entries");

  // emfrac
  BOOK1D(Emfrac, "EM fraction", 100, 0., 1.01);
  SETAXES(Emfrac, "em fraction", "Entries");

  BOOK1D(Emfrac_barrel, "EM fraction barrel", 100, 0., 1.01);
  SETAXES(Emfrac_barrel, "em fraction barrel", "Entries");

  BOOK1D(Emfrac_endcap, "EM fraction endcap", 100, 0., 1.01);
  SETAXES(Emfrac_endcap, "em fraction endcap", "Entries");

  BOOK1D(Emfrac_preid0, "EM fraction (preid=0)", 100, 0., 1.01);
  SETAXES(Emfrac_preid0, "em fraction", "Entries");

  BOOK1D(Emfrac_preid1, "EM fraction (preid=1)", 100, 0., 1.01);
  SETAXES(Emfrac_preid1, "em fraction", "Entries");

  // PreID
  BOOK1D(ElecPreID, "PFElectron PreID decision", 6, 0., 1.01);
  SETAXES(ElecPreID, "PFElectron PreID decision", "Entries");

  // MVA
  BOOK1D(ElecMVA, "PFElectron MVA", 100, -1.01, 1.01);
  SETAXES(ElecMVA, "PFElectron MVA", "Entries");

  // Discriminant
  BOOK1D(TauElecDiscriminant, "PFTau-Electron Discriminant", 6, 0., 1.01);
  SETAXES(TauElecDiscriminant, "PFTau-Electron Discriminant", "Entries");

  // PFCand clusters
  BOOK1D(pfcand_deltaEta, "PFCand cluster dEta", 100, 0., 0.8);
  SETAXES(pfcand_deltaEta, "PFCand cluster #Delta(#eta)", "Entries");

  BOOK1D(pfcand_deltaEta_weightE, "PFCand cluster dEta, energy weighted", 100, 0., 0.8);
  SETAXES(pfcand_deltaEta_weightE, "PFCand cluster #Delta(#eta)", "Entries");

  BOOK1D(pfcand_deltaPhiOverQ, "PFCand cluster dPhi/q", 100, -0.8, 0.8);
  SETAXES(pfcand_deltaPhiOverQ, "PFCand cluster #Delta(#phi)/q", "Entries");

  BOOK1D(pfcand_deltaPhiOverQ_weightE, "PFCand cluster dEta/q, energy weighted", 100, -0.8, 0.8);
  SETAXES(pfcand_deltaPhiOverQ_weightE, "PFCand cluster #Delta(#phi)/q", "Entries");

  // Leading KF track
  BOOK1D(leadTk_pt, "leading KF track pt", 100, 0., 80.);
  SETAXES(leadTk_pt, "leading KF track p_{T} (GeV)", "Entries");

  BOOK1D(leadTk_eta, "leading KF track eta", 100, -4., 4.);
  SETAXES(leadTk_eta, "leading KF track #eta", "Entries");

  BOOK1D(leadTk_phi, "leading KF track phi", 100, -3.2, 3.2);
  SETAXES(leadTk_phi, "leading KF track #phi", "Entries");

  // Leading Gsf track
  BOOK1D(leadGsfTk_pt, "leading Gsf track pt", 100, 0., 80.);
  SETAXES(leadGsfTk_pt, "leading Gsf track p_{T} (GeV)", "Entries");

  BOOK1D(leadGsfTk_eta, "leading Gsf track eta", 100, -4., 4.);
  SETAXES(leadGsfTk_eta, "leading Gsf track #eta", "Entries");

  BOOK1D(leadGsfTk_phi, "leading Gsf track phi", 100, -3.2, 3.2);
  SETAXES(leadGsfTk_phi, "leading Gsf track #phi", "Entries");

  // H/p vs E/p
  BOOK2D(HoPvsEoP, "H/p vs. E/p", 100, 0., 2., 100, 0., 2.);
  SETAXES(HoPvsEoP, "E/p", "H/p");

  BOOK2D(HoPvsEoP_preid0, "H/p vs. E/p (preid=0)", 100, 0., 2., 100, 0., 2.);
  SETAXES(HoPvsEoP_preid0, "E/p", "H/p");

  BOOK2D(HoPvsEoP_preid1, "H/p vs. E/p (preid=0)", 100, 0., 2., 100, 0., 2.);
  SETAXES(HoPvsEoP_preid1, "E/p", "H/p");

  // em fraction vs E/p
  BOOK2D(EmfracvsEoP, "emfrac vs. E/p", 100, 0., 2., 100, 0., 1.01);
  SETAXES(EmfracvsEoP, "E/p", "em fraction");

  BOOK2D(EmfracvsEoP_preid0, "emfrac vs. E/p (preid=0)", 100, 0., 2., 100, 0., 1.01);
  SETAXES(EmfracvsEoP_preid0, "E/p", "em fraction");

  BOOK2D(EmfracvsEoP_preid1, "emfrac vs. E/p (preid=0)", 100, 0., 2., 100, 0., 1.01);
  SETAXES(EmfracvsEoP_preid1, "E/p", "em fraction");
}

void PFTauElecRejectionBenchmark::process(edm::Handle<edm::HepMCProduct> mcevt,
                                          edm::Handle<reco::PFTauCollection> pfTaus,
                                          edm::Handle<reco::PFTauDiscriminator> pfTauIsoDiscr,
                                          edm::Handle<reco::PFTauDiscriminator> pfTauElecDiscr) {
  // Find Gen Objects to be matched with
  HepMC::GenEvent* generated_event = new HepMC::GenEvent(*(mcevt->GetEvent()));
  _GenObjects.clear();

  TLorentzVector taunet;
  HepMC::GenEvent::particle_iterator p;
  for (p = generated_event->particles_begin(); p != generated_event->particles_end(); p++) {
    if (std::abs((*p)->pdg_id()) == 15 && (*p)->status() == 2) {
      bool lept_decay = false;
      TLorentzVector tau((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
      HepMC::GenVertex::particle_iterator z = (*p)->end_vertex()->particles_begin(HepMC::descendants);
      for (; z != (*p)->end_vertex()->particles_end(HepMC::descendants); z++) {
        if (std::abs((*z)->pdg_id()) == 11 || std::abs((*z)->pdg_id()) == 13)
          lept_decay = true;
        if (std::abs((*z)->pdg_id()) == 16)
          taunet.SetPxPyPzE((*z)->momentum().px(), (*z)->momentum().py(), (*z)->momentum().pz(), (*z)->momentum().e());
      }
      if (lept_decay == false) {
        TLorentzVector jetMom = tau - taunet;
        if (sGenMatchObjectLabel_ == "tau")
          _GenObjects.push_back(jetMom);
      }
    } else if (std::abs((*p)->pdg_id()) == 11 && (*p)->status() == 1) {
      TLorentzVector elec((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
      if (sGenMatchObjectLabel_ == "e")
        _GenObjects.push_back(elec);
    }
  }


  // Loop over all PFTaus
  math::XYZPointF myleadTkEcalPos;
  for (PFTauCollection::size_type iPFTau = 0; iPFTau < pfTaus->size(); iPFTau++) {
    PFTauRef thePFTau(pfTaus, iPFTau);
    if ((*pfTauIsoDiscr)[thePFTau] == 1) {
      if ((*thePFTau).et() > minRecoPt_ && std::abs((*thePFTau).eta()) < maxRecoAbsEta_) {
        // Check if track goes to Ecal crack
        reco::TrackRef myleadTk;
        const reco::PFCandidatePtr& pflch = thePFTau->leadPFChargedHadrCand();
        if (pflch.isNonnull()) {
          myleadTk = pflch->trackRef();
          myleadTkEcalPos = pflch->positionAtECALEntrance();

          if (myleadTk.isNonnull()) {
            if (applyEcalCrackCut_ && isInEcalCrack(std::abs((double)myleadTkEcalPos.eta()))) {
              continue;  // do nothing
            } else {
              // Match with gen object
              for (unsigned int i = 0; i < _GenObjects.size(); i++) {
                if (_GenObjects[i].Et() >= minMCPt_ && std::abs(_GenObjects[i].Eta()) < maxMCAbsEta_) {
                  TLorentzVector pftau((*thePFTau).px(), (*thePFTau).py(), (*thePFTau).pz(), (*thePFTau).energy());
                  double GenDeltaR = pftau.DeltaR(_GenObjects[i]);
                  if (GenDeltaR < maxDeltaR_) {
                    hleadTk_pt->Fill((float)myleadTk->pt());
                    hleadTk_eta->Fill((float)myleadTk->eta());
                    hleadTk_phi->Fill((float)myleadTk->phi());

                    hEoverP->Fill((*thePFTau).ecalStripSumEOverPLead());
                    hHoverP->Fill((*thePFTau).hcal3x3OverPLead());
                    hEmfrac->Fill((*thePFTau).emFraction());

                    if (std::abs(myleadTk->eta()) < 1.5) {
                      hEoverP_barrel->Fill((*thePFTau).ecalStripSumEOverPLead());
                      hHoverP_barrel->Fill((*thePFTau).hcal3x3OverPLead());
                      hEmfrac_barrel->Fill((*thePFTau).emFraction());
                    } else if (std::abs(myleadTk->eta()) > 1.5 && std::abs(myleadTk->eta()) < 2.5) {
                      hEoverP_endcap->Fill((*thePFTau).ecalStripSumEOverPLead());
                      hHoverP_endcap->Fill((*thePFTau).hcal3x3OverPLead());
                      hEmfrac_endcap->Fill((*thePFTau).emFraction());
                    }

                    // if -999 fill in -1 bin!
                    if ((*thePFTau).electronPreIDOutput() < -1)
                      hElecMVA->Fill(-1);
                    else
                      hElecMVA->Fill((*thePFTau).electronPreIDOutput());

                    hTauElecDiscriminant->Fill((*pfTauElecDiscr)[thePFTau]);

                    hHoPvsEoP->Fill((*thePFTau).ecalStripSumEOverPLead(), (*thePFTau).hcal3x3OverPLead());
                    hEmfracvsEoP->Fill((*thePFTau).emFraction(), (*thePFTau).hcal3x3OverPLead());

                    if ((*thePFTau).electronPreIDDecision() == 1) {
                      hEoverP_preid1->Fill((*thePFTau).ecalStripSumEOverPLead());
                      hHoverP_preid1->Fill((*thePFTau).hcal3x3OverPLead());
                      hEmfrac_preid1->Fill((*thePFTau).emFraction());
                      hHoPvsEoP_preid1->Fill((*thePFTau).ecalStripSumEOverPLead(), (*thePFTau).hcal3x3OverPLead());
                      hEmfracvsEoP_preid1->Fill((*thePFTau).emFraction(), (*thePFTau).hcal3x3OverPLead());
                    } else {
                      hEoverP_preid0->Fill((*thePFTau).ecalStripSumEOverPLead());
                      hHoverP_preid0->Fill((*thePFTau).hcal3x3OverPLead());
                      hEmfrac_preid0->Fill((*thePFTau).emFraction());
                      hHoPvsEoP_preid0->Fill((*thePFTau).ecalStripSumEOverPLead(), (*thePFTau).hcal3x3OverPLead());
                      hEmfracvsEoP_preid0->Fill((*thePFTau).emFraction(), (*thePFTau).hcal3x3OverPLead());
                    }
                  }
                }

                // Loop over all PFCands for cluster plots
                std::vector<CandidatePtr> myPFCands = (*thePFTau).pfTauTagInfoRef()->PFCands();
                for (int i = 0; i < (int)myPFCands.size(); i++) {
                  const reco::PFCandidate* pfCand = dynamic_cast<const reco::PFCandidate*>(myPFCands[i].get());
                  if (pfCand == nullptr)
                    continue;

                  math::XYZPointF candPos;
                  if (std::abs(pfCand->pdgId()) == 211 ||
                      std::abs(pfCand->pdgId()) == 11)  // if charged hadron or electron
                    candPos = pfCand->positionAtECALEntrance();
                  else
                    candPos = math::XYZPointF(pfCand->px(), pfCand->py(), pfCand->pz());

                  //double deltaR   = ROOT::Math::VectorUtil::DeltaR(myleadTkEcalPos,candPos);
                  double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myleadTkEcalPos, candPos);
                  double deltaEta = std::abs(myleadTkEcalPos.eta() - myPFCands[i]->eta());
                  double deltaPhiOverQ = deltaPhi / (double)myleadTk->charge();

                  hpfcand_deltaEta->Fill(deltaEta);
                  hpfcand_deltaEta_weightE->Fill(deltaEta * pfCand->ecalEnergy());
                  hpfcand_deltaPhiOverQ->Fill(deltaPhiOverQ);
                  hpfcand_deltaPhiOverQ_weightE->Fill(deltaPhiOverQ * pfCand->ecalEnergy());
                }
              }
            }
          }
        }
      }
    }
  }
}

// Ecal crack  map from Egamma POG
bool PFTauElecRejectionBenchmark::isInEcalCrack(double eta) const {
  return (eta < 0.018 || (eta > 0.423 && eta < 0.461) || (eta > 0.770 && eta < 0.806) || (eta > 1.127 && eta < 1.163) ||
          (eta > 1.460 && eta < 1.558));
}