NanoAOD/plugins/BJetEnergyRegressionVarProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PhysicsTools/NanoAOD
0004 // Class:      BJetEnergyRegressionVarProducer
0005 //
0006 /**\class BJetEnergyRegressionVarProducer BJetEnergyRegressionVarProducer.cc PhysicsTools/NanoAOD/plugins/BJetEnergyRegressionVarProducer.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 #include <vector>
0022
0023 // user include files
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/global/EDProducer.h"
0026
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/StreamID.h"
0032
0033 #include "DataFormats/Common/interface/View.h"
0034 #include "DataFormats/PatCandidates/interface/Jet.h"
0035 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0036
0037 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0038 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0039 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
0040 #include "DataFormats/VertexReco/interface/Vertex.h"
0041 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0042
0043 //
0044 // class declaration
0045 //
0046
0047 template <typename T>
0048 class BJetEnergyRegressionVarProducer : public edm::global::EDProducer<> {
0049 public:
0050   explicit BJetEnergyRegressionVarProducer(const edm::ParameterSet& iConfig)
0051       : srcJet_(consumes<edm::View<pat::Jet>>(iConfig.getParameter<edm::InputTag>("src"))),
0052         srcVtx_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("pvsrc"))),
0053         srcSV_(consumes<edm::View<reco::VertexCompositePtrCandidate>>(iConfig.getParameter<edm::InputTag>("svsrc"))) {
0054     //un prodotto da copiare
0055     produces<edm::ValueMap<float>>("leptonPtRelv0");
0056     produces<edm::ValueMap<float>>("leptonPtRelInvv0");  //v0 ~ heppy?
0057     produces<edm::ValueMap<int>>("leptonPdgId");
0058     produces<edm::ValueMap<float>>("leptonDeltaR");
0059     produces<edm::ValueMap<float>>("leadTrackPt");
0060     produces<edm::ValueMap<float>>("vtxPt");
0061     produces<edm::ValueMap<float>>("vtxMass");
0062     produces<edm::ValueMap<float>>("vtx3dL");
0063     produces<edm::ValueMap<float>>("vtx3deL");
0064     produces<edm::ValueMap<int>>("vtxNtrk");
0065     produces<edm::ValueMap<float>>("ptD");
0066   }
0067   ~BJetEnergyRegressionVarProducer() override{};
0068
0069   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0070
0071 private:
0072   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0073
0074   std::tuple<float, float, float> calculatePtRatioRelSimple(edm::Ptr<reco::Candidate> lep,
0075                                                             edm::Ptr<pat::Jet> jet) const;  //old version?
0076
0077   // ----------member data ---------------------------
0078
0079   edm::EDGetTokenT<edm::View<pat::Jet>> srcJet_;
0080   edm::EDGetTokenT<std::vector<reco::Vertex>> srcVtx_;
0081   edm::EDGetTokenT<edm::View<reco::VertexCompositePtrCandidate>> srcSV_;
0082 };
0083
0084 //
0085 // constants, enums and typedefs
0086 //
0087
0088 //
0089 // static data member definitions
0090 //
0091
0092 //
0093 // member functions
0094 //
0095
0096 // ------------ method called to produce the data  ------------
0097 template <typename T>
0098 void BJetEnergyRegressionVarProducer<T>::produce(edm::StreamID streamID,
0099                                                  edm::Event& iEvent,
0100                                                  const edm::EventSetup& iSetup) const {
0101   const auto& srcJet = iEvent.getHandle(srcJet_);
0102   const auto& vtxProd = iEvent.get(srcVtx_);
0103   const auto& svProd = iEvent.get(srcSV_);
0104
0105   auto nJet = srcJet->size();
0106
0107   std::vector<float> leptonPtRel_v0(nJet, 0);
0108   std::vector<float> leptonPtRelInv_v0(nJet, 0);
0109   std::vector<int> leptonPdgId(nJet, 0);
0110   std::vector<float> leptonDeltaR(nJet, 0);
0111   std::vector<float> leadTrackPt(nJet, 0);
0112   std::vector<float> vtxPt(nJet, 0);
0113   std::vector<float> vtxMass(nJet, 0);
0114   std::vector<float> vtx3dL(nJet, 0);
0115   std::vector<float> vtx3deL(nJet, 0);
0116   std::vector<int> vtxNtrk(nJet, 0);
0117   std::vector<float> ptD(nJet, 0);
0118
0119   const auto& pv = vtxProd.at(0);
0120   for (unsigned int ij = 0; ij < nJet; ij++) {
0121     auto jet = srcJet->ptrAt(ij);
0122
0123     float ptMax = 0;
0124     float sumWeight = 0;
0125     float sumPt = 0;
0126
0127     for (const auto& d : jet->daughterPtrVector()) {
0128       sumWeight += (d->pt()) * (d->pt());
0129       sumPt += d->pt();
0130       if (d->pt() > ptMax)
0131         ptMax = d->pt();
0132     }
0133     leadTrackPt[ij] = ptMax;
0134     ptD[ij] = (sumWeight > 0 ? sqrt(sumWeight) / sumPt : 0);
0135
0136     //lepton properties
0137     float maxLepPt = 0;
0138
0139     for (const auto& d : jet->daughterPtrVector()) {
0140       if (abs(d->pdgId()) == 11 || abs(d->pdgId()) == 13) {
0141         if (d->pt() < maxLepPt)
0142           continue;
0143         auto res2 = calculatePtRatioRelSimple(d, jet);
0144         leptonPtRel_v0[ij] = std::get<1>(res2);
0145         leptonPtRelInv_v0[ij] = std::get<2>(res2);
0146         leptonPdgId[ij] = d->pdgId();
0147         leptonDeltaR[ij] = reco::deltaR(jet->p4(), d->p4());
0148         maxLepPt = d->pt();
0149       }
0150     }
0151
0152     //Fill vertex properties
0153     VertexDistance3D vdist;
0154     float maxFoundSignificance = 0;
0155
0156     vtxPt[ij] = 0;
0157     vtxMass[ij] = 0;
0158     vtx3dL[ij] = 0;
0159     vtx3deL[ij] = 0;
0160     vtxNtrk[ij] = 0;
0161
0162     for (const auto& sv : svProd) {
0163       GlobalVector flightDir(sv.vertex().x() - pv.x(), sv.vertex().y() - pv.y(), sv.vertex().z() - pv.z());
0164       GlobalVector jetDir(jet->px(), jet->py(), jet->pz());
0165       if (reco::deltaR2(flightDir, jetDir) < 0.09) {
0166         Measurement1D dl = vdist.distance(
0167             pv, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0168         if (dl.significance() > maxFoundSignificance) {
0169           maxFoundSignificance = dl.significance();
0170           vtxPt[ij] = sv.pt();
0171           vtxMass[ij] = sv.p4().M();
0172           vtx3dL[ij] = dl.value();
0173           vtx3deL[ij] = dl.error();
0174           vtxNtrk[ij] = sv.numberOfSourceCandidatePtrs();
0175         }
0176       }
0177     }
0178   }
0179
0180   std::unique_ptr<edm::ValueMap<float>> leptonPtRelV_v0(new edm::ValueMap<float>());
0181   edm::ValueMap<float>::Filler fillerRel_v0(*leptonPtRelV_v0);
0182   fillerRel_v0.insert(srcJet, leptonPtRel_v0.begin(), leptonPtRel_v0.end());
0183   fillerRel_v0.fill();
0184   iEvent.put(std::move(leptonPtRelV_v0), "leptonPtRelv0");
0185
0186   std::unique_ptr<edm::ValueMap<float>> leptonPtRelInvV_v0(new edm::ValueMap<float>());
0187   edm::ValueMap<float>::Filler fillerRelInv_v0(*leptonPtRelInvV_v0);
0188   fillerRelInv_v0.insert(srcJet, leptonPtRelInv_v0.begin(), leptonPtRelInv_v0.end());
0189   fillerRelInv_v0.fill();
0190   iEvent.put(std::move(leptonPtRelInvV_v0), "leptonPtRelInvv0");
0191
0192   std::unique_ptr<edm::ValueMap<float>> leptonDeltaRV(new edm::ValueMap<float>());
0193   edm::ValueMap<float>::Filler fillerLdR(*leptonDeltaRV);
0194   fillerLdR.insert(srcJet, leptonDeltaR.begin(), leptonDeltaR.end());
0195   fillerLdR.fill();
0196   iEvent.put(std::move(leptonDeltaRV), "leptonDeltaR");
0197
0198   std::unique_ptr<edm::ValueMap<int>> leptonPtPdgIdV(new edm::ValueMap<int>());
0199   edm::ValueMap<int>::Filler fillerId(*leptonPtPdgIdV);
0200   fillerId.insert(srcJet, leptonPdgId.begin(), leptonPdgId.end());
0201   fillerId.fill();
0202   iEvent.put(std::move(leptonPtPdgIdV), "leptonPdgId");
0203
0204   std::unique_ptr<edm::ValueMap<float>> leadTrackPtV(new edm::ValueMap<float>());
0205   edm::ValueMap<float>::Filler fillerLT(*leadTrackPtV);
0206   fillerLT.insert(srcJet, leadTrackPt.begin(), leadTrackPt.end());
0207   fillerLT.fill();
0208   iEvent.put(std::move(leadTrackPtV), "leadTrackPt");
0209
0210   std::unique_ptr<edm::ValueMap<float>> vtxPtV(new edm::ValueMap<float>());
0211   edm::ValueMap<float>::Filler fillerVtxPt(*vtxPtV);
0212   fillerVtxPt.insert(srcJet, vtxPt.begin(), vtxPt.end());
0213   fillerVtxPt.fill();
0214   iEvent.put(std::move(vtxPtV), "vtxPt");
0215
0216   std::unique_ptr<edm::ValueMap<float>> vtxMassV(new edm::ValueMap<float>());
0217   edm::ValueMap<float>::Filler fillerVtxMass(*vtxMassV);
0218   fillerVtxMass.insert(srcJet, vtxMass.begin(), vtxMass.end());
0219   fillerVtxMass.fill();
0220   iEvent.put(std::move(vtxMassV), "vtxMass");
0221
0222   std::unique_ptr<edm::ValueMap<float>> vtx3dLV(new edm::ValueMap<float>());
0223   edm::ValueMap<float>::Filler fillerVtx3dL(*vtx3dLV);
0224   fillerVtx3dL.insert(srcJet, vtx3dL.begin(), vtx3dL.end());
0225   fillerVtx3dL.fill();
0226   iEvent.put(std::move(vtx3dLV), "vtx3dL");
0227
0228   std::unique_ptr<edm::ValueMap<float>> vtx3deLV(new edm::ValueMap<float>());
0229   edm::ValueMap<float>::Filler fillerVtx3deL(*vtx3deLV);
0230   fillerVtx3deL.insert(srcJet, vtx3deL.begin(), vtx3deL.end());
0231   fillerVtx3deL.fill();
0232   iEvent.put(std::move(vtx3deLV), "vtx3deL");
0233
0234   std::unique_ptr<edm::ValueMap<int>> vtxNtrkV(new edm::ValueMap<int>());
0235   edm::ValueMap<int>::Filler fillerVtxNT(*vtxNtrkV);
0236   fillerVtxNT.insert(srcJet, vtxNtrk.begin(), vtxNtrk.end());
0237   fillerVtxNT.fill();
0238   iEvent.put(std::move(vtxNtrkV), "vtxNtrk");
0239
0240   std::unique_ptr<edm::ValueMap<float>> ptDV(new edm::ValueMap<float>());
0241   edm::ValueMap<float>::Filler fillerPtD(*ptDV);
0242   fillerPtD.insert(srcJet, ptD.begin(), ptD.end());
0243   fillerPtD.fill();
0244   iEvent.put(std::move(ptDV), "ptD");
0245 }
0246
0247 template <typename T>
0248 std::tuple<float, float, float> BJetEnergyRegressionVarProducer<T>::calculatePtRatioRelSimple(
0249     edm::Ptr<reco::Candidate> lep, edm::Ptr<pat::Jet> jet) const {
0250   auto lepp4 = lep->p4();
0251   auto rawp4 = jet->correctedP4("Uncorrected");
0252
0253   if ((rawp4 - lepp4).R() < 1e-4)
0254     return std::tuple<float, float, float>(1.0, 0.0, 0.0);
0255
0256   auto jetp4 = jet->p4();  //(rawp4 - lepp4*(1.0/jet->jecFactor("L1FastJet")))*(jet->pt()/rawp4.pt())+lepp4;
0257   auto ptratio = lepp4.pt() / jetp4.pt();
0258   auto ptrel = lepp4.Vect().Cross((jetp4).Vect().Unit()).R();
0259   auto ptrelinv = jetp4.Vect().Cross((lepp4).Vect().Unit()).R();
0260
0261   return std::tuple<float, float, float>(ptratio, ptrel, ptrelinv);
0262 }
0263
0264 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0265 template <typename T>
0266 void BJetEnergyRegressionVarProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0267   //The following says we do not know what parameters are allowed so do no validation
0268   // Please change this to state exactly what you do use, even if it is no parameters
0269   edm::ParameterSetDescription desc;
0270   desc.add<edm::InputTag>("src")->setComment("jet input collection");
0271   desc.add<edm::InputTag>("pvsrc")->setComment("primary vertex input collection");
0272   desc.add<edm::InputTag>("svsrc")->setComment("secondary vertex input collection");
0273   std::string modname;
0274   if (typeid(T) == typeid(pat::Jet))
0275     modname += "Jet";
0276   modname += "RegressionVarProducer";
0277   descriptions.add(modname, desc);
0278 }
0279
0280 typedef BJetEnergyRegressionVarProducer<pat::Jet> JetRegressionVarProducer;
0281
0282 //define this as a plug-in
0283 DEFINE_FWK_MODULE(JetRegressionVarProducer);