Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-01-13 02:35:28

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