JetExamples/test/printJetFlavourInfo.cc

0001 // user include files
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009
0010 #include "DataFormats/JetReco/interface/Jet.h"
0011 #include "DataFormats/JetMatching/interface/JetFlavourInfo.h"
0012 #include "DataFormats/JetMatching/interface/JetFlavourInfoMatching.h"
0013
0014 #include "DataFormats/Math/interface/deltaR.h"
0015
0016 // system include files
0017 #include <memory>
0018
0019 class printJetFlavourInfo : public edm::one::EDAnalyzer<> {
0020 public:
0021   explicit printJetFlavourInfo(const edm::ParameterSet&);
0022   ~printJetFlavourInfo() {}
0023   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup);
0024
0025 private:
0026   edm::EDGetTokenT<reco::JetFlavourInfoMatchingCollection> jetFlavourInfosToken_;
0027   edm::EDGetTokenT<reco::JetFlavourInfoMatchingCollection> subjetFlavourInfosToken_;
0028   edm::EDGetTokenT<edm::View<reco::Jet> > groomedJetsToken_;
0029   bool useSubjets_;
0030 };
0031
0032 printJetFlavourInfo::printJetFlavourInfo(const edm::ParameterSet& iConfig) {
0033   jetFlavourInfosToken_ =
0034       consumes<reco::JetFlavourInfoMatchingCollection>(iConfig.getParameter<edm::InputTag>("jetFlavourInfos"));
0035   subjetFlavourInfosToken_ = mayConsume<reco::JetFlavourInfoMatchingCollection>(
0036       iConfig.exists("subjetFlavourInfos") ? iConfig.getParameter<edm::InputTag>("subjetFlavourInfos")
0037                                            : edm::InputTag());
0038   groomedJetsToken_ = mayConsume<edm::View<reco::Jet> >(
0039       iConfig.exists("groomedJets") ? iConfig.getParameter<edm::InputTag>("groomedJets") : edm::InputTag());
0040   useSubjets_ = (iConfig.exists("subjetFlavourInfos") && iConfig.exists("groomedJets"));
0041 }
0042
0043 void printJetFlavourInfo::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0044   edm::Handle<reco::JetFlavourInfoMatchingCollection> theJetFlavourInfos;
0045   iEvent.getByToken(jetFlavourInfosToken_, theJetFlavourInfos);
0046
0047   edm::Handle<reco::JetFlavourInfoMatchingCollection> theSubjetFlavourInfos;
0048   edm::Handle<edm::View<reco::Jet> > groomedJets;
0049
0050   std::vector<int> matchedIndices;
0051   if (useSubjets_) {
0052     iEvent.getByToken(subjetFlavourInfosToken_, theSubjetFlavourInfos);
0053     iEvent.getByToken(groomedJetsToken_, groomedJets);
0054
0055     // match groomed and original jet
0056     std::vector<bool> jetLocks(theJetFlavourInfos->size(), false);
0057     std::vector<int> jetIndices;
0058
0059     for (size_t gj = 0; gj < groomedJets->size(); ++gj) {
0060       double matchedDR2 = 1e9;
0061       int matchedIdx = -1;
0062
0063       if (groomedJets->at(gj).pt() > 0.)  // skips pathological cases of groomed jets with Pt=0
0064       {
0065         for (reco::JetFlavourInfoMatchingCollection::const_iterator j = theJetFlavourInfos->begin();
0066              j != theJetFlavourInfos->end();
0067              ++j) {
0068           if (jetLocks.at(j - theJetFlavourInfos->begin()))
0069             continue;  // skip jets that have already been matched
0070
0071           double tempDR2 = reco::deltaR2(
0072               j->first->rapidity(), j->first->phi(), groomedJets->at(gj).rapidity(), groomedJets->at(gj).phi());
0073           if (tempDR2 < matchedDR2) {
0074             matchedDR2 = tempDR2;
0075             matchedIdx = (j - theJetFlavourInfos->begin());
0076           }
0077         }
0078       }
0079
0080       if (matchedIdx >= 0)
0081         jetLocks.at(matchedIdx) = true;
0082       jetIndices.push_back(matchedIdx);
0083     }
0084
0085     for (size_t j = 0; j < theJetFlavourInfos->size(); ++j) {
0086       std::vector<int>::iterator matchedIndex = std::find(jetIndices.begin(), jetIndices.end(), j);
0087
0088       matchedIndices.push_back(matchedIndex != jetIndices.end() ? std::distance(jetIndices.begin(), matchedIndex) : -1);
0089     }
0090   }
0091
0092   for (reco::JetFlavourInfoMatchingCollection::const_iterator j = theJetFlavourInfos->begin();
0093        j != theJetFlavourInfos->end();
0094        ++j) {
0095     std::cout << "-------------------- Jet Flavour Info --------------------" << std::endl;
0096
0097     const reco::Jet* aJet = (*j).first.get();
0098     reco::JetFlavourInfo aInfo = (*j).second;
0099     std::cout << std::setprecision(2) << std::setw(6) << std::fixed << "[printJetFlavourInfo] Jet "
0100               << (j - theJetFlavourInfos->begin()) << " pt, eta, rapidity, phi = " << aJet->pt() << ", " << aJet->eta()
0101               << ", " << aJet->rapidity() << ", " << aJet->phi() << std::endl;
0102     // ----------------------- Hadrons -------------------------------
0103     std::cout << "                      Hadron-based flavour: " << aInfo.getHadronFlavour() << std::endl;
0104
0105     const reco::GenParticleRefVector& bHadrons = aInfo.getbHadrons();
0106     std::cout << "                      # of clustered b hadrons: " << bHadrons.size() << std::endl;
0107     for (reco::GenParticleRefVector::const_iterator it = bHadrons.begin(); it != bHadrons.end(); ++it) {
0108       float dist = reco::deltaR(aJet->eta(), aJet->phi(), (*it)->eta(), (*it)->phi());
0109       float dist2 = reco::deltaR(aJet->rapidity(), aJet->phi(), (*it)->rapidity(), (*it)->phi());
0110       std::cout << "                        b hadron " << (it - bHadrons.begin())
0111                 << " PdgID, status, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId() << ", "
0112                 << (*it)->status() << ", (" << (*it)->pt() << "," << (*it)->eta() << "," << (*it)->rapidity() << ","
0113                 << (*it)->phi() << "), " << dist << ", " << dist2 << std::endl;
0114     }
0115
0116     const reco::GenParticleRefVector& cHadrons = aInfo.getcHadrons();
0117     std::cout << "                      # of clustered c hadrons: " << cHadrons.size() << std::endl;
0118     for (reco::GenParticleRefVector::const_iterator it = cHadrons.begin(); it != cHadrons.end(); ++it) {
0119       float dist = reco::deltaR(aJet->eta(), aJet->phi(), (*it)->eta(), (*it)->phi());
0120       float dist2 = reco::deltaR(aJet->rapidity(), aJet->phi(), (*it)->rapidity(), (*it)->phi());
0121       std::cout << "                        c hadron " << (it - cHadrons.begin())
0122                 << " PdgID, status, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId() << ", "
0123                 << (*it)->status() << ", (" << (*it)->pt() << "," << (*it)->eta() << "," << (*it)->rapidity() << ","
0124                 << (*it)->phi() << "), " << dist << ", " << dist2 << std::endl;
0125     }
0126     // ----------------------- Partons -------------------------------
0127     std::cout << "                      Parton-based flavour: " << aInfo.getPartonFlavour() << std::endl;
0128
0129     const reco::GenParticleRefVector& partons = aInfo.getPartons();
0130     std::cout << "                      # of clustered partons: " << partons.size() << std::endl;
0131     for (reco::GenParticleRefVector::const_iterator it = partons.begin(); it != partons.end(); ++it) {
0132       float dist = reco::deltaR(aJet->eta(), aJet->phi(), (*it)->eta(), (*it)->phi());
0133       float dist2 = reco::deltaR(aJet->rapidity(), aJet->phi(), (*it)->rapidity(), (*it)->phi());
0134       std::cout << "                        Parton " << (it - partons.begin())
0135                 << " PdgID, status, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId() << ", "
0136                 << (*it)->status() << ", (" << (*it)->pt() << "," << (*it)->eta() << "," << (*it)->rapidity() << ","
0137                 << (*it)->phi() << "), " << dist << ", " << dist2 << std::endl;
0138     }
0139     // ----------------------- Leptons -------------------------------
0140     const reco::GenParticleRefVector& leptons = aInfo.getLeptons();
0141     std::cout << "                      # of clustered leptons: " << leptons.size() << std::endl;
0142     for (reco::GenParticleRefVector::const_iterator it = leptons.begin(); it != leptons.end(); ++it) {
0143       float dist = reco::deltaR(aJet->eta(), aJet->phi(), (*it)->eta(), (*it)->phi());
0144       float dist2 = reco::deltaR(aJet->rapidity(), aJet->phi(), (*it)->rapidity(), (*it)->phi());
0145       std::cout << "                        Lepton " << (it - leptons.begin())
0146                 << " PdgID, status, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId() << ", "
0147                 << (*it)->status() << ", (" << (*it)->pt() << "," << (*it)->eta() << "," << (*it)->rapidity() << ","
0148                 << (*it)->phi() << "), " << dist << ", " << dist2 << std::endl;
0149     }
0150
0151     if (useSubjets_) {
0152       if (matchedIndices.at(j - theJetFlavourInfos->begin()) < 0) {
0153         std::cout << "  ----------------------- Subjet Flavour Info -----------------------" << std::endl;
0154         std::cout << "  No subjets assigned to this jet" << std::endl;
0155         continue;
0156       }
0157
0158       // loop over subjets
0159       std::cout << "  ----------------------- Subjet Flavour Info -----------------------" << std::endl;
0160       for (size_t s = 0; s < groomedJets->at(matchedIndices.at(j - theJetFlavourInfos->begin())).numberOfDaughters();
0161            ++s) {
0162         const edm::Ptr<reco::Candidate>& subjet =
0163             groomedJets->at(matchedIndices.at(j - theJetFlavourInfos->begin())).daughterPtr(s);
0164
0165         for (reco::JetFlavourInfoMatchingCollection::const_iterator sj = theSubjetFlavourInfos->begin();
0166              sj != theSubjetFlavourInfos->end();
0167              ++sj) {
0168           if (subjet != edm::Ptr<reco::Candidate>((*sj).first.id(), (*sj).first.get(), (*sj).first.key()))
0169             continue;
0170
0171           const reco::Jet* aSubjet = (*sj).first.get();
0172           aInfo = (*sj).second;
0173           std::cout << std::setprecision(2) << std::setw(6) << std::fixed << "  [printSubjetFlavourInfo] Subjet " << s
0174                     << " pt, eta, rapidity, phi, dR(eta-phi), dR(rap-phi) = " << aSubjet->pt() << ", " << aSubjet->eta()
0175                     << ", " << aSubjet->rapidity() << ", " << aSubjet->phi() << ", "
0176                     << reco::deltaR(aSubjet->eta(), aSubjet->phi(), aJet->eta(), aJet->phi()) << ", "
0177                     << reco::deltaR(aSubjet->rapidity(), aSubjet->phi(), aJet->rapidity(), aJet->phi()) << std::endl;
0178           // ----------------------- Hadrons -------------------------------
0179           std::cout << "                           Hadron-based flavour: " << aInfo.getHadronFlavour() << std::endl;
0180
0181           const reco::GenParticleRefVector& bHadrons = aInfo.getbHadrons();
0182           std::cout << "                           # of assigned b hadrons: " << bHadrons.size() << std::endl;
0183           for (reco::GenParticleRefVector::const_iterator it = bHadrons.begin(); it != bHadrons.end(); ++it) {
0184             float dist = reco::deltaR(aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi());
0185             float dist2 = reco::deltaR(aSubjet->rapidity(), aSubjet->phi(), (*it)->rapidity(), (*it)->phi());
0186             std::cout << "                             b hadron " << (it - bHadrons.begin())
0187                       << " PdgID, status, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId() << ", "
0188                       << (*it)->status() << ", (" << (*it)->pt() << "," << (*it)->eta() << "," << (*it)->rapidity()
0189                       << "," << (*it)->phi() << "), " << dist << ", " << dist2 << std::endl;
0190           }
0191
0192           const reco::GenParticleRefVector& cHadrons = aInfo.getcHadrons();
0193           std::cout << "                           # of assigned c hadrons: " << cHadrons.size() << std::endl;
0194           for (reco::GenParticleRefVector::const_iterator it = cHadrons.begin(); it != cHadrons.end(); ++it) {
0195             float dist = reco::deltaR(aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi());
0196             float dist2 = reco::deltaR(aSubjet->rapidity(), aSubjet->phi(), (*it)->rapidity(), (*it)->phi());
0197             std::cout << "                             c hadron " << (it - cHadrons.begin())
0198                       << " PdgID, status, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId() << ", "
0199                       << (*it)->status() << ", (" << (*it)->pt() << "," << (*it)->eta() << "," << (*it)->rapidity()
0200                       << "," << (*it)->phi() << "), " << dist << ", " << dist2 << std::endl;
0201           }
0202           // ----------------------- Partons -------------------------------
0203           std::cout << "                           Parton-based flavour: " << aInfo.getPartonFlavour() << std::endl;
0204
0205           const reco::GenParticleRefVector& partons = aInfo.getPartons();
0206           std::cout << "                           # of assigned partons: " << partons.size() << std::endl;
0207           for (reco::GenParticleRefVector::const_iterator it = partons.begin(); it != partons.end(); ++it) {
0208             float dist = reco::deltaR(aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi());
0209             float dist2 = reco::deltaR(aSubjet->rapidity(), aSubjet->phi(), (*it)->rapidity(), (*it)->phi());
0210             std::cout << "                             Parton " << (it - partons.begin())
0211                       << " PdgID, status, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId() << ", "
0212                       << (*it)->status() << ", (" << (*it)->pt() << "," << (*it)->eta() << "," << (*it)->rapidity()
0213                       << "," << (*it)->phi() << "), " << dist << ", " << dist2 << std::endl;
0214           }
0215           // ----------------------- Leptons -------------------------------
0216           const reco::GenParticleRefVector& leptons = aInfo.getLeptons();
0217           std::cout << "                           # of assigned leptons: " << leptons.size() << std::endl;
0218           for (reco::GenParticleRefVector::const_iterator it = leptons.begin(); it != leptons.end(); ++it) {
0219             float dist = reco::deltaR(aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi());
0220             float dist2 = reco::deltaR(aSubjet->rapidity(), aSubjet->phi(), (*it)->rapidity(), (*it)->phi());
0221             std::cout << "                             Lepton " << (it - leptons.begin())
0222                       << " PdgID, status, (pt,eta,rapidity,phi), dR(eta-phi), dR(rap-phi) = " << (*it)->pdgId() << ", "
0223                       << (*it)->status() << ", (" << (*it)->pt() << "," << (*it)->eta() << "," << (*it)->rapidity()
0224                       << "," << (*it)->phi() << "), " << dist << ", " << dist2 << std::endl;
0225           }
0226         }
0227       }
0228     }
0229   }
0230 }
0231
0232 DEFINE_FWK_MODULE(printJetFlavourInfo);