File indexing completed on 2024-04-06 12:24:27
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 edm::EDGetTokenT<reco::VertexCollection> token_primaryVertex;
0144 edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> token_trackBuilder;
0145 edm::ESGetToken<TrackProbabilityCalibration, BTagTrackProbability2DRcd> token_calib2D;
0146 edm::ESGetToken<TrackProbabilityCalibration, BTagTrackProbability3DRcd> token_calib3D;
0147
0148 bool m_computeProbabilities;
0149 bool m_computeGhostTrack;
0150 double m_ghostTrackPriorDeltaR;
0151 std::unique_ptr<HistogramProbabilityEstimator> m_probabilityEstimator;
0152 unsigned long long m_calibrationCacheId2D;
0153 unsigned long long m_calibrationCacheId3D;
0154 bool m_useDB;
0155
0156 int m_cutPixelHits;
0157 int m_cutTotalHits;
0158 double m_cutMaxTIP;
0159 double m_cutMinPt;
0160 double m_cutMaxChiSquared;
0161 double m_cutMaxLIP;
0162 bool m_directionWithTracks;
0163 bool m_directionWithGhostTrack;
0164 bool m_useTrackQuality;
0165 Helper m_helper;
0166 };
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189 template <class Container, class Base, class Helper>
0190 IPProducer<Container, Base, Helper>::IPProducer(const edm::ParameterSet& iConfig)
0191 : m_helper(iConfig, consumesCollector()) {
0192 m_calibrationCacheId3D = 0;
0193 m_calibrationCacheId2D = 0;
0194
0195 token_primaryVertex = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryVertex"));
0196 token_trackBuilder =
0197 esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"));
0198 token_calib2D = esConsumes<TrackProbabilityCalibration, BTagTrackProbability2DRcd>();
0199 token_calib3D = esConsumes<TrackProbabilityCalibration, BTagTrackProbability3DRcd>();
0200
0201 m_computeProbabilities = iConfig.getParameter<bool>("computeProbabilities");
0202 m_computeGhostTrack = iConfig.getParameter<bool>("computeGhostTrack");
0203 m_ghostTrackPriorDeltaR = iConfig.getParameter<double>("ghostTrackPriorDeltaR");
0204 m_cutPixelHits = iConfig.getParameter<int>("minimumNumberOfPixelHits");
0205 m_cutTotalHits = iConfig.getParameter<int>("minimumNumberOfHits");
0206 m_cutMaxTIP = iConfig.getParameter<double>("maximumTransverseImpactParameter");
0207 m_cutMinPt = iConfig.getParameter<double>("minimumTransverseMomentum");
0208 m_cutMaxChiSquared = iConfig.getParameter<double>("maximumChiSquared");
0209 m_cutMaxLIP = iConfig.getParameter<double>("maximumLongitudinalImpactParameter");
0210 m_directionWithTracks = iConfig.getParameter<bool>("jetDirectionUsingTracks");
0211 m_directionWithGhostTrack = iConfig.getParameter<bool>("jetDirectionUsingGhostTrack");
0212 m_useTrackQuality = iConfig.getParameter<bool>("useTrackQuality");
0213
0214 if (m_computeGhostTrack)
0215 produces<reco::TrackCollection>("ghostTracks");
0216 produces<Product>();
0217 }
0218
0219 template <class Container, class Base, class Helper>
0220 IPProducer<Container, Base, Helper>::~IPProducer() {}
0221
0222
0223
0224
0225
0226 template <class Container, class Base, class Helper>
0227 void IPProducer<Container, Base, Helper>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0228
0229 if (m_computeProbabilities)
0230 checkEventSetup(iSetup);
0231
0232 edm::Handle<reco::VertexCollection> primaryVertex;
0233 iEvent.getByToken(token_primaryVertex, primaryVertex);
0234
0235 edm::ESHandle<TransientTrackBuilder> builder = iSetup.getHandle(token_trackBuilder);
0236
0237
0238
0239 auto result = std::make_unique<Product>();
0240
0241 std::unique_ptr<reco::TrackCollection> ghostTracks;
0242 reco::TrackRefProd ghostTrackRefProd;
0243 if (m_computeGhostTrack) {
0244 ghostTracks = std::make_unique<reco::TrackCollection>();
0245 ghostTrackRefProd = iEvent.getRefBeforePut<reco::TrackCollection>("ghostTracks");
0246 }
0247
0248
0249 reco::Vertex dummy;
0250 const reco::Vertex* pv = &dummy;
0251 edm::Ref<reco::VertexCollection> pvRef;
0252 if (!primaryVertex->empty()) {
0253 pv = &*primaryVertex->begin();
0254
0255 pvRef = edm::Ref<reco::VertexCollection>(primaryVertex, 0);
0256 } else {
0257 reco::Vertex::Error e;
0258 e(0, 0) = 0.0015 * 0.0015;
0259 e(1, 1) = 0.0015 * 0.0015;
0260 e(2, 2) = 15. * 15.;
0261 reco::Vertex::Point p(0, 0, 0);
0262 dummy = reco::Vertex(p, e, 0, 0, 0);
0263 }
0264
0265 std::vector<Base> baseTagInfos = m_helper.makeBaseVector(iEvent);
0266 for (typename std::vector<Base>::const_iterator it = baseTagInfos.begin(); it != baseTagInfos.end(); it++) {
0267 Container tracks = m_helper.tracks(*it);
0268 math::XYZVector jetMomentum = it->jet()->momentum();
0269
0270 if (m_directionWithTracks) {
0271 jetMomentum *= 0.5;
0272 for (typename Container::const_iterator itTrack = tracks.begin(); itTrack != tracks.end(); ++itTrack)
0273 if (reco::btag::toTrack(*itTrack)->numberOfValidHits() >= m_cutTotalHits)
0274 jetMomentum += (*itTrack)->momentum();
0275 }
0276
0277 Container selectedTracks;
0278 std::vector<reco::TransientTrack> transientTracks;
0279
0280 for (typename Container::const_iterator itTrack = tracks.begin(); itTrack != tracks.end(); ++itTrack) {
0281 reco::TransientTrack transientTrack = builder->build(*itTrack);
0282 const reco::Track& track = transientTrack.track();
0283
0284
0285
0286
0287
0288
0289
0290 if (track.pt() > m_cutMinPt &&
0291 track.hitPattern().numberOfValidHits() >= m_cutTotalHits &&
0292 track.hitPattern().numberOfValidPixelHits() >= m_cutPixelHits &&
0293 track.normalizedChi2() < m_cutMaxChiSquared && std::abs(track.dxy(pv->position())) < m_cutMaxTIP &&
0294 std::abs(track.dz(pv->position())) < m_cutMaxLIP) {
0295
0296 selectedTracks.push_back(*itTrack);
0297 transientTracks.push_back(transientTrack);
0298 }
0299 }
0300
0301 GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());
0302
0303 std::unique_ptr<reco::GhostTrack> ghostTrack;
0304 reco::TrackRef ghostTrackRef;
0305 if (m_computeGhostTrack) {
0306 reco::GhostTrackFitter fitter;
0307 GlobalPoint origin = RecoVertex::convertPos(pv->position());
0308 GlobalError error = RecoVertex::convertError(pv->error());
0309 ghostTrack = std::make_unique<reco::GhostTrack>(
0310 fitter.fit(origin, error, direction, m_ghostTrackPriorDeltaR, transientTracks));
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328 ghostTrackRef = reco::TrackRef(ghostTrackRefProd, ghostTracks->size());
0329 ghostTracks->push_back(*ghostTrack);
0330
0331 if (m_directionWithGhostTrack) {
0332 const reco::GhostTrackPrediction& pred = ghostTrack->prediction();
0333 double lambda = pred.lambda(origin);
0334 dummy = reco::Vertex(RecoVertex::convertPos(pred.position(lambda)),
0335 RecoVertex::convertError(pred.positionError(lambda)),
0336 0,
0337 0,
0338 0);
0339 pv = &dummy;
0340 direction = pred.direction();
0341 }
0342 }
0343
0344 std::vector<float> prob2D, prob3D;
0345 std::vector<reco::btag::TrackIPData> ipData;
0346
0347 for (unsigned int ind = 0; ind < transientTracks.size(); ind++) {
0348 const reco::TransientTrack& transientTrack = transientTracks[ind];
0349 const reco::Track& track = transientTrack.track();
0350
0351 reco::btag::TrackIPData trackIP;
0352 trackIP.ip3d = IPTools::signedImpactParameter3D(transientTrack, direction, *pv).second;
0353 trackIP.ip2d = IPTools::signedTransverseImpactParameter(transientTrack, direction, *pv).second;
0354
0355 TrajectoryStateOnSurface closest =
0356 IPTools::closestApproachToJet(transientTrack.impactPointState(), *pv, direction, transientTrack.field());
0357 if (closest.isValid())
0358 trackIP.closestToJetAxis = closest.globalPosition();
0359
0360
0361 trackIP.distanceToJetAxis = IPTools::jetTrackDistance(transientTrack, direction, *pv).second;
0362
0363 if (ghostTrack.get()) {
0364 const std::vector<reco::GhostTrackState>& states = ghostTrack->states();
0365 std::vector<reco::GhostTrackState>::const_iterator pos =
0366 std::find_if(states.begin(), states.end(), [&](auto& arg) { return arg.track() == transientTrack; });
0367
0368 if (pos != states.end() && pos->isValid()) {
0369 VertexDistance3D dist;
0370 const reco::GhostTrackPrediction& pred = ghostTrack->prediction();
0371 GlobalPoint p1 = pos->tsos().globalPosition();
0372 GlobalError e1 = pos->tsos().cartesianError().position();
0373 GlobalPoint p2 = pred.position(pos->lambda());
0374 GlobalError e2 = pred.positionError(pos->lambda());
0375 trackIP.closestToGhostTrack = p1;
0376 trackIP.distanceToGhostTrack = dist.distance(VertexState(p1, e1), VertexState(p2, e2));
0377 trackIP.ghostTrackWeight = pos->weight();
0378 } else {
0379 trackIP.distanceToGhostTrack = Measurement1D(0., -1.);
0380 trackIP.ghostTrackWeight = 0.;
0381 }
0382 } else {
0383 trackIP.distanceToGhostTrack = Measurement1D(0., -1.);
0384 trackIP.ghostTrackWeight = 1.;
0385 }
0386
0387 ipData.push_back(trackIP);
0388
0389 if (m_computeProbabilities) {
0390
0391 std::pair<bool, double> probability = m_probabilityEstimator->probability(
0392 m_useTrackQuality, 0, ipData.back().ip3d.significance(), track, *(it->jet()), *pv);
0393 prob3D.push_back(probability.first ? probability.second : -1.);
0394
0395
0396 probability = m_probabilityEstimator->probability(
0397 m_useTrackQuality, 1, ipData.back().ip2d.significance(), track, *(it->jet()), *pv);
0398 prob2D.push_back(probability.first ? probability.second : -1.);
0399 }
0400 }
0401
0402 result->push_back(
0403 typename Product::value_type(ipData, prob2D, prob3D, selectedTracks, *it, pvRef, direction, ghostTrackRef));
0404 }
0405
0406 if (m_computeGhostTrack)
0407 iEvent.put(std::move(ghostTracks), "ghostTracks");
0408 iEvent.put(std::move(result));
0409 }
0410
0411 #include "CondFormats/BTauObjects/interface/TrackProbabilityCalibration.h"
0412 #include "CondFormats/DataRecord/interface/BTagTrackProbability2DRcd.h"
0413 #include "CondFormats/DataRecord/interface/BTagTrackProbability3DRcd.h"
0414 #include "FWCore/Framework/interface/EventSetupRecord.h"
0415 #include "FWCore/Framework/interface/EventSetupRecordImplementation.h"
0416 #include "FWCore/Framework/interface/EventSetupRecordKey.h"
0417
0418 template <class Container, class Base, class Helper>
0419 void IPProducer<Container, Base, Helper>::checkEventSetup(const edm::EventSetup& iSetup) {
0420 const edm::eventsetup::EventSetupRecord& re2D = iSetup.get<BTagTrackProbability2DRcd>();
0421 const edm::eventsetup::EventSetupRecord& re3D = iSetup.get<BTagTrackProbability3DRcd>();
0422 unsigned long long cacheId2D = re2D.cacheIdentifier();
0423 unsigned long long cacheId3D = re3D.cacheIdentifier();
0424
0425 if (cacheId2D != m_calibrationCacheId2D || cacheId3D != m_calibrationCacheId3D)
0426 {
0427
0428 edm::ESHandle<TrackProbabilityCalibration> calib2DHandle = iSetup.getHandle(token_calib2D);
0429 edm::ESHandle<TrackProbabilityCalibration> calib3DHandle = iSetup.getHandle(token_calib3D);
0430
0431 const TrackProbabilityCalibration* ca2D = calib2DHandle.product();
0432 const TrackProbabilityCalibration* ca3D = calib3DHandle.product();
0433
0434 m_probabilityEstimator = std::make_unique<HistogramProbabilityEstimator>(ca3D, ca2D);
0435 }
0436 m_calibrationCacheId3D = cacheId3D;
0437 m_calibrationCacheId2D = cacheId2D;
0438 }
0439
0440
0441
0442 template <>
0443 inline void IPProducer<reco::TrackRefVector, reco::JTATagInfo, IPProducerHelpers::FromJTA>::fillDescriptions(
0444 edm::ConfigurationDescriptions& descriptions) {
0445 edm::ParameterSetDescription desc;
0446 desc.add<double>("maximumTransverseImpactParameter", 0.2);
0447 desc.add<int>("minimumNumberOfHits", 8);
0448 desc.add<double>("minimumTransverseMomentum", 1.0);
0449 desc.add<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
0450 desc.add<double>("maximumLongitudinalImpactParameter", 17.0);
0451 desc.add<bool>("computeGhostTrack", true);
0452 desc.add<double>("ghostTrackPriorDeltaR", 0.03);
0453 desc.add<edm::InputTag>("jetTracks", edm::InputTag("ak4JetTracksAssociatorAtVertexPF"));
0454 desc.add<bool>("jetDirectionUsingGhostTrack", false);
0455 desc.add<int>("minimumNumberOfPixelHits", 2);
0456 desc.add<bool>("jetDirectionUsingTracks", false);
0457 desc.add<bool>("computeProbabilities", true);
0458 desc.add<bool>("useTrackQuality", false);
0459 desc.add<double>("maximumChiSquared", 5.0);
0460 descriptions.addDefault(desc);
0461 }
0462
0463 template <>
0464 inline void
0465 IPProducer<std::vector<reco::CandidatePtr>, reco::JetTagInfo, IPProducerHelpers::FromJetAndCands>::fillDescriptions(
0466 edm::ConfigurationDescriptions& descriptions) {
0467 edm::ParameterSetDescription desc;
0468 desc.add<double>("maximumTransverseImpactParameter", 0.2);
0469 desc.add<int>("minimumNumberOfHits", 8);
0470 desc.add<double>("minimumTransverseMomentum", 1.0);
0471 desc.add<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
0472 desc.add<double>("maximumLongitudinalImpactParameter", 17.0);
0473 desc.add<bool>("computeGhostTrack", true);
0474 desc.add<double>("maxDeltaR", 0.4);
0475 desc.add<edm::InputTag>("candidates", edm::InputTag("particleFlow"));
0476 desc.add<bool>("jetDirectionUsingGhostTrack", false);
0477 desc.add<int>("minimumNumberOfPixelHits", 2);
0478 desc.add<bool>("jetDirectionUsingTracks", false);
0479 desc.add<bool>("computeProbabilities", true);
0480 desc.add<bool>("useTrackQuality", false);
0481 desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
0482 desc.add<double>("ghostTrackPriorDeltaR", 0.03);
0483 desc.add<double>("maximumChiSquared", 5.0);
0484 desc.addOptional<bool>("explicitJTA", false);
0485 descriptions.addDefault(desc);
0486 }
0487
0488 #endif