Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:31

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