PFCandidateChecker.cc

Go to the documentation of this file.

#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
#include "DataFormats/JetReco/interface/PFJet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/stream/EDAnalyzer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "DataFormats/JetReco/interface/PFJetCollection.h"

namespace edm {
  class EventSetup;
  class Run;
}  // namespace edm

class PFCandidateChecker : public edm::stream::EDAnalyzer<> {
public:
  explicit PFCandidateChecker(const edm::ParameterSet&);

  void analyze(const edm::Event&, const edm::EventSetup&) override;

private:
  void printJets(const reco::PFJetCollection& pfJetsReco, const reco::PFJetCollection& pfJetsReReco) const;

  void printMet(const reco::PFCandidateCollection& pfReco, const reco::PFCandidateCollection& pfReReco) const;

  void printElementsInBlocks(const reco::PFCandidate& cand, std::ostream& out = std::cout) const;

  edm::InputTag inputTagPFCandidatesReco_;
  edm::InputTag inputTagPFCandidatesReReco_;
  edm::InputTag inputTagPFJetsReco_;
  edm::InputTag inputTagPFJetsReReco_;

  double deltaEMax_;
  double deltaEtaMax_;
  double deltaPhiMax_;

  bool verbose_;

  bool printBlocks_;

  bool rankByPt_;

  unsigned entry_;

  static bool greaterPt(const reco::PFCandidate& a, const reco::PFCandidate& b) { return (a.pt() > b.pt()); }
};

DEFINE_FWK_MODULE(PFCandidateChecker);

using namespace std;
using namespace edm;
using namespace reco;

PFCandidateChecker::PFCandidateChecker(const edm::ParameterSet& iConfig) {
  inputTagPFCandidatesReco_ = iConfig.getParameter<InputTag>("pfCandidatesReco");

  inputTagPFCandidatesReReco_ = iConfig.getParameter<InputTag>("pfCandidatesReReco");

  inputTagPFJetsReco_ = iConfig.getParameter<InputTag>("pfJetsReco");

  inputTagPFJetsReReco_ = iConfig.getParameter<InputTag>("pfJetsReReco");

  deltaEMax_ = iConfig.getParameter<double>("deltaEMax");

  deltaEtaMax_ = iConfig.getParameter<double>("deltaEtaMax");

  deltaPhiMax_ = iConfig.getParameter<double>("deltaPhiMax");

  verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);

  printBlocks_ = iConfig.getUntrackedParameter<bool>("printBlocks", false);

  rankByPt_ = iConfig.getUntrackedParameter<bool>("rankByPt", false);

  entry_ = 0;

  LogDebug("PFCandidateChecker") << " input collections : " << inputTagPFCandidatesReco_ << " "
                                 << inputTagPFCandidatesReReco_;
}

void PFCandidateChecker::analyze(const Event& iEvent, const EventSetup& iSetup) {
  LogDebug("PFCandidateChecker") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run() << endl;

  // get PFCandidates

  Handle<PFCandidateCollection> pfCandidatesReco;
  iEvent.getByLabel(inputTagPFCandidatesReco_, pfCandidatesReco);

  Handle<PFCandidateCollection> pfCandidatesReReco;
  iEvent.getByLabel(inputTagPFCandidatesReReco_, pfCandidatesReReco);

  Handle<PFJetCollection> pfJetsReco;
  iEvent.getByLabel(inputTagPFJetsReco_, pfJetsReco);

  Handle<PFJetCollection> pfJetsReReco;
  iEvent.getByLabel(inputTagPFJetsReReco_, pfJetsReReco);

  reco::PFCandidateCollection pfReco, pfReReco;

  // to sort, one needs to copy
  if (rankByPt_) {
    pfReco = *pfCandidatesReco;
    pfReReco = *pfCandidatesReReco;
    sort(pfReco.begin(), pfReco.end(), greaterPt);
    sort(pfReReco.begin(), pfReReco.end(), greaterPt);
  }

  unsigned minSize = pfReco.size() < pfReReco.size() ? pfReco.size() : pfReReco.size();
  bool differentCand = false;
  bool differentSize = pfReco.size() != pfReReco.size();
  if (differentSize)
    std::cout << "+++WARNING+++ PFCandidate size changed for entry " << entry_ << " !" << endl
              << " - RECO    size : " << pfReco.size() << endl
              << " - Re-RECO size : " << pfReReco.size() << endl;

  unsigned npr = 0;
  for (unsigned i = 0; i < minSize; i++) {
    const reco::PFCandidate& candReco = (rankByPt_) ? pfReco[i] : (*pfCandidatesReco)[i];
    const reco::PFCandidate& candReReco = (rankByPt_) ? pfReReco[i] : (*pfCandidatesReReco)[i];

    double deltaE = (candReReco.energy() - candReco.energy()) / (candReReco.energy() + candReco.energy());
    double deltaEta = candReReco.eta() - candReco.eta();
    double deltaPhi = candReReco.phi() - candReco.phi();
    if (fabs(deltaE) > deltaEMax_ || fabs(deltaEta) > deltaEtaMax_ || fabs(deltaPhi) > deltaPhiMax_) {
      differentCand = true;
      std::cout << "+++WARNING+++ PFCandidate " << i << " changed  for entry " << entry_ << " ! " << std::endl
                << " - RECO     : " << candReco << std::endl
                << " - Re-RECO  : " << candReReco << std::endl
                << " DeltaE   = : " << deltaE << std::endl
                << " DeltaEta = : " << deltaEta << std::endl
                << " DeltaPhi = : " << deltaPhi << std::endl
                << std::endl;
      if (printBlocks_) {
        std::cout << "Elements in Block for RECO: " << std::endl;
        printElementsInBlocks(candReco);
        std::cout << "Elements in Block for Re-RECO: " << std::endl;
        printElementsInBlocks(candReReco);
      }
      if (++npr == 5)
        break;
    }
  }

  if (differentSize || differentCand) {
    printJets(*pfJetsReco, *pfJetsReReco);
    printMet(pfReco, pfReReco);
  }

  ++entry_;
  LogDebug("PFCandidateChecker") << "STOP event: " << iEvent.id().event() << " in run " << iEvent.id().run()
                                 << std::endl;
}

void PFCandidateChecker::printMet(const PFCandidateCollection& pfReco, const PFCandidateCollection& pfReReco) const {
  double metX = 0.;
  double metY = 0.;
  for (unsigned i = 0; i < pfReco.size(); i++) {
    metX += pfReco[i].px();
    metY += pfReco[i].py();
  }
  double met = std::sqrt(metX * metX + metY * metY);
  std::cout << "MET RECO    = " << metX << " " << metY << " " << met << std::endl;

  metX = 0.;
  metY = 0.;
  for (unsigned i = 0; i < pfReReco.size(); i++) {
    metX += pfReReco[i].px();
    metY += pfReReco[i].py();
  }
  met = std::sqrt(metX * metX + metY * metY);
  std::cout << "MET Re-RECO = " << metX << " " << metY << " " << met << std::endl;
}

void PFCandidateChecker::printJets(const PFJetCollection& pfJetsReco, const PFJetCollection& pfJetsReReco) const {
  bool differentSize = pfJetsReco.size() != pfJetsReReco.size();
  if (differentSize)
    std::cout << "+++WARNING+++ PFJet size changed for entry " << entry_ << " !" << endl
              << " - RECO    size : " << pfJetsReco.size() << endl
              << " - Re-RECO size : " << pfJetsReReco.size() << endl;
  unsigned minSize = pfJetsReco.size() < pfJetsReReco.size() ? pfJetsReco.size() : pfJetsReReco.size();
  unsigned npr = 0;
  for (unsigned i = 0; i < minSize; ++i) {
    const reco::PFJet& candReco = pfJetsReco[i];
    const reco::PFJet& candReReco = pfJetsReReco[i];
    if (candReco.et() < 20. && candReReco.et() < 20.)
      break;
    double deltaE = (candReReco.et() - candReco.et()) / (candReReco.et() + candReco.et());
    double deltaEta = candReReco.eta() - candReco.eta();
    double deltaPhi = candReReco.phi() - candReco.phi();
    if (fabs(deltaE) > deltaEMax_ || fabs(deltaEta) > deltaEtaMax_ || fabs(deltaPhi) > deltaPhiMax_) {
      std::cout << "+++WARNING+++ PFJet " << i << " changed  for entry " << entry_ << " ! " << std::endl
                << " - RECO     : " << candReco.et() << " " << candReco.eta() << " " << candReco.phi() << std::endl
                << " - Re-RECO  : " << candReReco.et() << " " << candReReco.eta() << " " << candReReco.phi()
                << std::endl
                << " DeltaE   = : " << deltaE << std::endl
                << " DeltaEta = : " << deltaEta << std::endl
                << " DeltaPhi = : " << deltaPhi << std::endl
                << std::endl;
      if (++npr == 5)
        break;
    } else {
      std::cout << "Jet " << i << " " << candReco.et() << std::endl;
    }
  }
}

void PFCandidateChecker::printElementsInBlocks(const PFCandidate& cand, ostream& out) const {
  if (!out)
    return;

  PFBlockRef firstRef;

  assert(!cand.elementsInBlocks().empty());
  for (unsigned i = 0; i < cand.elementsInBlocks().size(); i++) {
    PFBlockRef blockRef = cand.elementsInBlocks()[i].first;

    if (blockRef.isNull()) {
      cerr << "ERROR! no block ref!";
      continue;
    }

    if (!i) {
      out << (*blockRef);
      firstRef = blockRef;
    } else if (blockRef != firstRef) {
      cerr << "WARNING! This PFCandidate is not made from a single block" << endl;
    }

    out << "\t" << cand.elementsInBlocks()[i].second << endl;
  }
}