File indexing completed on 2023-03-17 11:15:54
0001
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
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
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.)
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;
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
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
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
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
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
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
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
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);