PixelJetPuId.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    PixelJetPuId
// Class:      PixelJetPuId
//
//
// Original Author:  Silvio DONATO
//         Created:  Wed Dec 18 10:05:40 CET 2013
//
//

// system include files
#include <memory>
#include <algorithm>

// 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/Framework/interface/makeRefToBaseProdFrom.h"

#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"

#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"

#include "TrackingTools/IPTools/interface/IPTools.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"

#include "DataFormats/Math/interface/deltaR.h"
#include "DataFormats/BTauReco/interface/JetTag.h"
#include "DataFormats/Common/interface/RefToBase.h"
//
// class declaration
//

class PixelJetPuId : public edm::global::EDProducer<> {
public:
  PixelJetPuId(const edm::ParameterSet&);
  ~PixelJetPuId() override;

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

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

  // ----------member data ---------------------------
  edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> const transientTrackRecordToken_;

  edm::InputTag m_primaryVertex;
  edm::InputTag m_tracks;
  edm::InputTag m_jets;
  edm::EDGetTokenT<std::vector<reco::Track> > tracksToken;
  edm::EDGetTokenT<edm::View<reco::CaloJet> > jetsToken;
  edm::EDGetTokenT<edm::View<reco::Jet> > generaljetsToken;
  edm::EDGetTokenT<reco::VertexCollection> primaryVertexToken;

  double m_MinTrackPt;
  double m_MaxTrackChi2;
  double m_MaxTrackDistanceToJet;

  bool m_fwjets;
  double m_mineta_fwjets;
  double m_minet_fwjets;

  double m_MinGoodJetTrackPt;
  double m_MinGoodJetTrackPtRatio;
};

//
// constructors and destructor
//
PixelJetPuId::PixelJetPuId(const edm::ParameterSet& iConfig)
    : transientTrackRecordToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) {
  //InputTag
  m_tracks = iConfig.getParameter<edm::InputTag>("tracks");
  tracksToken = consumes<std::vector<reco::Track> >(m_tracks);
  m_jets = iConfig.getParameter<edm::InputTag>("jets");
  jetsToken = consumes<edm::View<reco::CaloJet> >(m_jets);
  generaljetsToken = consumes<edm::View<reco::Jet> >(m_jets);
  m_primaryVertex = iConfig.getParameter<edm::InputTag>("primaryVertex");
  primaryVertexToken = consumes<reco::VertexCollection>(m_primaryVertex);

  //Tracks Selection
  m_MinTrackPt = iConfig.getParameter<double>("MinTrackPt");
  m_MaxTrackDistanceToJet = iConfig.getParameter<double>("MaxTrackDistanceToJet");
  m_MaxTrackChi2 = iConfig.getParameter<double>("MaxTrackChi2");

  //A jet is defined as a signal jet if Sum(trackPt) > minPt or Sum(comp.trackPt)/CaloJetPt > minPtRatio
  m_MinGoodJetTrackPt = iConfig.getParameter<double>("MinGoodJetTrackPt");
  m_MinGoodJetTrackPtRatio = iConfig.getParameter<double>("MinGoodJetTrackPtRatio");

  m_fwjets = iConfig.getParameter<bool>("UseForwardJetsAsNoPU");
  m_mineta_fwjets = iConfig.getParameter<double>("MinEtaForwardJets");
  m_minet_fwjets = iConfig.getParameter<double>("MinEtForwardJets");

  produces<std::vector<reco::CaloJet> >();
  produces<std::vector<reco::CaloJet> >("PUjets");
  produces<reco::JetTagCollection>();
}

PixelJetPuId::~PixelJetPuId() = default;

void PixelJetPuId::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("jets", edm::InputTag("hltCaloJetL1FastJetCorrected"));
  desc.add<edm::InputTag>("tracks", edm::InputTag("hltPixelTracksNoPU"));
  desc.add<edm::InputTag>("primaryVertex", edm::InputTag("hltFastPVPixelVertices"));
  desc.add<double>("MinGoodJetTrackPtRatio", 0.045);
  desc.add<double>("MinGoodJetTrackPt", 1.8);
  desc.add<double>("MaxTrackDistanceToJet", 0.04);
  desc.add<double>("MinTrackPt", 0.6);
  desc.add<double>("MaxTrackChi2", 20.);
  desc.add<bool>("UseForwardJetsAsNoPU", true);
  desc.add<double>("MinEtaForwardJets", 2.4);
  desc.add<double>("MinEtForwardJets", 40.);
  descriptions.add("pixelJetPuId", desc);
}

//
// member functions
//

// ------------ method called on each new Event  ------------
void PixelJetPuId::produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
  using namespace edm;
  std::unique_ptr<std::vector<reco::CaloJet> > pOut(new std::vector<reco::CaloJet>);
  std::unique_ptr<std::vector<reco::CaloJet> > pOut_PUjets(new std::vector<reco::CaloJet>);
  std::unique_ptr<reco::JetTagCollection> pOut_jetTagCollection(new reco::JetTagCollection);

  //get tracks
  Handle<std::vector<reco::Track> > tracks;
  iEvent.getByToken(tracksToken, tracks);
  uint const asize = std::max(1u, (uint)tracks->size());
  float teta[asize], tphi[asize];
  uint tsize = 0;
  for (auto const& tr : *tracks) {
    teta[tsize] = tr.eta();
    tphi[tsize] = tr.phi();
    ++tsize;
  }

  //get jets
  Handle<edm::View<reco::CaloJet> > jets;
  iEvent.getByToken(jetsToken, jets);

  Handle<edm::View<reco::Jet> > generaljets;
  iEvent.getByToken(generaljetsToken, generaljets);

  //get primary vertices
  Handle<reco::VertexCollection> primaryVertex;
  iEvent.getByToken(primaryVertexToken, primaryVertex);

  //get Transient Track Builder
  auto const& builder = iSetup.getHandle(transientTrackRecordToken_);

  //init JetTagCollection
  if (!generaljets.product()->empty()) {
    edm::RefToBase<reco::Jet> jj = edm::RefToBase<reco::Jet>(generaljets, 0);
    pOut_jetTagCollection = std::make_unique<reco::JetTagCollection>(edm::makeRefToBaseProdFrom(jj, iEvent));
  }

  //loop on trackIPTagInfos
  if (!primaryVertex->empty()) {
    const reco::Vertex* pv = &*primaryVertex->begin();
    //loop on jets
    for (edm::View<reco::CaloJet>::const_iterator itJet = jets->begin(); itJet != jets->end(); itJet++) {
      math::XYZVector jetMomentum = itJet->momentum();
      GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());

      math::XYZVector trMomentum;

      if (fabs(itJet->eta()) > m_mineta_fwjets) {
        if ((m_fwjets) && (itJet->et() > m_minet_fwjets))
          pOut->push_back(*itJet);  // fill forward jet as signal jet
      } else {
        //loop on tracks
        auto itTrack = tracks->begin();
        for (unsigned int i = 0; i < tsize; ++i) {
          float deltaR2 = reco::deltaR2(itJet->eta(), itJet->phi(), teta[i], tphi[i]);
          if (deltaR2 < 0.25) {
            reco::TransientTrack transientTrack = builder->build(*itTrack);
            float jetTrackDistance = -((IPTools::jetTrackDistance(transientTrack, direction, *pv)).second).value();

            //select the tracks compabible with the jet
            if ((itTrack->pt() > m_MinTrackPt) && (itTrack->normalizedChi2() < m_MaxTrackChi2) &&
                (jetTrackDistance < m_MaxTrackDistanceToJet)) {
              trMomentum += itTrack->momentum();  //calculate the Sum(trackPt)
            }
          }
          itTrack++;
        }
        //if Sum(comp.trackPt)/CaloJetPt > minPtRatio or Sum(trackPt) > minPt  the jet is a signal jet
        if (trMomentum.rho() / jetMomentum.rho() > m_MinGoodJetTrackPtRatio || trMomentum.rho() > m_MinGoodJetTrackPt) {
          pOut->push_back(*itJet);  // fill it as signal jet
        } else                      //else it is a PUjet
        {
          pOut_PUjets->push_back(*itJet);  // fill it as PUjets
        }
      }
      RefToBase<reco::Jet> jRef(generaljets, itJet - jets->begin());
      (*pOut_jetTagCollection)[jRef] = trMomentum.rho();  // fill jetTagCollection
    }
  }
  iEvent.put(std::move(pOut));
  iEvent.put(std::move(pOut_PUjets), "PUjets");
  iEvent.put(std::move(pOut_jetTagCollection));
}

// declare this class as a framework plugin
DEFINE_FWK_MODULE(PixelJetPuId);