File indexing completed on 2022-11-09 11:29:54
0001
0002 #include <memory>
0003
0004
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/global/EDProducer.h"
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/StreamID.h"
0013
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016
0017 #include "DataFormats/PatCandidates/interface/Muon.h"
0018 #include "DataFormats/PatCandidates/interface/Jet.h"
0019 #include "DataFormats/PatCandidates/interface/Electron.h"
0020 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0021
0022 #include "DataFormats/NanoAOD/interface/FlatTable.h"
0023 #include "fastjet/PseudoJet.hh"
0024 #include <fastjet/JetDefinition.hh>
0025 #include <TLorentzVector.h>
0026 #include <TMath.h>
0027
0028 #include "PhysicsTools/NanoAOD/interface/MatchingUtils.h"
0029
0030 template <typename T>
0031 class LeptonInJetProducer : public edm::global::EDProducer<> {
0032 public:
0033 explicit LeptonInJetProducer(const edm::ParameterSet &iConfig)
0034 : srcJet_(consumes<edm::View<pat::Jet>>(iConfig.getParameter<edm::InputTag>("src"))),
0035 srcEle_(consumes<edm::View<pat::Electron>>(iConfig.getParameter<edm::InputTag>("srcEle"))),
0036 srcMu_(consumes<edm::View<pat::Muon>>(iConfig.getParameter<edm::InputTag>("srcMu"))) {
0037 produces<edm::ValueMap<float>>("lsf3");
0038 produces<edm::ValueMap<int>>("muIdx3SJ");
0039 produces<edm::ValueMap<int>>("eleIdx3SJ");
0040 }
0041 ~LeptonInJetProducer() override{};
0042
0043 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0044
0045 private:
0046 void produce(edm::StreamID, edm::Event &, edm::EventSetup const &) const override;
0047
0048 static bool orderPseudoJet(fastjet::PseudoJet j1, fastjet::PseudoJet j2);
0049 std::tuple<float, float> calculateLSF(std::vector<fastjet::PseudoJet> iCParticles,
0050 std::vector<fastjet::PseudoJet> &ljets,
0051 float ilPt,
0052 float ilEta,
0053 float ilPhi,
0054 int ilId,
0055 double dr,
0056 int nsj) const;
0057
0058 edm::EDGetTokenT<edm::View<pat::Jet>> srcJet_;
0059 edm::EDGetTokenT<edm::View<pat::Electron>> srcEle_;
0060 edm::EDGetTokenT<edm::View<pat::Muon>> srcMu_;
0061 };
0062
0063
0064 template <typename T>
0065 void LeptonInJetProducer<T>::produce(edm::StreamID streamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0066
0067 auto srcJet = iEvent.getHandle(srcJet_);
0068 const auto &eleProd = iEvent.get(srcEle_);
0069 const auto &muProd = iEvent.get(srcMu_);
0070
0071 unsigned int nJet = srcJet->size();
0072 unsigned int nEle = eleProd.size();
0073 unsigned int nMu = muProd.size();
0074
0075 std::vector<float> vlsf3;
0076 std::vector<int> vmuIdx3SJ;
0077 std::vector<int> veleIdx3SJ;
0078
0079 int ele_pfmatch_index = -1;
0080 int mu_pfmatch_index = -1;
0081
0082
0083 for (unsigned int ij = 0; ij < nJet; ij++) {
0084 const pat::Jet &itJet = (*srcJet)[ij];
0085 if (itJet.pt() <= 10)
0086 continue;
0087 std::vector<fastjet::PseudoJet> lClusterParticles;
0088 float lepPt(-1), lepEta(-1), lepPhi(-1);
0089 int lepId(-1);
0090 for (auto const &d : itJet.daughterPtrVector()) {
0091 fastjet::PseudoJet p(d->px(), d->py(), d->pz(), d->energy());
0092 lClusterParticles.emplace_back(p);
0093 }
0094
0095
0096 double dRmin(0.8), dRele(999), dRmu(999), dRtmp(999);
0097 for (unsigned int il(0); il < nEle; il++) {
0098 const auto &lep = eleProd.at(il);
0099 if (matchByCommonSourceCandidatePtr(lep, itJet)) {
0100 dRtmp = reco::deltaR(itJet.eta(), itJet.phi(), lep.eta(), lep.phi());
0101 if (dRtmp < dRmin && dRtmp < dRele && lep.pt() > lepPt) {
0102 lepPt = lep.pt();
0103 lepEta = lep.eta();
0104 lepPhi = lep.phi();
0105 lepId = 11;
0106 ele_pfmatch_index = il;
0107 dRele = dRtmp;
0108 break;
0109 }
0110 }
0111 }
0112 for (unsigned int il(0); il < nMu; il++) {
0113 const auto &lep = muProd.at(il);
0114 if (matchByCommonSourceCandidatePtr(lep, itJet)) {
0115 dRtmp = reco::deltaR(itJet.eta(), itJet.phi(), lep.eta(), lep.phi());
0116 if (dRtmp < dRmin && dRtmp < dRele && dRtmp < dRmu && lep.pt() > lepPt) {
0117 lepPt = lep.pt();
0118 lepEta = lep.eta();
0119 lepPhi = lep.phi();
0120 lepId = 13;
0121 ele_pfmatch_index = -1;
0122 mu_pfmatch_index = il;
0123 dRmu = dRtmp;
0124 break;
0125 }
0126 }
0127 }
0128
0129 std::vector<fastjet::PseudoJet> psub_3;
0130 std::sort(lClusterParticles.begin(), lClusterParticles.end(), orderPseudoJet);
0131 auto lsf_3 = calculateLSF(lClusterParticles, psub_3, lepPt, lepEta, lepPhi, lepId, 2.0, 3);
0132 vlsf3.push_back(std::get<0>(lsf_3));
0133 veleIdx3SJ.push_back(ele_pfmatch_index);
0134 vmuIdx3SJ.push_back(mu_pfmatch_index);
0135 }
0136
0137
0138 auto lsf3V = std::make_unique<edm::ValueMap<float>>();
0139 edm::ValueMap<float>::Filler fillerlsf3(*lsf3V);
0140 fillerlsf3.insert(srcJet, vlsf3.begin(), vlsf3.end());
0141 fillerlsf3.fill();
0142 iEvent.put(std::move(lsf3V), "lsf3");
0143
0144 auto muIdx3SJV = std::make_unique<edm::ValueMap<int>>();
0145 edm::ValueMap<int>::Filler fillermuIdx3SJ(*muIdx3SJV);
0146 fillermuIdx3SJ.insert(srcJet, vmuIdx3SJ.begin(), vmuIdx3SJ.end());
0147 fillermuIdx3SJ.fill();
0148 iEvent.put(std::move(muIdx3SJV), "muIdx3SJ");
0149
0150 auto eleIdx3SJV = std::make_unique<edm::ValueMap<int>>();
0151 edm::ValueMap<int>::Filler fillereleIdx3SJ(*eleIdx3SJV);
0152 fillereleIdx3SJ.insert(srcJet, veleIdx3SJ.begin(), veleIdx3SJ.end());
0153 fillereleIdx3SJ.fill();
0154 iEvent.put(std::move(eleIdx3SJV), "eleIdx3SJ");
0155 }
0156
0157 template <typename T>
0158 bool LeptonInJetProducer<T>::orderPseudoJet(fastjet::PseudoJet j1, fastjet::PseudoJet j2) {
0159 return j1.perp2() > j2.perp2();
0160 }
0161
0162 template <typename T>
0163 std::tuple<float, float> LeptonInJetProducer<T>::calculateLSF(std::vector<fastjet::PseudoJet> iCParticles,
0164 std::vector<fastjet::PseudoJet> &lsubjets,
0165 float ilPt,
0166 float ilEta,
0167 float ilPhi,
0168 int ilId,
0169 double dr,
0170 int nsj) const {
0171 float lsf(-1), lmd(-1);
0172 if (ilPt > 0 && (ilId == 11 || ilId == 13)) {
0173 TLorentzVector ilep;
0174 if (ilId == 11)
0175 ilep.SetPtEtaPhiM(ilPt, ilEta, ilPhi, 0.000511);
0176 if (ilId == 13)
0177 ilep.SetPtEtaPhiM(ilPt, ilEta, ilPhi, 0.105658);
0178 fastjet::JetDefinition lCJet_def(fastjet::kt_algorithm, dr);
0179 fastjet::ClusterSequence lCClust_seq(iCParticles, lCJet_def);
0180 if (dr > 0.5) {
0181 lsubjets = sorted_by_pt(lCClust_seq.exclusive_jets_up_to(nsj));
0182 } else {
0183 lsubjets = sorted_by_pt(lCClust_seq.inclusive_jets());
0184 }
0185 int lId(-1);
0186 double dRmin = 999.;
0187 for (unsigned int i0 = 0; i0 < lsubjets.size(); i0++) {
0188 double dR = reco::deltaR(lsubjets[i0].eta(), lsubjets[i0].phi(), ilep.Eta(), ilep.Phi());
0189 if (dR < dRmin) {
0190 dRmin = dR;
0191 lId = i0;
0192 }
0193 }
0194 if (lId != -1) {
0195 TLorentzVector pVec;
0196 pVec.SetPtEtaPhiM(lsubjets[lId].pt(), lsubjets[lId].eta(), lsubjets[lId].phi(), lsubjets[lId].m());
0197 lsf = ilep.Pt() / pVec.Pt();
0198 lmd = (ilep - pVec).M() / pVec.M();
0199 }
0200 }
0201 return std::tuple<float, float>(lsf, lmd);
0202 }
0203
0204 template <typename T>
0205 void LeptonInJetProducer<T>::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0206 edm::ParameterSetDescription desc;
0207 desc.add<edm::InputTag>("src")->setComment("jet input collection");
0208 desc.add<edm::InputTag>("srcEle")->setComment("electron input collection");
0209 desc.add<edm::InputTag>("srcMu")->setComment("muon input collection");
0210 descriptions.addWithDefaultLabel(desc);
0211 }
0212
0213 typedef LeptonInJetProducer<pat::Jet> LepInJetProducer;
0214
0215
0216 DEFINE_FWK_MODULE(LepInJetProducer);
