BJetEnergyRegressionVarProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    PhysicsTools/NanoAOD
// Class:      BJetEnergyRegressionVarProducer
//
//
// Original Author:  Marco Peruzzi
//         Created:  Tue, 05 Sep 2017 12:24:38 GMT
//
//

// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/StreamID.h"

#include "DataFormats/Math/interface/deltaR.h"

#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/PatCandidates/interface/Jet.h"

#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"

#include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
#include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
#include "RecoVertex/VertexPrimitives/interface/VertexState.h"

#include <vector>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/StreamID.h"

#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/VertexReco/interface/Vertex.h"

#include "DataFormats/Common/interface/View.h"
#include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"

#include "PhysicsTools/NanoAOD/interface/MatchingUtils.h"

//
// class declaration
//

template <typename T>
class BJetEnergyRegressionVarProducer : public edm::global::EDProducer<> {
public:
  explicit BJetEnergyRegressionVarProducer(const edm::ParameterSet& iConfig)
      : srcJet_(consumes<edm::View<pat::Jet>>(iConfig.getParameter<edm::InputTag>("src"))),
        srcVtx_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("pvsrc"))),
        srcSV_(consumes<edm::View<reco::VertexCompositePtrCandidate>>(iConfig.getParameter<edm::InputTag>("svsrc"))),
        srcGP_(consumes<std::vector<reco::GenParticle>>(iConfig.getParameter<edm::InputTag>("gpsrc"))) {
    //un prodotto da copiare
    produces<edm::ValueMap<float>>("leptonPtRel");
    produces<edm::ValueMap<float>>("leptonPtRatio");
    produces<edm::ValueMap<float>>("leptonPtRelInv");  //wrong variable?
    produces<edm::ValueMap<float>>("leptonPtRelv0");
    produces<edm::ValueMap<float>>("leptonPtRatiov0");
    produces<edm::ValueMap<float>>("leptonPtRelInvv0");  //v0 ~ heppy?
    produces<edm::ValueMap<float>>("leptonPt");
    produces<edm::ValueMap<int>>("leptonPdgId");
    produces<edm::ValueMap<float>>("leptonDeltaR");
    produces<edm::ValueMap<float>>("leadTrackPt");
    produces<edm::ValueMap<float>>("vtxPt");
    produces<edm::ValueMap<float>>("vtxMass");
    produces<edm::ValueMap<float>>("vtx3dL");
    produces<edm::ValueMap<float>>("vtx3deL");
    produces<edm::ValueMap<int>>("vtxNtrk");
    produces<edm::ValueMap<float>>("ptD");
    produces<edm::ValueMap<float>>("genPtwNu");
  }
  ~BJetEnergyRegressionVarProducer() override{};

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;

  std::tuple<float, float, float> calculatePtRatioRel(edm::Ptr<reco::Candidate> lep, edm::Ptr<pat::Jet> jet) const;
  std::tuple<float, float, float> calculatePtRatioRelSimple(edm::Ptr<reco::Candidate> lep,
                                                            edm::Ptr<pat::Jet> jet) const;  //old version?

  // ----------member data ---------------------------

  edm::EDGetTokenT<edm::View<pat::Jet>> srcJet_;
  edm::EDGetTokenT<std::vector<reco::Vertex>> srcVtx_;
  edm::EDGetTokenT<edm::View<reco::VertexCompositePtrCandidate>> srcSV_;
  edm::EDGetTokenT<std::vector<reco::GenParticle>> srcGP_;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// member functions
//

// ------------ method called to produce the data  ------------
template <typename T>
void BJetEnergyRegressionVarProducer<T>::produce(edm::StreamID streamID,
                                                 edm::Event& iEvent,
                                                 const edm::EventSetup& iSetup) const {
  edm::Handle<edm::View<pat::Jet>> srcJet;
  iEvent.getByToken(srcJet_, srcJet);
  edm::Handle<std::vector<reco::Vertex>> srcVtx;
  iEvent.getByToken(srcVtx_, srcVtx);
  edm::Handle<edm::View<reco::VertexCompositePtrCandidate>> srcSV;
  iEvent.getByToken(srcSV_, srcSV);
  edm::Handle<std::vector<reco::GenParticle>> srcGP;
  iEvent.getByToken(srcGP_, srcGP);

  unsigned int nJet = srcJet->size();
  //   unsigned int nLep = srcLep->size();

  std::vector<float> leptonPtRel(nJet, 0);
  std::vector<float> leptonPtRatio(nJet, 0);
  std::vector<float> leptonPtRelInv(nJet, 0);
  std::vector<float> leptonPtRel_v0(nJet, 0);
  std::vector<float> leptonPtRatio_v0(nJet, 0);
  std::vector<float> leptonPtRelInv_v0(nJet, 0);
  std::vector<int> leptonPdgId(nJet, 0);
  std::vector<float> leptonPt(nJet, 0);
  std::vector<float> leptonDeltaR(nJet, 0);
  std::vector<float> leadTrackPt(nJet, 0);
  std::vector<float> vtxPt(nJet, 0);
  std::vector<float> vtxMass(nJet, 0);
  std::vector<float> vtx3dL(nJet, 0);
  std::vector<float> vtx3deL(nJet, 0);
  std::vector<int> vtxNtrk(nJet, 0);
  std::vector<float> ptD(nJet, 0);
  std::vector<float> genPtwNu(nJet, 0);

  const auto& pv = (*srcVtx)[0];
  for (unsigned int ij = 0; ij < nJet; ij++) {
    auto jet = srcJet->ptrAt(ij);

    if (jet->genJet() != nullptr) {
      auto genp4 = jet->genJet()->p4();
      auto gep4wNu = genp4;
      for (const auto& gp : *srcGP) {
        if ((abs(gp.pdgId()) == 12 || abs(gp.pdgId()) == 14 || abs(gp.pdgId()) == 16) && gp.status() == 1) {
          if (reco::deltaR(genp4, gp.p4()) < 0.4) {
            //                    std::cout<<" from "<<gep4wNu.pt()<<std::endl;
            gep4wNu = gep4wNu + gp.p4();
            //                    std::cout<<" to "<<gep4wNu.pt()<<std::endl;
          }
        }
      }

      genPtwNu[ij] = gep4wNu.pt();
    }

    float ptMax = 0;
    float sumWeight = 0;
    float sumPt = 0;

    for (const auto& d : jet->daughterPtrVector()) {
      sumWeight += (d->pt()) * (d->pt());
      sumPt += d->pt();
      if (d->pt() > ptMax)
        ptMax = d->pt();
    }
    leadTrackPt[ij] = ptMax;
    ptD[ij] = (sumWeight > 0 ? sqrt(sumWeight) / sumPt : 0);

    //lepton properties
    float maxLepPt = 0;
    leptonPtRel[ij] = 0;

    for (const auto& d : jet->daughterPtrVector()) {
      if (abs(d->pdgId()) == 11 || abs(d->pdgId()) == 13) {
        if (d->pt() < maxLepPt)
          continue;
        auto res = calculatePtRatioRel(d, jet);
        leptonPtRatio[ij] = std::get<0>(res);
        leptonPtRel[ij] = std::get<1>(res);
        leptonPtRelInv[ij] = std::get<2>(res);
        auto res2 = calculatePtRatioRelSimple(d, jet);
        leptonPtRatio_v0[ij] = std::get<0>(res2);
        leptonPtRel_v0[ij] = std::get<1>(res2);
        leptonPtRelInv_v0[ij] = std::get<2>(res2);
        leptonPdgId[ij] = d->pdgId();
        leptonDeltaR[ij] = reco::deltaR(jet->p4(), d->p4());
        leptonPt[ij] = d->pt();
        maxLepPt = d->pt();
      }
    }

    //Fill vertex properties
    VertexDistance3D vdist;
    float maxFoundSignificance = 0;

    vtxPt[ij] = 0;
    vtxMass[ij] = 0;
    vtx3dL[ij] = 0;
    vtx3deL[ij] = 0;
    vtxNtrk[ij] = 0;

    for (const auto& sv : *srcSV) {
      GlobalVector flightDir(sv.vertex().x() - pv.x(), sv.vertex().y() - pv.y(), sv.vertex().z() - pv.z());
      GlobalVector jetDir(jet->px(), jet->py(), jet->pz());
      if (reco::deltaR2(flightDir, jetDir) < 0.09) {
        Measurement1D dl = vdist.distance(
            pv, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
        if (dl.significance() > maxFoundSignificance) {
          maxFoundSignificance = dl.significance();
          vtxPt[ij] = sv.pt();
          vtxMass[ij] = sv.p4().M();
          vtx3dL[ij] = dl.value();
          vtx3deL[ij] = dl.error();
          vtxNtrk[ij] = sv.numberOfSourceCandidatePtrs();
        }
      }
    }
  }

  std::unique_ptr<edm::ValueMap<float>> leptonPtRelV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerRel(*leptonPtRelV);
  fillerRel.insert(srcJet, leptonPtRel.begin(), leptonPtRel.end());
  fillerRel.fill();
  iEvent.put(std::move(leptonPtRelV), "leptonPtRel");

  std::unique_ptr<edm::ValueMap<float>> leptonPtRatioV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerRatio(*leptonPtRatioV);
  fillerRatio.insert(srcJet, leptonPtRatio.begin(), leptonPtRatio.end());
  fillerRatio.fill();
  iEvent.put(std::move(leptonPtRatioV), "leptonPtRatio");

  std::unique_ptr<edm::ValueMap<float>> leptonPtRelInvV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerRelInv(*leptonPtRelInvV);
  fillerRelInv.insert(srcJet, leptonPtRelInv.begin(), leptonPtRelInv.end());
  fillerRelInv.fill();
  iEvent.put(std::move(leptonPtRelInvV), "leptonPtRelInv");

  std::unique_ptr<edm::ValueMap<float>> leptonPtRelV_v0(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerRel_v0(*leptonPtRelV_v0);
  fillerRel_v0.insert(srcJet, leptonPtRel_v0.begin(), leptonPtRel_v0.end());
  fillerRel_v0.fill();
  iEvent.put(std::move(leptonPtRelV_v0), "leptonPtRelv0");

  std::unique_ptr<edm::ValueMap<float>> leptonPtRatioV_v0(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerRatio_v0(*leptonPtRatioV_v0);
  fillerRatio_v0.insert(srcJet, leptonPtRatio_v0.begin(), leptonPtRatio_v0.end());
  fillerRatio_v0.fill();
  iEvent.put(std::move(leptonPtRatioV_v0), "leptonPtRatiov0");

  std::unique_ptr<edm::ValueMap<float>> leptonPtRelInvV_v0(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerRelInv_v0(*leptonPtRelInvV_v0);
  fillerRelInv_v0.insert(srcJet, leptonPtRelInv_v0.begin(), leptonPtRelInv_v0.end());
  fillerRelInv_v0.fill();
  iEvent.put(std::move(leptonPtRelInvV_v0), "leptonPtRelInvv0");

  std::unique_ptr<edm::ValueMap<float>> leptonPtV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerLpt(*leptonPtV);
  fillerLpt.insert(srcJet, leptonPt.begin(), leptonPt.end());
  fillerLpt.fill();
  iEvent.put(std::move(leptonPtV), "leptonPt");

  std::unique_ptr<edm::ValueMap<float>> leptonDeltaRV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerLdR(*leptonDeltaRV);
  fillerLdR.insert(srcJet, leptonDeltaR.begin(), leptonDeltaR.end());
  fillerLdR.fill();
  iEvent.put(std::move(leptonDeltaRV), "leptonDeltaR");

  std::unique_ptr<edm::ValueMap<int>> leptonPtPdgIdV(new edm::ValueMap<int>());
  edm::ValueMap<int>::Filler fillerId(*leptonPtPdgIdV);
  fillerId.insert(srcJet, leptonPdgId.begin(), leptonPdgId.end());
  fillerId.fill();
  iEvent.put(std::move(leptonPtPdgIdV), "leptonPdgId");

  std::unique_ptr<edm::ValueMap<float>> leadTrackPtV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerLT(*leadTrackPtV);
  fillerLT.insert(srcJet, leadTrackPt.begin(), leadTrackPt.end());
  fillerLT.fill();
  iEvent.put(std::move(leadTrackPtV), "leadTrackPt");

  std::unique_ptr<edm::ValueMap<float>> vtxPtV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerVtxPt(*vtxPtV);
  fillerVtxPt.insert(srcJet, vtxPt.begin(), vtxPt.end());
  fillerVtxPt.fill();
  iEvent.put(std::move(vtxPtV), "vtxPt");

  std::unique_ptr<edm::ValueMap<float>> vtxMassV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerVtxMass(*vtxMassV);
  fillerVtxMass.insert(srcJet, vtxMass.begin(), vtxMass.end());
  fillerVtxMass.fill();
  iEvent.put(std::move(vtxMassV), "vtxMass");

  std::unique_ptr<edm::ValueMap<float>> vtx3dLV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerVtx3dL(*vtx3dLV);
  fillerVtx3dL.insert(srcJet, vtx3dL.begin(), vtx3dL.end());
  fillerVtx3dL.fill();
  iEvent.put(std::move(vtx3dLV), "vtx3dL");

  std::unique_ptr<edm::ValueMap<float>> vtx3deLV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerVtx3deL(*vtx3deLV);
  fillerVtx3deL.insert(srcJet, vtx3deL.begin(), vtx3deL.end());
  fillerVtx3deL.fill();
  iEvent.put(std::move(vtx3deLV), "vtx3deL");

  std::unique_ptr<edm::ValueMap<int>> vtxNtrkV(new edm::ValueMap<int>());
  edm::ValueMap<int>::Filler fillerVtxNT(*vtxNtrkV);
  fillerVtxNT.insert(srcJet, vtxNtrk.begin(), vtxNtrk.end());
  fillerVtxNT.fill();
  iEvent.put(std::move(vtxNtrkV), "vtxNtrk");

  std::unique_ptr<edm::ValueMap<float>> ptDV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerPtD(*ptDV);
  fillerPtD.insert(srcJet, ptD.begin(), ptD.end());
  fillerPtD.fill();
  iEvent.put(std::move(ptDV), "ptD");

  std::unique_ptr<edm::ValueMap<float>> genptV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillergenpt(*genptV);
  fillergenpt.insert(srcJet, genPtwNu.begin(), genPtwNu.end());
  fillergenpt.fill();
  iEvent.put(std::move(genptV), "genPtwNu");
}

template <typename T>
std::tuple<float, float, float> BJetEnergyRegressionVarProducer<T>::calculatePtRatioRel(edm::Ptr<reco::Candidate> lep,
                                                                                        edm::Ptr<pat::Jet> jet) const {
  auto rawp4 = jet->correctedP4("Uncorrected");
  auto lepp4 = lep->p4();

  if ((rawp4 - lepp4).R() < 1e-4)
    return std::tuple<float, float, float>(1.0, 0.0, 0.0);

  auto jetp4 = (rawp4 - lepp4 * (1.0 / jet->jecFactor("L1FastJet"))) * (jet->pt() / rawp4.pt()) + lepp4;
  auto ptratio = lepp4.pt() / jetp4.pt();
  auto ptrel = lepp4.Vect().Cross((jetp4 - lepp4).Vect().Unit()).R();
  auto ptrelinv = (jetp4 - lepp4).Vect().Cross((lepp4).Vect().Unit()).R();

  return std::tuple<float, float, float>(ptratio, ptrel, ptrelinv);
}

template <typename T>
std::tuple<float, float, float> BJetEnergyRegressionVarProducer<T>::calculatePtRatioRelSimple(
    edm::Ptr<reco::Candidate> lep, edm::Ptr<pat::Jet> jet) const {
  auto lepp4 = lep->p4();
  auto rawp4 = jet->correctedP4("Uncorrected");

  if ((rawp4 - lepp4).R() < 1e-4)
    return std::tuple<float, float, float>(1.0, 0.0, 0.0);

  auto jetp4 = jet->p4();  //(rawp4 - lepp4*(1.0/jet->jecFactor("L1FastJet")))*(jet->pt()/rawp4.pt())+lepp4;
  auto ptratio = lepp4.pt() / jetp4.pt();
  auto ptrel = lepp4.Vect().Cross((jetp4).Vect().Unit()).R();
  auto ptrelinv = jetp4.Vect().Cross((lepp4).Vect().Unit()).R();

  return std::tuple<float, float, float>(ptratio, ptrel, ptrelinv);
}

// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
template <typename T>
void BJetEnergyRegressionVarProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  //The following says we do not know what parameters are allowed so do no validation
  // Please change this to state exactly what you do use, even if it is no parameters
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("src")->setComment("jet input collection");
  //   desc.add<edm::InputTag>("musrc")->setComment("muons input collection");
  //   desc.add<edm::InputTag>("elesrc")->setComment("electrons input collection");
  desc.add<edm::InputTag>("pvsrc")->setComment("primary vertex input collection");
  desc.add<edm::InputTag>("svsrc")->setComment("secondary vertex input collection");
  desc.add<edm::InputTag>("gpsrc")->setComment("genparticles for nu recovery");
  std::string modname;
  if (typeid(T) == typeid(pat::Jet))
    modname += "Jet";
  modname += "RegressionVarProducer";
  descriptions.add(modname, desc);
}

typedef BJetEnergyRegressionVarProducer<pat::Jet> JetRegressionVarProducer;

//define this as a plug-in
DEFINE_FWK_MODULE(JetRegressionVarProducer);