PixelJetPuId

Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
// -*- C++ -*-
//
// Package:    PixelJetPuId
// Class:      PixelJetPuId
//
/**\class PixelJetPuId PixelJetPuId.cc RecoBTag/PixelJetPuId/src/PixelJetPuId.cc

Description:
The PixelJetPuId module select all the pixel tracks compatible with a jet.
If the sum of the tracks momentum is under a threshold the jet is tagged as "PUjets".

Implementation:
[Notes on implementation]
 */
//
// 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);