CMS 3D CMS Logo
DuplicateTrackMerger.cc
Go to the documentation of this file.
1 
8 #include "FWCore/Framework/interface/stream/EDProducer.h"
9 
10 #include "FWCore/ParameterSet/interface/ParameterSet.h"
11 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
12 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
13 #include "FWCore/MessageLogger/interface/MessageLogger.h"
14 #include "FWCore/Utilities/interface/ESGetToken.h"
15 
16 #include "DataFormats/TrackReco/interface/TrackFwd.h"
17 #include "DataFormats/TrackReco/interface/Track.h"
18 
19 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
20 
21 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
22 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
23 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
24 #include "TrackMerger.h"
25 
26 #include "CommonTools/Utils/interface/DynArray.h"
27 #include <vector>
28 #include <algorithm>
29 #include <string>
30 #include <iostream>
31 #include <atomic>
32 
33 #include "CondFormats/EgammaObjects/interface/GBRForest.h"
34 #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
35 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
36 
37 using namespace reco;
38 namespace {
39  class DuplicateTrackMerger final : public edm::stream::EDProducer<> {
40  public:
42  explicit DuplicateTrackMerger(const edm::ParameterSet &iPara);
44  ~DuplicateTrackMerger() override;
45 
46  using CandidateToDuplicate = std::vector<std::pair<int, int>>;
47 
48  static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
49 
50  private:
52  void produce(edm::Event &, const edm::EventSetup &) override;
53 
54  bool checkForDisjointTracks(const reco::Track *t1, const reco::Track *t2, TSCPBuilderNoMaterial &tscpBuilder) const;
55  bool checkForOverlappingTracks(
56  const reco::Track *t1, const reco::Track *t2, unsigned int nvh1, unsigned int nvh2, double cosT) const;
57 
58  private:
60  const GBRForest *forest_;
61 
63  std::string dbFileName_;
64  bool useForestFromDB_;
65  std::string forestLabel_;
66 
67  std::string propagatorName_;
68  std::string chi2EstimatorName_;
69 
71  edm::EDGetTokenT<reco::TrackCollection> trackSource_;
73  double minDeltaR3d2_;
75  double minBDTG_;
77  double minpT2_;
79  double minP_;
81  float maxDCA2_;
83  float maxDPhi_;
85  float maxDLambda_;
87  float maxDdxy_;
89  float maxDdsz_;
91  float maxDQoP_;
93  unsigned int overlapCheckMaxHits_;
95  unsigned int overlapCheckMaxMissingLayers_;
97  double overlapCheckMinCosT_;
98 
99  const MagneticField *magfield_;
100  const TrackerTopology *ttopo_;
101  const TrackerGeometry *geom_;
102  const Propagator *propagator_;
103  const Chi2MeasurementEstimatorBase *chi2Estimator_;
104 
105  edm::ESGetToken<GBRForest, GBRWrapperRcd> forestToken_;
106  edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
107  edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken_;
108  edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geometryToken_;
109  edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorToken_;
110  edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> estimatorToken_;
111 
113  TrackMerger merger_;
114 
115 #ifdef EDM_ML_DEBUG
116  bool debug_;
117 #endif
118  };
119 } // namespace
120 
121 #include "FWCore/Framework/interface/Event.h"
122 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
123 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
124 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
125 #include "FWCore/Framework/interface/EventSetup.h"
126 #include "FWCore/Framework/interface/ESHandle.h"
127 #include "TFile.h"
128 
129 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
130 #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
131 
132 #include "DuplicateTrackType.h"
133 
134 namespace {
135 
136  void DuplicateTrackMerger::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
137  edm::ParameterSetDescription desc;
138  desc.add<edm::InputTag>("source", edm::InputTag());
139  desc.add<double>("minDeltaR3d", -4.0);
140  desc.add<double>("minBDTG", -0.1);
141  desc.add<double>("minpT", 0.2);
142  desc.add<double>("minP", 0.4);
143  desc.add<double>("maxDCA", 30.0);
144  desc.add<double>("maxDPhi", 0.30);
145  desc.add<double>("maxDLambda", 0.30);
146  desc.add<double>("maxDdsz", 10.0);
147  desc.add<double>("maxDdxy", 10.0);
148  desc.add<double>("maxDQoP", 0.25);
149  desc.add<unsigned>("overlapCheckMaxHits", 4);
150  desc.add<unsigned>("overlapCheckMaxMissingLayers", 1);
151  desc.add<double>("overlapCheckMinCosT", 0.99);
152  desc.add<std::string>("forestLabel", "MVADuplicate");
153  desc.add<std::string>("GBRForestFileName", "");
154  desc.add<bool>("useInnermostState", true);
155  desc.add<std::string>("ttrhBuilderName", "WithAngleAndTemplate");
156  desc.add<std::string>("propagatorName", "PropagatorWithMaterial");
157  desc.add<std::string>("chi2EstimatorName", "DuplicateTrackMergerChi2Est");
158  descriptions.add("DuplicateTrackMerger", desc);
159  }
160 
161  DuplicateTrackMerger::DuplicateTrackMerger(const edm::ParameterSet &iPara)
162  : forest_(nullptr), merger_(iPara, consumesCollector()) {
163  trackSource_ = consumes<reco::TrackCollection>(iPara.getParameter<edm::InputTag>("source"));
164  minDeltaR3d2_ = iPara.getParameter<double>("minDeltaR3d");
165  minDeltaR3d2_ *= std::abs(minDeltaR3d2_);
166  minBDTG_ = iPara.getParameter<double>("minBDTG");
167  minpT2_ = iPara.getParameter<double>("minpT");
168  minpT2_ *= minpT2_;
169  minP_ = iPara.getParameter<double>("minP");
170  maxDCA2_ = iPara.getParameter<double>("maxDCA");
171  maxDCA2_ *= maxDCA2_;
172  maxDPhi_ = iPara.getParameter<double>("maxDPhi");
173  maxDLambda_ = iPara.getParameter<double>("maxDLambda");
174  maxDdsz_ = iPara.getParameter<double>("maxDdsz");
175  maxDdxy_ = iPara.getParameter<double>("maxDdxy");
176  maxDQoP_ = iPara.getParameter<double>("maxDQoP");
177  overlapCheckMaxHits_ = iPara.getParameter<unsigned>("overlapCheckMaxHits");
178  overlapCheckMaxMissingLayers_ = iPara.getParameter<unsigned>("overlapCheckMaxMissingLayers");
179  overlapCheckMinCosT_ = iPara.getParameter<double>("overlapCheckMinCosT");
180 
181  produces<std::vector<TrackCandidate>>("candidates");
182  produces<CandidateToDuplicate>("candidateMap");
183 
184  forestLabel_ = iPara.getParameter<std::string>("forestLabel");
185 
186  dbFileName_ = iPara.getParameter<std::string>("GBRForestFileName");
187  useForestFromDB_ = dbFileName_.empty();
188 
189  propagatorName_ = iPara.getParameter<std::string>("propagatorName");
190  chi2EstimatorName_ = iPara.getParameter<std::string>("chi2EstimatorName");
191  if (useForestFromDB_) {
192  forestToken_ = esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", forestLabel_));
193  }
194  magFieldToken_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
195  trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
196  geometryToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
197  propagatorToken_ = esConsumes<Propagator, TrackingComponentsRecord>(edm::ESInputTag("", propagatorName_));
198  estimatorToken_ =
199  esConsumes<Chi2MeasurementEstimatorBase, TrackingComponentsRecord>(edm::ESInputTag("", chi2EstimatorName_));
200 
201  /*
202  tmvaReader_ = new TMVA::Reader("!Color:Silent");
203  tmvaReader_->AddVariable("ddsz",&tmva_ddsz_);
204  tmvaReader_->AddVariable("ddxy",&tmva_ddxy_);
205  tmvaReader_->AddVariable("dphi",&tmva_dphi_);
206  tmvaReader_->AddVariable("dlambda",&tmva_dlambda_);
207  tmvaReader_->AddVariable("dqoverp",&tmva_dqoverp_);
208  tmvaReader_->AddVariable("delta3d_r",&tmva_d3dr_);
209  tmvaReader_->AddVariable("delta3d_z",&tmva_d3dz_);
210  tmvaReader_->AddVariable("outer_nMissingInner",&tmva_outer_nMissingInner_);
211  tmvaReader_->AddVariable("inner_nMissingOuter",&tmva_inner_nMissingOuter_);
212  tmvaReader_->BookMVA("BDTG",mvaFilePath);
213  */
214  }
215 
216  DuplicateTrackMerger::~DuplicateTrackMerger() {
217  if (!useForestFromDB_)
218  delete forest_;
219  }
220 
221 #ifdef VI_STAT
222  struct Stat {
223  Stat() : maxCos(1.1), nCand(0), nLoop0(0) {}
224  ~Stat() { std::cout << "Stats " << nCand << ' ' << nLoop0 << ' ' << maxCos << std::endl; }
225  std::atomic<float> maxCos;
226  std::atomic<int> nCand, nLoop0;
227  };
228  Stat stat;
229 #endif
230 
231  template <typename T>
232  void update_maximum(std::atomic<T> &maximum_value, T const &value) noexcept {
233  T prev_value = maximum_value;
234  while (prev_value < value && !maximum_value.compare_exchange_weak(prev_value, value))
235  ;
236  }
237 
238  template <typename T>
239  void update_minimum(std::atomic<T> &minimum_value, T const &value) noexcept {
240  T prev_value = minimum_value;
241  while (prev_value > value && !minimum_value.compare_exchange_weak(prev_value, value))
242  ;
243  }
244 
245  void DuplicateTrackMerger::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
246  merger_.init(iSetup);
247 
248  if (!forest_) {
249  if (useForestFromDB_) {
250  edm::ESHandle<GBRForest> forestHandle = iSetup.getHandle(forestToken_);
251  forest_ = forestHandle.product();
252  } else {
253  TFile gbrfile(dbFileName_.c_str());
254  forest_ = dynamic_cast<const GBRForest *>(gbrfile.Get(forestLabel_.c_str()));
255  }
256  }
257 
258  //edm::Handle<edm::View<reco::Track> >handle;
259  edm::Handle<reco::TrackCollection> handle;
260  iEvent.getByToken(trackSource_, handle);
261  auto const &tracks = *handle;
262 
263  edm::ESHandle<MagneticField> hmagfield = iSetup.getHandle(magFieldToken_);
264  magfield_ = hmagfield.product();
265 
266  edm::ESHandle<TrackerTopology> httopo = iSetup.getHandle(trackerTopoToken_);
267  ttopo_ = httopo.product();
268 
269  edm::ESHandle<TrackerGeometry> hgeom = iSetup.getHandle(geometryToken_);
270  geom_ = hgeom.product();
271 
272  edm::ESHandle<Propagator> hpropagator = iSetup.getHandle(propagatorToken_);
273  propagator_ = hpropagator.product();
274 
275  edm::ESHandle<Chi2MeasurementEstimatorBase> hestimator = iSetup.getHandle(estimatorToken_);
276  chi2Estimator_ = hestimator.product();
277 
278  TSCPBuilderNoMaterial tscpBuilder;
279  auto out_duplicateCandidates = std::make_unique<std::vector<TrackCandidate>>();
280 
281  auto out_candidateMap = std::make_unique<CandidateToDuplicate>();
282  LogDebug("DuplicateTrackMerger") << "Number of tracks to be checked for merging: " << tracks.size();
283 
284 #ifdef EDM_ML_DEBUG
285  auto test = [&](const reco::Track *a, const reco::Track *b) {
286  const auto ev = iEvent.id().event();
287  const auto aOriAlgo = a->originalAlgo();
288  const auto bOriAlgo = b->originalAlgo();
289  const auto aSeed = a->seedRef().key();
290  const auto bSeed = b->seedRef().key();
291  return ((ev == 6903 && ((aOriAlgo == 23 && aSeed == 695 && bOriAlgo == 5 && bSeed == 652) ||
292  (aOriAlgo == 23 && aSeed == 400 && bOriAlgo == 7 && bSeed == 156) ||
293  (aOriAlgo == 4 && aSeed == 914 && bOriAlgo == 22 && bSeed == 503) ||
294  (aOriAlgo == 5 && aSeed == 809 && bOriAlgo == 4 && bSeed == 1030) ||
295  (aOriAlgo == 23 && aSeed == 749 && bOriAlgo == 5 && bSeed == 659) ||
296  (aOriAlgo == 4 && aSeed == 1053 && bOriAlgo == 23 && bSeed == 1035) ||
297  (aOriAlgo == 4 && aSeed == 810 && bOriAlgo == 5 && bSeed == 666) ||
298  (aOriAlgo == 4 && aSeed == 974 && bOriAlgo == 5 && bSeed == 778))) ||
299  (ev == 6904 && ((aOriAlgo == 23 && aSeed == 526 && bOriAlgo == 5 && bSeed == 307) ||
300  (aOriAlgo == 4 && aSeed == 559 && bOriAlgo == 22 && bSeed == 472))) ||
301  (ev == 6902 && ((aOriAlgo == 4 && aSeed == 750 && bOriAlgo == 22 && bSeed == 340) ||
302  (aOriAlgo == 4 && aSeed == 906 && bOriAlgo == 5 && bSeed == 609) ||
303  (aOriAlgo == 4 && aSeed == 724 && bOriAlgo == 5 && bSeed == 528) ||
304  (aOriAlgo == 4 && aSeed == 943 && bOriAlgo == 23 && bSeed == 739) ||
305  (aOriAlgo == 8 && aSeed == 2 && bOriAlgo == 9 && bSeed == 2282) ||
306  (aOriAlgo == 23 && aSeed == 827 && bOriAlgo == 5 && bSeed == 656) ||
307  (aOriAlgo == 22 && aSeed == 667 && bOriAlgo == 7 && bSeed == 516))));
308  };
309 #endif
310 
311  // cache few "heavy to compute quantities
312  int nTracks = 0;
313  declareDynArray(const reco::Track *, tracks.size(), selTracks);
314  declareDynArray(unsigned int, tracks.size(), nValidHits);
315  declareDynArray(unsigned int, tracks.size(), oriIndex);
316  for (auto i = 0U; i < tracks.size(); i++) {
317  const reco::Track *rt1 = &tracks[i];
318  if (rt1->innerMomentum().perp2() < minpT2_)
319  continue;
320  selTracks[nTracks] = rt1;
321  nValidHits[nTracks] = rt1->numberOfValidHits(); // yes it is extremely heavy!
322  oriIndex[nTracks] = i;
323  ++nTracks;
324  }
325 
326  for (int i = 0; i < nTracks; i++) {
327  const reco::Track *rt1 = selTracks[i];
328  for (int j = i + 1; j < nTracks; j++) {
329  const reco::Track *rt2 = selTracks[j];
330 
331 #ifdef EDM_ML_DEBUG
332  debug_ = false;
333  if (test(rt1, rt2) || test(rt2, rt1)) {
334  debug_ = true;
335  LogTrace("DuplicateTrackMerger")
336  << "Track1 " << i << " originalAlgo " << rt1->originalAlgo() << " seed " << rt1->seedRef().key() << " pT "
337  << std::sqrt(rt1->innerMomentum().perp2()) << " charge " << rt1->charge() << " outerPosition2 "
338  << rt1->outerPosition().perp2() << "\n"
339  << "Track2 " << j << " originalAlgo " << rt2->originalAlgo() << " seed " << rt2->seedRef().key() << " pT "
340  << std::sqrt(rt2->innerMomentum().perp2()) << " charge " << rt2->charge() << " outerPosition2 "
341  << rt2->outerPosition().perp2();
342  }
343 #endif
344 
345  if (rt1->charge() != rt2->charge())
346  continue;
347  auto cosT = (*rt1).momentum().Dot((*rt2).momentum()); // not normalized!
348  IfLogTrace(debug_, "DuplicateTrackMerger") << " cosT " << cosT;
349  if (cosT < 0.)
350  continue;
351  cosT /= std::sqrt((*rt1).momentum().Mag2() * (*rt2).momentum().Mag2());
352 
353  const reco::Track *t1, *t2;
354  unsigned int nhv1, nhv2;
355  if (rt1->outerPosition().perp2() < rt2->outerPosition().perp2()) {
356  t1 = rt1;
357  nhv1 = nValidHits[i];
358  t2 = rt2;
359  nhv2 = nValidHits[j];
360  } else {
361  t1 = rt2;
362  nhv1 = nValidHits[j];
363  t2 = rt1;
364  nhv2 = nValidHits[i];
365  }
366  auto deltaR3d2 = (t1->outerPosition() - t2->innerPosition()).mag2();
367 
368  if (t1->outerPosition().perp2() > t2->innerPosition().perp2())
369  deltaR3d2 *= -1.0;
370  IfLogTrace(debug_, "DuplicateTrackMerger")
371  << " deltaR3d2 " << deltaR3d2 << " t1.outerPos2 " << t1->outerPosition().perp2() << " t2.innerPos2 "
372  << t2->innerPosition().perp2();
373 
374  bool compatible = false;
375  DuplicateTrackType duplType;
376  if (deltaR3d2 >= minDeltaR3d2_) {
377  compatible = checkForDisjointTracks(t1, t2, tscpBuilder);
378  duplType = DuplicateTrackType::Disjoint;
379  } else {
380  compatible = checkForOverlappingTracks(t1, t2, nhv1, nhv2, cosT);
381  duplType = DuplicateTrackType::Overlapping;
382  }
383  if (!compatible)
384  continue;
385 
386  IfLogTrace(debug_, "DuplicateTrackMerger") << " marking as duplicates" << oriIndex[i] << ',' << oriIndex[j];
387  out_duplicateCandidates->push_back(merger_.merge(*t1, *t2, duplType));
388  out_candidateMap->emplace_back(oriIndex[i], oriIndex[j]);
389 
390 #ifdef VI_STAT
391  ++stat.nCand;
392  // auto cosT = float((*t1).momentum().unit().Dot((*t2).momentum().unit()));
393  if (cosT > 0)
394  update_minimum(stat.maxCos, float(cosT));
395  else
396  ++stat.nLoop0;
397 #endif
398  }
399  }
400  iEvent.put(std::move(out_duplicateCandidates), "candidates");
401  iEvent.put(std::move(out_candidateMap), "candidateMap");
402  }
403 
404  bool DuplicateTrackMerger::checkForDisjointTracks(const reco::Track *t1,
405  const reco::Track *t2,
406  TSCPBuilderNoMaterial &tscpBuilder) const {
407  IfLogTrace(debug_, "DuplicateTrackMerger") << " Checking for disjoint duplicates";
408 
409  FreeTrajectoryState fts1 = trajectoryStateTransform::outerFreeState(*t1, &*magfield_, false);
410  FreeTrajectoryState fts2 = trajectoryStateTransform::innerFreeState(*t2, &*magfield_, false);
411  GlobalPoint avgPoint((t1->outerPosition().x() + t2->innerPosition().x()) * 0.5,
412  (t1->outerPosition().y() + t2->innerPosition().y()) * 0.5,
413  (t1->outerPosition().z() + t2->innerPosition().z()) * 0.5);
414  TrajectoryStateClosestToPoint TSCP1 = tscpBuilder(fts1, avgPoint);
415  TrajectoryStateClosestToPoint TSCP2 = tscpBuilder(fts2, avgPoint);
416  IfLogTrace(debug_, "DuplicateTrackMerger")
417  << " TSCP1.isValid " << TSCP1.isValid() << " TSCP2.isValid " << TSCP2.isValid();
418  if (!TSCP1.isValid())
419  return false;
420  if (!TSCP2.isValid())
421  return false;
422 
423  const FreeTrajectoryState &ftsn1 = TSCP1.theState();
424  const FreeTrajectoryState &ftsn2 = TSCP2.theState();
425 
426  IfLogTrace(debug_, "DuplicateTrackMerger") << " DCA2 " << (ftsn2.position() - ftsn1.position()).mag2();
427  if ((ftsn2.position() - ftsn1.position()).mag2() > maxDCA2_)
428  return false;
429 
430  auto qoverp1 = ftsn1.signedInverseMomentum();
431  auto qoverp2 = ftsn2.signedInverseMomentum();
432  float tmva_dqoverp_ = qoverp1 - qoverp2;
433  IfLogTrace(debug_, "DuplicateTrackMerger") << " dqoverp " << tmva_dqoverp_;
434  if (std::abs(tmva_dqoverp_) > maxDQoP_)
435  return false;
436 
437  //auto pp = [&](TrajectoryStateClosestToPoint const & ts) { std::cout << ' ' << ts.perigeeParameters().vector()[0] << '/' << ts.perigeeError().transverseCurvatureError();};
438  //if(qoverp1*qoverp2 <0) { std::cout << "charge different " << qoverp1 <<',' << qoverp2; pp(TSCP1); pp(TSCP2); std::cout << std::endl;}
439 
440  // lambda = pi/2 - theta.... so l1-l2 == t2-t1
441  float tmva_dlambda_ = ftsn2.momentum().theta() - ftsn1.momentum().theta();
442  IfLogTrace(debug_, "DuplicateTrackMerger") << " dlambda " << tmva_dlambda_;
443  if (std::abs(tmva_dlambda_) > maxDLambda_)
444  return false;
445 
446  auto phi1 = ftsn1.momentum().phi();
447  auto phi2 = ftsn2.momentum().phi();
448  float tmva_dphi_ = phi1 - phi2;
449  if (std::abs(tmva_dphi_) > float(M_PI))
450  tmva_dphi_ = 2.f * float(M_PI) - std::abs(tmva_dphi_);
451  IfLogTrace(debug_, "DuplicateTrackMerger") << " dphi " << tmva_dphi_;
452  if (std::abs(tmva_dphi_) > maxDPhi_)
453  return false;
454 
455  auto dxy1 =
456  (-ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x()) / TSCP1.pt();
457  auto dxy2 =
458  (-ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x()) / TSCP2.pt();
459  float tmva_ddxy_ = dxy1 - dxy2;
460  IfLogTrace(debug_, "DuplicateTrackMerger") << " ddxy " << tmva_ddxy_;
461  if (std::abs(tmva_ddxy_) > maxDdxy_)
462  return false;
463 
464  auto dsz1 = ftsn1.position().z() * TSCP1.pt() / TSCP1.momentum().mag() -
465  (ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x()) /
466  TSCP1.pt() * ftsn1.momentum().z() / ftsn1.momentum().mag();
467  auto dsz2 = ftsn2.position().z() * TSCP2.pt() / TSCP2.momentum().mag() -
468  (ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x()) /
469  TSCP2.pt() * ftsn2.momentum().z() / ftsn2.momentum().mag();
470  float tmva_ddsz_ = dsz1 - dsz2;
471  IfLogTrace(debug_, "DuplicateTrackMerger") << " ddsz " << tmva_ddsz_;
472  if (std::abs(tmva_ddsz_) > maxDdsz_)
473  return false;
474 
475  float tmva_d3dr_ = avgPoint.perp();
476  float tmva_d3dz_ = avgPoint.z();
477  float tmva_outer_nMissingInner_ = t2->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
478  float tmva_inner_nMissingOuter_ = t1->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
479 
480  float gbrVals_[9];
481  gbrVals_[0] = tmva_ddsz_;
482  gbrVals_[1] = tmva_ddxy_;
483  gbrVals_[2] = tmva_dphi_;
484  gbrVals_[3] = tmva_dlambda_;
485  gbrVals_[4] = tmva_dqoverp_;
486  gbrVals_[5] = tmva_d3dr_;
487  gbrVals_[6] = tmva_d3dz_;
488  gbrVals_[7] = tmva_outer_nMissingInner_;
489  gbrVals_[8] = tmva_inner_nMissingOuter_;
490 
491  auto mvaBDTG = forest_->GetClassifier(gbrVals_);
492  IfLogTrace(debug_, "DuplicateTrackMerger") << " mvaBDTG " << mvaBDTG;
493  if (mvaBDTG < minBDTG_)
494  return false;
495 
496  // std::cout << "to merge " << mvaBDTG << ' ' << std::copysign(std::sqrt(std::abs(deltaR3d2)),deltaR3d2) << ' ' << tmva_dphi_ << ' ' << TSCP1.pt() <<'/'<<TSCP2.pt() << std::endl;
497  return true;
498  }
499 
500  bool DuplicateTrackMerger::checkForOverlappingTracks(
501  const reco::Track *t1, const reco::Track *t2, unsigned int nvh1, unsigned int nvh2, double cosT) const {
502  // ensure t1 is the shorter track
503  if (nvh2 < nvh1) {
504  std::swap(t1, t2);
505  std::swap(nvh1, nvh2);
506  }
507 
508  IfLogTrace(debug_, "DuplicateTrackMerger")
509  << " Checking for overlapping duplicates, cosT " << cosT << " t1 hits " << nvh1;
510  if (cosT < overlapCheckMinCosT_)
511  return false;
512  if (nvh1 > overlapCheckMaxHits_)
513  return false;
514 
515  // find the hit on the longer track on layer of the first hit of the shorter track
516  auto findHitOnT2 = [&](const TrackingRecHit *hit1) {
517  const auto hitDet = hit1->geographicalId().det();
518  const auto hitSubdet = hit1->geographicalId().subdetId();
519  const auto hitLayer = ttopo_->layer(hit1->geographicalId());
520  return std::find_if(t2->recHitsBegin(), t2->recHitsEnd(), [&](const TrackingRecHit *hit2) {
521  const auto &detId = hit2->geographicalId();
522  return (detId.det() == hitDet && detId.subdetId() == hitSubdet && ttopo_->layer(detId) == hitLayer);
523  });
524  };
525 
526  auto t1HitIter = t1->recHitsBegin();
527  if (!(*t1HitIter)->isValid()) {
528  IfLogTrace(debug_, "DuplicateTrackMerger") << " first t1 hit invalid";
529  return false;
530  }
531  auto t2HitIter = findHitOnT2(*t1HitIter);
532  if (t2HitIter == t2->recHitsEnd()) {
533  // if first hit not found, try with second
534  // if that fails, then reject
535  ++t1HitIter;
536  assert(t1HitIter != t1->recHitsEnd());
537 
538  if (!(*t1HitIter)->isValid()) {
539  IfLogTrace(debug_, "DuplicateTrackMerger") << " second t1 hit invalid";
540  return false;
541  }
542  t2HitIter = findHitOnT2(*t1HitIter);
543  if (t2HitIter == t2->recHitsEnd())
544  return false;
545  }
546  IfLogTrace(debug_, "DuplicateTrackMerger")
547  << " starting overlap check from t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
548  << std::distance(t2->recHitsBegin(), t2HitIter);
549 
550  auto prevSubdet = (*t1HitIter)->geographicalId().subdetId();
551  const TrajectoryStateOnSurface tsosInner = trajectoryStateTransform::innerStateOnSurface(*t2, *geom_, magfield_);
552 
553  ++t1HitIter;
554  ++t2HitIter;
555  unsigned int missedLayers = 0;
556  while (t1HitIter != t1->recHitsEnd() && t2HitIter != t2->recHitsEnd()) {
557  // in case of invalid hits, reject immediately
558  if ((*t1HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t1HitIter) ||
559  (*t2HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t2HitIter)) {
560  IfLogTrace(debug_, "DuplicateTrackMerger")
561  << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
562  << std::distance(t2->recHitsBegin(), t2HitIter) << " either is invalid, types t1 "
563  << (*t1HitIter)->getType() << " t2 " << (*t2HitIter)->getType();
564  return false;
565  }
566 
567  const auto &t1DetId = (*t1HitIter)->geographicalId();
568  const auto &t2DetId = (*t2HitIter)->geographicalId();
569 
570  const auto t1Det = t1DetId.det();
571  const auto t2Det = t2DetId.det();
572  if (t1Det != DetId::Tracker || t2Det != DetId::Tracker) {
573  IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
574  << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
575  << " either not from Tracker, dets t1 " << t1Det << " t2 " << t2Det;
576  return false;
577  }
578 
579  const auto t1Subdet = t1DetId.subdetId();
580  const auto t1Layer = ttopo_->layer(t1DetId);
581 
582  // reject if hits have the same DetId but are different
583  if (t1DetId == t2DetId) {
584  if (!(*t1HitIter)->sharesInput(*t2HitIter, TrackingRecHit::all)) {
585  IfLogTrace(debug_, "DuplicateTrackMerger")
586  << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
587  << std::distance(t2->recHitsBegin(), t2HitIter) << " same DetId (" << t1DetId.rawId()
588  << ") but do not share all input";
589  return false;
590  }
591  } else {
592  const auto t2Subdet = t2DetId.subdetId();
593  const auto t2Layer = ttopo_->layer(t2DetId);
594 
595  // reject if hits are on different layers
596  if (t1Subdet != t2Subdet || t1Layer != t2Layer) {
597  bool recovered = false;
598  // but try to recover first by checking if either one has skipped over a layer
599  if (t1Subdet == prevSubdet && t2Subdet != prevSubdet) {
600  // t1 has a layer t2 doesn't
601  ++t1HitIter;
602  recovered = true;
603  } else if (t1Subdet != prevSubdet && t2Subdet == prevSubdet) {
604  // t2 has a layer t1 doesn't
605  ++t2HitIter;
606  recovered = true;
607  } else if (t1Subdet == t2Subdet) {
608  prevSubdet = t1Subdet;
609  // same subdet, so layer must be different
610  if (t2Layer > t1Layer) {
611  // t1 has a layer t2 doesn't
612  ++t1HitIter;
613  recovered = true;
614  } else if (t1Layer > t2Layer) {
615  // t2 has a layer t1 doesn't
616  ++t2HitIter;
617  recovered = true;
618  }
619  }
620  if (recovered) {
621  ++missedLayers;
622  if (missedLayers > overlapCheckMaxMissingLayers_) {
623  IfLogTrace(debug_, "DuplicateTrackMerger") << " max number of missed layers exceeded";
624  return false;
625  }
626  continue;
627  }
628 
629  IfLogTrace(debug_, "DuplicateTrackMerger")
630  << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
631  << std::distance(t2->recHitsBegin(), t2HitIter) << " are on different layers (subdet, layer) t1 "
632  << t1Subdet << "," << t1Layer << " t2 " << t2Subdet << "," << t2Layer;
633  return false;
634  }
635  // reject if same layer (but not same hit) in non-pixel detector
636  else if (t1Subdet != PixelSubdetector::PixelBarrel && t1Subdet != PixelSubdetector::PixelEndcap) {
637  IfLogTrace(debug_, "DuplicateTrackMerger")
638  << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
639  << std::distance(t2->recHitsBegin(), t2HitIter) << " are on same layer, but in non-pixel detector (det "
640  << t1Det << " subdet " << t1Subdet << " layer " << t1Layer << ")";
641  return false;
642  }
643  }
644 
645  // Propagate longer track to the shorter track hit surface, check compatibility
646  TrajectoryStateOnSurface tsosPropagated = propagator_->propagate(tsosInner, (*t1HitIter)->det()->surface());
647  if (!tsosPropagated.isValid()) { // reject immediately if TSOS is not valid
648  IfLogTrace(debug_, "DuplicateTrackMerger")
649  << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
650  << std::distance(t2->recHitsBegin(), t2HitIter) << " TSOS not valid";
651  return false;
652  }
653  auto passChi2Pair = chi2Estimator_->estimate(tsosPropagated, **t1HitIter);
654  if (!passChi2Pair.first) {
655  IfLogTrace(debug_, "DuplicateTrackMerger")
656  << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
657  << std::distance(t2->recHitsBegin(), t2HitIter) << " hit chi2 compatibility failed with chi2 "
658  << passChi2Pair.second;
659  return false;
660  }
661 
662  prevSubdet = t1Subdet;
663  ++t1HitIter;
664  ++t2HitIter;
665  }
666  if (t1HitIter != t1->recHitsEnd()) {
667  IfLogTrace(debug_, "DuplicateTrackMerger") << " hits on t2 ended before hits on t1";
668  return false;
669  }
670 
671  IfLogTrace(debug_, "DuplicateTrackMerger") << " all hits on t2 are on layers whre t1 has also a hit";
672  return true;
673  }
674 } // namespace
675 
676 #include "FWCore/PluginManager/interface/ModuleDef.h"
677 #include "FWCore/Framework/interface/MakerMacros.h"
678 
679 DEFINE_FWK_MODULE(DuplicateTrackMerger);
reco::Track::outerPosition
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:62
ConfigurationDescriptions.h
Chi2MeasurementEstimatorBase.h
Propagator.h
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
RandomServiceHelper.t2
t2
Definition: RandomServiceHelper.py:257
reco::HitPattern::MISSING_OUTER_HITS
Definition: HitPattern.h:155
FreeTrajectoryState::momentum
GlobalVector momentum() const
Definition: FreeTrajectoryState.h:68
TrajectoryStateOnSurface.h
PDWG_EXOHSCP_cff.tracks
tracks
Definition: PDWG_EXOHSCP_cff.py:28
mps_fire.i
i
Definition: mps_fire.py:428
edm::ESInputTag
Definition: ESInputTag.h:87
PixelSubdetector.h
FreeTrajectoryState.h
MessageLogger.h
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
PixelSubdetector::PixelEndcap
Definition: PixelSubdetector.h:11
trajectoryStateTransform::outerFreeState
FreeTrajectoryState outerFreeState(const reco::Track &tk, const MagneticField *field, bool withErr=true)
Definition: TrajectoryStateTransform.cc:98
PixelSubdetector::PixelBarrel
Definition: PixelSubdetector.h:11
ESHandle.h
edm::RefToBase::key
size_t key() const
Definition: RefToBase.h:219
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
DuplicateTrackType.h
patZpeak.handle
handle
Definition: patZpeak.py:23
edm::EDGetTokenT< reco::TrackCollection >
TrackerTopology
Definition: TrackerTopology.h:16
GBRWrapperRcd.h
PV3DBase::theta
Geom::Theta< T > theta() const
Definition: PV3DBase.h:72
GBRForest
Definition: GBRForest.h:25
gather_cfg.cout
cout
Definition: gather_cfg.py:144
HLT_FULL_cff.InputTag
InputTag
Definition: HLT_FULL_cff.py:89287
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
cms::cuda::assert
assert(be >=bs)
EDProducer.h
trackerHitRTTI::isUndef
bool isUndef(TrackingRecHit const &hit)
Definition: trackerHitRTTI.h:24
FreeTrajectoryState::signedInverseMomentum
double signedInverseMomentum() const
Definition: FreeTrajectoryState.h:70
reco::TrackBase::originalAlgo
TrackAlgorithm originalAlgo() const
Definition: TrackBase.h:548
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
FreeTrajectoryState::position
GlobalPoint position() const
Definition: FreeTrajectoryState.h:67
BeamSpotPI::nTracks
Definition: BeamSpotPayloadInspectorHelper.h:41
GBRForest.h
edm::Handle< reco::TrackCollection >
TrackMerger
Definition: TrackMerger.h:19
reco::TrackBase::numberOfValidHits
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:798
ESGetToken.h
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
Propagator
Definition: Propagator.h:44
TrajectoryStateOnSurface
Definition: TrajectoryStateOnSurface.h:16
MakerMacros.h
reco::Track::innerMomentum
const math::XYZVector & innerMomentum() const
momentum vector at the innermost hit position
Definition: Track.h:59
test
Definition: SmallWORMDict.h:13
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
ctpps_dqm_sourceclient-live_cfg.test
test
Definition: ctpps_dqm_sourceclient-live_cfg.py:7
Track.h
TrackFwd.h
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
hgcalPlots.stat
stat
Definition: hgcalPlots.py:1119
TrackMerger.h
std::swap
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
Definition: DataFrameContainer.h:209
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
reco::Track
Definition: Track.h:27
edm::ESHandle< GBRForest >
declareDynArray
#define declareDynArray(T, n, x)
Definition: DynArray.h:91
reco::TrackBase::charge
int charge() const
track electric charge
Definition: TrackBase.h:596
IfLogTrace
#define IfLogTrace(cond, cat)
Definition: MessageLogger.h:250
Point3DBase< float, GlobalTag >
ParameterSetDescription.h
b
double b
Definition: hdecay.h:118
mitigatedMETSequence_cff.U
U
Definition: mitigatedMETSequence_cff.py:36
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:223
edm::ParameterSet
Definition: ParameterSet.h:47
a
double a
Definition: hdecay.h:119
DetId::Tracker
Definition: DetId.h:25
TrackingRecHit::all
Definition: TrackingRecHit.h:59
Event.h
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
TrajectoryStateClosestToPoint::pt
double pt() const
Definition: TrajectoryStateClosestToPoint.h:78
ModuleDef.h
DuplicateTrackType
DuplicateTrackType
Definition: DuplicateTrackType.h:4
iEvent
int iEvent
Definition: GenABIO.cc:224
value
Definition: value.py:1
mag2
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
Definition: Basic3DVectorLD.h:124
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:49
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:148
edm::stream::EDProducer
Definition: EDProducer.h:38
edm::EventSetup
Definition: EventSetup.h:57
TrajectoryStateClosestToPoint
Definition: TrajectoryStateClosestToPoint.h:18
TSCPBuilderNoMaterial.h
edm::ESGetToken< GBRForest, GBRWrapperRcd >
TrackingRecHit::valid
Definition: TrackingRecHit.h:46
DynArray.h
Chi2MeasurementEstimatorBase
Definition: Chi2MeasurementEstimatorBase.h:14
TrackingRecHit
Definition: TrackingRecHit.h:21
BaseTrackerRecHit.h
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:64
trajectoryStateTransform::innerFreeState
FreeTrajectoryState innerFreeState(const reco::Track &tk, const MagneticField *field, bool withErr=true)
Definition: TrajectoryStateTransform.cc:86
submitPVResolutionJobs.desc
string desc
Definition: submitPVResolutionJobs.py:251
eostools.move
def move(src, dest)
Definition: eostools.py:511
TSCPBuilderNoMaterial
Definition: TSCPBuilderNoMaterial.h:17
FreeTrajectoryState
Definition: FreeTrajectoryState.h:27
reco::Track::seedRef
const edm::RefToBase< TrajectorySeed > & seedRef() const
Definition: Track.h:155
T
long double T
Definition: Basic3DVectorLD.h:48
DuplicateTrackType::Overlapping
TrackingComponentsRecord.h
ev
bool ev
Definition: Hydjet2Hadronizer.cc:95
reco::HitPattern::MISSING_INNER_HITS
Definition: HitPattern.h:155
EventSetup.h
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
TrajectoryStateTransform.h
DuplicateTrackType::Disjoint
TrajectoryStateClosestToPoint::isValid
bool isValid() const
Definition: TrajectoryStateClosestToPoint.h:111
DuplicateTrackMerger
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:224
ParameterSet.h
TrackCandidate.h
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
trajectoryStateTransform::innerStateOnSurface
TrajectoryStateOnSurface innerStateOnSurface(const reco::Track &tk, const TrackingGeometry &geom, const MagneticField *field, bool withErr=true)
Definition: TrajectoryStateTransform.cc:110
TrajectoryStateClosestToPoint::theState
const FreeTrajectoryState & theState() const
Definition: TrajectoryStateClosestToPoint.h:96
TrajectoryStateClosestToPoint::momentum
GlobalVector momentum() const
Definition: TrajectoryStateClosestToPoint.h:92
edm::Event
Definition: Event.h:73
HLT_FULL_cff.distance
distance
Definition: HLT_FULL_cff.py:7733
MagneticField
Definition: MagneticField.h:19
edm::InputTag
Definition: InputTag.h:15
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
TrajectoryStateOnSurface::isValid
bool isValid() const
Definition: TrajectoryStateOnSurface.h:54
TrackerGeometry
Definition: TrackerGeometry.h:14
Generated for CMSSW Reference Manual by doxygen 1.8.16