CMS 3D CMS Logo

TrackCollectionMerger.cc
Go to the documentation of this file.
2 
5 
12 
21 
24 
25 #include <vector>
26 #include <memory>
27 #include <cassert>
28 namespace {
29  class TrackCollectionMerger final : public edm::global::EDProducer<> {
30  public:
31  explicit TrackCollectionMerger(const edm::ParameterSet& conf)
32  : collectionCloner(producesCollector(), conf, true),
33  priorityName_(conf.getParameter<std::string>("trackAlgoPriorityOrder")),
34  m_priorityToken(esConsumes<TrackAlgoPriorityOrder, CkfComponentsRecord>(edm::ESInputTag("", priorityName_))),
35  m_foundHitBonus(conf.getParameter<double>("foundHitBonus")),
36  m_lostHitPenalty(conf.getParameter<double>("lostHitPenalty")),
37  m_shareFrac(conf.getParameter<double>("shareFrac")),
38  m_minShareHits(conf.getParameter<unsigned int>("minShareHits")),
39  m_minQuality(reco::TrackBase::qualityByName(conf.getParameter<std::string>("minQuality"))),
40  m_allowFirstHitShare(conf.getParameter<bool>("allowFirstHitShare")),
41  m_enableMerging(conf.getParameter<bool>("enableMerging")) {
42  for (auto const& it : conf.getParameter<std::vector<edm::InputTag>>("trackProducers"))
43  srcColls.emplace_back(it, consumesCollector());
44  for (auto const& it : conf.getParameter<std::vector<std::string>>("inputClassifiers")) {
45  srcMVAs.push_back(consumes<MVACollection>(edm::InputTag(it, "MVAValues")));
46  srcQuals.push_back(consumes<QualityMaskCollection>(edm::InputTag(it, "QualityMasks")));
47  }
48 
49  assert(srcColls.size() == srcQuals.size());
50 
51  produces<MVACollection>("MVAValues");
52  produces<QualityMaskCollection>("QualityMasks");
53  }
54 
55  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
57  desc.add<std::vector<edm::InputTag>>("trackProducers", std::vector<edm::InputTag>());
58  desc.add<std::vector<std::string>>("inputClassifiers", std::vector<std::string>());
59  desc.add<std::string>("trackAlgoPriorityOrder", "trackAlgoPriorityOrder");
60  desc.add<double>("shareFrac", .19);
61  desc.add<double>("foundHitBonus", 10.);
62  desc.add<double>("lostHitPenalty", 5.);
63  desc.add<unsigned int>("minShareHits", 2);
64  desc.add<bool>("allowFirstHitShare", true);
65  desc.add<bool>("enableMerging", true);
66  desc.add<std::string>("minQuality", "loose");
68  descriptions.add("TrackCollectionMerger", desc);
69  }
70 
71  private:
72  TrackCollectionCloner collectionCloner;
73  std::vector<TrackCollectionCloner::Tokens> srcColls;
74 
75  using MVACollection = std::vector<float>;
76  using QualityMaskCollection = std::vector<unsigned char>;
77 
78  using IHit = std::pair<unsigned int, TrackingRecHit const*>;
79  using IHitV = std::vector<IHit>;
80 
81  std::vector<edm::EDGetTokenT<MVACollection>> srcMVAs;
82  std::vector<edm::EDGetTokenT<QualityMaskCollection>> srcQuals;
83 
84  std::string priorityName_;
86 
87  float m_foundHitBonus;
88  float m_lostHitPenalty;
89  float m_shareFrac;
90  unsigned int m_minShareHits;
91  reco::TrackBase::TrackQuality m_minQuality;
92  bool m_allowFirstHitShare;
93  bool m_enableMerging;
94 
95  void produce(edm::StreamID, edm::Event& evt, const edm::EventSetup&) const override;
96 
97  bool areDuplicate(IHitV const& rh1, IHitV const& rh2) const;
98  };
99 
100 } // namespace
101 
103 
104 namespace {
105 
106  void TrackCollectionMerger::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& es) const {
107  TrackCollectionCloner::Producer producer(evt, collectionCloner);
108 
109  edm::ESHandle<TrackAlgoPriorityOrder> priorityH = es.getHandle(m_priorityToken);
110  auto const& trackAlgoPriorityOrder = *priorityH;
111 
112  // load collections
113  auto collsSize = srcColls.size();
114  auto rSize = 0U;
115  declareDynArray(reco::TrackCollection const*, collsSize, trackColls);
116  for (auto i = 0U; i < collsSize; ++i) {
117  trackColls[i] = &srcColls[i].tracks(evt);
118  rSize += (*trackColls[i]).size();
119  }
120 
121  unsigned char qualMask = ~0;
122  const bool acceptAll = m_minQuality == reco::TrackBase::undefQuality;
123  if (!acceptAll)
124  qualMask = 1 << m_minQuality;
125 
126  // load tracks
127  initDynArray(unsigned int, collsSize, nGoods, 0);
128  declareDynArray(float, rSize, mvas);
129  declareDynArray(unsigned char, rSize, quals);
130  declareDynArray(unsigned int, rSize, tkInds);
131  auto k = 0U;
132  for (auto i = 0U; i < collsSize; ++i) {
133  auto const& tkColl = *trackColls[i];
134  auto size = tkColl.size();
136  evt.getByToken(srcMVAs[i], hmva);
137  assert((*hmva).size() == size);
139  evt.getByToken(srcQuals[i], hqual);
140  for (auto j = 0U; j < size; ++j) {
141  if (!(acceptAll || (qualMask & (*hqual)[j])))
142  continue;
143  mvas[k] = (*hmva)[j];
144  quals[k] = (*hqual)[j];
145  tkInds[k] = j;
146  ++k;
147  ++nGoods[i];
148  }
149  }
150 
151  auto ntotTk = k;
152  std::vector<bool> selected(ntotTk, true);
153 
156  declareDynArray(reco::TrackBase::AlgoMask, ntotTk, algoMask);
157 
158  // merge (if more than one collection...)
159 
160  auto merger = [&]() -> void {
161  // load momentum, hits and score
163  declareDynArray(float, ntotTk, score);
164  declareDynArray(IHitV, ntotTk, rh1);
165 
166  k = 0U;
167  for (auto i = 0U; i < collsSize; ++i) {
168  auto const& tkColl = *trackColls[i];
169  for (auto j = 0U; j < nGoods[i]; ++j) {
170  auto const& track = tkColl[tkInds[k]];
171  algo[k] = track.algo();
172  oriAlgo[k] = track.originalAlgo();
173  algoMask[k] = track.algoMask();
174  mom[k] = track.isLooper() ? reco::TrackBase::Vector() : track.momentum();
175  auto validHits = track.numberOfValidHits();
176  auto validPixelHits = track.hitPattern().numberOfValidPixelHits();
177  auto lostHits = track.numberOfLostHits();
178  score[k] = m_foundHitBonus * validPixelHits + m_foundHitBonus * validHits - m_lostHitPenalty * lostHits -
179  track.chi2();
180 
181  auto& rhv = rh1[k];
182  rhv.reserve(validHits);
183  auto compById = [](IHit const& h1, IHit const& h2) { return h1.first < h2.first; };
184  for (auto it = track.recHitsBegin(); it != track.recHitsEnd(); ++it) {
185  auto const& hit = *(*it);
186  auto id = hit.rawId();
187  if LIKELY (hit.isValid()) {
188  rhv.emplace_back(id, &hit);
189  std::push_heap(rhv.begin(), rhv.end(), compById);
190  }
191  }
192  std::sort_heap(rhv.begin(), rhv.end(), compById);
193 
194  ++k;
195  }
196  }
197  assert(ntotTk == k);
198 
199  auto seti = [&](unsigned int ii, unsigned int jj) {
200  selected[jj] = false;
201  selected[ii] = true;
202  if (trackAlgoPriorityOrder.priority(oriAlgo[jj]) < trackAlgoPriorityOrder.priority(oriAlgo[ii]))
203  oriAlgo[ii] = oriAlgo[jj];
204  algoMask[ii] |= algoMask[jj];
205  quals[ii] |= (1 << reco::TrackBase::confirmed);
206  algoMask[jj] = algoMask[ii]; // in case we keep discarded
207  quals[jj] |= (1 << reco::TrackBase::discarded);
208  };
209 
210  auto iStart2 = 0U;
211  for (auto i = 0U; i < collsSize - 1; ++i) {
212  auto iStart1 = iStart2;
213  iStart2 = iStart1 + nGoods[i];
214  for (auto t1 = iStart1; t1 < iStart2; ++t1) {
215  if (!selected[t1])
216  continue;
217  auto score1 = score[t1];
218  for (auto t2 = iStart2; t2 < ntotTk; ++t2) {
219  if (!selected[t1])
220  break;
221  if (!selected[t2])
222  continue;
223  if (mom[t1].Dot(mom[t2]) < 0)
224  continue; // do not bother if in opposite hemespheres...
225  if (!areDuplicate(rh1[t1], rh1[t2]))
226  continue;
227  auto score2 = score[t2];
228 
229  // difference rather than ratio due to possible negative values for score
230  constexpr float almostSame = 0.01f;
231  if (score1 - score2 > almostSame) {
232  seti(t1, t2);
233  } else if (score2 - score1 > almostSame) {
234  seti(t2, t1);
235  } else { // take best
236  constexpr unsigned int qmask =
238  if ((quals[t1] & qmask) == (quals[t2] & qmask)) {
239  // take first
240  if (trackAlgoPriorityOrder.priority(algo[t1]) <= trackAlgoPriorityOrder.priority(algo[t2])) {
241  seti(t1, t2);
242  } else {
243  seti(t2, t1);
244  }
245  } else if ((quals[t1] & qmask) > (quals[t2] & qmask))
246  seti(t1, t2);
247  else
248  seti(t2, t1);
249  } // end ifs...
250 
251  } // end t2
252  } // end t1
253  } // end colls
254  }; // end merger;
255 
256  const bool doMerging = m_enableMerging && collsSize > 1;
257  if (doMerging)
258  merger();
259 
260  // products
261  auto pmvas = std::make_unique<MVACollection>();
262  auto pquals = std::make_unique<QualityMaskCollection>();
263 
264  // clone selected tracks...
265  auto nsel = 0U;
266  auto iStart2 = 0U;
267  auto isel = 0U;
268  for (auto i = 0U; i < collsSize; ++i) {
269  std::vector<unsigned int> selId;
270  std::vector<unsigned int> tid;
271  auto iStart1 = iStart2;
272  iStart2 = iStart1 + nGoods[i];
273  assert(producer.selTracks_->size() == isel);
274  for (auto t1 = iStart1; t1 < iStart2; ++t1) {
275  if (!selected[t1])
276  continue;
277  ++nsel;
278  tid.push_back(t1);
279  selId.push_back(tkInds[t1]);
280  pmvas->push_back(mvas[t1]);
281  pquals->push_back(quals[t1]);
282  }
283  producer(srcColls[i], selId);
284  assert(producer.selTracks_->size() == nsel);
285  assert(tid.size() == nsel - isel);
286  auto k = 0U;
287  for (; isel < nsel; ++isel) {
288  auto& otk = (*producer.selTracks_)[isel];
289  otk.setQualityMask((*pquals)[isel]); // needed also without merging
290  if (doMerging) {
291  otk.setOriginalAlgorithm(oriAlgo[tid[k]]);
292  otk.setAlgoMask(algoMask[tid[k++]]);
293  }
294  }
295  if (doMerging)
296  assert(tid.size() == k);
297  }
298 
299  assert(producer.selTracks_->size() == pmvas->size());
300 
301  // std::cout << "TrackCollectionMerger: sel tracks " << producer.selTracks_->size() << std::endl;
302 
303  evt.put(std::move(pmvas), "MVAValues");
304  evt.put(std::move(pquals), "QualityMasks");
305  }
306 
307  bool TrackCollectionMerger::areDuplicate(IHitV const& rh1, IHitV const& rh2) const {
308  auto nh1 = rh1.size();
309  auto nh2 = rh2.size();
310 
311  auto share = [](const TrackingRecHit* it, const TrackingRecHit* jt) -> bool {
312  return it->sharesInput(jt, TrackingRecHit::some);
313  };
314 
315  //loop over rechits
316  int noverlap = 0;
317  int firstoverlap = 0;
318  // check first hit (should use REAL first hit?)
319  if (m_allowFirstHitShare && rh1[0].first == rh2[0].first) {
320  if (share(rh1[0].second, rh2[0].second))
321  firstoverlap = 1;
322  }
323 
324  // exploit sorting
325  unsigned int jh = 0;
326  unsigned int ih = 0;
327  while (ih != nh1 && jh != nh2) {
328  auto const id1 = rh1[ih].first;
329  auto const id2 = rh2[jh].first;
330  if (id1 < id2)
331  ++ih;
332  else if (id2 < id1)
333  ++jh;
334  else {
335  if (share(rh1[ih].second, rh2[jh].second))
336  noverlap++;
337  ++jh;
338  ++ih;
339  } // equal ids
340  } //loop over ih & jh
341 
342  return noverlap >= int(m_minShareHits) &&
343  (noverlap - firstoverlap) > (std::min(nh1, nh2) - firstoverlap) * m_shareFrac;
344  }
345 
346 } // namespace
347 
350 
351 DEFINE_FWK_MODULE(TrackCollectionMerger);
ConfigurationDescriptions.h
trackAlgoPriorityOrder
constexpr std::array< unsigned int, reco::TrackBase::algoSize > trackAlgoPriorityOrder
Definition: trackAlgoPriorityOrder.h:157
edm::StreamID
Definition: StreamID.h:30
TrackExtra.h
RandomServiceHelper.t2
t2
Definition: RandomServiceHelper.py:257
electrons_cff.bool
bool
Definition: electrons_cff.py:366
mps_fire.i
i
Definition: mps_fire.py:428
HLT_FULL_cff.track
track
Definition: HLT_FULL_cff.py:11724
TrackCollectionCloner::Producer
Definition: TrackCollectionCloner.h:52
reco::TrackBase::tight
Definition: TrackBase.h:153
reco::TrackBase::loose
Definition: TrackBase.h:152
ESHandle.h
ESInputTag
min
T min(T a, T b)
Definition: MathUtil.h:58
reco::TrackBase::confirmed
Definition: TrackBase.h:155
edm
HLT enums.
Definition: AlignableModifier.h:19
reco::TrackBase::undefQuality
Definition: TrackBase.h:151
globals_cff.id1
id1
Definition: globals_cff.py:33
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
reco::TrackBase::TrackAlgorithm
TrackAlgorithm
track algorithm
Definition: TrackBase.h:89
TrackAlgoPriorityOrder.h
RawToDigi_Repacked_cff.producer
producer
Definition: RawToDigi_Repacked_cff.py:6
cms::cuda::assert
assert(be >=bs)
reco::TrackBase::TrackQuality
TrackQuality
track quality
Definition: TrackBase.h:150
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:222
edm::EDConsumerBase::consumesCollector
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
Definition: EDConsumerBase.cc:47
edm::global::EDProducerBase::produce
virtual void produce(StreamID, Event &, EventSetup const &) const =0
initDynArray
#define initDynArray(T, n, x, i)
Definition: DynArray.h:94
CkfComponentsRecord.h
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:46
edm::global::EDProducerBase::fillDescriptions
static void fillDescriptions(ConfigurationDescriptions &descriptions)
Definition: EDProducerBase.cc:294
edm::Handle
Definition: AssociativeIterator.h:50
ESGetToken.h
TrackCollectionCloner
Definition: TrackCollectionCloner.h:25
cmsdt::algo
algo
Definition: constants.h:171
MakerMacros.h
TrackingRecHit.h
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
Track.h
TrackFwd.h
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
CkfComponentsRecord
Definition: CkfComponentsRecord.h:22
TrackAlgoPriorityOrder
Definition: TrackAlgoPriorityOrder.h:6
edm::ESHandle
Definition: DTSurvey.h:22
declareDynArray
#define declareDynArray(T, n, x)
Definition: DynArray.h:91
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
dqmdumpme.k
k
Definition: dqmdumpme.py:60
ParameterSetDescription.h
TrackingRecHit::some
Definition: TrackingRecHit.h:59
TrajTrackAssociation.h
first
auto first
Definition: CAHitNtupletGeneratorKernelsImpl.h:125
mitigatedMETSequence_cff.U
U
Definition: mitigatedMETSequence_cff.py:36
edm::global::EDProducer
Definition: EDProducer.h:32
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
TrackCollectionCloner::fill
static void fill(edm::ParameterSetDescription &desc)
Definition: TrackCollectionCloner.cc:3
funct::true
true
Definition: Factorize.h:173
edm::ParameterSet
Definition: ParameterSet.h:47
reco::TrackBase::discarded
Definition: TrackBase.h:159
Event.h
ModuleDef.h
createfilelist.int
int
Definition: createfilelist.py:10
edm::Event::put
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:133
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:155
edm::EventSetup
Definition: EventSetup.h:58
edm::ESGetToken< TrackAlgoPriorityOrder, CkfComponentsRecord >
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
DynArray.h
TrackCollectionCloner.h
InputTag.h
TrackingRecHit
Definition: TrackingRecHit.h:21
Trajectory.h
submitPVResolutionJobs.desc
string desc
Definition: submitPVResolutionJobs.py:251
eostools.move
def move(src, dest)
Definition: eostools.py:511
std
Definition: JetResolutionObject.h:76
Frameworkfwd.h
electrons_cff.lostHits
lostHits
Definition: electrons_cff.py:391
LIKELY
#define LIKELY(x)
Definition: Likely.h:20
findQualityFiles.jj
string jj
Definition: findQualityFiles.py:188
reco::TrackBase::AlgoMask
std::bitset< algoSize > AlgoMask
algo mask
Definition: TrackBase.h:145
EventSetup.h
TrackingRecHit::sharesInput
virtual bool sharesInput(const TrackingRecHit *other, SharedInputType what) const
Definition: TrackingRecHit.cc:12
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
pixelTrack::qualityByName
Quality qualityByName(std::string const &name)
Definition: TrackSoAHeterogeneousT.h:16
ParameterSet.h
offlineSlimmedPrimaryVertices_cfi.score
score
Definition: offlineSlimmedPrimaryVertices_cfi.py:6
EDProducer.h
globals_cff.id2
id2
Definition: globals_cff.py:34
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
DeDxTools::esConsumes
ESGetTokenH3DDVariant esConsumes(std::string const &Reccord, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
edm::Event
Definition: Event.h:73
cuy.ii
ii
Definition: cuy.py:589
edm::InputTag
Definition: InputTag.h:15
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
hit
Definition: SiStripHitEffFromCalibTree.cc:88
reco::TrackBase::Vector
math::XYZVector Vector
spatial vector
Definition: TrackBase.h:77
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
reco::TrackBase::highPurity
Definition: TrackBase.h:154