Belle II Software  release-05-02-19
MCRecoTracksMatcherModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Oliver Frost *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include <tracking/modules/mcMatcher/MCRecoTracksMatcherModule.h>
11 
12 // datastore types
13 #include <framework/datastore/StoreArray.h>
14 
15 #include <framework/gearbox/Const.h>
16 
17 #include <tracking/dataobjects/RecoTrack.h>
18 #include <mdst/dataobjects/MCParticle.h>
19 
20 // hit types
21 #include <pxd/dataobjects/PXDCluster.h>
22 #include <svd/dataobjects/SVDCluster.h>
23 #include <cdc/dataobjects/CDCHit.h>
24 
25 #include <map>
26 
27 #include <Eigen/Dense>
28 
29 /* disables the false positive warning of Intel icc about iter_pair_range::epmty() being unused
30 tracking/modules/mcMatcher/src/MCRecoTracksMatcherModule.cc:55: warning #177: function
31 "::iter_pair_range::empty [with Iter=std::_Rb_tree_const_iterator::DetId={Belle2::Const::EDetector},
32 ::HitId={int}>, ::WeightedRecoTrackId>>]" was declared but never referenced
33 bool empty() const
34 */
35 #ifdef __INTEL_COMPILER
36 #pragma warning disable 177
37 #endif
38 
39 using namespace Belle2;
40 
41 REG_MODULE(MCRecoTracksMatcher);
42 
43 namespace {
44  // small anonymous helper construct making converting a pair of iterators usable
45  // with range based for
46  template <class Iter>
47  struct iter_pair_range : std::pair<Iter, Iter> {
48  explicit iter_pair_range(std::pair<Iter, Iter> const& x)
49  : std::pair<Iter, Iter>(x)
50  {
51  }
52 
53  Iter begin() const
54  {
55  return this->first;
56  }
57 
58  Iter end() const
59  {
60  return this->second;
61  }
62 
63  bool empty() const
64  {
65  return begin() == end();
66  }
67  };
68 
69  template <class Iter>
70  inline iter_pair_range<Iter> as_range(std::pair<Iter, Iter> const& x)
71  {
72  return iter_pair_range<Iter>(x);
73  }
74 
75  using DetId = Const::EDetector;
76  using HitId = int;
77  using RecoTrackId = int;
78  struct WeightedRecoTrackId {
79  operator int() const { return id; }
80  RecoTrackId id;
81  double weight;
82  };
83 
84  using DetHitIdPair = std::pair<DetId, HitId>;
85 
86  struct CompDetHitIdPair {
87  bool operator()(const std::pair<DetId, HitId>& lhs,
88  const std::pair<std::pair<DetId, HitId>, WeightedRecoTrackId>& rhs)
89  {
90  return lhs < rhs.first;
91  }
92 
93  bool operator()(const std::pair<std::pair<DetId, HitId>, WeightedRecoTrackId>& lhs,
94  const std::pair<DetId, HitId>& rhs)
95  {
96  return lhs.first < rhs;
97  }
98  };
99 
100  // anonymous helper function to fill a set or a map with the hit IDs and det IDs (we need both a set or a map in the following).
101  template <class AMapOrSet>
102  void fillIDsFromStoreArray(AMapOrSet& recoTrackID_by_hitID,
103  const StoreArray<RecoTrack>& storedRecoTracks)
104  {
105  RecoTrackId recoTrackId = -1;
106  for (const RecoTrack& recoTrack : storedRecoTracks) {
107  ++recoTrackId;
108  std::vector<std::pair<DetHitIdPair, WeightedRecoTrackId> > hitIDsInTrack;
109  double totalWeight = 0;
110  using OriginTrackFinder = RecoHitInformation::OriginTrackFinder;
111  const OriginTrackFinder c_MCTrackFinderAuxiliaryHit =
112  OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
113 
114  for (const RecoHitInformation::UsedCDCHit* cdcHit : recoTrack.getCDCHitList()) {
115  OriginTrackFinder originFinder = recoTrack.getFoundByTrackFinder(cdcHit);
116  double weight = originFinder == c_MCTrackFinderAuxiliaryHit ? 0 : 1;
117  totalWeight += weight;
118  hitIDsInTrack.push_back({{Const::CDC, cdcHit->getArrayIndex()}, {recoTrackId, weight}});
119  }
120  for (const RecoHitInformation::UsedSVDHit* svdHit : recoTrack.getSVDHitList()) {
121  OriginTrackFinder originFinder = recoTrack.getFoundByTrackFinder(svdHit);
122  double weight = originFinder == c_MCTrackFinderAuxiliaryHit ? 0 : 1;
123  totalWeight += weight;
124  hitIDsInTrack.push_back({{Const::SVD, svdHit->getArrayIndex()}, {recoTrackId, weight}});
125  }
126  for (const RecoHitInformation::UsedPXDHit* pxdHit : recoTrack.getPXDHitList()) {
127  OriginTrackFinder originFinder = recoTrack.getFoundByTrackFinder(pxdHit);
128  double weight = originFinder == c_MCTrackFinderAuxiliaryHit ? 0 : 1;
129  totalWeight += weight;
130  hitIDsInTrack.push_back({{Const::PXD, pxdHit->getArrayIndex()}, {recoTrackId, weight}});
131  }
132 
133  // In case all hits are auxiliary for a track - reset all weights to 1
134  if (totalWeight == 0) {
135  for (std::pair<DetHitIdPair, WeightedRecoTrackId>& recoTrack_for_hitID : hitIDsInTrack) {
136  recoTrack_for_hitID.second.weight = 1;
137  }
138  }
139 
140  // Commit to output
141  typename AMapOrSet::iterator itInsertHint = recoTrackID_by_hitID.end();
142  for (std::pair<DetHitIdPair, WeightedRecoTrackId>& recoTrack_for_hitID : hitIDsInTrack) {
143  itInsertHint = recoTrackID_by_hitID.insert(itInsertHint, recoTrack_for_hitID);
144  }
145  }
146  }
147 }
148 
150  : Module()
151 {
152  setDescription("This module compares reconstructed tracks generated by some pattern recognition "
153  "algorithm for PXD, SVD and/or CDC to ideal Monte Carlo tracks and performs a "
154  "matching from the former to the underlying MCParticles.");
155 
157 
158  //Parameter definition
159  // Inputs
160  addParam("prRecoTracksStoreArrayName",
162  "Name of the collection containing the tracks as generate a patter recognition algorithm to be evaluated ",
163  std::string(""));
164 
165  addParam("mcRecoTracksStoreArrayName",
167  "Name of the collection containing the reference tracks as generate by a Monte-Carlo-Tracker (e.g. MCTrackFinder)",
168  std::string("MCGFTrackCands"));
169 
170  // Hit content to be evaluated
171  addParam("UsePXDHits",
173  "Set true if PXDHits or PXDClusters should be used in the matching in case they are present",
174  true);
175 
176  addParam("UseSVDHits",
178  "Set true if SVDHits or SVDClusters should be used in the matching in case they are present",
179  true);
180 
181  addParam("UseCDCHits",
183  "Set true if CDCHits should be used in the matching in case they are present",
184  true);
185 
186  addParam("UseOnlyAxialCDCHits",
188  "Set true if only the axial CDCHits should be used",
189  false);
190 
191  addParam("MinimalPurity",
193  "Minimal purity of a PRTrack to be considered matchable to a MCTrack. "
194  "This number encodes how many correct hits are minimally need to compensate for a false hits. "
195  "The default 2.0 / 3.0 suggests that for each background hit can be compensated by two correct hits.",
196  2.0 / 3.0);
197 
198  addParam("MinimalEfficiency",
200  "Minimal efficiency of a MCTrack to be considered matchable to a PRTrack. "
201  "This number encodes which fraction of the true hits must at least be in the reconstructed track. "
202  "The default 0.05 suggests that at least 5% of the true hits should have been picked up.",
203  0.05);
204 }
205 
207 {
208  // Check if there are MC Particles
209  StoreArray<MCParticle> storeMCParticles;
210 
211  if (storeMCParticles.isOptional()) {
212  m_mcParticlesPresent = true;
213 
214  // Actually retrieve the StoreArrays
217 
218  // Require both RecoTrack arrays and the MCParticles to be present in the DataStore
219  storePRRecoTracks.isRequired();
220  storeMCRecoTracks.isRequired();
221  storeMCParticles.isRequired();
222 
223  // Extract the default names for the case empty stings were given
224  m_param_prRecoTracksStoreArrayName = storePRRecoTracks.getName();
225  m_param_mcRecoTracksStoreArrayName = storeMCRecoTracks.getName();
226 
227  // Purity relation - for each PRTrack to store the purest MCTrack
228  storePRRecoTracks.registerRelationTo(storeMCRecoTracks);
229 
230  // Efficiency relation - for each MCTrack to store the most efficient PRTrack
231  storeMCRecoTracks.registerRelationTo(storePRRecoTracks);
232 
233  // MCParticle matching relation - purity
234  storePRRecoTracks.registerRelationTo(storeMCParticles);
235 
236  // MCParticle matching relation - efficiency
237  storeMCParticles.registerRelationTo(storePRRecoTracks);
238 
239  // Announce optional store arrays to the hits or clusters in case they should be used
240  // We make them optional in case of limited detector setup.
241  // PXD
242  if (m_param_usePXDHits) {
243  StoreArray<PXDCluster> pxdClusters;
244  pxdClusters.isOptional();
245  }
246 
247  // SVD
248  if (m_param_useSVDHits) {
249  StoreArray<SVDCluster> svdClusters;
250  svdClusters.isOptional();
251  }
252 
253  // CDC
254  if (m_param_useCDCHits) {
255  StoreArray<CDCHit> cdcHits;
256  cdcHits.isOptional();
257  }
258  }
259 }
260 
262 {
263  // Skip in the case there are no MC particles present.
264  if (not m_mcParticlesPresent) {
265  B2DEBUG(100, "Skipping MC Track Matcher as there are no MC Particles registered in the DataStore.");
266  return;
267  }
268 
269  B2DEBUG(100, "########## MCRecoTracksMatcherModule ############");
270 
271  // Fetch store array
274  StoreArray<MCParticle> mcParticles("");
275 
276  int nMCRecoTracks = mcRecoTracks.getEntries();
277  int nPRRecoTracks = prRecoTracks.getEntries();
278 
279  B2DEBUG(100, "Number patter recognition tracks is " << nPRRecoTracks);
280  B2DEBUG(100, "Number Monte-Carlo tracks is " << nMCRecoTracks);
281 
282  if (not nMCRecoTracks or not nPRRecoTracks) {
283  // Neither no pattern recognition tracks
284  // or no Monte Carlo tracks are present in this event
285  // Cannot perform matching.
286  return;
287  }
288 
289  // ### Build a detector_id hit_id to reco track map for easier lookup later ###
290  std::multimap<DetHitIdPair, WeightedRecoTrackId > mcId_by_hitId;
291  fillIDsFromStoreArray(mcId_by_hitId, mcRecoTracks);
292 
293  // Use set instead of multimap to handel to following situation
294  // * One hit may be assigned to multiple tracks should contribute to the efficiency of both tracks
295  // * One hit assigned twice or more to the same track should not contribute to the purity multiple times
296  // The first part is handled well by the multimap. But to enforce that one hit is also only assigned
297  // once to a track we use a set.
298  std::set<std::pair<DetHitIdPair, WeightedRecoTrackId>> prId_by_hitId;
299  fillIDsFromStoreArray(prId_by_hitId, prRecoTracks);
300 
301  // ### Get the number of relevant hits for each detector ###
302  // Since we are mostly dealing with indices in this module, this is all we need from the StoreArray
303  // Silently skip store arrays that are not present in reduced detector setups.
304  std::map<DetId, int> nHits_by_detId;
305 
306  // PXD
307  if (m_param_usePXDHits) {
308  StoreArray<PXDCluster> pxdClusters;
309  if (pxdClusters.isOptional()) {
310  int nHits = pxdClusters.getEntries();
311  nHits_by_detId[Const::PXD] = nHits;
312  }
313  }
314 
315  // SVD
316  if (m_param_useSVDHits) {
317  StoreArray<SVDCluster> svdClusters;
318  if (svdClusters.isOptional()) {
319  int nHits = svdClusters.getEntries();
320  nHits_by_detId[Const::SVD] = nHits;
321  }
322  }
323 
324  // CDC
325  if (m_param_useCDCHits) {
326  StoreArray<CDCHit> cdcHits;
327  if (cdcHits.isOptional()) {
328  nHits_by_detId[Const::CDC] = cdcHits.getEntries();
329  }
330  }
331 
332  //### Build the confusion matrix ###
333 
334  // Reserve enough space for the confusion matrix
335  // The last row is meant for hits not assigned to a mcRecoTrack (aka background hits)
336  Eigen::MatrixXd confusionMatrix = Eigen::MatrixXd::Zero(nPRRecoTracks, nMCRecoTracks + 1);
337  Eigen::MatrixXd weightedConfusionMatrix = Eigen::MatrixXd::Zero(nPRRecoTracks, nMCRecoTracks + 1);
338 
339  // Accumulated the total number of hits/ndf for each Monte-Carlo track separately to avoid double counting,
340  // in case patter recognition tracks share hits.
341  Eigen::RowVectorXd totalNDF_by_mcId = Eigen::RowVectorXd::Zero(nMCRecoTracks + 1);
342  Eigen::RowVectorXd totalWeight_by_mcId = Eigen::RowVectorXd::Zero(nMCRecoTracks + 1);
343 
344  // Accumulated the total number of hits/ndf for each pattern recognition track separately to avoid double counting,
345  // in case Monte-Carlo tracks share hits.
346  Eigen::VectorXd totalNDF_by_prId = Eigen::VectorXd::Zero(nPRRecoTracks);
347 
348  // Column index for the hits not assigned to any MCRecoTrack
349  const int mcBkgId = nMCRecoTracks;
350 
351  // for each detector examine every hit to which mcRecoTrack and prRecoTrack it belongs
352  // if the hit is not part of any mcRecoTrack put the hit in the background column.
353  for (const std::pair<const DetId, NDF>& detId_nHits_pair : nHits_by_detId) {
354 
355  DetId detId = detId_nHits_pair.first;
356  int nHits = detId_nHits_pair.second;
357  NDF ndfForOneHit = m_ndf_by_detId[detId];
358 
359  for (HitId hitId = 0; hitId < nHits; ++hitId) {
360  DetHitIdPair detId_hitId_pair(detId, hitId);
361 
362  if (m_param_useOnlyAxialCDCHits and detId == Const::CDC) {
363  StoreArray<CDCHit> cdcHits;
364  const CDCHit* cdcHit = cdcHits[hitId];
365  if (cdcHit->getISuperLayer() % 2 != 0) {
366  // Skip stereo hits
367  continue;
368  }
369  }
370 
371  // Seek all Monte Carlo tracks with the given hit
372  const auto mcIds_for_detId_hitId_pair =
373  as_range(mcId_by_hitId.equal_range(detId_hitId_pair));
374 
375  // Seek all pattern recognition tracks with the given hit
376  const auto prIds_for_detId_hitId_pair =
377  as_range(std::equal_range(prId_by_hitId.begin(),
378  prId_by_hitId.end(),
379  detId_hitId_pair,
380  CompDetHitIdPair()));
381 
382  // Assign the hits/ndf to the total ndf vector separately to avoid double counting,
383  // if patter recognition track share hits.
384  if (mcIds_for_detId_hitId_pair.empty()) {
385  // If the hit is not assigned to any mcRecoTrack
386  // The hit is assigned to the background column
387  RecoTrackId mcId = mcBkgId;
388  double mcWeight = 1;
389  totalNDF_by_mcId(mcId) += ndfForOneHit;
390  totalWeight_by_mcId(mcId) += ndfForOneHit * mcWeight;
391  } else {
392  for (const auto& detId_hitId_pair_and_mcId : mcIds_for_detId_hitId_pair) {
393  WeightedRecoTrackId mcId = detId_hitId_pair_and_mcId.second;
394  double mcWeight = mcId.weight;
395  totalNDF_by_mcId(mcId) += ndfForOneHit;
396  totalWeight_by_mcId(mcId) += ndfForOneHit * mcWeight;
397  }
398  }
399 
400  // Assign the hits/ndf to the total ndf vector separately here to avoid double counting,
401  // if Monte-Carlo track share hits.
402  for (const auto& detId_hitId_pair_and_prId : prIds_for_detId_hitId_pair) {
403  RecoTrackId prId = detId_hitId_pair_and_prId.second;
404  totalNDF_by_prId(prId) += ndfForOneHit;
405  }
406 
407  // Fill the confusion matrix
408  for (const auto& detId_hitId_pair_and_prId : prIds_for_detId_hitId_pair) {
409  RecoTrackId prId = detId_hitId_pair_and_prId.second;
410  if (mcIds_for_detId_hitId_pair.empty()) {
411  RecoTrackId mcId = mcBkgId;
412  double mcWeight = 1;
413  confusionMatrix(prId, mcId) += ndfForOneHit;
414  weightedConfusionMatrix(prId, mcId) += ndfForOneHit * mcWeight;
415  } else {
416  for (const auto& detId_hitId_pair_and_mcId : mcIds_for_detId_hitId_pair) {
417  WeightedRecoTrackId mcId = detId_hitId_pair_and_mcId.second;
418  double mcWeight = mcId.weight;
419  confusionMatrix(prId, mcId) += ndfForOneHit;
420  weightedConfusionMatrix(prId, mcId) += ndfForOneHit * mcWeight;
421  }
422  }
423  } // end for
424  } // end for hitId
425  } // end for detId
426 
427  B2DEBUG(200, "Confusion matrix of the event : " << std::endl << confusionMatrix);
428  B2DEBUG(200, "Weighted confusion matrix of the event : " << std::endl << weightedConfusionMatrix);
429 
430  B2DEBUG(200, "totalNDF_by_mcId : " << std::endl << totalNDF_by_mcId);
431  B2DEBUG(200, "totalWeight_by_mcId : " << std::endl << totalWeight_by_mcId);
432 
433  B2DEBUG(200, "totalNDF_by_prId : " << std::endl << totalNDF_by_prId);
434 
435  Eigen::MatrixXd purityMatrix = confusionMatrix.array().colwise() / totalNDF_by_prId.array();
436  Eigen::MatrixXd efficiencyMatrix = confusionMatrix.array().rowwise() / totalNDF_by_mcId.array();
437  Eigen::MatrixXd weightedEfficiencyMatrix = weightedConfusionMatrix.array().rowwise() / totalWeight_by_mcId.array();
438 
439  B2DEBUG(100, "Purities");
440  B2DEBUG(100, purityMatrix);
441 
442  B2DEBUG(100, "Efficiencies");
443  B2DEBUG(100, efficiencyMatrix);
444 
445  B2DEBUG(100, "Weighted efficiencies");
446  B2DEBUG(100, weightedEfficiencyMatrix);
447 
448  // ### Building the Monte-Carlo track to highest efficiency patter recognition track relation ###
449  // Weighted efficiency
450  using Efficiency = float;
451  using Purity = float;
452 
453  struct MostWeightEfficientPRId {
454  RecoTrackId id;
455  Efficiency weightedEfficiency;
456  Efficiency efficiency;
457  };
458  std::vector<MostWeightEfficientPRId> mostWeightEfficientPRId_by_mcId(nMCRecoTracks);
459  for (RecoTrackId mcId = 0; mcId < nMCRecoTracks; ++mcId) {
460  Eigen::VectorXd efficiencyCol = efficiencyMatrix.col(mcId);
461  Eigen::VectorXd weightedEfficiencyCol = weightedEfficiencyMatrix.col(mcId);
462 
463  RecoTrackId bestPrId = 0;
464  Efficiency bestWeightedEfficiency = weightedEfficiencyCol(0);
465  Efficiency bestEfficiency = efficiencyCol(0);
466  Purity bestPurity = purityMatrix.row(0)(mcId);
467 
468  // Reject efficiency smaller than the minimal one
469  if (bestWeightedEfficiency < m_param_minimalEfficiency) {
470  bestWeightedEfficiency = 0;
471  }
472 
473  // In case of a tie in the weighted efficiency we use the regular efficiency to break it.
474  for (RecoTrackId prId = 1; prId < nPRRecoTracks; ++prId) {
475  Eigen::RowVectorXd purityRow = purityMatrix.row(prId);
476 
477  Efficiency currentWeightedEfficiency = weightedEfficiencyCol(prId);
478  Efficiency currentEfficiency = efficiencyCol(prId);
479  Purity currentPurity = purityRow(mcId);
480 
481  // Reject efficiency smaller than the minimal one
482  if (currentWeightedEfficiency < m_param_minimalEfficiency) {
483  currentWeightedEfficiency = 0;
484  }
485 
486  if (std::tie(currentWeightedEfficiency, currentEfficiency, currentPurity) >
487  std::tie(bestWeightedEfficiency, bestEfficiency, bestPurity)) {
488  bestPrId = prId;
489  bestEfficiency = currentEfficiency;
490  bestWeightedEfficiency = currentWeightedEfficiency;
491  bestPurity = currentPurity;
492  }
493  }
494 
495  bestWeightedEfficiency = weightedEfficiencyCol(bestPrId);
496  bestEfficiency = efficiencyCol(bestPrId);
497  mostWeightEfficientPRId_by_mcId[mcId] = {bestPrId, bestWeightedEfficiency, bestEfficiency};
498  }
499 
500  // ### Building the patter recognition track to highest purity Monte-Carlo track relation ###
501  // Unweighted purity
502  struct MostPureMCId {
503  RecoTrackId id;
504  Purity purity;
505  };
506 
507  std::vector<MostPureMCId> mostPureMCId_by_prId(nPRRecoTracks);
508  for (int prId = 0; prId < nPRRecoTracks; ++prId) {
509  Eigen::RowVectorXd purityRow = purityMatrix.row(prId);
510 
511  int mcId;
512  Purity highestPurity = purityRow.maxCoeff(&mcId);
513 
514  mostPureMCId_by_prId[prId] = {mcId, highestPurity};
515  }
516 
517  // Log the Monte-Carlo track to highest efficiency patter recognition track relation
518  // Weighted efficiency
519  {
520  RecoTrackId mcId = -1;
521  B2DEBUG(200, "MCTrack to highest weighted efficiency PRTrack relation");
522  for (const auto& mostWeightEfficientPRId_for_mcId : mostWeightEfficientPRId_by_mcId) {
523  ++mcId;
524  const Efficiency& weightedEfficiency = mostWeightEfficientPRId_for_mcId.weightedEfficiency;
525  const RecoTrackId& prId = mostWeightEfficientPRId_for_mcId.id;
526  B2DEBUG(200,
527  "mcId : " << mcId << " -> prId : " << prId << " with weighted efficiency "
528  << weightedEfficiency);
529  }
530  }
531 
532  // Log the patter recognition track to highest purity Monte-Carlo track relation
533  // Unweighted purity
534  {
535  RecoTrackId prId = -1;
536  B2DEBUG(200, "PRTrack to highest purity MCTrack relation");
537  for (const auto& mostPureMCId_for_prId : mostPureMCId_by_prId) {
538  ++prId;
539  const RecoTrackId& mcId = mostPureMCId_for_prId.id;
540  const Purity& purity = mostPureMCId_for_prId.purity;
541  B2DEBUG(200, "prId : " << prId << " -> mcId : " << mcId << " with purity " << purity);
542  }
543  }
544 
545  // Count the categories
546  int nMatched{}, nBackground{}, nClones{}, nGhost{};
547 
548  // ### Classify the patter recognition tracks ###
549  // Means saving the highest purity relation to the data store
550  // + setup the PRTrack to MCParticle relation
551  // + save the set the MatchingStatus
552  for (RecoTrackId prId = 0; prId < nPRRecoTracks; ++prId) {
553  RecoTrack* prRecoTrack = prRecoTracks[prId];
554 
555  const MostPureMCId& mostPureMCId = mostPureMCId_by_prId[prId];
556 
557  const RecoTrackId& mcId = mostPureMCId.id;
558  const Purity& purity = mostPureMCId.purity;
559 
560  // GHOST
561  if (not(purity > 0) or not(purity >= m_param_minimalPurity)) {
562  prRecoTrack->setMatchingStatus(RecoTrack::MatchingStatus::c_ghost);
563 
564  B2DEBUG(100, "Stored PRTrack " << prId << " as ghost because of too low purity");
565  continue;
566  }
567 
568  // BACKGROUND
569  if (mcId == mcBkgId) {
570  prRecoTrack->setMatchingStatus(RecoTrack::MatchingStatus::c_background);
571 
572  B2DEBUG(100, "Stored PRTrack " << prId << " as background because of too low purity.");
573  ++nBackground;
574  continue;
575  }
576 
577  // After the classification for bad purity and background we examine,
578  // whether the highest purity Monte-Carlo track has in turn the highest efficiency
579  // patter recognition track equal to this track.
580  // All extra patter recognition tracks are marked as clones.
581 
582  RecoTrack* mcRecoTrack = mcRecoTracks[mcId];
583  MCParticle* mcParticle = mcRecoTrack->getRelated<MCParticle>();
584  B2ASSERT("No relation from MCRecoTrack to MCParticle.", mcParticle);
585 
586  const MostWeightEfficientPRId& mostWeightEfficientPRId_for_mcId =
587  mostWeightEfficientPRId_by_mcId[mcId];
588 
589  const RecoTrackId& mostWeightEfficientPRId = mostWeightEfficientPRId_for_mcId.id;
590  const Efficiency& weightedEfficiency = mostWeightEfficientPRId_for_mcId.weightedEfficiency;
591  // const Efficiency& efficiency = mostWeightEfficientPRId_for_mcId.efficiency;
592 
593  // Note : The matched category may also contain higher order clones recognisable by their low
594  // absolute efficiency
595 
596  // MATCHED
597  if (prId == mostWeightEfficientPRId) {
598  // Setup the relation purity relation
599  prRecoTrack->setMatchingStatus(RecoTrack::MatchingStatus::c_matched);
600  prRecoTrack->addRelationTo(mcRecoTrack, purity);
601 
602  // Add the mc matching relation
603  prRecoTrack->addRelationTo(mcParticle, purity);
604 
605  B2DEBUG(100, "Stored PRTrack " << prId << " as matched.");
606  B2DEBUG(100, "MC Match prId " << prId << " to mcPartId " << mcParticle->getArrayIndex());
607  B2DEBUG(100, "Purity rel: prId " << prId << " -> mcId " << mcId << " : " << purity);
608  ++nMatched;
609  continue;
610  }
611 
612  // GHOST
613  // Pattern recognition track fails the minimal efficiency requirement to be matched.
614  // We might want to introduce a different classification here, if we see problems
615  // with too many ghosts and want to investigate the specific source of the mismatch.
616  if (not(weightedEfficiency >= m_param_minimalEfficiency)) {
617  prRecoTrack->setMatchingStatus(RecoTrack::MatchingStatus::c_ghost);
618  B2DEBUG(100, "Stored PRTrack " << prId << " as ghost because of too low efficiency.");
619  ++nGhost;
620  continue;
621  }
622 
623  // Final category
624  // CLONE
625  // Setup the relation purity relation
626  prRecoTrack->setMatchingStatus(RecoTrack::MatchingStatus::c_clone);
627  prRecoTrack->addRelationTo(mcRecoTrack, -purity);
628  prRecoTrack->addRelationTo(mcParticle, -purity);
629  ++nClones;
630 
631  // Add the mc matching relation to the mc particle
632  B2DEBUG(100, "Stored PRTrack " << prId << " as clone.");
633  B2DEBUG(100, "MC Match prId " << prId << " to mcPartId " << mcParticle->getArrayIndex());
634  B2DEBUG(100, "Purity rel: prId " << prId << " -> mcId " << mcId << " : " << -purity);
635  } // end for prId
636 
637 
638  B2DEBUG(100, "Number of matches " << nMatched);
639  B2DEBUG(100, "Number of clones " << nClones);
640  B2DEBUG(100, "Number of bkg " << nBackground);
641  B2DEBUG(100, "Number of ghost " << nGhost);
642 
643  // ### Classify the Monte-Carlo tracks ###
644  // Meaning save the highest weighted efficiency relation to the data store.
645  for (RecoTrackId mcId = 0; mcId < nMCRecoTracks; ++mcId) {
646  RecoTrack* mcRecoTrack = mcRecoTracks[mcId];
647  MCParticle* mcParticle = mcRecoTrack->getRelated<MCParticle>();
648 
649  const MostWeightEfficientPRId& mostWeightEfficiencyPRId = mostWeightEfficientPRId_by_mcId[mcId];
650 
651  const RecoTrackId& prId = mostWeightEfficiencyPRId.id;
652  const Efficiency& weightedEfficiency = mostWeightEfficiencyPRId.weightedEfficiency;
653  // const Efficiency& efficiency = mostWeightEfficiencyPRId.efficiency;
654 
655  B2ASSERT("Index of pattern recognition tracks out of range.", prId < nPRRecoTracks and prId >= 0);
656 
657  RecoTrack* prRecoTrack = prRecoTracks[prId];
658 
659  const MostPureMCId& mostPureMCId_for_prId = mostPureMCId_by_prId[prId];
660  const RecoTrackId& mostPureMCId = mostPureMCId_for_prId.id;
661 
662  // MATCHED
663  if (mcId == mostPureMCId and
664  (prRecoTrack->getMatchingStatus() == RecoTrack::MatchingStatus::c_matched or
665  prRecoTrack->getMatchingStatus() == RecoTrack::MatchingStatus::c_clone)) {
666  // Setup the relation with positive weighted efficiency for this case.
667  mcRecoTrack->addRelationTo(prRecoTrack, weightedEfficiency);
668  mcParticle->addRelationTo(prRecoTrack, weightedEfficiency);
669  B2DEBUG(100, "Efficiency rel: mcId " << mcId << " -> prId " << prId << " : " << weightedEfficiency);
670  continue;
671  }
672 
673  // MERGED
674  // This MCTrack has a significant portion of hits in a PRTrack
675  // which in turn better describes a MCTrack different form this.
676  // Setup the relation with negative weighted efficiency for this case.
677  bool isMergedMCRecoTrack =
678  weightedEfficiency >= m_param_minimalEfficiency and
679  (prRecoTrack->getMatchingStatus() == RecoTrack::MatchingStatus::c_matched or
680  prRecoTrack->getMatchingStatus() == RecoTrack::MatchingStatus::c_clone);
681 
682  if (isMergedMCRecoTrack) {
683  mcRecoTrack->addRelationTo(prRecoTrack, -weightedEfficiency);
684  mcParticle->addRelationTo(prRecoTrack, -weightedEfficiency);
685  B2DEBUG(100, "Efficiency rel: mcId " << mcId << " -> prId " << prId << " : " << -weightedEfficiency);
686  continue;
687  }
688 
689  // MISSING
690  // No related pattern recognition track
691  // Do not create a relation.
692  B2DEBUG(100, "mcId " << mcId << " is missing. No relation created.");
693  B2DEBUG(100, "is Primary" << mcRecoTracks[mcId]->getRelatedTo<MCParticle>()->isPrimaryParticle());
694  B2DEBUG(100, "best prId " << prId << " with purity " << mostPureMCId_for_prId.purity << " -> " << mostPureMCId);
695  B2DEBUG(100, "MC Total ndf " << totalNDF_by_mcId[mcId]);
696  B2DEBUG(100, "MC Total weight" << totalWeight_by_mcId[mcId]);
697  B2DEBUG(100, "MC Overlap ndf\n " << confusionMatrix.col(mcId).transpose());
698  B2DEBUG(100, "MC Overlap weight\n " << weightedConfusionMatrix.col(mcId).transpose());
699  B2DEBUG(100, "MC Efficiencies for the track\n" << efficiencyMatrix.col(mcId).transpose());
700  B2DEBUG(100, "MC Weighted efficiencies for the track\n" << weightedEfficiencyMatrix.col(mcId).transpose());
701  } // end for mcId
702 
703  B2DEBUG(100, "########## End MCRecoTracksMatcherModule ############");
704 
705 } //end event()
Belle2::CDCHit::getISuperLayer
unsigned short getISuperLayer() const
Getter for iSuperLayer.
Definition: CDCHit.h:195
Belle2::MCRecoTracksMatcherModule::m_ndf_by_detId
std::map< int, NDF > m_ndf_by_detId
Map storing the standard number degrees of freedom for a single hit by detector *‍/.
Definition: MCRecoTracksMatcherModule.h:187
Belle2::MCRecoTracksMatcherModule::event
void event() final
Process the event.
Definition: MCRecoTracksMatcherModule.cc:261
Belle2::MCRecoTracksMatcherModule::m_param_minimalPurity
double m_param_minimalPurity
Parameter : Minimal purity of a PRTrack to be considered matchable to a MCTrack.
Definition: MCRecoTracksMatcherModule.h:172
Belle2::StoreArray::registerRelationTo
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:150
Belle2::Module::setDescription
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:216
Belle2::MCRecoTracksMatcherModule::m_mcParticlesPresent
bool m_mcParticlesPresent
Flag to indicated whether the Monte Carlo track are on the DataStore.
Definition: MCRecoTracksMatcherModule.h:190
Belle2::MCRecoTracksMatcherModule::initialize
void initialize() final
Signal the beginning of the event processing.
Definition: MCRecoTracksMatcherModule.cc:206
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::Module::c_ParallelProcessingCertified
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition: Module.h:82
Belle2::MCRecoTracksMatcherModule::m_param_useSVDHits
bool m_param_useSVDHits
Parameter : Switch whether SVDHits should be used in the matching.
Definition: MCRecoTracksMatcherModule.h:158
Belle2::CDCHit
Class containing the result of the unpacker in raw data and the result of the digitizer in simulation...
Definition: CDCHit.h:51
Belle2::MCRecoTracksMatcherModule::MCRecoTracksMatcherModule
MCRecoTracksMatcherModule()
Constructor setting up the parameter.
Definition: MCRecoTracksMatcherModule.cc:149
Belle2::RelationsInterface::addRelationTo
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
Definition: RelationsObject.h:144
Belle2::MCRecoTracksMatcherModule::m_param_usePXDHits
bool m_param_usePXDHits
Parameter : Switch whether PXDHits should be used in the matching.
Definition: MCRecoTracksMatcherModule.h:155
Belle2::RelationsInterface::getRelated
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
Definition: RelationsObject.h:280
Belle2::Const::EDetector
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition: Const.h:44
Belle2::RecoTrack::getMatchingStatus
MatchingStatus getMatchingStatus() const
Return the matching status set by the TrackMatcher module.
Definition: RecoTrack.h:718
Belle2::MCParticle::getArrayIndex
int getArrayIndex() const
Get 0-based index of the particle in the corresponding MCParticle list.
Definition: MCParticle.h:255
Belle2::MCRecoTracksMatcherModule::m_param_mcRecoTracksStoreArrayName
std::string m_param_mcRecoTracksStoreArrayName
Parameter : RecoTracks store array name from the mc recognition.
Definition: MCRecoTracksMatcherModule.h:152
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::RecoHitInformation::OriginTrackFinder
OriginTrackFinder
The TrackFinder which has added the hit to the track.
Definition: RecoHitInformation.h:80
Belle2::Module::setPropertyFlags
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition: Module.cc:210
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2::MCRecoTracksMatcherModule::m_param_minimalEfficiency
double m_param_minimalEfficiency
Parameter : Minimal efficiency for a MCTrack to be considered matchable to a PRTrack.
Definition: MCRecoTracksMatcherModule.h:180
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::RecoTrack::setMatchingStatus
void setMatchingStatus(MatchingStatus matchingStatus)
Set the matching status (used by the TrackMatcher module)
Definition: RecoTrack.h:724
Belle2::MCRecoTracksMatcherModule::NDF
int NDF
Descriptive type defintion for a number of degrees of freedom.
Definition: MCRecoTracksMatcherModule.h:184
Belle2::MCRecoTracksMatcherModule::m_param_useCDCHits
bool m_param_useCDCHits
Parameter : Switch whether CDCHits should be used in the matching.
Definition: MCRecoTracksMatcherModule.h:161
Belle2::PXDCluster
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition: PXDCluster.h:41
Belle2::MCRecoTracksMatcherModule::m_param_prRecoTracksStoreArrayName
std::string m_param_prRecoTracksStoreArrayName
Parameter : RecoTracks store array name from the patter recognition.
Definition: MCRecoTracksMatcherModule.h:149
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:38
Belle2::Module::addParam
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:562
Belle2::MCRecoTracksMatcherModule::m_param_useOnlyAxialCDCHits
bool m_param_useOnlyAxialCDCHits
Parameter : Switch whether only axial CDCHits should be used.
Definition: MCRecoTracksMatcherModule.h:164
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226