NanoAOD/plugins/LeptonJetVarProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PhysicsTools/NanoAOD
0004 // Class:      LeptonJetVarProducer
0005 //
0006 /**\class LeptonJetVarProducer LeptonJetVarProducer.cc PhysicsTools/NanoAOD/plugins/LeptonJetVarProducer.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Marco Peruzzi
0015 //         Created:  Tue, 05 Sep 2017 12:24:38 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/global/EDProducer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/StreamID.h"
0031
0032 #include "DataFormats/PatCandidates/interface/Muon.h"
0033 #include "DataFormats/PatCandidates/interface/Jet.h"
0034 #include "DataFormats/PatCandidates/interface/Electron.h"
0035 #include "DataFormats/VertexReco/interface/Vertex.h"
0036
0037 #include "DataFormats/Common/interface/View.h"
0038 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0039
0040 #include "PhysicsTools/NanoAOD/interface/MatchingUtils.h"
0041
0042 //
0043 // class declaration
0044 //
0045
0046 template <typename T>
0047 class LeptonJetVarProducer : public edm::global::EDProducer<> {
0048 public:
0049   explicit LeptonJetVarProducer(const edm::ParameterSet& iConfig)
0050       : srcJet_(consumes<edm::View<pat::Jet>>(iConfig.getParameter<edm::InputTag>("srcJet"))),
0051         srcLep_(consumes<edm::View<T>>(iConfig.getParameter<edm::InputTag>("srcLep"))),
0052         srcVtx_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("srcVtx"))) {
0053     produces<edm::ValueMap<float>>("ptRatio");
0054     produces<edm::ValueMap<float>>("ptRel");
0055     produces<edm::ValueMap<float>>("jetNDauChargedMVASel");
0056     produces<edm::ValueMap<reco::CandidatePtr>>("jetForLepJetVar");
0057   }
0058   ~LeptonJetVarProducer() override {}
0059
0060   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0061
0062 private:
0063   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0064
0065   std::tuple<float, float, float> calculatePtRatioRel(edm::Ptr<reco::Candidate> lep,
0066                                                       edm::Ptr<pat::Jet> jet,
0067                                                       const reco::Vertex& vtx) const;
0068
0069   // ----------member data ---------------------------
0070
0071   edm::EDGetTokenT<edm::View<pat::Jet>> srcJet_;
0072   edm::EDGetTokenT<edm::View<T>> srcLep_;
0073   edm::EDGetTokenT<std::vector<reco::Vertex>> srcVtx_;
0074 };
0075
0076 //
0077 // constants, enums and typedefs
0078 //
0079
0080 //
0081 // static data member definitions
0082 //
0083
0084 //
0085 // member functions
0086 //
0087
0088 // ------------ method called to produce the data  ------------
0089 template <typename T>
0090 void LeptonJetVarProducer<T>::produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0091   auto srcJet = iEvent.getHandle(srcJet_);
0092   auto srcLep = iEvent.getHandle(srcLep_);
0093   const auto& vtxProd = iEvent.get(srcVtx_);
0094
0095   unsigned int nJet = srcJet->size();
0096   unsigned int nLep = srcLep->size();
0097
0098   std::vector<float> ptRatio(nLep, -1);
0099   std::vector<float> ptRel(nLep, -1);
0100   std::vector<float> jetNDauChargedMVASel(nLep, 0);
0101   std::vector<reco::CandidatePtr> jetForLepJetVar(nLep, reco::CandidatePtr());
0102
0103   const auto& pv = vtxProd.at(0);
0104
0105   for (unsigned int il = 0; il < nLep; il++) {
0106     for (unsigned int ij = 0; ij < nJet; ij++) {
0107       auto lep = srcLep->ptrAt(il);
0108       auto jet = srcJet->ptrAt(ij);
0109       if (matchByCommonSourceCandidatePtr(*lep, *jet)) {
0110         auto res = calculatePtRatioRel(lep, jet, pv);
0111         ptRatio[il] = std::get<0>(res);
0112         ptRel[il] = std::get<1>(res);
0113         jetNDauChargedMVASel[il] = std::get<2>(res);
0114         jetForLepJetVar[il] = jet;
0115         break;  // take leading jet with shared source candidates
0116       }
0117     }
0118   }
0119
0120   auto ptRatioV = std::make_unique<edm::ValueMap<float>>();
0121   edm::ValueMap<float>::Filler fillerRatio(*ptRatioV);
0122   fillerRatio.insert(srcLep, ptRatio.begin(), ptRatio.end());
0123   fillerRatio.fill();
0124   iEvent.put(std::move(ptRatioV), "ptRatio");
0125
0126   auto ptRelV = std::make_unique<edm::ValueMap<float>>();
0127   edm::ValueMap<float>::Filler fillerRel(*ptRelV);
0128   fillerRel.insert(srcLep, ptRel.begin(), ptRel.end());
0129   fillerRel.fill();
0130   iEvent.put(std::move(ptRelV), "ptRel");
0131
0132   auto jetNDauChargedMVASelV = std::make_unique<edm::ValueMap<float>>();
0133   edm::ValueMap<float>::Filler fillerNDau(*jetNDauChargedMVASelV);
0134   fillerNDau.insert(srcLep, jetNDauChargedMVASel.begin(), jetNDauChargedMVASel.end());
0135   fillerNDau.fill();
0136   iEvent.put(std::move(jetNDauChargedMVASelV), "jetNDauChargedMVASel");
0137
0138   auto jetForLepJetVarV = std::make_unique<edm::ValueMap<reco::CandidatePtr>>();
0139   edm::ValueMap<reco::CandidatePtr>::Filler fillerjetForLepJetVar(*jetForLepJetVarV);
0140   fillerjetForLepJetVar.insert(srcLep, jetForLepJetVar.begin(), jetForLepJetVar.end());
0141   fillerjetForLepJetVar.fill();
0142   iEvent.put(std::move(jetForLepJetVarV), "jetForLepJetVar");
0143 }
0144
0145 template <typename T>
0146 std::tuple<float, float, float> LeptonJetVarProducer<T>::calculatePtRatioRel(edm::Ptr<reco::Candidate> lep,
0147                                                                              edm::Ptr<pat::Jet> jet,
0148                                                                              const reco::Vertex& vtx) const {
0149   auto rawp4 = jet->correctedP4("Uncorrected");
0150   auto lepp4 = lep->p4();
0151
0152   if ((rawp4 - lepp4).R() < 1e-4)
0153     return std::tuple<float, float, float>(1.0, 0.0, 0.0);
0154
0155   auto l1corrFactor = jet->jecFactor("L1FastJet") / jet->jecFactor("Uncorrected");
0156
0157   auto jetp4 = (rawp4 - lepp4 * (1.0 / l1corrFactor)) * (jet->pt() / rawp4.pt()) + lepp4;
0158   auto ptratio = lepp4.pt() / jetp4.pt();
0159   auto ptrel = lepp4.Vect().Cross((jetp4 - lepp4).Vect().Unit()).R();
0160
0161   unsigned int jndau = 0;
0162   for (const auto& _d : jet->daughterPtrVector()) {
0163     const auto d = dynamic_cast<const pat::PackedCandidate*>(_d.get());
0164     if (d->charge() == 0)
0165       continue;
0166     if (d->fromPV() <= 1)
0167       continue;
0168     if (deltaR(*d, *lep) > 0.4)
0169       continue;
0170     if (!(d->hasTrackDetails()))
0171       continue;
0172     auto tk = d->pseudoTrack();
0173     if (tk.pt() > 1 && tk.hitPattern().numberOfValidHits() >= 8 && tk.hitPattern().numberOfValidPixelHits() >= 2 &&
0174         tk.normalizedChi2() < 5 && fabs(tk.dxy(vtx.position())) < 0.2 && fabs(tk.dz(vtx.position())) < 17)
0175       jndau++;
0176   }
0177
0178   return std::tuple<float, float, float>(ptratio, ptrel, float(jndau));
0179 }
0180
0181 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0182 template <typename T>
0183 void LeptonJetVarProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0184   //The following says we do not know what parameters are allowed so do no validation
0185   // Please change this to state exactly what you do use, even if it is no parameters
0186   edm::ParameterSetDescription desc;
0187   desc.add<edm::InputTag>("srcJet")->setComment("jet input collection");
0188   desc.add<edm::InputTag>("srcLep")->setComment("lepton input collection");
0189   desc.add<edm::InputTag>("srcVtx")->setComment("primary vertex input collection");
0190   std::string modname;
0191   if (typeid(T) == typeid(pat::Muon))
0192     modname += "Muon";
0193   else if (typeid(T) == typeid(pat::Electron))
0194     modname += "Electron";
0195   modname += "JetVarProducer";
0196   descriptions.add(modname, desc);
0197 }
0198
0199 typedef LeptonJetVarProducer<pat::Muon> MuonJetVarProducer;
0200 typedef LeptonJetVarProducer<pat::Electron> ElectronJetVarProducer;
0201
0202 //define this as a plug-in
0203 DEFINE_FWK_MODULE(MuonJetVarProducer);
0204 DEFINE_FWK_MODULE(ElectronJetVarProducer);