File indexing completed on 2025-02-28 23:34:30
0001
0002
0003
0004
0005
0006
0007
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/Utilities/interface/ESGetToken.h"
0015
0016 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0017 #include "DataFormats/TrackReco/interface/Track.h"
0018
0019 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
0020
0021 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0022 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0023 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0024 #include "TrackMerger.h"
0025
0026 #include "CommonTools/Utils/interface/DynArray.h"
0027 #include <vector>
0028 #include <algorithm>
0029 #include <string>
0030 #include <iostream>
0031 #include <atomic>
0032
0033 #include "CondFormats/GBRForest/interface/GBRForest.h"
0034 #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
0035 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0036
0037 using namespace reco;
0038 namespace {
0039 class DuplicateTrackMerger final : public edm::stream::EDProducer<> {
0040 public:
0041
0042 explicit DuplicateTrackMerger(const edm::ParameterSet &iPara);
0043
0044 ~DuplicateTrackMerger() override;
0045
0046 using CandidateToDuplicate = std::vector<std::pair<int, int>>;
0047
0048 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0049
0050 private:
0051
0052 void produce(edm::Event &, const edm::EventSetup &) override;
0053
0054 bool checkForDisjointTracks(const reco::Track *t1, const reco::Track *t2, TSCPBuilderNoMaterial &tscpBuilder) const;
0055 bool checkForOverlappingTracks(
0056 const reco::Track *t1, const reco::Track *t2, unsigned int nvh1, unsigned int nvh2, double cosT) const;
0057
0058 private:
0059
0060 const GBRForest *forest_;
0061
0062
0063 std::string dbFileName_;
0064 bool useForestFromDB_;
0065 std::string forestLabel_;
0066
0067 std::string propagatorName_;
0068 std::string chi2EstimatorName_;
0069
0070
0071 edm::EDGetTokenT<reco::TrackCollection> trackSource_;
0072
0073 double minDeltaR3d2_;
0074
0075 double minBDTG_;
0076
0077 double minpT2_;
0078
0079 double minP_;
0080
0081 float maxDCA2_;
0082
0083 float maxDPhi_;
0084
0085 float maxDLambda_;
0086
0087 float maxDdxy_;
0088
0089 float maxDdsz_;
0090
0091 float maxDQoP_;
0092
0093 unsigned int overlapCheckMaxHits_;
0094
0095 unsigned int overlapCheckMaxMissingLayers_;
0096
0097 double overlapCheckMinCosT_;
0098
0099 const MagneticField *magfield_;
0100 const TrackerTopology *ttopo_;
0101 const TrackerGeometry *geom_;
0102 const Propagator *propagator_;
0103 const Chi2MeasurementEstimatorBase *chi2Estimator_;
0104
0105 edm::ESGetToken<GBRForest, GBRWrapperRcd> forestToken_;
0106 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0107 edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken_;
0108 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geometryToken_;
0109 edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorToken_;
0110 edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> estimatorToken_;
0111
0112
0113 TrackMerger merger_;
0114
0115 #ifdef EDM_ML_DEBUG
0116 bool debug_;
0117 #endif
0118 };
0119 }
0120
0121 #include "FWCore/Framework/interface/Event.h"
0122 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0123 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0124 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0125 #include "FWCore/Framework/interface/EventSetup.h"
0126 #include "FWCore/Framework/interface/ESHandle.h"
0127 #include "TFile.h"
0128
0129 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0130 #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
0131
0132 #include "DuplicateTrackType.h"
0133
0134 namespace {
0135
0136 void DuplicateTrackMerger::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0137 edm::ParameterSetDescription desc;
0138 desc.add<edm::InputTag>("source", edm::InputTag());
0139 desc.add<double>("minDeltaR3d", -4.0);
0140 desc.add<double>("minBDTG", -0.1);
0141 desc.add<double>("minpT", 0.2);
0142 desc.add<double>("minP", 0.4);
0143 desc.add<double>("maxDCA", 30.0);
0144 desc.add<double>("maxDPhi", 0.30);
0145 desc.add<double>("maxDLambda", 0.30);
0146 desc.add<double>("maxDdsz", 10.0);
0147 desc.add<double>("maxDdxy", 10.0);
0148 desc.add<double>("maxDQoP", 0.25);
0149 desc.add<unsigned>("overlapCheckMaxHits", 4);
0150 desc.add<unsigned>("overlapCheckMaxMissingLayers", 1);
0151 desc.add<double>("overlapCheckMinCosT", 0.99);
0152 desc.add<std::string>("forestLabel", "MVADuplicate");
0153 desc.add<std::string>("GBRForestFileName", "");
0154 desc.add<bool>("useInnermostState", true);
0155 desc.add<std::string>("ttrhBuilderName", "WithAngleAndTemplate");
0156 desc.add<std::string>("propagatorName", "PropagatorWithMaterial");
0157 desc.add<std::string>("chi2EstimatorName", "DuplicateTrackMergerChi2Est");
0158 descriptions.add("DuplicateTrackMerger", desc);
0159 }
0160
0161 DuplicateTrackMerger::DuplicateTrackMerger(const edm::ParameterSet &iPara)
0162 : forest_(nullptr), merger_(iPara, consumesCollector()) {
0163 trackSource_ = consumes<reco::TrackCollection>(iPara.getParameter<edm::InputTag>("source"));
0164 minDeltaR3d2_ = iPara.getParameter<double>("minDeltaR3d");
0165 minDeltaR3d2_ *= std::abs(minDeltaR3d2_);
0166 minBDTG_ = iPara.getParameter<double>("minBDTG");
0167 minpT2_ = iPara.getParameter<double>("minpT");
0168 minpT2_ *= minpT2_;
0169 minP_ = iPara.getParameter<double>("minP");
0170 maxDCA2_ = iPara.getParameter<double>("maxDCA");
0171 maxDCA2_ *= maxDCA2_;
0172 maxDPhi_ = iPara.getParameter<double>("maxDPhi");
0173 maxDLambda_ = iPara.getParameter<double>("maxDLambda");
0174 maxDdsz_ = iPara.getParameter<double>("maxDdsz");
0175 maxDdxy_ = iPara.getParameter<double>("maxDdxy");
0176 maxDQoP_ = iPara.getParameter<double>("maxDQoP");
0177 overlapCheckMaxHits_ = iPara.getParameter<unsigned>("overlapCheckMaxHits");
0178 overlapCheckMaxMissingLayers_ = iPara.getParameter<unsigned>("overlapCheckMaxMissingLayers");
0179 overlapCheckMinCosT_ = iPara.getParameter<double>("overlapCheckMinCosT");
0180
0181 produces<std::vector<TrackCandidate>>("candidates");
0182 produces<CandidateToDuplicate>("candidateMap");
0183
0184 forestLabel_ = iPara.getParameter<std::string>("forestLabel");
0185
0186 dbFileName_ = iPara.getParameter<std::string>("GBRForestFileName");
0187 useForestFromDB_ = dbFileName_.empty();
0188
0189 propagatorName_ = iPara.getParameter<std::string>("propagatorName");
0190 chi2EstimatorName_ = iPara.getParameter<std::string>("chi2EstimatorName");
0191 if (useForestFromDB_) {
0192 forestToken_ = esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", forestLabel_));
0193 }
0194 magFieldToken_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0195 trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0196 geometryToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0197 propagatorToken_ = esConsumes<Propagator, TrackingComponentsRecord>(edm::ESInputTag("", propagatorName_));
0198 estimatorToken_ =
0199 esConsumes<Chi2MeasurementEstimatorBase, TrackingComponentsRecord>(edm::ESInputTag("", chi2EstimatorName_));
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214 }
0215
0216 DuplicateTrackMerger::~DuplicateTrackMerger() {
0217 if (!useForestFromDB_)
0218 delete forest_;
0219 }
0220
0221 #ifdef VI_STAT
0222 struct Stat {
0223 Stat() : maxCos(1.1), nCand(0), nLoop0(0) {}
0224 ~Stat() { std::cout << "Stats " << nCand << ' ' << nLoop0 << ' ' << maxCos << std::endl; }
0225 std::atomic<float> maxCos;
0226 std::atomic<int> nCand, nLoop0;
0227 };
0228 Stat stat;
0229 #endif
0230
0231 template <typename T>
0232 void update_maximum(std::atomic<T> &maximum_value, T const &value) noexcept {
0233 T prev_value = maximum_value;
0234 while (prev_value < value && !maximum_value.compare_exchange_weak(prev_value, value))
0235 ;
0236 }
0237
0238 template <typename T>
0239 void update_minimum(std::atomic<T> &minimum_value, T const &value) noexcept {
0240 T prev_value = minimum_value;
0241 while (prev_value > value && !minimum_value.compare_exchange_weak(prev_value, value))
0242 ;
0243 }
0244
0245 void DuplicateTrackMerger::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0246 merger_.init(iSetup);
0247
0248 if (!forest_) {
0249 if (useForestFromDB_) {
0250 edm::ESHandle<GBRForest> forestHandle = iSetup.getHandle(forestToken_);
0251 forest_ = forestHandle.product();
0252 } else {
0253 TFile gbrfile(dbFileName_.c_str());
0254 forest_ = dynamic_cast<const GBRForest *>(gbrfile.Get(forestLabel_.c_str()));
0255 }
0256 }
0257
0258
0259 edm::Handle<reco::TrackCollection> handle;
0260 iEvent.getByToken(trackSource_, handle);
0261 auto const &tracks = *handle;
0262
0263 edm::ESHandle<MagneticField> hmagfield = iSetup.getHandle(magFieldToken_);
0264 magfield_ = hmagfield.product();
0265
0266 edm::ESHandle<TrackerTopology> httopo = iSetup.getHandle(trackerTopoToken_);
0267 ttopo_ = httopo.product();
0268
0269 edm::ESHandle<TrackerGeometry> hgeom = iSetup.getHandle(geometryToken_);
0270 geom_ = hgeom.product();
0271
0272 edm::ESHandle<Propagator> hpropagator = iSetup.getHandle(propagatorToken_);
0273 propagator_ = hpropagator.product();
0274
0275 edm::ESHandle<Chi2MeasurementEstimatorBase> hestimator = iSetup.getHandle(estimatorToken_);
0276 chi2Estimator_ = hestimator.product();
0277
0278 TSCPBuilderNoMaterial tscpBuilder;
0279 auto out_duplicateCandidates = std::make_unique<std::vector<TrackCandidate>>();
0280
0281 auto out_candidateMap = std::make_unique<CandidateToDuplicate>();
0282 LogDebug("DuplicateTrackMerger") << "Number of tracks to be checked for merging: " << tracks.size();
0283
0284 #ifdef EDM_ML_DEBUG
0285 auto test = [&](const reco::Track *a, const reco::Track *b) {
0286 const auto ev = iEvent.id().event();
0287 const auto aOriAlgo = a->originalAlgo();
0288 const auto bOriAlgo = b->originalAlgo();
0289 const auto aSeedRef = a->seedRef();
0290 if (!aSeedRef) {
0291 LogDebug("DuplicateTrackMerger") << "aSeedRef is null!";
0292 return false;
0293 }
0294 const auto aSeed = a->seedRef().key();
0295 const auto bSeedRef = b->seedRef();
0296 if (!bSeedRef) {
0297 LogDebug("DuplicateTrackMerger") << "bSeedRef is null!";
0298 return false;
0299 }
0300 const auto bSeed = b->seedRef().key();
0301 return ((ev == 6903 && ((aOriAlgo == 23 && aSeed == 695 && bOriAlgo == 5 && bSeed == 652) ||
0302 (aOriAlgo == 23 && aSeed == 400 && bOriAlgo == 7 && bSeed == 156) ||
0303 (aOriAlgo == 4 && aSeed == 914 && bOriAlgo == 22 && bSeed == 503) ||
0304 (aOriAlgo == 5 && aSeed == 809 && bOriAlgo == 4 && bSeed == 1030) ||
0305 (aOriAlgo == 23 && aSeed == 749 && bOriAlgo == 5 && bSeed == 659) ||
0306 (aOriAlgo == 4 && aSeed == 1053 && bOriAlgo == 23 && bSeed == 1035) ||
0307 (aOriAlgo == 4 && aSeed == 810 && bOriAlgo == 5 && bSeed == 666) ||
0308 (aOriAlgo == 4 && aSeed == 974 && bOriAlgo == 5 && bSeed == 778))) ||
0309 (ev == 6904 && ((aOriAlgo == 23 && aSeed == 526 && bOriAlgo == 5 && bSeed == 307) ||
0310 (aOriAlgo == 4 && aSeed == 559 && bOriAlgo == 22 && bSeed == 472))) ||
0311 (ev == 6902 && ((aOriAlgo == 4 && aSeed == 750 && bOriAlgo == 22 && bSeed == 340) ||
0312 (aOriAlgo == 4 && aSeed == 906 && bOriAlgo == 5 && bSeed == 609) ||
0313 (aOriAlgo == 4 && aSeed == 724 && bOriAlgo == 5 && bSeed == 528) ||
0314 (aOriAlgo == 4 && aSeed == 943 && bOriAlgo == 23 && bSeed == 739) ||
0315 (aOriAlgo == 8 && aSeed == 2 && bOriAlgo == 9 && bSeed == 2282) ||
0316 (aOriAlgo == 23 && aSeed == 827 && bOriAlgo == 5 && bSeed == 656) ||
0317 (aOriAlgo == 22 && aSeed == 667 && bOriAlgo == 7 && bSeed == 516))));
0318 };
0319 #endif
0320
0321
0322 int nTracks = 0;
0323 declareDynArray(const reco::Track *, tracks.size(), selTracks);
0324 declareDynArray(unsigned int, tracks.size(), nValidHits);
0325 declareDynArray(unsigned int, tracks.size(), oriIndex);
0326 for (auto i = 0U; i < tracks.size(); i++) {
0327 const reco::Track *rt1 = &tracks[i];
0328 if (rt1->innerMomentum().perp2() < minpT2_)
0329 continue;
0330 selTracks[nTracks] = rt1;
0331 nValidHits[nTracks] = rt1->numberOfValidHits();
0332 oriIndex[nTracks] = i;
0333 ++nTracks;
0334 }
0335
0336 for (int i = 0; i < nTracks; i++) {
0337 const reco::Track *rt1 = selTracks[i];
0338 for (int j = i + 1; j < nTracks; j++) {
0339 const reco::Track *rt2 = selTracks[j];
0340
0341 #ifdef EDM_ML_DEBUG
0342 debug_ = false;
0343 if (test(rt1, rt2) || test(rt2, rt1)) {
0344 debug_ = true;
0345 LogTrace("DuplicateTrackMerger")
0346 << "Track1 " << i << " originalAlgo " << rt1->originalAlgo() << " seed " << rt1->seedRef().key() << " pT "
0347 << std::sqrt(rt1->innerMomentum().perp2()) << " charge " << rt1->charge() << " outerPosition2 "
0348 << rt1->outerPosition().perp2() << "\n"
0349 << "Track2 " << j << " originalAlgo " << rt2->originalAlgo() << " seed " << rt2->seedRef().key() << " pT "
0350 << std::sqrt(rt2->innerMomentum().perp2()) << " charge " << rt2->charge() << " outerPosition2 "
0351 << rt2->outerPosition().perp2();
0352 }
0353 #endif
0354
0355 if (rt1->charge() != rt2->charge())
0356 continue;
0357 auto cosT = (*rt1).momentum().Dot((*rt2).momentum());
0358 IfLogTrace(debug_, "DuplicateTrackMerger") << " cosT " << cosT;
0359 if (cosT < 0.)
0360 continue;
0361 cosT /= std::sqrt((*rt1).momentum().Mag2() * (*rt2).momentum().Mag2());
0362
0363 const reco::Track *t1, *t2;
0364 unsigned int nhv1, nhv2;
0365 if (rt1->outerPosition().perp2() < rt2->outerPosition().perp2()) {
0366 t1 = rt1;
0367 nhv1 = nValidHits[i];
0368 t2 = rt2;
0369 nhv2 = nValidHits[j];
0370 } else {
0371 t1 = rt2;
0372 nhv1 = nValidHits[j];
0373 t2 = rt1;
0374 nhv2 = nValidHits[i];
0375 }
0376 auto deltaR3d2 = (t1->outerPosition() - t2->innerPosition()).mag2();
0377
0378 if (t1->outerPosition().perp2() > t2->innerPosition().perp2())
0379 deltaR3d2 *= -1.0;
0380 IfLogTrace(debug_, "DuplicateTrackMerger")
0381 << " deltaR3d2 " << deltaR3d2 << " t1.outerPos2 " << t1->outerPosition().perp2() << " t2.innerPos2 "
0382 << t2->innerPosition().perp2();
0383
0384 bool compatible = false;
0385 DuplicateTrackType duplType;
0386 if (deltaR3d2 >= minDeltaR3d2_) {
0387 compatible = checkForDisjointTracks(t1, t2, tscpBuilder);
0388 duplType = DuplicateTrackType::Disjoint;
0389 } else {
0390 compatible = checkForOverlappingTracks(t1, t2, nhv1, nhv2, cosT);
0391 duplType = DuplicateTrackType::Overlapping;
0392 }
0393 if (!compatible)
0394 continue;
0395
0396 IfLogTrace(debug_, "DuplicateTrackMerger") << " marking as duplicates" << oriIndex[i] << ',' << oriIndex[j];
0397 out_duplicateCandidates->push_back(merger_.merge(*t1, *t2, duplType));
0398 out_candidateMap->emplace_back(oriIndex[i], oriIndex[j]);
0399
0400 #ifdef VI_STAT
0401 ++stat.nCand;
0402
0403 if (cosT > 0)
0404 update_minimum(stat.maxCos, float(cosT));
0405 else
0406 ++stat.nLoop0;
0407 #endif
0408 }
0409 }
0410 iEvent.put(std::move(out_duplicateCandidates), "candidates");
0411 iEvent.put(std::move(out_candidateMap), "candidateMap");
0412 }
0413
0414 bool DuplicateTrackMerger::checkForDisjointTracks(const reco::Track *t1,
0415 const reco::Track *t2,
0416 TSCPBuilderNoMaterial &tscpBuilder) const {
0417 IfLogTrace(debug_, "DuplicateTrackMerger") << " Checking for disjoint duplicates";
0418
0419 FreeTrajectoryState fts1 = trajectoryStateTransform::outerFreeState(*t1, &*magfield_, false);
0420 FreeTrajectoryState fts2 = trajectoryStateTransform::innerFreeState(*t2, &*magfield_, false);
0421 GlobalPoint avgPoint((t1->outerPosition().x() + t2->innerPosition().x()) * 0.5,
0422 (t1->outerPosition().y() + t2->innerPosition().y()) * 0.5,
0423 (t1->outerPosition().z() + t2->innerPosition().z()) * 0.5);
0424 TrajectoryStateClosestToPoint TSCP1 = tscpBuilder(fts1, avgPoint);
0425 TrajectoryStateClosestToPoint TSCP2 = tscpBuilder(fts2, avgPoint);
0426 IfLogTrace(debug_, "DuplicateTrackMerger")
0427 << " TSCP1.isValid " << TSCP1.isValid() << " TSCP2.isValid " << TSCP2.isValid();
0428 if (!TSCP1.isValid())
0429 return false;
0430 if (!TSCP2.isValid())
0431 return false;
0432
0433 const FreeTrajectoryState &ftsn1 = TSCP1.theState();
0434 const FreeTrajectoryState &ftsn2 = TSCP2.theState();
0435
0436 IfLogTrace(debug_, "DuplicateTrackMerger") << " DCA2 " << (ftsn2.position() - ftsn1.position()).mag2();
0437 if ((ftsn2.position() - ftsn1.position()).mag2() > maxDCA2_)
0438 return false;
0439
0440 auto qoverp1 = ftsn1.signedInverseMomentum();
0441 auto qoverp2 = ftsn2.signedInverseMomentum();
0442 float tmva_dqoverp_ = qoverp1 - qoverp2;
0443 IfLogTrace(debug_, "DuplicateTrackMerger") << " dqoverp " << tmva_dqoverp_;
0444 if (std::abs(tmva_dqoverp_) > maxDQoP_)
0445 return false;
0446
0447
0448
0449
0450
0451 float tmva_dlambda_ = ftsn2.momentum().theta() - ftsn1.momentum().theta();
0452 IfLogTrace(debug_, "DuplicateTrackMerger") << " dlambda " << tmva_dlambda_;
0453 if (std::abs(tmva_dlambda_) > maxDLambda_)
0454 return false;
0455
0456 auto phi1 = ftsn1.momentum().phi();
0457 auto phi2 = ftsn2.momentum().phi();
0458 float tmva_dphi_ = phi1 - phi2;
0459 if (std::abs(tmva_dphi_) > float(M_PI))
0460 tmva_dphi_ = 2.f * float(M_PI) - std::abs(tmva_dphi_);
0461 IfLogTrace(debug_, "DuplicateTrackMerger") << " dphi " << tmva_dphi_;
0462 if (std::abs(tmva_dphi_) > maxDPhi_)
0463 return false;
0464
0465 auto dxy1 =
0466 (-ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x()) / TSCP1.pt();
0467 auto dxy2 =
0468 (-ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x()) / TSCP2.pt();
0469 float tmva_ddxy_ = dxy1 - dxy2;
0470 IfLogTrace(debug_, "DuplicateTrackMerger") << " ddxy " << tmva_ddxy_;
0471 if (std::abs(tmva_ddxy_) > maxDdxy_)
0472 return false;
0473
0474 auto dsz1 = ftsn1.position().z() * TSCP1.pt() / TSCP1.momentum().mag() -
0475 (ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x()) /
0476 TSCP1.pt() * ftsn1.momentum().z() / ftsn1.momentum().mag();
0477 auto dsz2 = ftsn2.position().z() * TSCP2.pt() / TSCP2.momentum().mag() -
0478 (ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x()) /
0479 TSCP2.pt() * ftsn2.momentum().z() / ftsn2.momentum().mag();
0480 float tmva_ddsz_ = dsz1 - dsz2;
0481 IfLogTrace(debug_, "DuplicateTrackMerger") << " ddsz " << tmva_ddsz_;
0482 if (std::abs(tmva_ddsz_) > maxDdsz_)
0483 return false;
0484
0485 float tmva_d3dr_ = avgPoint.perp();
0486 float tmva_d3dz_ = avgPoint.z();
0487 float tmva_outer_nMissingInner_ = t2->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0488 float tmva_inner_nMissingOuter_ = t1->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
0489
0490 float gbrVals_[9];
0491 gbrVals_[0] = tmva_ddsz_;
0492 gbrVals_[1] = tmva_ddxy_;
0493 gbrVals_[2] = tmva_dphi_;
0494 gbrVals_[3] = tmva_dlambda_;
0495 gbrVals_[4] = tmva_dqoverp_;
0496 gbrVals_[5] = tmva_d3dr_;
0497 gbrVals_[6] = tmva_d3dz_;
0498 gbrVals_[7] = tmva_outer_nMissingInner_;
0499 gbrVals_[8] = tmva_inner_nMissingOuter_;
0500
0501 auto mvaBDTG = forest_->GetClassifier(gbrVals_);
0502 IfLogTrace(debug_, "DuplicateTrackMerger") << " mvaBDTG " << mvaBDTG;
0503 if (mvaBDTG < minBDTG_)
0504 return false;
0505
0506
0507 return true;
0508 }
0509
0510 bool DuplicateTrackMerger::checkForOverlappingTracks(
0511 const reco::Track *t1, const reco::Track *t2, unsigned int nvh1, unsigned int nvh2, double cosT) const {
0512
0513 if (nvh2 < nvh1) {
0514 std::swap(t1, t2);
0515 std::swap(nvh1, nvh2);
0516 }
0517
0518 IfLogTrace(debug_, "DuplicateTrackMerger")
0519 << " Checking for overlapping duplicates, cosT " << cosT << " t1 hits " << nvh1;
0520 if (cosT < overlapCheckMinCosT_)
0521 return false;
0522 if (nvh1 > overlapCheckMaxHits_)
0523 return false;
0524
0525
0526 auto findHitOnT2 = [&](const TrackingRecHit *hit1) {
0527 const auto hitDet = hit1->geographicalId().det();
0528 const auto hitSubdet = hit1->geographicalId().subdetId();
0529 const auto hitLayer = ttopo_->layer(hit1->geographicalId());
0530 return std::find_if(t2->recHitsBegin(), t2->recHitsEnd(), [&](const TrackingRecHit *hit2) {
0531 const auto &detId = hit2->geographicalId();
0532 return (detId.det() == hitDet && detId.subdetId() == hitSubdet && ttopo_->layer(detId) == hitLayer);
0533 });
0534 };
0535
0536 auto t1HitIter = t1->recHitsBegin();
0537 if (!(*t1HitIter)->isValid()) {
0538 IfLogTrace(debug_, "DuplicateTrackMerger") << " first t1 hit invalid";
0539 return false;
0540 }
0541 auto t2HitIter = findHitOnT2(*t1HitIter);
0542 if (t2HitIter == t2->recHitsEnd()) {
0543
0544
0545 ++t1HitIter;
0546 assert(t1HitIter != t1->recHitsEnd());
0547
0548 if (!(*t1HitIter)->isValid()) {
0549 IfLogTrace(debug_, "DuplicateTrackMerger") << " second t1 hit invalid";
0550 return false;
0551 }
0552 t2HitIter = findHitOnT2(*t1HitIter);
0553 if (t2HitIter == t2->recHitsEnd())
0554 return false;
0555 }
0556 IfLogTrace(debug_, "DuplicateTrackMerger")
0557 << " starting overlap check from t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0558 << std::distance(t2->recHitsBegin(), t2HitIter);
0559
0560 auto prevSubdet = (*t1HitIter)->geographicalId().subdetId();
0561 const TrajectoryStateOnSurface tsosInner = trajectoryStateTransform::innerStateOnSurface(*t2, *geom_, magfield_);
0562
0563 ++t1HitIter;
0564 ++t2HitIter;
0565 unsigned int missedLayers = 0;
0566 while (t1HitIter != t1->recHitsEnd() && t2HitIter != t2->recHitsEnd()) {
0567
0568 if ((*t1HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t1HitIter) ||
0569 (*t2HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t2HitIter)) {
0570 IfLogTrace(debug_, "DuplicateTrackMerger")
0571 << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0572 << std::distance(t2->recHitsBegin(), t2HitIter) << " either is invalid, types t1 "
0573 << (*t1HitIter)->getType() << " t2 " << (*t2HitIter)->getType();
0574 return false;
0575 }
0576
0577 const auto &t1DetId = (*t1HitIter)->geographicalId();
0578 const auto &t2DetId = (*t2HitIter)->geographicalId();
0579
0580 const auto t1Det = t1DetId.det();
0581 const auto t2Det = t2DetId.det();
0582 if (t1Det != DetId::Tracker || t2Det != DetId::Tracker) {
0583 IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
0584 << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
0585 << " either not from Tracker, dets t1 " << t1Det << " t2 " << t2Det;
0586 return false;
0587 }
0588
0589 const auto t1Subdet = t1DetId.subdetId();
0590 const auto t1Layer = ttopo_->layer(t1DetId);
0591
0592
0593 if (t1DetId == t2DetId) {
0594 if (!(*t1HitIter)->sharesInput(*t2HitIter, TrackingRecHit::all)) {
0595 IfLogTrace(debug_, "DuplicateTrackMerger")
0596 << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0597 << std::distance(t2->recHitsBegin(), t2HitIter) << " same DetId (" << t1DetId.rawId()
0598 << ") but do not share all input";
0599 return false;
0600 }
0601 } else {
0602 const auto t2Subdet = t2DetId.subdetId();
0603 const auto t2Layer = ttopo_->layer(t2DetId);
0604
0605
0606 if (t1Subdet != t2Subdet || t1Layer != t2Layer) {
0607 bool recovered = false;
0608
0609 if (t1Subdet == prevSubdet && t2Subdet != prevSubdet) {
0610
0611 ++t1HitIter;
0612 recovered = true;
0613 } else if (t1Subdet != prevSubdet && t2Subdet == prevSubdet) {
0614
0615 ++t2HitIter;
0616 recovered = true;
0617 } else if (t1Subdet == t2Subdet) {
0618 prevSubdet = t1Subdet;
0619
0620 if (t2Layer > t1Layer) {
0621
0622 ++t1HitIter;
0623 recovered = true;
0624 } else if (t1Layer > t2Layer) {
0625
0626 ++t2HitIter;
0627 recovered = true;
0628 }
0629 }
0630 if (recovered) {
0631 ++missedLayers;
0632 if (missedLayers > overlapCheckMaxMissingLayers_) {
0633 IfLogTrace(debug_, "DuplicateTrackMerger") << " max number of missed layers exceeded";
0634 return false;
0635 }
0636 continue;
0637 }
0638
0639 IfLogTrace(debug_, "DuplicateTrackMerger")
0640 << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0641 << std::distance(t2->recHitsBegin(), t2HitIter) << " are on different layers (subdet, layer) t1 "
0642 << t1Subdet << "," << t1Layer << " t2 " << t2Subdet << "," << t2Layer;
0643 return false;
0644 }
0645
0646 else if (t1Subdet != PixelSubdetector::PixelBarrel && t1Subdet != PixelSubdetector::PixelEndcap) {
0647 IfLogTrace(debug_, "DuplicateTrackMerger")
0648 << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0649 << std::distance(t2->recHitsBegin(), t2HitIter) << " are on same layer, but in non-pixel detector (det "
0650 << t1Det << " subdet " << t1Subdet << " layer " << t1Layer << ")";
0651 return false;
0652 }
0653 }
0654
0655
0656 TrajectoryStateOnSurface tsosPropagated = propagator_->propagate(tsosInner, (*t1HitIter)->det()->surface());
0657 if (!tsosPropagated.isValid()) {
0658 IfLogTrace(debug_, "DuplicateTrackMerger")
0659 << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0660 << std::distance(t2->recHitsBegin(), t2HitIter) << " TSOS not valid";
0661 return false;
0662 }
0663 auto passChi2Pair = chi2Estimator_->estimate(tsosPropagated, **t1HitIter);
0664 if (!passChi2Pair.first) {
0665 IfLogTrace(debug_, "DuplicateTrackMerger")
0666 << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0667 << std::distance(t2->recHitsBegin(), t2HitIter) << " hit chi2 compatibility failed with chi2 "
0668 << passChi2Pair.second;
0669 return false;
0670 }
0671
0672 prevSubdet = t1Subdet;
0673 ++t1HitIter;
0674 ++t2HitIter;
0675 }
0676 if (t1HitIter != t1->recHitsEnd()) {
0677 IfLogTrace(debug_, "DuplicateTrackMerger") << " hits on t2 ended before hits on t1";
0678 return false;
0679 }
0680
0681 IfLogTrace(debug_, "DuplicateTrackMerger") << " all hits on t2 are on layers whre t1 has also a hit";
0682 return true;
0683 }
0684 }
0685
0686 #include "FWCore/PluginManager/interface/ModuleDef.h"
0687 #include "FWCore/Framework/interface/MakerMacros.h"
0688
0689 DEFINE_FWK_MODULE(DuplicateTrackMerger);