Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:32:45

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