Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-11-09 11:29:54

0001 // system include files
0002 #include <memory>
0003 
0004 // user include files
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 // ------------ method called to produce the data  ------------
0064 template <typename T>
0065 void LeptonInJetProducer<T>::produce(edm::StreamID streamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0066   // needs jet collection (srcJet), leptons collection
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   // Find leptons in jets
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     // match to leading and closest electron or muon
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   // Filling table
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 //define this as a plug-in
0216 DEFINE_FWK_MODULE(LepInJetProducer);