JetMET/plugins/PixelJetPuId.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PixelJetPuId
0004 // Class:      PixelJetPuId
0005 //
0006 /**\class PixelJetPuId PixelJetPuId.cc RecoBTag/PixelJetPuId/src/PixelJetPuId.cc
0007
0008 Description:
0009 The PixelJetPuId module select all the pixel tracks compatible with a jet.
0010 If the sum of the tracks momentum is under a threshold the jet is tagged as "PUjets".
0011
0012 Implementation:
0013 [Notes on implementation]
0014  */
0015 //
0016 // Original Author:  Silvio DONATO
0017 //         Created:  Wed Dec 18 10:05:40 CET 2013
0018 //
0019 //
0020
0021 // system include files
0022 #include <memory>
0023 #include <algorithm>
0024
0025 // user include files
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/global/EDProducer.h"
0028
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 #include "FWCore/Framework/interface/makeRefToBaseProdFrom.h"
0032
0033 #include "FWCore/Utilities/interface/InputTag.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0037
0038 #include "DataFormats/TrackReco/interface/Track.h"
0039 #include "DataFormats/JetReco/interface/CaloJet.h"
0040 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0041 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0042 #include "DataFormats/VertexReco/interface/Vertex.h"
0043
0044 #include "TrackingTools/IPTools/interface/IPTools.h"
0045 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0046 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0047 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0048
0049 #include "DataFormats/Math/interface/deltaR.h"
0050 #include "DataFormats/BTauReco/interface/JetTag.h"
0051 #include "DataFormats/Common/interface/RefToBase.h"
0052 //
0053 // class declaration
0054 //
0055
0056 class PixelJetPuId : public edm::global::EDProducer<> {
0057 public:
0058   PixelJetPuId(const edm::ParameterSet&);
0059   ~PixelJetPuId() override;
0060
0061   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0062
0063 private:
0064   void produce(edm::StreamID sid, edm::Event&, const edm::EventSetup&) const override;
0065
0066   // ----------member data ---------------------------
0067   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> const transientTrackRecordToken_;
0068
0069   edm::InputTag m_primaryVertex;
0070   edm::InputTag m_tracks;
0071   edm::InputTag m_jets;
0072   edm::EDGetTokenT<std::vector<reco::Track> > tracksToken;
0073   edm::EDGetTokenT<edm::View<reco::CaloJet> > jetsToken;
0074   edm::EDGetTokenT<edm::View<reco::Jet> > generaljetsToken;
0075   edm::EDGetTokenT<reco::VertexCollection> primaryVertexToken;
0076
0077   double m_MinTrackPt;
0078   double m_MaxTrackChi2;
0079   double m_MaxTrackDistanceToJet;
0080
0081   bool m_fwjets;
0082   double m_mineta_fwjets;
0083   double m_minet_fwjets;
0084
0085   double m_MinGoodJetTrackPt;
0086   double m_MinGoodJetTrackPtRatio;
0087 };
0088
0089 //
0090 // constructors and destructor
0091 //
0092 PixelJetPuId::PixelJetPuId(const edm::ParameterSet& iConfig)
0093     : transientTrackRecordToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) {
0094   //InputTag
0095   m_tracks = iConfig.getParameter<edm::InputTag>("tracks");
0096   tracksToken = consumes<std::vector<reco::Track> >(m_tracks);
0097   m_jets = iConfig.getParameter<edm::InputTag>("jets");
0098   jetsToken = consumes<edm::View<reco::CaloJet> >(m_jets);
0099   generaljetsToken = consumes<edm::View<reco::Jet> >(m_jets);
0100   m_primaryVertex = iConfig.getParameter<edm::InputTag>("primaryVertex");
0101   primaryVertexToken = consumes<reco::VertexCollection>(m_primaryVertex);
0102
0103   //Tracks Selection
0104   m_MinTrackPt = iConfig.getParameter<double>("MinTrackPt");
0105   m_MaxTrackDistanceToJet = iConfig.getParameter<double>("MaxTrackDistanceToJet");
0106   m_MaxTrackChi2 = iConfig.getParameter<double>("MaxTrackChi2");
0107
0108   //A jet is defined as a signal jet if Sum(trackPt) > minPt or Sum(comp.trackPt)/CaloJetPt > minPtRatio
0109   m_MinGoodJetTrackPt = iConfig.getParameter<double>("MinGoodJetTrackPt");
0110   m_MinGoodJetTrackPtRatio = iConfig.getParameter<double>("MinGoodJetTrackPtRatio");
0111
0112   m_fwjets = iConfig.getParameter<bool>("UseForwardJetsAsNoPU");
0113   m_mineta_fwjets = iConfig.getParameter<double>("MinEtaForwardJets");
0114   m_minet_fwjets = iConfig.getParameter<double>("MinEtForwardJets");
0115
0116   produces<std::vector<reco::CaloJet> >();
0117   produces<std::vector<reco::CaloJet> >("PUjets");
0118   produces<reco::JetTagCollection>();
0119 }
0120
0121 PixelJetPuId::~PixelJetPuId() = default;
0122
0123 void PixelJetPuId::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0124   edm::ParameterSetDescription desc;
0125   desc.add<edm::InputTag>("jets", edm::InputTag("hltCaloJetL1FastJetCorrected"));
0126   desc.add<edm::InputTag>("tracks", edm::InputTag("hltPixelTracksNoPU"));
0127   desc.add<edm::InputTag>("primaryVertex", edm::InputTag("hltFastPVPixelVertices"));
0128   desc.add<double>("MinGoodJetTrackPtRatio", 0.045);
0129   desc.add<double>("MinGoodJetTrackPt", 1.8);
0130   desc.add<double>("MaxTrackDistanceToJet", 0.04);
0131   desc.add<double>("MinTrackPt", 0.6);
0132   desc.add<double>("MaxTrackChi2", 20.);
0133   desc.add<bool>("UseForwardJetsAsNoPU", true);
0134   desc.add<double>("MinEtaForwardJets", 2.4);
0135   desc.add<double>("MinEtForwardJets", 40.);
0136   descriptions.add("pixelJetPuId", desc);
0137 }
0138
0139 //
0140 // member functions
0141 //
0142
0143 // ------------ method called on each new Event  ------------
0144 void PixelJetPuId::produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0145   using namespace edm;
0146   std::unique_ptr<std::vector<reco::CaloJet> > pOut(new std::vector<reco::CaloJet>);
0147   std::unique_ptr<std::vector<reco::CaloJet> > pOut_PUjets(new std::vector<reco::CaloJet>);
0148   std::unique_ptr<reco::JetTagCollection> pOut_jetTagCollection(new reco::JetTagCollection);
0149
0150   //get tracks
0151   Handle<std::vector<reco::Track> > tracks;
0152   iEvent.getByToken(tracksToken, tracks);
0153   uint const asize = std::max(1u, (uint)tracks->size());
0154   float teta[asize], tphi[asize];
0155   uint tsize = 0;
0156   for (auto const& tr : *tracks) {
0157     teta[tsize] = tr.eta();
0158     tphi[tsize] = tr.phi();
0159     ++tsize;
0160   }
0161
0162   //get jets
0163   Handle<edm::View<reco::CaloJet> > jets;
0164   iEvent.getByToken(jetsToken, jets);
0165
0166   Handle<edm::View<reco::Jet> > generaljets;
0167   iEvent.getByToken(generaljetsToken, generaljets);
0168
0169   //get primary vertices
0170   Handle<reco::VertexCollection> primaryVertex;
0171   iEvent.getByToken(primaryVertexToken, primaryVertex);
0172
0173   //get Transient Track Builder
0174   auto const& builder = iSetup.getHandle(transientTrackRecordToken_);
0175
0176   //init JetTagCollection
0177   if (!generaljets.product()->empty()) {
0178     edm::RefToBase<reco::Jet> jj = edm::RefToBase<reco::Jet>(generaljets, 0);
0179     pOut_jetTagCollection = std::make_unique<reco::JetTagCollection>(edm::makeRefToBaseProdFrom(jj, iEvent));
0180   }
0181
0182   //loop on trackIPTagInfos
0183   if (!primaryVertex->empty()) {
0184     const reco::Vertex* pv = &*primaryVertex->begin();
0185     //loop on jets
0186     for (edm::View<reco::CaloJet>::const_iterator itJet = jets->begin(); itJet != jets->end(); itJet++) {
0187       math::XYZVector jetMomentum = itJet->momentum();
0188       GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());
0189
0190       math::XYZVector trMomentum;
0191
0192       if (fabs(itJet->eta()) > m_mineta_fwjets) {
0193         if ((m_fwjets) && (itJet->et() > m_minet_fwjets))
0194           pOut->push_back(*itJet);  // fill forward jet as signal jet
0195       } else {
0196         //loop on tracks
0197         auto itTrack = tracks->begin();
0198         for (unsigned int i = 0; i < tsize; ++i) {
0199           float deltaR2 = reco::deltaR2(itJet->eta(), itJet->phi(), teta[i], tphi[i]);
0200           if (deltaR2 < 0.25) {
0201             reco::TransientTrack transientTrack = builder->build(*itTrack);
0202             float jetTrackDistance = -((IPTools::jetTrackDistance(transientTrack, direction, *pv)).second).value();
0203
0204             //select the tracks compabible with the jet
0205             if ((itTrack->pt() > m_MinTrackPt) && (itTrack->normalizedChi2() < m_MaxTrackChi2) &&
0206                 (jetTrackDistance < m_MaxTrackDistanceToJet)) {
0207               trMomentum += itTrack->momentum();  //calculate the Sum(trackPt)
0208             }
0209           }
0210           itTrack++;
0211         }
0212         //if Sum(comp.trackPt)/CaloJetPt > minPtRatio or Sum(trackPt) > minPt  the jet is a signal jet
0213         if (trMomentum.rho() / jetMomentum.rho() > m_MinGoodJetTrackPtRatio || trMomentum.rho() > m_MinGoodJetTrackPt) {
0214           pOut->push_back(*itJet);  // fill it as signal jet
0215         } else                      //else it is a PUjet
0216         {
0217           pOut_PUjets->push_back(*itJet);  // fill it as PUjets
0218         }
0219       }
0220       RefToBase<reco::Jet> jRef(generaljets, itJet - jets->begin());
0221       (*pOut_jetTagCollection)[jRef] = trMomentum.rho();  // fill jetTagCollection
0222     }
0223   }
0224   iEvent.put(std::move(pOut));
0225   iEvent.put(std::move(pOut_PUjets), "PUjets");
0226   iEvent.put(std::move(pOut_jetTagCollection));
0227 }
0228
0229 // declare this class as a framework plugin
0230 DEFINE_FWK_MODULE(PixelJetPuId);