File indexing completed on 2024-04-06 12:23:42
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
0080 for (unsigned int ij = 0; ij < nJet; ij++) {
0081 const pat::Jet &itJet = (*srcJet)[ij];
0082 if (itJet.pt() <= 10)
0083 continue;
0084 std::vector<fastjet::PseudoJet> lClusterParticles;
0085 float lepPt(-1), lepEta(-1), lepPhi(-1);
0086 int lepId(-1);
0087 for (auto const &d : itJet.daughterPtrVector()) {
0088 fastjet::PseudoJet p(d->px(), d->py(), d->pz(), d->energy());
0089 lClusterParticles.emplace_back(p);
0090 }
0091 int ele_pfmatch_index = -1;
0092 int mu_pfmatch_index = -1;
0093
0094
0095 double dRmin(0.8), dRele(999), dRmu(999), dRtmp(999);
0096 for (unsigned int il(0); il < nEle; il++) {
0097 const auto &lep = eleProd.at(il);
0098 if (matchByCommonSourceCandidatePtr(lep, itJet)) {
0099 dRtmp = reco::deltaR(itJet.eta(), itJet.phi(), lep.eta(), lep.phi());
0100 if (dRtmp < dRmin && dRtmp < dRele && lep.pt() > lepPt) {
0101 lepPt = lep.pt();
0102 lepEta = lep.eta();
0103 lepPhi = lep.phi();
0104 lepId = 11;
0105 ele_pfmatch_index = il;
0106 dRele = dRtmp;
0107 break;
0108 }
0109 }
0110 }
0111 for (unsigned int il(0); il < nMu; il++) {
0112 const auto &lep = muProd.at(il);
0113 if (matchByCommonSourceCandidatePtr(lep, itJet)) {
0114 dRtmp = reco::deltaR(itJet.eta(), itJet.phi(), lep.eta(), lep.phi());
0115 if (dRtmp < dRmin && dRtmp < dRele && dRtmp < dRmu && lep.pt() > lepPt) {
0116 lepPt = lep.pt();
0117 lepEta = lep.eta();
0118 lepPhi = lep.phi();
0119 lepId = 13;
0120 ele_pfmatch_index = -1;
0121 mu_pfmatch_index = il;
0122 dRmu = dRtmp;
0123 break;
0124 }
0125 }
0126 }
0127
0128 std::vector<fastjet::PseudoJet> psub_3;
0129 std::sort(lClusterParticles.begin(), lClusterParticles.end(), orderPseudoJet);
0130 auto lsf_3 = calculateLSF(lClusterParticles, psub_3, lepPt, lepEta, lepPhi, lepId, 2.0, 3);
0131 vlsf3.push_back(std::get<0>(lsf_3));
0132 veleIdx3SJ.push_back(ele_pfmatch_index);
0133 vmuIdx3SJ.push_back(mu_pfmatch_index);
0134 }
0135
0136
0137 auto lsf3V = std::make_unique<edm::ValueMap<float>>();
0138 edm::ValueMap<float>::Filler fillerlsf3(*lsf3V);
0139 fillerlsf3.insert(srcJet, vlsf3.begin(), vlsf3.end());
0140 fillerlsf3.fill();
0141 iEvent.put(std::move(lsf3V), "lsf3");
0142
0143 auto muIdx3SJV = std::make_unique<edm::ValueMap<int>>();
0144 edm::ValueMap<int>::Filler fillermuIdx3SJ(*muIdx3SJV);
0145 fillermuIdx3SJ.insert(srcJet, vmuIdx3SJ.begin(), vmuIdx3SJ.end());
0146 fillermuIdx3SJ.fill();
0147 iEvent.put(std::move(muIdx3SJV), "muIdx3SJ");
0148
0149 auto eleIdx3SJV = std::make_unique<edm::ValueMap<int>>();
0150 edm::ValueMap<int>::Filler fillereleIdx3SJ(*eleIdx3SJV);
0151 fillereleIdx3SJ.insert(srcJet, veleIdx3SJ.begin(), veleIdx3SJ.end());
0152 fillereleIdx3SJ.fill();
0153 iEvent.put(std::move(eleIdx3SJV), "eleIdx3SJ");
0154 }
0155
0156 template <typename T>
0157 bool LeptonInJetProducer<T>::orderPseudoJet(fastjet::PseudoJet j1, fastjet::PseudoJet j2) {
0158 return j1.perp2() > j2.perp2();
0159 }
0160
0161 template <typename T>
0162 std::tuple<float, float> LeptonInJetProducer<T>::calculateLSF(std::vector<fastjet::PseudoJet> iCParticles,
0163 std::vector<fastjet::PseudoJet> &lsubjets,
0164 float ilPt,
0165 float ilEta,
0166 float ilPhi,
0167 int ilId,
0168 double dr,
0169 int nsj) const {
0170 float lsf(-1), lmd(-1);
0171 if (ilPt > 0 && (ilId == 11 || ilId == 13)) {
0172 TLorentzVector ilep;
0173 if (ilId == 11)
0174 ilep.SetPtEtaPhiM(ilPt, ilEta, ilPhi, 0.000511);
0175 if (ilId == 13)
0176 ilep.SetPtEtaPhiM(ilPt, ilEta, ilPhi, 0.105658);
0177 fastjet::JetDefinition lCJet_def(fastjet::kt_algorithm, dr);
0178 fastjet::ClusterSequence lCClust_seq(iCParticles, lCJet_def);
0179 if (dr > 0.5) {
0180 lsubjets = sorted_by_pt(lCClust_seq.exclusive_jets_up_to(nsj));
0181 } else {
0182 lsubjets = sorted_by_pt(lCClust_seq.inclusive_jets());
0183 }
0184 int lId(-1);
0185 double dRmin = 999.;
0186 for (unsigned int i0 = 0; i0 < lsubjets.size(); i0++) {
0187 double dR = reco::deltaR(lsubjets[i0].eta(), lsubjets[i0].phi(), ilep.Eta(), ilep.Phi());
0188 if (dR < dRmin) {
0189 dRmin = dR;
0190 lId = i0;
0191 }
0192 }
0193 if (lId != -1) {
0194 TLorentzVector pVec;
0195 pVec.SetPtEtaPhiM(lsubjets[lId].pt(), lsubjets[lId].eta(), lsubjets[lId].phi(), lsubjets[lId].m());
0196 lsf = ilep.Pt() / pVec.Pt();
0197 lmd = (ilep - pVec).M() / pVec.M();
0198 }
0199 }
0200 return std::tuple<float, float>(lsf, lmd);
0201 }
0202
0203 template <typename T>
0204 void LeptonInJetProducer<T>::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0205 edm::ParameterSetDescription desc;
0206 desc.add<edm::InputTag>("src")->setComment("jet input collection");
0207 desc.add<edm::InputTag>("srcEle")->setComment("electron input collection");
0208 desc.add<edm::InputTag>("srcMu")->setComment("muon input collection");
0209 descriptions.addWithDefaultLabel(desc);
0210 }
0211
0212 typedef LeptonInJetProducer<pat::Jet> LepInJetProducer;
0213
0214
0215 DEFINE_FWK_MODULE(LepInJetProducer);
