HLTProducers/src/L2TauPixelIsoTagProducer.cc

0001 //
0002 // Original Author:  Vadim Khotilovich
0003 //         Created:  2012-03-07
0004 //
0005
0006 #include <memory>
0007
0008 #include "RecoTauTag/HLTProducers/interface/L2TauPixelIsoTagProducer.h"
0009
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 //#include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 #include "DataFormats/Math/interface/deltaR.h"
0016 #include "DataFormats/BTauReco/interface/JetTag.h"
0017
0018 L2TauPixelIsoTagProducer::L2TauPixelIsoTagProducer(const edm::ParameterSet& conf)
0019     : m_jetSrc_token(consumes<edm::View<reco::Jet> >(conf.getParameter<edm::InputTag>("JetSrc"))),
0020       m_vertexSrc_token(consumes<reco::VertexCollection>(conf.getParameter<edm::InputTag>("VertexSrc"))),
0021       m_trackSrc_token(consumes<reco::TrackCollection>(
0022           conf.getParameter<edm::InputTag>("TrackSrc"))),  // for future use (now tracks are taken directly from PV)
0023       m_beamSpotSrc_token(consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("BeamSpotSrc"))),
0024       m_maxNumberPV(conf.getParameter<int>("MaxNumberPV")),  // for future use, now is assumed to be = 1
0025       m_trackMinPt(conf.getParameter<double>("TrackMinPt")),
0026       m_trackMaxDxy(conf.getParameter<double>("TrackMaxDxy")),
0027       m_trackMaxNChi2(conf.getParameter<double>("TrackMaxNChi2")),
0028       m_trackMinNHits(conf.getParameter<int>("TrackMinNHits")),
0029       m_trackPVMaxDZ(conf.getParameter<double>("TrackPVMaxDZ")),  // for future use with tracks not from PV
0030       m_isoCone2Min(std::pow(conf.getParameter<double>("IsoConeMin"), 2)),
0031       m_isoCone2Max(std::pow(conf.getParameter<double>("IsoConeMax"), 2)) {
0032   produces<reco::JetTagCollection>();
0033 }
0034
0035 void L2TauPixelIsoTagProducer::produce(edm::StreamID sid, edm::Event& ev, const edm::EventSetup& es) const {
0036   using namespace reco;
0037   using namespace std;
0038   using namespace edm;
0039   //m_trackSrc.encode();
0040
0041   edm::Handle<BeamSpot> bs;
0042   ev.getByToken(m_beamSpotSrc_token, bs);
0043
0044   // Get jets
0045   Handle<View<Jet> > jets_h;
0046   ev.getByToken(m_jetSrc_token, jets_h);
0047   vector<RefToBase<Jet> > jets;
0048   jets.reserve(jets_h->size());
0049   for (size_t i = 0; i < jets_h->size(); ++i)
0050     jets.push_back(jets_h->refAt(i));
0051
0052   // define the product to store
0053   unique_ptr<JetTagCollection> jetTagCollection;
0054   if (jets.empty()) {
0055     jetTagCollection = std::make_unique<JetTagCollection>();
0056   } else {
0057     jetTagCollection = std::make_unique<JetTagCollection>(RefToBaseProd<Jet>(jets_h));
0058   }
0059   // by default, initialize all the jets as isolated:
0060   for (const auto& jet : jets)
0061     (*jetTagCollection)[jet] = 0.f;
0062
0063   // Get pixel vertices (their x,y positions are already supposed to be defined from the BeamSpot)
0064   Handle<VertexCollection> vertices;
0065   ev.getByToken(m_vertexSrc_token, vertices);
0066
0067   // find the primary vertex (the 1st valid non-fake vertex in the collection)
0068   const Vertex* pv = nullptr;
0069   for (const auto& v : *(vertices.product())) {
0070     if (!v.isValid() || v.isFake())
0071       continue;
0072     pv = &v;
0073     break;
0074   }
0075
0076   // If primary vertex exists, calculate jets' isolation:
0077   if (pv && !jets.empty()) {
0078     for (const auto& jet : jets) {
0079       // re-calculate jet eta in PV:
0080       float jet_eta = Jet::physicsEta(pv->z(), jet->eta());
0081       float jet_phi = jet->phi();
0082
0083       // to calculate isolation, use only tracks that were assigned to the vertex
0084       float iso = 0.f;
0085       for (vector<TrackBaseRef>::const_iterator tr = pv->tracks_begin(); tr != pv->tracks_end(); ++tr) {
0086         if ((*tr)->pt() < m_trackMinPt)
0087           continue;
0088         if ((*tr)->numberOfValidHits() < m_trackMinNHits)
0089           continue;
0090         if ((*tr)->normalizedChi2() > m_trackMaxNChi2)
0091           continue;
0092         if (std::abs((*tr)->dxy(*bs)) > m_trackMaxDxy)
0093           continue;
0094
0095         float dr2 = deltaR2(jet_eta, jet_phi, (*tr)->eta(), (*tr)->phi());
0096
0097         // sum pT based isolation
0098         if (dr2 >= m_isoCone2Min && dr2 <= m_isoCone2Max)
0099           iso += (*tr)->pt();
0100       }
0101
0102       (*jetTagCollection)[jet] = iso;
0103     }
0104   }
0105
0106   ev.put(std::move(jetTagCollection));
0107 }
0108
0109 void L2TauPixelIsoTagProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0110   edm::ParameterSetDescription desc;
0111   desc.add<edm::InputTag>("JetSrc", edm::InputTag("hltL2DiTauCaloJets"))->setComment("Jet source collection");
0112   desc.add<edm::InputTag>("BeamSpotSrc", edm::InputTag("hltOnlineBeamSpot"));
0113   desc.add<edm::InputTag>("VertexSrc", edm::InputTag("hltPixelVertices"))
0114       ->setComment("Collection of vertices where isolation tracks come from");
0115   desc.add<int>("MaxNumberPV", 1)->setComment("No. of considered vertices (not used yet)");
0116   desc.add<double>("IsoConeMax", 0.4)->setComment("Outer radius of isolation annulus");
0117   desc.add<double>("IsoConeMin", 0.2)->setComment("Inner radius of isolation annulus");
0118   desc.add<double>("TrackMinPt", 1.6)->setComment("Isolation track quality: min. pT");
0119   desc.add<int>("TrackMinNHits", 3)->setComment("Isolation track quality: min. no. of hits");
0120   desc.add<double>("TrackMaxNChi2", 100.0)->setComment("Isolation track quality: max. chi2/ndof");
0121   desc.add<double>("TrackPVMaxDZ", 0.1)->setComment("Isolation track quality: max. dz");
0122   ;
0123   desc.add<double>("TrackMaxDxy", 0.2)->setComment("Isolation track quality: max. dxy");
0124   ;
0125   desc.add<edm::InputTag>("TrackSrc", edm::InputTag(""))->setComment("Not used yet");
0126   descriptions.setComment("Produces isolation tag for caloJets/L2Taus");
0127   descriptions.add("L2TauPixelIsoTagProducer", desc);
0128 }
0129
0130 #include "FWCore/Framework/interface/MakerMacros.h"
0131 DEFINE_FWK_MODULE(L2TauPixelIsoTagProducer);