Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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 
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/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/Math/interface/deltaR.h"
0033 
0034 #include "DataFormats/PatCandidates/interface/Muon.h"
0035 #include "DataFormats/PatCandidates/interface/Electron.h"
0036 #include "DataFormats/PatCandidates/interface/Jet.h"
0037 
0038 #include "DataFormats/VertexReco/interface/Vertex.h"
0039 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0040 
0041 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0042 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0043 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
0044 
0045 #include <vector>
0046 
0047 // user include files
0048 #include "FWCore/Framework/interface/Frameworkfwd.h"
0049 #include "FWCore/Framework/interface/global/EDProducer.h"
0050 
0051 #include "FWCore/Framework/interface/Event.h"
0052 #include "FWCore/Framework/interface/MakerMacros.h"
0053 
0054 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0055 #include "FWCore/Utilities/interface/StreamID.h"
0056 
0057 #include "DataFormats/PatCandidates/interface/Muon.h"
0058 #include "DataFormats/PatCandidates/interface/Jet.h"
0059 #include "DataFormats/PatCandidates/interface/Electron.h"
0060 #include "DataFormats/VertexReco/interface/Vertex.h"
0061 
0062 #include "DataFormats/Common/interface/View.h"
0063 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0064 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0065 
0066 #include "PhysicsTools/NanoAOD/interface/MatchingUtils.h"
0067 
0068 //
0069 // class declaration
0070 //
0071 
0072 template <typename T>
0073 class BJetEnergyRegressionVarProducer : public edm::global::EDProducer<> {
0074 public:
0075   explicit BJetEnergyRegressionVarProducer(const edm::ParameterSet& iConfig)
0076       : srcJet_(consumes<edm::View<pat::Jet>>(iConfig.getParameter<edm::InputTag>("src"))),
0077         srcVtx_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("pvsrc"))),
0078         srcSV_(consumes<edm::View<reco::VertexCompositePtrCandidate>>(iConfig.getParameter<edm::InputTag>("svsrc"))),
0079         srcGP_(consumes<std::vector<reco::GenParticle>>(iConfig.getParameter<edm::InputTag>("gpsrc"))) {
0080     //un prodotto da copiare
0081     produces<edm::ValueMap<float>>("leptonPtRel");
0082     produces<edm::ValueMap<float>>("leptonPtRatio");
0083     produces<edm::ValueMap<float>>("leptonPtRelInv");  //wrong variable?
0084     produces<edm::ValueMap<float>>("leptonPtRelv0");
0085     produces<edm::ValueMap<float>>("leptonPtRatiov0");
0086     produces<edm::ValueMap<float>>("leptonPtRelInvv0");  //v0 ~ heppy?
0087     produces<edm::ValueMap<float>>("leptonPt");
0088     produces<edm::ValueMap<int>>("leptonPdgId");
0089     produces<edm::ValueMap<float>>("leptonDeltaR");
0090     produces<edm::ValueMap<float>>("leadTrackPt");
0091     produces<edm::ValueMap<float>>("vtxPt");
0092     produces<edm::ValueMap<float>>("vtxMass");
0093     produces<edm::ValueMap<float>>("vtx3dL");
0094     produces<edm::ValueMap<float>>("vtx3deL");
0095     produces<edm::ValueMap<int>>("vtxNtrk");
0096     produces<edm::ValueMap<float>>("ptD");
0097     produces<edm::ValueMap<float>>("genPtwNu");
0098   }
0099   ~BJetEnergyRegressionVarProducer() override{};
0100 
0101   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0102 
0103 private:
0104   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0105 
0106   std::tuple<float, float, float> calculatePtRatioRel(edm::Ptr<reco::Candidate> lep, edm::Ptr<pat::Jet> jet) const;
0107   std::tuple<float, float, float> calculatePtRatioRelSimple(edm::Ptr<reco::Candidate> lep,
0108                                                             edm::Ptr<pat::Jet> jet) const;  //old version?
0109 
0110   // ----------member data ---------------------------
0111 
0112   edm::EDGetTokenT<edm::View<pat::Jet>> srcJet_;
0113   edm::EDGetTokenT<std::vector<reco::Vertex>> srcVtx_;
0114   edm::EDGetTokenT<edm::View<reco::VertexCompositePtrCandidate>> srcSV_;
0115   edm::EDGetTokenT<std::vector<reco::GenParticle>> srcGP_;
0116 };
0117 
0118 //
0119 // constants, enums and typedefs
0120 //
0121 
0122 //
0123 // static data member definitions
0124 //
0125 
0126 //
0127 // member functions
0128 //
0129 
0130 // ------------ method called to produce the data  ------------
0131 template <typename T>
0132 void BJetEnergyRegressionVarProducer<T>::produce(edm::StreamID streamID,
0133                                                  edm::Event& iEvent,
0134                                                  const edm::EventSetup& iSetup) const {
0135   edm::Handle<edm::View<pat::Jet>> srcJet;
0136   iEvent.getByToken(srcJet_, srcJet);
0137   edm::Handle<std::vector<reco::Vertex>> srcVtx;
0138   iEvent.getByToken(srcVtx_, srcVtx);
0139   edm::Handle<edm::View<reco::VertexCompositePtrCandidate>> srcSV;
0140   iEvent.getByToken(srcSV_, srcSV);
0141   edm::Handle<std::vector<reco::GenParticle>> srcGP;
0142   iEvent.getByToken(srcGP_, srcGP);
0143 
0144   unsigned int nJet = srcJet->size();
0145   //   unsigned int nLep = srcLep->size();
0146 
0147   std::vector<float> leptonPtRel(nJet, 0);
0148   std::vector<float> leptonPtRatio(nJet, 0);
0149   std::vector<float> leptonPtRelInv(nJet, 0);
0150   std::vector<float> leptonPtRel_v0(nJet, 0);
0151   std::vector<float> leptonPtRatio_v0(nJet, 0);
0152   std::vector<float> leptonPtRelInv_v0(nJet, 0);
0153   std::vector<int> leptonPdgId(nJet, 0);
0154   std::vector<float> leptonPt(nJet, 0);
0155   std::vector<float> leptonDeltaR(nJet, 0);
0156   std::vector<float> leadTrackPt(nJet, 0);
0157   std::vector<float> vtxPt(nJet, 0);
0158   std::vector<float> vtxMass(nJet, 0);
0159   std::vector<float> vtx3dL(nJet, 0);
0160   std::vector<float> vtx3deL(nJet, 0);
0161   std::vector<int> vtxNtrk(nJet, 0);
0162   std::vector<float> ptD(nJet, 0);
0163   std::vector<float> genPtwNu(nJet, 0);
0164 
0165   const auto& pv = (*srcVtx)[0];
0166   for (unsigned int ij = 0; ij < nJet; ij++) {
0167     auto jet = srcJet->ptrAt(ij);
0168 
0169     if (jet->genJet() != nullptr) {
0170       auto genp4 = jet->genJet()->p4();
0171       auto gep4wNu = genp4;
0172       for (const auto& gp : *srcGP) {
0173         if ((abs(gp.pdgId()) == 12 || abs(gp.pdgId()) == 14 || abs(gp.pdgId()) == 16) && gp.status() == 1) {
0174           if (reco::deltaR(genp4, gp.p4()) < 0.4) {
0175             //                    std::cout<<" from "<<gep4wNu.pt()<<std::endl;
0176             gep4wNu = gep4wNu + gp.p4();
0177             //                    std::cout<<" to "<<gep4wNu.pt()<<std::endl;
0178           }
0179         }
0180       }
0181 
0182       genPtwNu[ij] = gep4wNu.pt();
0183     }
0184 
0185     float ptMax = 0;
0186     float sumWeight = 0;
0187     float sumPt = 0;
0188 
0189     for (const auto& d : jet->daughterPtrVector()) {
0190       sumWeight += (d->pt()) * (d->pt());
0191       sumPt += d->pt();
0192       if (d->pt() > ptMax)
0193         ptMax = d->pt();
0194     }
0195     leadTrackPt[ij] = ptMax;
0196     ptD[ij] = (sumWeight > 0 ? sqrt(sumWeight) / sumPt : 0);
0197 
0198     //lepton properties
0199     float maxLepPt = 0;
0200     leptonPtRel[ij] = 0;
0201 
0202     for (const auto& d : jet->daughterPtrVector()) {
0203       if (abs(d->pdgId()) == 11 || abs(d->pdgId()) == 13) {
0204         if (d->pt() < maxLepPt)
0205           continue;
0206         auto res = calculatePtRatioRel(d, jet);
0207         leptonPtRatio[ij] = std::get<0>(res);
0208         leptonPtRel[ij] = std::get<1>(res);
0209         leptonPtRelInv[ij] = std::get<2>(res);
0210         auto res2 = calculatePtRatioRelSimple(d, jet);
0211         leptonPtRatio_v0[ij] = std::get<0>(res2);
0212         leptonPtRel_v0[ij] = std::get<1>(res2);
0213         leptonPtRelInv_v0[ij] = std::get<2>(res2);
0214         leptonPdgId[ij] = d->pdgId();
0215         leptonDeltaR[ij] = reco::deltaR(jet->p4(), d->p4());
0216         leptonPt[ij] = d->pt();
0217         maxLepPt = d->pt();
0218       }
0219     }
0220 
0221     //Fill vertex properties
0222     VertexDistance3D vdist;
0223     float maxFoundSignificance = 0;
0224 
0225     vtxPt[ij] = 0;
0226     vtxMass[ij] = 0;
0227     vtx3dL[ij] = 0;
0228     vtx3deL[ij] = 0;
0229     vtxNtrk[ij] = 0;
0230 
0231     for (const auto& sv : *srcSV) {
0232       GlobalVector flightDir(sv.vertex().x() - pv.x(), sv.vertex().y() - pv.y(), sv.vertex().z() - pv.z());
0233       GlobalVector jetDir(jet->px(), jet->py(), jet->pz());
0234       if (reco::deltaR2(flightDir, jetDir) < 0.09) {
0235         Measurement1D dl = vdist.distance(
0236             pv, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0237         if (dl.significance() > maxFoundSignificance) {
0238           maxFoundSignificance = dl.significance();
0239           vtxPt[ij] = sv.pt();
0240           vtxMass[ij] = sv.p4().M();
0241           vtx3dL[ij] = dl.value();
0242           vtx3deL[ij] = dl.error();
0243           vtxNtrk[ij] = sv.numberOfSourceCandidatePtrs();
0244         }
0245       }
0246     }
0247   }
0248 
0249   std::unique_ptr<edm::ValueMap<float>> leptonPtRelV(new edm::ValueMap<float>());
0250   edm::ValueMap<float>::Filler fillerRel(*leptonPtRelV);
0251   fillerRel.insert(srcJet, leptonPtRel.begin(), leptonPtRel.end());
0252   fillerRel.fill();
0253   iEvent.put(std::move(leptonPtRelV), "leptonPtRel");
0254 
0255   std::unique_ptr<edm::ValueMap<float>> leptonPtRatioV(new edm::ValueMap<float>());
0256   edm::ValueMap<float>::Filler fillerRatio(*leptonPtRatioV);
0257   fillerRatio.insert(srcJet, leptonPtRatio.begin(), leptonPtRatio.end());
0258   fillerRatio.fill();
0259   iEvent.put(std::move(leptonPtRatioV), "leptonPtRatio");
0260 
0261   std::unique_ptr<edm::ValueMap<float>> leptonPtRelInvV(new edm::ValueMap<float>());
0262   edm::ValueMap<float>::Filler fillerRelInv(*leptonPtRelInvV);
0263   fillerRelInv.insert(srcJet, leptonPtRelInv.begin(), leptonPtRelInv.end());
0264   fillerRelInv.fill();
0265   iEvent.put(std::move(leptonPtRelInvV), "leptonPtRelInv");
0266 
0267   std::unique_ptr<edm::ValueMap<float>> leptonPtRelV_v0(new edm::ValueMap<float>());
0268   edm::ValueMap<float>::Filler fillerRel_v0(*leptonPtRelV_v0);
0269   fillerRel_v0.insert(srcJet, leptonPtRel_v0.begin(), leptonPtRel_v0.end());
0270   fillerRel_v0.fill();
0271   iEvent.put(std::move(leptonPtRelV_v0), "leptonPtRelv0");
0272 
0273   std::unique_ptr<edm::ValueMap<float>> leptonPtRatioV_v0(new edm::ValueMap<float>());
0274   edm::ValueMap<float>::Filler fillerRatio_v0(*leptonPtRatioV_v0);
0275   fillerRatio_v0.insert(srcJet, leptonPtRatio_v0.begin(), leptonPtRatio_v0.end());
0276   fillerRatio_v0.fill();
0277   iEvent.put(std::move(leptonPtRatioV_v0), "leptonPtRatiov0");
0278 
0279   std::unique_ptr<edm::ValueMap<float>> leptonPtRelInvV_v0(new edm::ValueMap<float>());
0280   edm::ValueMap<float>::Filler fillerRelInv_v0(*leptonPtRelInvV_v0);
0281   fillerRelInv_v0.insert(srcJet, leptonPtRelInv_v0.begin(), leptonPtRelInv_v0.end());
0282   fillerRelInv_v0.fill();
0283   iEvent.put(std::move(leptonPtRelInvV_v0), "leptonPtRelInvv0");
0284 
0285   std::unique_ptr<edm::ValueMap<float>> leptonPtV(new edm::ValueMap<float>());
0286   edm::ValueMap<float>::Filler fillerLpt(*leptonPtV);
0287   fillerLpt.insert(srcJet, leptonPt.begin(), leptonPt.end());
0288   fillerLpt.fill();
0289   iEvent.put(std::move(leptonPtV), "leptonPt");
0290 
0291   std::unique_ptr<edm::ValueMap<float>> leptonDeltaRV(new edm::ValueMap<float>());
0292   edm::ValueMap<float>::Filler fillerLdR(*leptonDeltaRV);
0293   fillerLdR.insert(srcJet, leptonDeltaR.begin(), leptonDeltaR.end());
0294   fillerLdR.fill();
0295   iEvent.put(std::move(leptonDeltaRV), "leptonDeltaR");
0296 
0297   std::unique_ptr<edm::ValueMap<int>> leptonPtPdgIdV(new edm::ValueMap<int>());
0298   edm::ValueMap<int>::Filler fillerId(*leptonPtPdgIdV);
0299   fillerId.insert(srcJet, leptonPdgId.begin(), leptonPdgId.end());
0300   fillerId.fill();
0301   iEvent.put(std::move(leptonPtPdgIdV), "leptonPdgId");
0302 
0303   std::unique_ptr<edm::ValueMap<float>> leadTrackPtV(new edm::ValueMap<float>());
0304   edm::ValueMap<float>::Filler fillerLT(*leadTrackPtV);
0305   fillerLT.insert(srcJet, leadTrackPt.begin(), leadTrackPt.end());
0306   fillerLT.fill();
0307   iEvent.put(std::move(leadTrackPtV), "leadTrackPt");
0308 
0309   std::unique_ptr<edm::ValueMap<float>> vtxPtV(new edm::ValueMap<float>());
0310   edm::ValueMap<float>::Filler fillerVtxPt(*vtxPtV);
0311   fillerVtxPt.insert(srcJet, vtxPt.begin(), vtxPt.end());
0312   fillerVtxPt.fill();
0313   iEvent.put(std::move(vtxPtV), "vtxPt");
0314 
0315   std::unique_ptr<edm::ValueMap<float>> vtxMassV(new edm::ValueMap<float>());
0316   edm::ValueMap<float>::Filler fillerVtxMass(*vtxMassV);
0317   fillerVtxMass.insert(srcJet, vtxMass.begin(), vtxMass.end());
0318   fillerVtxMass.fill();
0319   iEvent.put(std::move(vtxMassV), "vtxMass");
0320 
0321   std::unique_ptr<edm::ValueMap<float>> vtx3dLV(new edm::ValueMap<float>());
0322   edm::ValueMap<float>::Filler fillerVtx3dL(*vtx3dLV);
0323   fillerVtx3dL.insert(srcJet, vtx3dL.begin(), vtx3dL.end());
0324   fillerVtx3dL.fill();
0325   iEvent.put(std::move(vtx3dLV), "vtx3dL");
0326 
0327   std::unique_ptr<edm::ValueMap<float>> vtx3deLV(new edm::ValueMap<float>());
0328   edm::ValueMap<float>::Filler fillerVtx3deL(*vtx3deLV);
0329   fillerVtx3deL.insert(srcJet, vtx3deL.begin(), vtx3deL.end());
0330   fillerVtx3deL.fill();
0331   iEvent.put(std::move(vtx3deLV), "vtx3deL");
0332 
0333   std::unique_ptr<edm::ValueMap<int>> vtxNtrkV(new edm::ValueMap<int>());
0334   edm::ValueMap<int>::Filler fillerVtxNT(*vtxNtrkV);
0335   fillerVtxNT.insert(srcJet, vtxNtrk.begin(), vtxNtrk.end());
0336   fillerVtxNT.fill();
0337   iEvent.put(std::move(vtxNtrkV), "vtxNtrk");
0338 
0339   std::unique_ptr<edm::ValueMap<float>> ptDV(new edm::ValueMap<float>());
0340   edm::ValueMap<float>::Filler fillerPtD(*ptDV);
0341   fillerPtD.insert(srcJet, ptD.begin(), ptD.end());
0342   fillerPtD.fill();
0343   iEvent.put(std::move(ptDV), "ptD");
0344 
0345   std::unique_ptr<edm::ValueMap<float>> genptV(new edm::ValueMap<float>());
0346   edm::ValueMap<float>::Filler fillergenpt(*genptV);
0347   fillergenpt.insert(srcJet, genPtwNu.begin(), genPtwNu.end());
0348   fillergenpt.fill();
0349   iEvent.put(std::move(genptV), "genPtwNu");
0350 }
0351 
0352 template <typename T>
0353 std::tuple<float, float, float> BJetEnergyRegressionVarProducer<T>::calculatePtRatioRel(edm::Ptr<reco::Candidate> lep,
0354                                                                                         edm::Ptr<pat::Jet> jet) const {
0355   auto rawp4 = jet->correctedP4("Uncorrected");
0356   auto lepp4 = lep->p4();
0357 
0358   if ((rawp4 - lepp4).R() < 1e-4)
0359     return std::tuple<float, float, float>(1.0, 0.0, 0.0);
0360 
0361   auto jetp4 = (rawp4 - lepp4 * (1.0 / jet->jecFactor("L1FastJet"))) * (jet->pt() / rawp4.pt()) + lepp4;
0362   auto ptratio = lepp4.pt() / jetp4.pt();
0363   auto ptrel = lepp4.Vect().Cross((jetp4 - lepp4).Vect().Unit()).R();
0364   auto ptrelinv = (jetp4 - lepp4).Vect().Cross((lepp4).Vect().Unit()).R();
0365 
0366   return std::tuple<float, float, float>(ptratio, ptrel, ptrelinv);
0367 }
0368 
0369 template <typename T>
0370 std::tuple<float, float, float> BJetEnergyRegressionVarProducer<T>::calculatePtRatioRelSimple(
0371     edm::Ptr<reco::Candidate> lep, edm::Ptr<pat::Jet> jet) const {
0372   auto lepp4 = lep->p4();
0373   auto rawp4 = jet->correctedP4("Uncorrected");
0374 
0375   if ((rawp4 - lepp4).R() < 1e-4)
0376     return std::tuple<float, float, float>(1.0, 0.0, 0.0);
0377 
0378   auto jetp4 = jet->p4();  //(rawp4 - lepp4*(1.0/jet->jecFactor("L1FastJet")))*(jet->pt()/rawp4.pt())+lepp4;
0379   auto ptratio = lepp4.pt() / jetp4.pt();
0380   auto ptrel = lepp4.Vect().Cross((jetp4).Vect().Unit()).R();
0381   auto ptrelinv = jetp4.Vect().Cross((lepp4).Vect().Unit()).R();
0382 
0383   return std::tuple<float, float, float>(ptratio, ptrel, ptrelinv);
0384 }
0385 
0386 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0387 template <typename T>
0388 void BJetEnergyRegressionVarProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0389   //The following says we do not know what parameters are allowed so do no validation
0390   // Please change this to state exactly what you do use, even if it is no parameters
0391   edm::ParameterSetDescription desc;
0392   desc.add<edm::InputTag>("src")->setComment("jet input collection");
0393   //   desc.add<edm::InputTag>("musrc")->setComment("muons input collection");
0394   //   desc.add<edm::InputTag>("elesrc")->setComment("electrons input collection");
0395   desc.add<edm::InputTag>("pvsrc")->setComment("primary vertex input collection");
0396   desc.add<edm::InputTag>("svsrc")->setComment("secondary vertex input collection");
0397   desc.add<edm::InputTag>("gpsrc")->setComment("genparticles for nu recovery");
0398   std::string modname;
0399   if (typeid(T) == typeid(pat::Jet))
0400     modname += "Jet";
0401   modname += "RegressionVarProducer";
0402   descriptions.add(modname, desc);
0403 }
0404 
0405 typedef BJetEnergyRegressionVarProducer<pat::Jet> JetRegressionVarProducer;
0406 
0407 //define this as a plug-in
0408 DEFINE_FWK_MODULE(JetRegressionVarProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!