File indexing completed on 2021-06-24 02:11:12
0001 #ifndef RecoBTag_IPProducer
0002 #define RecoBTag_IPProducer
0003
0004
0005 #include <algorithm>
0006 #include <cmath>
0007 #include <functional>
0008 #include <iostream>
0009 #include <memory>
0010
0011
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "DataFormats/Common/interface/View.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0021 #include "FWCore/Utilities/interface/ESGetToken.h"
0022
0023 #include "DataFormats/TrackReco/interface/Track.h"
0024 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0025 #include "DataFormats/BTauReco/interface/JetTag.h"
0026 #include "DataFormats/BTauReco/interface/TrackIPTagInfo.h"
0027 #include "DataFormats/Candidate/interface/Candidate.h"
0028
0029 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0030 #include "TrackingTools/IPTools/interface/IPTools.h"
0031 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0032
0033 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
0034 #include "RecoVertex/VertexPrimitives/interface/ConvertError.h"
0035 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0036 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0037 #include "RecoVertex/GhostTrackFitter/interface/GhostTrack.h"
0038 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackState.h"
0039 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
0040 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackFitter.h"
0041
0042 #include "RecoBTag/TrackProbability/interface/HistogramProbabilityEstimator.h"
0043 #include "DataFormats/BTauReco/interface/JTATagInfo.h"
0044 #include "DataFormats/BTauReco/interface/JetTagInfo.h"
0045 #include "CondFormats/DataRecord/interface/BTagTrackProbability2DRcd.h"
0046 #include "CondFormats/DataRecord/interface/BTagTrackProbability3DRcd.h"
0047
0048
0049 #include <memory>
0050
0051
0052 #include "FWCore/Framework/interface/Frameworkfwd.h"
0053 #include "FWCore/Framework/interface/stream/EDProducer.h"
0054 #include "FWCore/Utilities/interface/InputTag.h"
0055 #include "FWCore/Framework/interface/ConsumesCollector.h"
0056
0057 class HistogramProbabilityEstimator;
0058
0059 namespace IPProducerHelpers {
0060 class FromJTA {
0061 public:
0062 FromJTA(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
0063 : token_associator(
0064 iC.consumes<reco::JetTracksAssociationCollection>(iConfig.getParameter<edm::InputTag>("jetTracks"))) {}
0065 reco::TrackRefVector tracks(const reco::JTATagInfo& it) { return it.tracks(); }
0066 std::vector<reco::JTATagInfo> makeBaseVector(const edm::Event& iEvent) {
0067 edm::Handle<reco::JetTracksAssociationCollection> jetTracksAssociation;
0068 iEvent.getByToken(token_associator, jetTracksAssociation);
0069 std::vector<reco::JTATagInfo> bases;
0070 size_t i = 0;
0071 for (reco::JetTracksAssociationCollection::const_iterator it = jetTracksAssociation->begin();
0072 it != jetTracksAssociation->end();
0073 it++, i++) {
0074 edm::Ref<reco::JetTracksAssociationCollection> jtaRef(jetTracksAssociation, i);
0075 bases.push_back(reco::JTATagInfo(jtaRef));
0076 }
0077 return bases;
0078 }
0079
0080 edm::EDGetTokenT<reco::JetTracksAssociationCollection> token_associator;
0081 };
0082 class FromJetAndCands {
0083 public:
0084 FromJetAndCands(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC, const std::string& jets = "jets")
0085 : token_jets(iC.consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>(jets))),
0086 token_cands(iC.consumes<edm::View<reco::Candidate> >(iConfig.getParameter<edm::InputTag>("candidates"))),
0087 maxDeltaR(iConfig.getParameter<double>("maxDeltaR")),
0088 explicitJTA(iConfig.existsAs<bool>("explicitJTA") ? iConfig.getParameter<bool>("explicitJTA") : false) {}
0089
0090 const std::vector<reco::CandidatePtr>& tracks(const reco::JetTagInfo& it) { return m_map[it.jet().key()]; }
0091 std::vector<reco::JetTagInfo> makeBaseVector(const edm::Event& iEvent) {
0092 edm::Handle<edm::View<reco::Jet> > jets;
0093 iEvent.getByToken(token_jets, jets);
0094 std::vector<reco::JetTagInfo> bases;
0095
0096 edm::Handle<edm::View<reco::Candidate> > cands;
0097 iEvent.getByToken(token_cands, cands);
0098 m_map.clear();
0099 m_map.resize(jets->size());
0100 double maxDeltaR2 = maxDeltaR * maxDeltaR;
0101 size_t i = 0;
0102 for (edm::View<reco::Jet>::const_iterator it = jets->begin(); it != jets->end(); it++, i++) {
0103 edm::RefToBase<reco::Jet> jRef(jets, i);
0104 bases.push_back(reco::JetTagInfo(jRef));
0105 if (explicitJTA) {
0106 for (size_t j = 0; j < it->numberOfDaughters(); ++j) {
0107 if (it->daughterPtr(j)->bestTrack() != nullptr && it->daughterPtr(j)->charge() != 0) {
0108 m_map[i].push_back(it->daughterPtr(j));
0109 }
0110 }
0111 } else {
0112 for (size_t j = 0; j < cands->size(); ++j) {
0113 if ((*cands)[j].bestTrack() != nullptr && (*cands)[j].charge() != 0 && (*cands)[j].pt() > 0 &&
0114 Geom::deltaR2((*cands)[j], (*jets)[i]) < maxDeltaR2) {
0115 m_map[i].push_back(cands->ptrAt(j));
0116 }
0117 }
0118 }
0119 }
0120 return bases;
0121 }
0122 std::vector<std::vector<reco::CandidatePtr> > m_map;
0123 edm::EDGetTokenT<edm::View<reco::Jet> > token_jets;
0124 edm::EDGetTokenT<edm::View<reco::Candidate> > token_cands;
0125 double maxDeltaR;
0126 bool explicitJTA;
0127 };
0128 }
0129 template <class Container, class Base, class Helper>
0130 class IPProducer : public edm::stream::EDProducer<> {
0131 public:
0132 typedef std::vector<reco::IPTagInfo<Container, Base> > Product;
0133
0134 explicit IPProducer(const edm::ParameterSet&);
0135 ~IPProducer() override;
0136 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0137
0138 void produce(edm::Event&, const edm::EventSetup&) override;
0139
0140 private:
0141 void checkEventSetup(const edm::EventSetup& iSetup);
0142
0143 const edm::ParameterSet& m_config;
0144 edm::EDGetTokenT<reco::VertexCollection> token_primaryVertex;
0145 edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> token_trackBuilder;
0146 edm::ESGetToken<TrackProbabilityCalibration, BTagTrackProbability2DRcd> token_calib2D;
0147 edm::ESGetToken<TrackProbabilityCalibration, BTagTrackProbability3DRcd> token_calib3D;
0148
0149 bool m_computeProbabilities;
0150 bool m_computeGhostTrack;
0151 double m_ghostTrackPriorDeltaR;
0152 std::unique_ptr<HistogramProbabilityEstimator> m_probabilityEstimator;
0153 unsigned long long m_calibrationCacheId2D;
0154 unsigned long long m_calibrationCacheId3D;
0155 bool m_useDB;
0156
0157 int m_cutPixelHits;
0158 int m_cutTotalHits;
0159 double m_cutMaxTIP;
0160 double m_cutMinPt;
0161 double m_cutMaxChiSquared;
0162 double m_cutMaxLIP;
0163 bool m_directionWithTracks;
0164 bool m_directionWithGhostTrack;
0165 bool m_useTrackQuality;
0166 Helper m_helper;
0167 };
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190 template <class Container, class Base, class Helper>
0191 IPProducer<Container, Base, Helper>::IPProducer(const edm::ParameterSet& iConfig)
0192 : m_config(iConfig), m_helper(iConfig, consumesCollector()) {
0193 m_calibrationCacheId3D = 0;
0194 m_calibrationCacheId2D = 0;
0195
0196 token_primaryVertex = consumes<reco::VertexCollection>(m_config.getParameter<edm::InputTag>("primaryVertex"));
0197 token_trackBuilder =
0198 esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"));
0199 token_calib2D = esConsumes<TrackProbabilityCalibration, BTagTrackProbability2DRcd>();
0200 token_calib3D = esConsumes<TrackProbabilityCalibration, BTagTrackProbability3DRcd>();
0201
0202 m_computeProbabilities = m_config.getParameter<bool>("computeProbabilities");
0203 m_computeGhostTrack = m_config.getParameter<bool>("computeGhostTrack");
0204 m_ghostTrackPriorDeltaR = m_config.getParameter<double>("ghostTrackPriorDeltaR");
0205 m_cutPixelHits = m_config.getParameter<int>("minimumNumberOfPixelHits");
0206 m_cutTotalHits = m_config.getParameter<int>("minimumNumberOfHits");
0207 m_cutMaxTIP = m_config.getParameter<double>("maximumTransverseImpactParameter");
0208 m_cutMinPt = m_config.getParameter<double>("minimumTransverseMomentum");
0209 m_cutMaxChiSquared = m_config.getParameter<double>("maximumChiSquared");
0210 m_cutMaxLIP = m_config.getParameter<double>("maximumLongitudinalImpactParameter");
0211 m_directionWithTracks = m_config.getParameter<bool>("jetDirectionUsingTracks");
0212 m_directionWithGhostTrack = m_config.getParameter<bool>("jetDirectionUsingGhostTrack");
0213 m_useTrackQuality = m_config.getParameter<bool>("useTrackQuality");
0214
0215 if (m_computeGhostTrack)
0216 produces<reco::TrackCollection>("ghostTracks");
0217 produces<Product>();
0218 }
0219
0220 template <class Container, class Base, class Helper>
0221 IPProducer<Container, Base, Helper>::~IPProducer() {}
0222
0223
0224
0225
0226
0227 template <class Container, class Base, class Helper>
0228 void IPProducer<Container, Base, Helper>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0229
0230 if (m_computeProbabilities)
0231 checkEventSetup(iSetup);
0232
0233 edm::Handle<reco::VertexCollection> primaryVertex;
0234 iEvent.getByToken(token_primaryVertex, primaryVertex);
0235
0236 edm::ESHandle<TransientTrackBuilder> builder = iSetup.getHandle(token_trackBuilder);
0237
0238
0239
0240 auto result = std::make_unique<Product>();
0241
0242 std::unique_ptr<reco::TrackCollection> ghostTracks;
0243 reco::TrackRefProd ghostTrackRefProd;
0244 if (m_computeGhostTrack) {
0245 ghostTracks = std::make_unique<reco::TrackCollection>();
0246 ghostTrackRefProd = iEvent.getRefBeforePut<reco::TrackCollection>("ghostTracks");
0247 }
0248
0249
0250 reco::Vertex dummy;
0251 const reco::Vertex* pv = &dummy;
0252 edm::Ref<reco::VertexCollection> pvRef;
0253 if (!primaryVertex->empty()) {
0254 pv = &*primaryVertex->begin();
0255
0256 pvRef = edm::Ref<reco::VertexCollection>(primaryVertex, 0);
0257 } else {
0258 reco::Vertex::Error e;
0259 e(0, 0) = 0.0015 * 0.0015;
0260 e(1, 1) = 0.0015 * 0.0015;
0261 e(2, 2) = 15. * 15.;
0262 reco::Vertex::Point p(0, 0, 0);
0263 dummy = reco::Vertex(p, e, 0, 0, 0);
0264 }
0265
0266 std::vector<Base> baseTagInfos = m_helper.makeBaseVector(iEvent);
0267 for (typename std::vector<Base>::const_iterator it = baseTagInfos.begin(); it != baseTagInfos.end(); it++) {
0268 Container tracks = m_helper.tracks(*it);
0269 math::XYZVector jetMomentum = it->jet()->momentum();
0270
0271 if (m_directionWithTracks) {
0272 jetMomentum *= 0.5;
0273 for (typename Container::const_iterator itTrack = tracks.begin(); itTrack != tracks.end(); ++itTrack)
0274 if (reco::btag::toTrack(*itTrack)->numberOfValidHits() >= m_cutTotalHits)
0275 jetMomentum += (*itTrack)->momentum();
0276 }
0277
0278 Container selectedTracks;
0279 std::vector<reco::TransientTrack> transientTracks;
0280
0281 for (typename Container::const_iterator itTrack = tracks.begin(); itTrack != tracks.end(); ++itTrack) {
0282 reco::TransientTrack transientTrack = builder->build(*itTrack);
0283 const reco::Track& track = transientTrack.track();
0284
0285
0286
0287
0288
0289
0290
0291 if (track.pt() > m_cutMinPt &&
0292 track.hitPattern().numberOfValidHits() >= m_cutTotalHits &&
0293 track.hitPattern().numberOfValidPixelHits() >= m_cutPixelHits &&
0294 track.normalizedChi2() < m_cutMaxChiSquared && std::abs(track.dxy(pv->position())) < m_cutMaxTIP &&
0295 std::abs(track.dz(pv->position())) < m_cutMaxLIP) {
0296
0297 selectedTracks.push_back(*itTrack);
0298 transientTracks.push_back(transientTrack);
0299 }
0300 }
0301
0302 GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());
0303
0304 std::unique_ptr<reco::GhostTrack> ghostTrack;
0305 reco::TrackRef ghostTrackRef;
0306 if (m_computeGhostTrack) {
0307 reco::GhostTrackFitter fitter;
0308 GlobalPoint origin = RecoVertex::convertPos(pv->position());
0309 GlobalError error = RecoVertex::convertError(pv->error());
0310 ghostTrack = std::make_unique<reco::GhostTrack>(
0311 fitter.fit(origin, error, direction, m_ghostTrackPriorDeltaR, transientTracks));
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329 ghostTrackRef = reco::TrackRef(ghostTrackRefProd, ghostTracks->size());
0330 ghostTracks->push_back(*ghostTrack);
0331
0332 if (m_directionWithGhostTrack) {
0333 const reco::GhostTrackPrediction& pred = ghostTrack->prediction();
0334 double lambda = pred.lambda(origin);
0335 dummy = reco::Vertex(RecoVertex::convertPos(pred.position(lambda)),
0336 RecoVertex::convertError(pred.positionError(lambda)),
0337 0,
0338 0,
0339 0);
0340 pv = &dummy;
0341 direction = pred.direction();
0342 }
0343 }
0344
0345 std::vector<float> prob2D, prob3D;
0346 std::vector<reco::btag::TrackIPData> ipData;
0347
0348 for (unsigned int ind = 0; ind < transientTracks.size(); ind++) {
0349 const reco::TransientTrack& transientTrack = transientTracks[ind];
0350 const reco::Track& track = transientTrack.track();
0351
0352 reco::btag::TrackIPData trackIP;
0353 trackIP.ip3d = IPTools::signedImpactParameter3D(transientTrack, direction, *pv).second;
0354 trackIP.ip2d = IPTools::signedTransverseImpactParameter(transientTrack, direction, *pv).second;
0355
0356 TrajectoryStateOnSurface closest =
0357 IPTools::closestApproachToJet(transientTrack.impactPointState(), *pv, direction, transientTrack.field());
0358 if (closest.isValid())
0359 trackIP.closestToJetAxis = closest.globalPosition();
0360
0361
0362 trackIP.distanceToJetAxis = IPTools::jetTrackDistance(transientTrack, direction, *pv).second;
0363
0364 if (ghostTrack.get()) {
0365 const std::vector<reco::GhostTrackState>& states = ghostTrack->states();
0366 std::vector<reco::GhostTrackState>::const_iterator pos =
0367 std::find_if(states.begin(), states.end(), [&](auto& arg) { return arg.track() == transientTrack; });
0368
0369 if (pos != states.end() && pos->isValid()) {
0370 VertexDistance3D dist;
0371 const reco::GhostTrackPrediction& pred = ghostTrack->prediction();
0372 GlobalPoint p1 = pos->tsos().globalPosition();
0373 GlobalError e1 = pos->tsos().cartesianError().position();
0374 GlobalPoint p2 = pred.position(pos->lambda());
0375 GlobalError e2 = pred.positionError(pos->lambda());
0376 trackIP.closestToGhostTrack = p1;
0377 trackIP.distanceToGhostTrack = dist.distance(VertexState(p1, e1), VertexState(p2, e2));
0378 trackIP.ghostTrackWeight = pos->weight();
0379 } else {
0380 trackIP.distanceToGhostTrack = Measurement1D(0., -1.);
0381 trackIP.ghostTrackWeight = 0.;
0382 }
0383 } else {
0384 trackIP.distanceToGhostTrack = Measurement1D(0., -1.);
0385 trackIP.ghostTrackWeight = 1.;
0386 }
0387
0388 ipData.push_back(trackIP);
0389
0390 if (m_computeProbabilities) {
0391
0392 std::pair<bool, double> probability = m_probabilityEstimator->probability(
0393 m_useTrackQuality, 0, ipData.back().ip3d.significance(), track, *(it->jet()), *pv);
0394 prob3D.push_back(probability.first ? probability.second : -1.);
0395
0396
0397 probability = m_probabilityEstimator->probability(
0398 m_useTrackQuality, 1, ipData.back().ip2d.significance(), track, *(it->jet()), *pv);
0399 prob2D.push_back(probability.first ? probability.second : -1.);
0400 }
0401 }
0402
0403 result->push_back(
0404 typename Product::value_type(ipData, prob2D, prob3D, selectedTracks, *it, pvRef, direction, ghostTrackRef));
0405 }
0406
0407 if (m_computeGhostTrack)
0408 iEvent.put(std::move(ghostTracks), "ghostTracks");
0409 iEvent.put(std::move(result));
0410 }
0411
0412 #include "CondFormats/BTauObjects/interface/TrackProbabilityCalibration.h"
0413 #include "CondFormats/DataRecord/interface/BTagTrackProbability2DRcd.h"
0414 #include "CondFormats/DataRecord/interface/BTagTrackProbability3DRcd.h"
0415 #include "FWCore/Framework/interface/EventSetupRecord.h"
0416 #include "FWCore/Framework/interface/EventSetupRecordImplementation.h"
0417 #include "FWCore/Framework/interface/EventSetupRecordKey.h"
0418
0419 template <class Container, class Base, class Helper>
0420 void IPProducer<Container, Base, Helper>::checkEventSetup(const edm::EventSetup& iSetup) {
0421 const edm::eventsetup::EventSetupRecord& re2D = iSetup.get<BTagTrackProbability2DRcd>();
0422 const edm::eventsetup::EventSetupRecord& re3D = iSetup.get<BTagTrackProbability3DRcd>();
0423 unsigned long long cacheId2D = re2D.cacheIdentifier();
0424 unsigned long long cacheId3D = re3D.cacheIdentifier();
0425
0426 if (cacheId2D != m_calibrationCacheId2D || cacheId3D != m_calibrationCacheId3D)
0427 {
0428
0429 edm::ESHandle<TrackProbabilityCalibration> calib2DHandle = iSetup.getHandle(token_calib2D);
0430 edm::ESHandle<TrackProbabilityCalibration> calib3DHandle = iSetup.getHandle(token_calib3D);
0431
0432 const TrackProbabilityCalibration* ca2D = calib2DHandle.product();
0433 const TrackProbabilityCalibration* ca3D = calib3DHandle.product();
0434
0435 m_probabilityEstimator = std::make_unique<HistogramProbabilityEstimator>(ca3D, ca2D);
0436 }
0437 m_calibrationCacheId3D = cacheId3D;
0438 m_calibrationCacheId2D = cacheId2D;
0439 }
0440
0441
0442
0443 template <>
0444 inline void IPProducer<reco::TrackRefVector, reco::JTATagInfo, IPProducerHelpers::FromJTA>::fillDescriptions(
0445 edm::ConfigurationDescriptions& descriptions) {
0446 edm::ParameterSetDescription desc;
0447 desc.add<double>("maximumTransverseImpactParameter", 0.2);
0448 desc.add<int>("minimumNumberOfHits", 8);
0449 desc.add<double>("minimumTransverseMomentum", 1.0);
0450 desc.add<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
0451 desc.add<double>("maximumLongitudinalImpactParameter", 17.0);
0452 desc.add<bool>("computeGhostTrack", true);
0453 desc.add<double>("ghostTrackPriorDeltaR", 0.03);
0454 desc.add<edm::InputTag>("jetTracks", edm::InputTag("ak4JetTracksAssociatorAtVertexPF"));
0455 desc.add<bool>("jetDirectionUsingGhostTrack", false);
0456 desc.add<int>("minimumNumberOfPixelHits", 2);
0457 desc.add<bool>("jetDirectionUsingTracks", false);
0458 desc.add<bool>("computeProbabilities", true);
0459 desc.add<bool>("useTrackQuality", false);
0460 desc.add<double>("maximumChiSquared", 5.0);
0461 descriptions.addDefault(desc);
0462 }
0463
0464 template <>
0465 inline void
0466 IPProducer<std::vector<reco::CandidatePtr>, reco::JetTagInfo, IPProducerHelpers::FromJetAndCands>::fillDescriptions(
0467 edm::ConfigurationDescriptions& descriptions) {
0468 edm::ParameterSetDescription desc;
0469 desc.add<double>("maximumTransverseImpactParameter", 0.2);
0470 desc.add<int>("minimumNumberOfHits", 8);
0471 desc.add<double>("minimumTransverseMomentum", 1.0);
0472 desc.add<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
0473 desc.add<double>("maximumLongitudinalImpactParameter", 17.0);
0474 desc.add<bool>("computeGhostTrack", true);
0475 desc.add<double>("maxDeltaR", 0.4);
0476 desc.add<edm::InputTag>("candidates", edm::InputTag("particleFlow"));
0477 desc.add<bool>("jetDirectionUsingGhostTrack", false);
0478 desc.add<int>("minimumNumberOfPixelHits", 2);
0479 desc.add<bool>("jetDirectionUsingTracks", false);
0480 desc.add<bool>("computeProbabilities", true);
0481 desc.add<bool>("useTrackQuality", false);
0482 desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
0483 desc.add<double>("ghostTrackPriorDeltaR", 0.03);
0484 desc.add<double>("maximumChiSquared", 5.0);
0485 desc.addOptional<bool>("explicitJTA", false);
0486 descriptions.addDefault(desc);
0487 }
0488
0489 #endif