PFProducer/plugins/PFCandidateChecker.cc

0001 /**\class PFCandidateChecker
0002 \brief Checks what a re-reco changes in PFCandidates.
0003
0004 \author Patrick Janot
0005 \date   August 2011
0006 */
0007
0008 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0009 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0010 #include "DataFormats/JetReco/interface/PFJet.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/Framework/interface/stream/EDAnalyzer.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0017
0018 namespace edm {
0019   class EventSetup;
0020   class Run;
0021 }  // namespace edm
0022
0023 class PFCandidateChecker : public edm::stream::EDAnalyzer<> {
0024 public:
0025   explicit PFCandidateChecker(const edm::ParameterSet&);
0026
0027   void analyze(const edm::Event&, const edm::EventSetup&) override;
0028
0029 private:
0030   void printJets(const reco::PFJetCollection& pfJetsReco, const reco::PFJetCollection& pfJetsReReco) const;
0031
0032   void printMet(const reco::PFCandidateCollection& pfReco, const reco::PFCandidateCollection& pfReReco) const;
0033
0034   void printElementsInBlocks(const reco::PFCandidate& cand, std::ostream& out = std::cout) const;
0035
0036   /// PFCandidates in which we'll look for pile up particles
0037   edm::EDGetTokenT<reco::PFCandidateCollection> inputTokenPFCandidatesReco_;
0038   edm::EDGetTokenT<reco::PFCandidateCollection> inputTokenPFCandidatesReReco_;
0039   edm::EDGetTokenT<reco::PFJetCollection> inputTokenPFJetsReco_;
0040   edm::EDGetTokenT<reco::PFJetCollection> inputTokenPFJetsReReco_;
0041
0042   /// Cuts for comparison
0043   double deltaEMax_;
0044   double deltaEtaMax_;
0045   double deltaPhiMax_;
0046
0047   /// verbose ?
0048   bool verbose_;
0049
0050   /// print the blocks associated to a given candidate ?
0051   bool printBlocks_;
0052
0053   /// rank the candidates by Pt
0054   bool rankByPt_;
0055
0056   /// Counter
0057   unsigned entry_;
0058
0059   static bool greaterPt(const reco::PFCandidate& a, const reco::PFCandidate& b) { return (a.pt() > b.pt()); }
0060 };
0061
0062 DEFINE_FWK_MODULE(PFCandidateChecker);
0063
0064 using namespace std;
0065 using namespace edm;
0066 using namespace reco;
0067
0068 PFCandidateChecker::PFCandidateChecker(const edm::ParameterSet& iConfig) {
0069   inputTokenPFCandidatesReco_ =
0070       consumes<reco::PFCandidateCollection>(iConfig.getParameter<InputTag>("pfCandidatesReco"));
0071
0072   inputTokenPFCandidatesReReco_ =
0073       consumes<reco::PFCandidateCollection>(iConfig.getParameter<InputTag>("pfCandidatesReReco"));
0074
0075   inputTokenPFJetsReco_ = consumes<reco::PFJetCollection>(iConfig.getParameter<InputTag>("pfJetsReco"));
0076
0077   inputTokenPFJetsReReco_ = consumes<reco::PFJetCollection>(iConfig.getParameter<InputTag>("pfJetsReReco"));
0078
0079   deltaEMax_ = iConfig.getParameter<double>("deltaEMax");
0080
0081   deltaEtaMax_ = iConfig.getParameter<double>("deltaEtaMax");
0082
0083   deltaPhiMax_ = iConfig.getParameter<double>("deltaPhiMax");
0084
0085   verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);
0086
0087   printBlocks_ = iConfig.getUntrackedParameter<bool>("printBlocks", false);
0088
0089   rankByPt_ = iConfig.getUntrackedParameter<bool>("rankByPt", false);
0090
0091   entry_ = 0;
0092
0093   LogDebug("PFCandidateChecker") << " input collections : " << iConfig.getParameter<InputTag>("pfCandidatesReco") << " "
0094                                  << iConfig.getParameter<InputTag>("pfCandidatesReReco");
0095 }
0096
0097 void PFCandidateChecker::analyze(const Event& iEvent, const EventSetup& iSetup) {
0098   LogDebug("PFCandidateChecker") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run() << endl;
0099
0100   // get PFCandidates
0101
0102   const auto& pfCandidatesReco = iEvent.getHandle(inputTokenPFCandidatesReco_);
0103   const auto& pfCandidatesReReco = iEvent.getHandle(inputTokenPFCandidatesReReco_);
0104   const auto& pfJetsReco = iEvent.getHandle(inputTokenPFJetsReco_);
0105   const auto& pfJetsReReco = iEvent.getHandle(inputTokenPFJetsReReco_);
0106
0107   reco::PFCandidateCollection pfReco, pfReReco;
0108
0109   // to sort, one needs to copy
0110   if (rankByPt_) {
0111     pfReco = *pfCandidatesReco;
0112     pfReReco = *pfCandidatesReReco;
0113     sort(pfReco.begin(), pfReco.end(), greaterPt);
0114     sort(pfReReco.begin(), pfReReco.end(), greaterPt);
0115   }
0116
0117   unsigned minSize = pfReco.size() < pfReReco.size() ? pfReco.size() : pfReReco.size();
0118   bool differentCand = false;
0119   bool differentSize = pfReco.size() != pfReReco.size();
0120   if (differentSize)
0121     std::cout << "+++WARNING+++ PFCandidate size changed for entry " << entry_ << " !" << endl
0122               << " - RECO    size : " << pfReco.size() << endl
0123               << " - Re-RECO size : " << pfReReco.size() << endl;
0124
0125   unsigned npr = 0;
0126   for (unsigned i = 0; i < minSize; i++) {
0127     const reco::PFCandidate& candReco = (rankByPt_) ? pfReco[i] : (*pfCandidatesReco)[i];
0128     const reco::PFCandidate& candReReco = (rankByPt_) ? pfReReco[i] : (*pfCandidatesReReco)[i];
0129
0130     double deltaE = (candReReco.energy() - candReco.energy()) / (candReReco.energy() + candReco.energy());
0131     double deltaEta = candReReco.eta() - candReco.eta();
0132     double deltaPhi = candReReco.phi() - candReco.phi();
0133     if (fabs(deltaE) > deltaEMax_ || fabs(deltaEta) > deltaEtaMax_ || fabs(deltaPhi) > deltaPhiMax_) {
0134       differentCand = true;
0135       std::cout << "+++WARNING+++ PFCandidate " << i << " changed  for entry " << entry_ << " ! " << std::endl
0136                 << " - RECO     : " << candReco << std::endl
0137                 << " - Re-RECO  : " << candReReco << std::endl
0138                 << " DeltaE   = : " << deltaE << std::endl
0139                 << " DeltaEta = : " << deltaEta << std::endl
0140                 << " DeltaPhi = : " << deltaPhi << std::endl
0141                 << std::endl;
0142       if (printBlocks_) {
0143         std::cout << "Elements in Block for RECO: " << std::endl;
0144         printElementsInBlocks(candReco);
0145         std::cout << "Elements in Block for Re-RECO: " << std::endl;
0146         printElementsInBlocks(candReReco);
0147       }
0148       if (++npr == 5)
0149         break;
0150     }
0151   }
0152
0153   if (differentSize || differentCand) {
0154     printJets(*pfJetsReco, *pfJetsReReco);
0155     printMet(pfReco, pfReReco);
0156   }
0157
0158   ++entry_;
0159   LogDebug("PFCandidateChecker") << "STOP event: " << iEvent.id().event() << " in run " << iEvent.id().run()
0160                                  << std::endl;
0161 }
0162
0163 void PFCandidateChecker::printMet(const PFCandidateCollection& pfReco, const PFCandidateCollection& pfReReco) const {
0164   double metX = 0.;
0165   double metY = 0.;
0166   for (unsigned i = 0; i < pfReco.size(); i++) {
0167     metX += pfReco[i].px();
0168     metY += pfReco[i].py();
0169   }
0170   double met = std::sqrt(metX * metX + metY * metY);
0171   std::cout << "MET RECO    = " << metX << " " << metY << " " << met << std::endl;
0172
0173   metX = 0.;
0174   metY = 0.;
0175   for (unsigned i = 0; i < pfReReco.size(); i++) {
0176     metX += pfReReco[i].px();
0177     metY += pfReReco[i].py();
0178   }
0179   met = std::sqrt(metX * metX + metY * metY);
0180   std::cout << "MET Re-RECO = " << metX << " " << metY << " " << met << std::endl;
0181 }
0182
0183 void PFCandidateChecker::printJets(const PFJetCollection& pfJetsReco, const PFJetCollection& pfJetsReReco) const {
0184   bool differentSize = pfJetsReco.size() != pfJetsReReco.size();
0185   if (differentSize)
0186     std::cout << "+++WARNING+++ PFJet size changed for entry " << entry_ << " !" << endl
0187               << " - RECO    size : " << pfJetsReco.size() << endl
0188               << " - Re-RECO size : " << pfJetsReReco.size() << endl;
0189   unsigned minSize = pfJetsReco.size() < pfJetsReReco.size() ? pfJetsReco.size() : pfJetsReReco.size();
0190   unsigned npr = 0;
0191   for (unsigned i = 0; i < minSize; ++i) {
0192     const reco::PFJet& candReco = pfJetsReco[i];
0193     const reco::PFJet& candReReco = pfJetsReReco[i];
0194     if (candReco.et() < 20. && candReReco.et() < 20.)
0195       break;
0196     double deltaE = (candReReco.et() - candReco.et()) / (candReReco.et() + candReco.et());
0197     double deltaEta = candReReco.eta() - candReco.eta();
0198     double deltaPhi = candReReco.phi() - candReco.phi();
0199     if (fabs(deltaE) > deltaEMax_ || fabs(deltaEta) > deltaEtaMax_ || fabs(deltaPhi) > deltaPhiMax_) {
0200       std::cout << "+++WARNING+++ PFJet " << i << " changed  for entry " << entry_ << " ! " << std::endl
0201                 << " - RECO     : " << candReco.et() << " " << candReco.eta() << " " << candReco.phi() << std::endl
0202                 << " - Re-RECO  : " << candReReco.et() << " " << candReReco.eta() << " " << candReReco.phi()
0203                 << std::endl
0204                 << " DeltaE   = : " << deltaE << std::endl
0205                 << " DeltaEta = : " << deltaEta << std::endl
0206                 << " DeltaPhi = : " << deltaPhi << std::endl
0207                 << std::endl;
0208       if (++npr == 5)
0209         break;
0210     } else {
0211       std::cout << "Jet " << i << " " << candReco.et() << std::endl;
0212     }
0213   }
0214 }
0215
0216 void PFCandidateChecker::printElementsInBlocks(const PFCandidate& cand, ostream& out) const {
0217   if (!out)
0218     return;
0219
0220   PFBlockRef firstRef;
0221
0222   assert(!cand.elementsInBlocks().empty());
0223   for (unsigned i = 0; i < cand.elementsInBlocks().size(); i++) {
0224     PFBlockRef blockRef = cand.elementsInBlocks()[i].first;
0225
0226     if (blockRef.isNull()) {
0227       cerr << "ERROR! no block ref!";
0228       continue;
0229     }
0230
0231     if (!i) {
0232       out << (*blockRef);
0233       firstRef = blockRef;
0234     } else if (blockRef != firstRef) {
0235       cerr << "WARNING! This PFCandidate is not made from a single block" << endl;
0236     }
0237
0238     out << "\t" << cand.elementsInBlocks()[i].second << endl;
0239   }
0240 }