LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArHierarchyHelper.cc
Go to the documentation of this file.
1 
9 #include "Pandora/PdgTable.h"
10 #include "Pandora/StatusCodes.h"
11 
14 
15 #include <numeric>
16 
17 namespace lar_content
18 {
19 
20 using namespace pandora;
21 
23  m_foldToLeadingShowers{false},
24  m_foldToTier{false},
25  m_foldDynamic{false},
26  m_cosAngleTolerance{0.9962f},
27  m_tier{1}
28 {
29 }
30 
31 //------------------------------------------------------------------------------------------------------------------------------------------
32 
33 LArHierarchyHelper::FoldingParameters::FoldingParameters(const bool foldDynamic, const float cosAngleTolerance) :
35  m_foldToTier{false},
36  m_foldDynamic{foldDynamic},
37  m_cosAngleTolerance{cosAngleTolerance},
38  m_tier{1}
39 {
40 }
41 
42 //------------------------------------------------------------------------------------------------------------------------------------------
43 
46  m_foldToTier{true},
47  m_foldDynamic{false},
48  m_cosAngleTolerance{0.9962f},
49  m_tier{foldingTier}
50 {
51  if (m_tier < 1)
52  {
53  std::cout << "LArHierarchyHelper: Error - attempting to fold to non-positive tier" << std::endl;
54  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
55  }
56 }
57 
58 //------------------------------------------------------------------------------------------------------------------------------------------
59 //------------------------------------------------------------------------------------------------------------------------------------------
60 
62  m_minPurity{0.8f},
63  m_minCompleteness{0.65f}
64 {
65 }
66 
67 //------------------------------------------------------------------------------------------------------------------------------------------
68 
69 LArHierarchyHelper::QualityCuts::QualityCuts(const float minPurity, const float minCompleteness) :
70  m_minPurity{minPurity},
71  m_minCompleteness{minCompleteness}
72 {
73 }
74 
75 //------------------------------------------------------------------------------------------------------------------------------------------
76 //------------------------------------------------------------------------------------------------------------------------------------------
77 
79  m_nextNodeId{1}
80 {
81 }
82 
83 //------------------------------------------------------------------------------------------------------------------------------------------
84 
86  m_recoCriteria(recoCriteria),
87  m_nextNodeId{1}
88 {
89 }
90 
91 //------------------------------------------------------------------------------------------------------------------------------------------
92 
94 {
95  for (const auto &[pRoot, nodeVector] : m_interactions)
96  {
97  (void)pRoot;
98  for (const Node *pNode : nodeVector)
99  delete pNode;
100  }
101  m_interactions.clear();
102 }
103 
104 //------------------------------------------------------------------------------------------------------------------------------------------
105 
106 void LArHierarchyHelper::MCHierarchy::FillHierarchy(const MCParticleList &mcParticleList, const CaloHitList &caloHitList, const FoldingParameters &foldParameters)
107 {
108  const auto predicate = [](const MCParticle *pMCParticle) { return std::abs(pMCParticle->GetParticleId()) == NEUTRON; };
109  m_mcToHitsMap.clear();
110  for (const CaloHit *pCaloHit : caloHitList)
111  {
112  try
113  {
114  const MCParticle *pMCParticle{MCParticleHelper::GetMainMCParticle(pCaloHit)};
115  m_mcToHitsMap[pMCParticle].emplace_back(pCaloHit);
116  }
117  catch (const StatusCodeException &)
118  {
119  }
120  }
121 
122  MCParticleList rootNodes;
123  for (const MCParticle *pMCParticle : mcParticleList)
124  {
125  const MCParticleList &parentList{pMCParticle->GetParentList()};
126  if (parentList.empty())
127  {
128  rootNodes.emplace_back(pMCParticle);
129  }
130  }
131 
132  for (const MCParticle *pRoot : rootNodes)
133  {
134  MCParticleSet primarySet;
135  LArHierarchyHelper::GetMCPrimaries(pRoot, primarySet);
136  MCParticleList primaries(primarySet.begin(), primarySet.end());
137  primaries.sort(LArMCParticleHelper::SortByMomentum);
139  primaries.erase(std::remove_if(primaries.begin(), primaries.end(), predicate), primaries.end());
140  if (foldParameters.m_foldToTier && foldParameters.m_tier == 1)
141  {
142  for (const MCParticle *pPrimary : primaries)
143  {
144  MCParticleList allParticles{pPrimary};
146  {
147  LArMCParticleHelper::GetAllDescendentMCParticles(pPrimary, allParticles);
148  }
149  else
150  {
151  // Collect track-like and shower-like particles together, but throw out neutrons and descendents
152  MCParticleList dummy;
153  LArMCParticleHelper::GetAllDescendentMCParticles(pPrimary, allParticles, allParticles, dummy);
154  }
155  CaloHitList allHits;
156  for (const MCParticle *pMCParticle : allParticles)
157  {
158  // Not all MC particles will have hits
159  if (m_mcToHitsMap.find(pMCParticle) != m_mcToHitsMap.end())
160  {
161  const CaloHitList &caloHits(m_mcToHitsMap.at(pMCParticle));
162  allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
163  }
164  }
165  m_interactions[pRoot].emplace_back(new Node(*this, allParticles, allHits));
166  }
167  }
168  else if (foldParameters.m_foldToLeadingShowers)
169  {
170  for (const MCParticle *pPrimary : primaries)
171  {
172  MCParticleList allParticles{pPrimary};
173  int pdg{std::abs(pPrimary->GetParticleId())};
174  const bool isShower{pdg == E_MINUS || pdg == PHOTON};
175  const bool isNeutron{pdg == NEUTRON};
176  if (isShower || (isNeutron && !m_recoCriteria.m_removeNeutrons))
177  LArMCParticleHelper::GetAllDescendentMCParticles(pPrimary, allParticles);
178  CaloHitList allHits;
179  for (const MCParticle *pMCParticle : allParticles)
180  {
181  // ATTN - Not all MC particles will have hits
182  if (m_mcToHitsMap.find(pMCParticle) != m_mcToHitsMap.end())
183  {
184  const CaloHitList &caloHits(m_mcToHitsMap.at(pMCParticle));
185  allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
186  }
187  }
188  Node *pNode{new Node(*this, allParticles, allHits)};
189  m_interactions[pRoot].emplace_back(pNode);
190  if (!(isShower || isNeutron))
191  {
192  // Find the children of this particle and recursively add them to the hierarchy
193  const MCParticleList &children{pPrimary->GetDaughterList()};
194  for (const MCParticle *pChild : children)
195  pNode->FillHierarchy(pChild, foldParameters);
196  }
197  }
198  }
199  else if (foldParameters.m_foldDynamic)
200  {
201  for (const MCParticle *pPrimary : primaries)
202  {
203  MCParticleList leadingParticles, childParticles;
204  this->InterpretHierarchy(pPrimary, leadingParticles, childParticles, foldParameters.m_cosAngleTolerance);
205  CaloHitList allHits;
206  for (const MCParticle *pMCParticle : leadingParticles)
207  {
208  // ATTN - Not all MC particles will have hits
209  if (m_mcToHitsMap.find(pMCParticle) != m_mcToHitsMap.end())
210  {
211  const CaloHitList &caloHits(m_mcToHitsMap.at(pMCParticle));
212  allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
213  }
214  }
215 
216  Node *pNode{new Node(*this, leadingParticles, allHits)};
217  m_interactions[pRoot].emplace_back(pNode);
218  for (const MCParticle *pChild : childParticles)
219  pNode->FillHierarchy(pChild, foldParameters);
220  }
221  }
222  else
223  {
224  // Unfolded and folded to tier > 1 have the same behaviour for primaries
225  for (const MCParticle *pPrimary : primaries)
226  {
227  MCParticleList allParticles{pPrimary};
228  CaloHitList allHits;
229  for (const MCParticle *pMCParticle : allParticles)
230  {
231  // ATTN - Not all MC particles will have hits
232  if (m_mcToHitsMap.find(pMCParticle) != m_mcToHitsMap.end())
233  {
234  const CaloHitList &caloHits(m_mcToHitsMap.at(pMCParticle));
235  allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
236  }
237  }
238  Node *pNode{new Node(*this, allParticles, allHits)};
239  m_interactions[pRoot].emplace_back(pNode);
240  // Find the children of this particle and recursively add them to the hierarchy
241  const MCParticleList &children{pPrimary->GetDaughterList()};
242  for (const MCParticle *pChild : children)
243  pNode->FillHierarchy(pChild, foldParameters);
244  }
245  }
246 
247  Node *pLeadingLepton{nullptr};
248  float leadingLeptonEnergy{-std::numeric_limits<float>::max()};
249  for (const Node *pNode : m_interactions[pRoot])
250  {
251  const MCParticle *pMC{pNode->GetLeadingMCParticle()};
252  if (pMC)
253  {
254  const int pdg{std::abs(pMC->GetParticleId())};
255  if ((pdg == MU_MINUS || pdg == E_MINUS || pdg == TAU_MINUS) && pMC->GetEnergy() > leadingLeptonEnergy)
256  {
257  pLeadingLepton = const_cast<Node *>(pNode);
258  leadingLeptonEnergy = pMC->GetEnergy();
259  }
260  }
261  }
262  if (pLeadingLepton)
263  pLeadingLepton->SetLeadingLepton();
264  }
265 }
266 
267 //------------------------------------------------------------------------------------------------------------------------------------------
268 
270 {
271  if (m_interactions.find(pRoot) == m_interactions.end())
272  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
273 
274  return m_interactions.at(pRoot);
275 }
276 
277 //------------------------------------------------------------------------------------------------------------------------------------------
278 
279 void LArHierarchyHelper::MCHierarchy::GetRootMCParticles(MCParticleList &rootMCParticles) const
280 {
281  for (auto iter = m_interactions.begin(); iter != m_interactions.end(); ++iter)
282  rootMCParticles.emplace_back(iter->first);
283 }
284 
285 //------------------------------------------------------------------------------------------------------------------------------------------
286 
288  const MCParticle *const pRoot, MCParticleList &leadingParticles, MCParticleList &childParticles, const float cosAngleTolerance) const
289 {
290  leadingParticles.emplace_back(pRoot);
291  MCParticleList foldCandidates, childCandidates;
292  const MCParticleList &children{pRoot->GetDaughterList()};
293  for (const MCParticle *pMCParticle : children)
294  {
295  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
296  if (!pLArMCParticle)
297  continue;
299  {
300  // Elastic and inelastic scattering can either lead to folding, distinct nodes or disposable hits, all other processes
301  // are either distinct nodes, or disposable
302  if (pMCParticle->GetParticleId() == pRoot->GetParticleId())
303  foldCandidates.emplace_back(pMCParticle);
304  else
305  childCandidates.emplace_back(pMCParticle);
306  }
307  else if (!m_recoCriteria.m_removeNeutrons || (m_recoCriteria.m_removeNeutrons && pLArMCParticle->GetProcess() != MC_PROC_N_CAPTURE))
308  {
309  // Non-scattering process particles become leading candidates unless it's neutron capture and we're removing neutrons
310  childCandidates.emplace_back(pMCParticle);
311  }
312  }
313  const MCParticle *pBestFoldCandidate{nullptr};
314  float bestDp{std::numeric_limits<float>::max()};
315  for (const MCParticle *pMCParticle : foldCandidates)
316  {
317  if (foldCandidates.size() == 1)
318  {
319  // No alternative options, so this is either the best folding option by default, or a sufficiently large scatter to
320  // treat as a new particle for reconstruction purposes
321  if (LArMCParticleHelper::AreTopologicallyContinuous(pRoot, pMCParticle, cosAngleTolerance))
322  pBestFoldCandidate = pMCParticle;
323  else
324  childCandidates.emplace_back(pMCParticle);
325  }
326  else
327  {
328  // Assess which, if any, of the children might be a continuation of the trajectory, otherwise move to child candidates
329  if (LArMCParticleHelper::AreTopologicallyContinuous(pRoot, pMCParticle, cosAngleTolerance))
330  {
331  const float dp{pRoot->GetMomentum().GetMagnitude() - pMCParticle->GetMomentum().GetMagnitude()};
332  if (dp < bestDp)
333  {
334  pBestFoldCandidate = pMCParticle;
335  bestDp = dp;
336  }
337  }
338  else
339  {
340  childCandidates.emplace_back(pMCParticle);
341  }
342  }
343  }
344  if (pBestFoldCandidate)
345  {
346  leadingParticles.emplace_back(pBestFoldCandidate);
347  // Having found a particle to fold back at this level, continue to explore its downstream hierarchy for further folding
348  // opportunities and make their respective children leading particles for the folded node we are creating
349  this->CollectContinuations(pBestFoldCandidate, leadingParticles, childCandidates, cosAngleTolerance);
350  }
351  for (const MCParticle *pMCParticle : childCandidates)
352  {
353  // Consider if the child particle will produce enough downstream hits to warrant inclusion
354  if (this->IsReconstructable(pMCParticle))
355  childParticles.emplace_back(pMCParticle);
356  else
357  {
358  MCParticleList localHierarchy{pMCParticle};
359  CaloHitList localHits;
360  LArMCParticleHelper::GetAllDescendentMCParticles(pMCParticle, localHierarchy);
361  for (const MCParticle *pLocalMCParticle : localHierarchy)
362  {
363  if (m_mcToHitsMap.find(pLocalMCParticle) != m_mcToHitsMap.end())
364  {
365  const CaloHitList &caloHits(m_mcToHitsMap.at(pLocalMCParticle));
366  localHits.insert(localHits.begin(), caloHits.begin(), caloHits.end());
367  }
368  }
369  if (this->IsReconstructable(localHits))
370  childParticles.emplace_back(pMCParticle);
371  }
372  }
373 }
374 
375 //------------------------------------------------------------------------------------------------------------------------------------------
376 
378  const MCParticle *pRoot, MCParticleList &continuingParticles, MCParticleList &childParticles, const float cosAngleTolerance) const
379 {
380  const MCParticleList &children{pRoot->GetDaughterList()};
381  MCParticleList foldCandidates;
382  for (const MCParticle *pMCParticle : children)
383  {
384  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
385  if (!pLArMCParticle)
386  continue;
387  // Only elastic and inelastic scattering can lead to folding
389  {
390  if (pMCParticle->GetParticleId() == pRoot->GetParticleId())
391  foldCandidates.emplace_back(pMCParticle);
392  }
393  else if (!m_recoCriteria.m_removeNeutrons || (m_recoCriteria.m_removeNeutrons && pLArMCParticle->GetProcess() != MC_PROC_N_CAPTURE))
394  {
395  // Non-scattering process particles become leading candidates unless it's neutron capture and we're removing neutrons
396  childParticles.emplace_back(pMCParticle);
397  }
398  }
399  const MCParticle *pBestFoldCandidate{nullptr};
400  float bestDp{std::numeric_limits<float>::max()};
401  for (const MCParticle *pMCParticle : foldCandidates)
402  {
403  if (foldCandidates.size() == 1)
404  {
405  // No alternative options, so this is either the best folding option by default, or a sufficiently large scatter to
406  // treat as a new particle for reconstruction purposes
407  if (LArMCParticleHelper::AreTopologicallyContinuous(pRoot, pMCParticle, cosAngleTolerance))
408  pBestFoldCandidate = pMCParticle;
409  }
410  else
411  {
412  // Assess which, if any, of the children might be a continuation of the trajectory, otherwise move to child candidates
413  if (LArMCParticleHelper::AreTopologicallyContinuous(pRoot, pMCParticle, cosAngleTolerance))
414  {
415  const float dp{pRoot->GetMomentum().GetMagnitude() - pMCParticle->GetMomentum().GetMagnitude()};
416  if (dp < bestDp)
417  {
418  pBestFoldCandidate = pMCParticle;
419  bestDp = dp;
420  }
421  }
422  }
423  }
424  if (pBestFoldCandidate)
425  {
426  continuingParticles.emplace_back(pBestFoldCandidate);
427  const MCParticleList &newLeadingParticles{pBestFoldCandidate->GetDaughterList()};
428  // We need to add the children as child particles to ensure these sub-hierarchies are explored...
429  childParticles.insert(childParticles.begin(), newLeadingParticles.begin(), newLeadingParticles.end());
430  // but this current best fold candidate may have been added to the child particles by previously, so remove it
431  const auto iter{std::find(childParticles.begin(), childParticles.end(), pBestFoldCandidate)};
432  if (iter != childParticles.end())
433  childParticles.erase(iter);
434  // Having found a particle to fold back at this level, continue to explore its downstream hierarchy for further folding
435  // opportunities and make their respective children child particles for the folded node we are creating
436  LArHierarchyHelper::MCHierarchy::CollectContinuations(pBestFoldCandidate, continuingParticles, childParticles, cosAngleTolerance);
437  }
438 }
439 
440 //------------------------------------------------------------------------------------------------------------------------------------------
441 
442 void LArHierarchyHelper::MCHierarchy::GetFlattenedNodes(const MCParticle *const pRoot, NodeVector &nodeVector) const
443 {
444  NodeList queue;
445  for (const Node *pNode : m_interactions.at(pRoot))
446  {
447  nodeVector.emplace_back(pNode);
448  queue.emplace_back(pNode);
449  }
450  while (!queue.empty())
451  {
452  const NodeVector &children{queue.front()->GetChildren()};
453  queue.pop_front();
454  for (const Node *pChild : children)
455  {
456  nodeVector.emplace_back(pChild);
457  queue.emplace_back(pChild);
458  }
459  }
460 }
461 
462 //------------------------------------------------------------------------------------------------------------------------------------------
463 
465 {
466  m_nodeToIdMap.insert(std::make_pair(pNode, m_nextNodeId));
467  ++m_nextNodeId;
468 }
469 
470 //------------------------------------------------------------------------------------------------------------------------------------------
471 
473 {
474  std::string str;
475  for (const auto &[pRoot, nodeVector] : m_interactions)
476  {
477  const LArMCParticle *const pLArRoot{dynamic_cast<const LArMCParticle *const>(pRoot)};
478  if (pLArRoot)
479  str += "=== MC Interaction : PDG " + std::to_string(pLArRoot->GetParticleId()) +
480  " Energy: " + std::to_string(pLArRoot->GetEnergy()) + " Nuance: " + std::to_string(pLArRoot->GetNuanceCode()) + "\n";
481  else
482  str += "=== MC Interaction : PDG " + std::to_string(pRoot->GetParticleId()) + " Energy: " + std::to_string(pRoot->GetEnergy()) + "\n";
483  for (const Node *pNode : nodeVector)
484  str += " " + pNode->ToString("") + "\n";
485  }
486 
487  return str;
488 }
489 
490 //------------------------------------------------------------------------------------------------------------------------------------------
491 
492 bool LArHierarchyHelper::MCHierarchy::IsReconstructable(const pandora::MCParticle *pMCParticle) const
493 {
494  if (m_mcToHitsMap.find(pMCParticle) != m_mcToHitsMap.end())
495  {
496  unsigned int nHitsU{0}, nHitsV{0}, nHitsW{0};
497  for (const CaloHit *pCaloHit : m_mcToHitsMap.at(pMCParticle))
498  {
499  const HitType view{pCaloHit->GetHitType()};
500  if (view == TPC_VIEW_U)
501  ++nHitsU;
502  else if (view == TPC_VIEW_V)
503  ++nHitsV;
504  else if (view == TPC_VIEW_W)
505  ++nHitsW;
506  }
507  const unsigned int nHits{nHitsU + nHitsV + nHitsW};
508  unsigned int nGoodViews{0};
509  nGoodViews += nHitsU >= m_recoCriteria.m_minHitsForGoodView ? 1 : 0;
510  nGoodViews += nHitsV >= m_recoCriteria.m_minHitsForGoodView ? 1 : 0;
511  nGoodViews += nHitsW >= m_recoCriteria.m_minHitsForGoodView ? 1 : 0;
512 
513  return nHits >= m_recoCriteria.m_minHits && nGoodViews >= m_recoCriteria.m_minGoodViews;
514  }
515 
516  return false;
517 }
518 
519 //------------------------------------------------------------------------------------------------------------------------------------------
520 
521 bool LArHierarchyHelper::MCHierarchy::IsReconstructable(const CaloHitList &caloHits) const
522 {
523  unsigned int nHitsU{0}, nHitsV{0}, nHitsW{0};
524  for (const CaloHit *pCaloHit : caloHits)
525  {
526  const HitType view{pCaloHit->GetHitType()};
527  if (view == TPC_VIEW_U)
528  ++nHitsU;
529  else if (view == TPC_VIEW_V)
530  ++nHitsV;
531  else if (view == TPC_VIEW_W)
532  ++nHitsW;
533  }
534  const unsigned int nHits{nHitsU + nHitsV + nHitsW};
535  unsigned int nGoodViews{0};
536  nGoodViews += nHitsU >= m_recoCriteria.m_minHitsForGoodView ? 1 : 0;
537  nGoodViews += nHitsV >= m_recoCriteria.m_minHitsForGoodView ? 1 : 0;
538  nGoodViews += nHitsW >= m_recoCriteria.m_minHitsForGoodView ? 1 : 0;
539 
540  return nHits >= m_recoCriteria.m_minHits && nGoodViews >= m_recoCriteria.m_minGoodViews;
541 }
542 
543 //------------------------------------------------------------------------------------------------------------------------------------------
544 //------------------------------------------------------------------------------------------------------------------------------------------
545 
546 LArHierarchyHelper::MCHierarchy::Node::Node(MCHierarchy &hierarchy, const MCParticle *pMCParticle, const int tier) :
547  m_hierarchy(hierarchy),
548  m_mainParticle(pMCParticle),
549  m_tier{tier},
550  m_pdg{0},
551  m_isLeadingLepton{false}
552 {
553  if (pMCParticle)
554  {
555  m_pdg = pMCParticle->GetParticleId();
556  m_mcParticles.emplace_back(pMCParticle);
557  }
558  m_hierarchy.RegisterNode(this);
559 }
560 
561 //------------------------------------------------------------------------------------------------------------------------------------------
562 
563 LArHierarchyHelper::MCHierarchy::Node::Node(MCHierarchy &hierarchy, const MCParticleList &mcParticleList, const CaloHitList &caloHitList, const int tier) :
564  m_hierarchy(hierarchy),
565  m_mcParticles(mcParticleList),
566  m_caloHits(caloHitList),
567  m_mainParticle(nullptr),
568  m_tier{tier},
569  m_pdg{0},
570  m_isLeadingLepton{false}
571 {
572  if (!mcParticleList.empty())
573  {
574  m_mainParticle = mcParticleList.front();
575  m_pdg = m_mainParticle->GetParticleId();
576  }
577  m_mcParticles.sort(LArMCParticleHelper::SortByMomentum);
578  m_caloHits.sort();
579  m_hierarchy.RegisterNode(this);
580 }
581 
582 //------------------------------------------------------------------------------------------------------------------------------------------
583 
585 {
586  m_mcParticles.clear();
587  m_caloHits.clear();
588  for (const Node *node : m_children)
589  delete node;
590  m_children.clear();
591 }
592 
593 //------------------------------------------------------------------------------------------------------------------------------------------
594 
595 void LArHierarchyHelper::MCHierarchy::Node::FillHierarchy(const MCParticle *pRoot, const FoldingParameters &foldParameters)
596 {
597  if (foldParameters.m_foldDynamic)
598  {
599  MCParticleList leadingParticles, childParticles;
600  m_hierarchy.InterpretHierarchy(pRoot, leadingParticles, childParticles, foldParameters.m_cosAngleTolerance);
601  CaloHitList allHits;
602  for (const MCParticle *pMCParticle : leadingParticles)
603  {
604  // ATTN - Not all MC particles will have hits
605  if (m_hierarchy.m_mcToHitsMap.find(pMCParticle) != m_hierarchy.m_mcToHitsMap.end())
606  {
607  const CaloHitList &caloHits(m_hierarchy.m_mcToHitsMap.at(pMCParticle));
608  allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
609  }
610  }
611 
612  Node *pNode{new Node(m_hierarchy, leadingParticles, allHits, this->m_tier + 1)};
613  m_children.emplace_back(pNode);
614  for (const MCParticle *pChild : childParticles)
615  pNode->FillHierarchy(pChild, foldParameters);
616  }
617  else
618  {
619  MCParticleList allParticles{pRoot};
620  const int pdg{std::abs(pRoot->GetParticleId())};
621  const bool isShower{pdg == E_MINUS || pdg == PHOTON};
622  const bool isNeutron{pdg == NEUTRON};
623 
624  if (foldParameters.m_foldToTier && LArMCParticleHelper::GetHierarchyTier(pRoot) >= foldParameters.m_tier)
626  else if (foldParameters.m_foldToLeadingShowers && (isShower || (isNeutron && !m_hierarchy.m_recoCriteria.m_removeNeutrons)))
628  else if (m_hierarchy.m_recoCriteria.m_removeNeutrons && isNeutron)
629  return;
630 
631  CaloHitList allHits;
632  for (const MCParticle *pMCParticle : allParticles)
633  {
634  // ATTN - Not all MC particles will have hits
635  if (m_hierarchy.m_mcToHitsMap.find(pMCParticle) != m_hierarchy.m_mcToHitsMap.end())
636  {
637  const CaloHitList &caloHits(m_hierarchy.m_mcToHitsMap.at(pMCParticle));
638  allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
639  }
640  }
641 
642  if (!allParticles.empty())
643  {
644  const bool hasChildren{(foldParameters.m_foldToTier && LArMCParticleHelper::GetHierarchyTier(pRoot) < foldParameters.m_tier) ||
645  (!foldParameters.m_foldToTier && !foldParameters.m_foldToLeadingShowers) ||
646  (foldParameters.m_foldToLeadingShowers && !(isShower || isNeutron))};
647  // Only add the node if it either has children, or is a leaf node with hits
648  if (hasChildren || (!hasChildren && !allHits.empty()))
649  {
650  Node *pNode{new Node(m_hierarchy, allParticles, allHits, this->m_tier + 1)};
651  m_children.emplace_back(pNode);
652  if (hasChildren)
653  {
654  // Find the children of this particle and recursively add them to the hierarchy
655  const MCParticleList &children{pRoot->GetDaughterList()};
656  for (const MCParticle *pChild : children)
657  pNode->FillHierarchy(pChild, foldParameters);
658  }
659  }
660  }
661  }
662 }
663 
664 //------------------------------------------------------------------------------------------------------------------------------------------
665 
667 {
668  MCParticleList allParticles{pRoot};
669  if (!m_hierarchy.m_recoCriteria.m_removeNeutrons)
670  {
672  }
673  else
674  {
675  MCParticleList neutrons;
676  LArMCParticleHelper::GetAllDescendentMCParticles(pRoot, allParticles, allParticles, neutrons);
677  }
678  CaloHitList allHits;
679  for (const MCParticle *pMCParticle : allParticles)
680  {
681  // ATTN - Not all MC particles will have hits
682  if (m_hierarchy.m_mcToHitsMap.find(pMCParticle) != m_hierarchy.m_mcToHitsMap.end())
683  {
684  const CaloHitList &caloHits(m_hierarchy.m_mcToHitsMap.at(pMCParticle));
685  allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
686  }
687  }
688  if (!allParticles.empty())
689  {
690  Node *pNode{new Node(m_hierarchy, allParticles, allHits, this->m_tier + 1)};
691  m_children.emplace_back(pNode);
692  }
693 }
694 
695 //------------------------------------------------------------------------------------------------------------------------------------------
696 
698 {
699  return m_hierarchy.m_nodeToIdMap.at(this);
700 }
701 
702 //------------------------------------------------------------------------------------------------------------------------------------------
703 
705 {
706  const bool enoughHits{m_caloHits.size() >= m_hierarchy.m_recoCriteria.m_minHits};
707  if (!enoughHits)
708  return false;
709  bool enoughGoodViews{false};
710  unsigned int nHitsU{0}, nHitsV{0}, nHitsW{0};
711  for (const CaloHit *pCaloHit : m_caloHits)
712  {
713  switch (pCaloHit->GetHitType())
714  {
715  case TPC_VIEW_U:
716  ++nHitsU;
717  break;
718  case TPC_VIEW_V:
719  ++nHitsV;
720  break;
721  case TPC_VIEW_W:
722  ++nHitsW;
723  break;
724  default:
725  break;
726  }
727  unsigned int nGoodViews{0};
728  if (nHitsU >= m_hierarchy.m_recoCriteria.m_minHitsForGoodView)
729  ++nGoodViews;
730  if (nHitsV >= m_hierarchy.m_recoCriteria.m_minHitsForGoodView)
731  ++nGoodViews;
732  if (nHitsW >= m_hierarchy.m_recoCriteria.m_minHitsForGoodView)
733  ++nGoodViews;
734  if (nGoodViews >= m_hierarchy.m_recoCriteria.m_minGoodViews)
735  {
736  enoughGoodViews = true;
737  break;
738  }
739  }
740 
741  return enoughGoodViews;
742 }
743 
744 //------------------------------------------------------------------------------------------------------------------------------------------
745 
747 {
748  if (m_mainParticle)
749  return LArMCParticleHelper::IsBeamParticle(m_mainParticle);
750  else
751  return false;
752 }
753 
754 //------------------------------------------------------------------------------------------------------------------------------------------
755 
757 {
758  if (m_mainParticle)
759  return LArMCParticleHelper::IsCosmicRay(m_mainParticle);
760  else
761  return false;
762 }
763 
764 //------------------------------------------------------------------------------------------------------------------------------------------
765 
766 const std::string LArHierarchyHelper::MCHierarchy::Node::ToString(const std::string &prefix) const
767 {
768  std::string str(prefix + "PDG: " + std::to_string(m_pdg) + " Energy: " + std::to_string(m_mainParticle ? m_mainParticle->GetEnergy() : 0) +
769  " Hits: " + std::to_string(m_caloHits.size()) + "\n");
770  for (const Node *pChild : m_children)
771  str += pChild->ToString(prefix + " ");
772 
773  return str;
774 }
775 
776 //------------------------------------------------------------------------------------------------------------------------------------------
777 //------------------------------------------------------------------------------------------------------------------------------------------
778 
780  m_minHits{30},
782  m_minGoodViews{2},
783  m_removeNeutrons{true}
784 {
785 }
786 
787 //------------------------------------------------------------------------------------------------------------------------------------------
788 
790  m_minHits{obj.m_minHits},
791  m_minHitsForGoodView{obj.m_minHitsForGoodView},
792  m_minGoodViews{obj.m_minGoodViews},
793  m_removeNeutrons{obj.m_removeNeutrons}
794 {
795 }
796 
797 //------------------------------------------------------------------------------------------------------------------------------------------
798 
800  const unsigned int minHits, const unsigned int minHitsForGoodView, const unsigned int minGoodViews, const bool removeNeutrons) :
801  m_minHits{minHits},
802  m_minHitsForGoodView{minHitsForGoodView},
803  m_minGoodViews{minGoodViews},
804  m_removeNeutrons{removeNeutrons}
805 {
806 }
807 
808 //------------------------------------------------------------------------------------------------------------------------------------------
809 //------------------------------------------------------------------------------------------------------------------------------------------
810 
812 {
813 }
814 
815 //------------------------------------------------------------------------------------------------------------------------------------------
816 
818 {
819  for (const auto &[pRoot, nodeVector] : m_interactions)
820  {
821  for (const Node *pNode : nodeVector)
822  delete pNode;
823  }
824  m_interactions.clear();
825 }
826 
827 //------------------------------------------------------------------------------------------------------------------------------------------
828 
829 void LArHierarchyHelper::RecoHierarchy::FillHierarchy(const PfoList &pfoList, const FoldingParameters &foldParameters)
830 {
831  PfoList rootNodes;
832  for (const ParticleFlowObject *pPfo : pfoList)
833  {
834  const PfoList &parentList{pPfo->GetParentPfoList()};
835  if (parentList.empty())
836  {
837  rootNodes.emplace_back(pPfo);
838  }
839  }
840 
841  for (const ParticleFlowObject *const pRoot : rootNodes)
842  {
843  PfoSet primarySet;
844  LArHierarchyHelper::GetRecoPrimaries(pRoot, primarySet);
845  PfoList primaries(primarySet.begin(), primarySet.end());
846  primaries.sort(LArPfoHelper::SortByNHits);
847  if (foldParameters.m_foldToTier && foldParameters.m_tier == 1)
848  {
849  for (const ParticleFlowObject *pPrimary : primaries)
850  {
851  PfoList allParticles;
852  // ATTN - pPrimary gets added to the list of downstream PFOs, not just the child PFOs
853  LArPfoHelper::GetAllDownstreamPfos(pPrimary, allParticles);
854  CaloHitList allHits;
855  for (const ParticleFlowObject *pPfo : allParticles)
856  LArPfoHelper::GetAllCaloHits(pPfo, allHits);
857  m_interactions[pRoot].emplace_back(new Node(*this, allParticles, allHits));
858  }
859  }
860  else if (foldParameters.m_foldToLeadingShowers)
861  {
862  for (const ParticleFlowObject *pPrimary : primaries)
863  {
864  PfoList allParticles;
865  int pdg{std::abs(pPrimary->GetParticleId())};
866  const bool isShower{pdg == E_MINUS};
867  // ATTN - pPrimary gets added to the list of downstream PFOs, not just the child PFOs
868  if (isShower)
869  LArPfoHelper::GetAllDownstreamPfos(pPrimary, allParticles);
870  else
871  allParticles.emplace_back(pPrimary);
872 
873  CaloHitList allHits;
874  for (const ParticleFlowObject *pPfo : allParticles)
875  LArPfoHelper::GetAllCaloHits(pPfo, allHits);
876  Node *pNode{new Node(*this, allParticles, allHits)};
877  m_interactions[pRoot].emplace_back(pNode);
878  if (!isShower)
879  {
880  // Find the children of this particle and recursively add them to the hierarchy
881  const PfoList &children{pPrimary->GetDaughterPfoList()};
882  for (const ParticleFlowObject *pChild : children)
883  pNode->FillHierarchy(pChild, foldParameters);
884  }
885  }
886  }
887  else
888  {
889  // Dynamic fold, Unfolded and fold to tier > 1 have the same behaviour for primaries
890  for (const ParticleFlowObject *pPrimary : primaries)
891  {
892  PfoList allParticles{pPrimary};
893  CaloHitList allHits;
894  for (const ParticleFlowObject *pPfo : allParticles)
895  LArPfoHelper::GetAllCaloHits(pPfo, allHits);
896  Node *pNode{new Node(*this, allParticles, allHits)};
897  m_interactions[pRoot].emplace_back(pNode);
898  // Find the children of this particle and recursively add them to the hierarchy
899  const PfoList &children{pPrimary->GetDaughterPfoList()};
900  for (const ParticleFlowObject *pChild : children)
901  pNode->FillHierarchy(pChild, foldParameters.m_foldToLeadingShowers);
902  }
903  }
904  }
905 }
906 
907 //------------------------------------------------------------------------------------------------------------------------------------------
908 
910 {
911  if (m_interactions.find(pRoot) == m_interactions.end())
912  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
913 
914  return m_interactions.at(pRoot);
915 }
916 
917 //------------------------------------------------------------------------------------------------------------------------------------------
918 
920 {
921  for (auto iter = m_interactions.begin(); iter != m_interactions.end(); ++iter)
922  rootPfos.emplace_back(iter->first);
923 }
924 
925 //------------------------------------------------------------------------------------------------------------------------------------------
926 
927 void LArHierarchyHelper::RecoHierarchy::GetFlattenedNodes(const ParticleFlowObject *pRoot, NodeVector &nodeVector) const
928 {
929  NodeList queue;
930  for (const Node *pNode : m_interactions.at(pRoot))
931  {
932  nodeVector.emplace_back(pNode);
933  queue.emplace_back(pNode);
934  }
935  while (!queue.empty())
936  {
937  const NodeVector &children{queue.front()->GetChildren()};
938  queue.pop_front();
939  for (const Node *pChild : children)
940  {
941  nodeVector.emplace_back(pChild);
942  queue.emplace_back(pChild);
943  }
944  }
945 }
946 
947 //------------------------------------------------------------------------------------------------------------------------------------------
948 
950 {
951  std::string str;
952  for (const auto &[pRoot, nodeVector] : m_interactions)
953  {
954  str += "=== Reco Interaction : PDG " + std::to_string(pRoot->GetParticleId()) + "\n";
955  for (const Node *pNode : nodeVector)
956  str += " " + pNode->ToString("") + "\n";
957  }
958 
959  return str;
960 }
961 
962 //------------------------------------------------------------------------------------------------------------------------------------------
963 //------------------------------------------------------------------------------------------------------------------------------------------
964 
965 LArHierarchyHelper::RecoHierarchy::Node::Node(const RecoHierarchy &hierarchy, const ParticleFlowObject *pPfo) :
966  m_hierarchy{hierarchy},
967  m_mainPfo{pPfo},
968  m_pdg{0}
969 {
970  if (pPfo)
971  {
972  m_pdg = pPfo->GetParticleId();
973  m_pfos.emplace_back(pPfo);
974  }
975 }
976 
977 //------------------------------------------------------------------------------------------------------------------------------------------
978 
979 LArHierarchyHelper::RecoHierarchy::Node::Node(const RecoHierarchy &hierarchy, const PfoList &pfoList, const CaloHitList &caloHitList) :
980  m_hierarchy(hierarchy),
981  m_mainPfo{nullptr},
982  m_pdg{0}
983 {
984  if (!pfoList.empty())
985  {
986  m_mainPfo = pfoList.front();
987  m_pdg = pfoList.front()->GetParticleId();
988  }
989  m_pfos = pfoList;
991  m_caloHits = caloHitList;
992  m_caloHits.sort();
993 }
994 
995 //------------------------------------------------------------------------------------------------------------------------------------------
996 
998 {
999  m_pfos.clear();
1000  m_caloHits.clear();
1001  for (const Node *node : m_children)
1002  delete node;
1003  m_children.clear();
1004 }
1005 
1006 //------------------------------------------------------------------------------------------------------------------------------------------
1007 
1008 void LArHierarchyHelper::RecoHierarchy::Node::FillHierarchy(const ParticleFlowObject *pRoot, const FoldingParameters &foldParameters)
1009 {
1010  PfoList allParticles;
1011  int pdg{std::abs(pRoot->GetParticleId())};
1012  const bool isShower{pdg == E_MINUS};
1013  if (foldParameters.m_foldToTier && LArPfoHelper::GetHierarchyTier(pRoot) >= foldParameters.m_tier)
1014  LArPfoHelper::GetAllDownstreamPfos(pRoot, allParticles);
1015  else if (foldParameters.m_foldToLeadingShowers && isShower)
1016  LArPfoHelper::GetAllDownstreamPfos(pRoot, allParticles);
1017  else
1018  allParticles.emplace_back(pRoot);
1019 
1020  CaloHitList allHits;
1021  for (const ParticleFlowObject *pPfo : allParticles)
1022  LArPfoHelper::GetAllCaloHits(pPfo, allHits);
1023  const bool hasChildren{(foldParameters.m_foldToTier && LArPfoHelper::GetHierarchyTier(pRoot) < foldParameters.m_tier) ||
1024  (!foldParameters.m_foldToTier && !foldParameters.m_foldToLeadingShowers) || (foldParameters.m_foldToLeadingShowers && !isShower)};
1025 
1026  if (hasChildren || (!hasChildren && !allHits.empty()))
1027  {
1028  Node *pNode{new Node(m_hierarchy, allParticles, allHits)};
1029  m_children.emplace_back(pNode);
1030 
1031  if (hasChildren)
1032  {
1033  const PfoList &children{pRoot->GetDaughterPfoList()};
1034  for (const ParticleFlowObject *pChild : children)
1035  pNode->FillHierarchy(pChild, foldParameters);
1036  }
1037  }
1038 }
1039 
1040 //------------------------------------------------------------------------------------------------------------------------------------------
1041 
1042 void LArHierarchyHelper::RecoHierarchy::Node::FillFlat(const ParticleFlowObject *pRoot)
1043 {
1044  PfoList allParticles;
1045  LArPfoHelper::GetAllDownstreamPfos(pRoot, allParticles);
1046  CaloHitList allHits;
1047  for (const ParticleFlowObject *pPfo : allParticles)
1048  LArPfoHelper::GetAllCaloHits(pPfo, allHits);
1049  Node *pNode{new Node(m_hierarchy, allParticles, allHits)};
1050  m_children.emplace_back(pNode);
1051 }
1052 
1053 //------------------------------------------------------------------------------------------------------------------------------------------
1054 
1056 {
1057  return m_pfos;
1058 }
1059 
1060 //------------------------------------------------------------------------------------------------------------------------------------------
1061 
1063 {
1064  return m_caloHits;
1065 }
1066 
1067 //------------------------------------------------------------------------------------------------------------------------------------------
1068 
1070 {
1071  return m_pdg;
1072 }
1073 
1074 //------------------------------------------------------------------------------------------------------------------------------------------
1075 
1076 const std::string LArHierarchyHelper::RecoHierarchy::Node::ToString(const std::string &prefix) const
1077 {
1078  std::string str(prefix + "PDG: " + std::to_string(m_pdg) + " Hits: " + std::to_string(m_caloHits.size()) + "\n");
1079  for (const Node *pChild : m_children)
1080  str += pChild->ToString(prefix + " ");
1081 
1082  return str;
1083 }
1084 
1085 //------------------------------------------------------------------------------------------------------------------------------------------
1086 //------------------------------------------------------------------------------------------------------------------------------------------
1087 
1089  m_pMCParticle{pMCParticle}
1090 {
1091 }
1092 
1093 //------------------------------------------------------------------------------------------------------------------------------------------
1094 
1095 void LArHierarchyHelper::MCMatches::AddRecoMatch(const RecoHierarchy::Node *pReco, const int nSharedHits)
1096 {
1097  m_recoNodes.emplace_back(pReco);
1098  m_sharedHits.emplace_back(nSharedHits);
1099 }
1100 
1101 //------------------------------------------------------------------------------------------------------------------------------------------
1102 
1104 {
1105  auto iter{std::find(m_recoNodes.begin(), m_recoNodes.end(), pReco)};
1106  if (iter == m_recoNodes.end())
1107  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1108  int index = iter - m_recoNodes.begin();
1109 
1110  return static_cast<int>(m_sharedHits[index]);
1111 }
1112 
1113 //------------------------------------------------------------------------------------------------------------------------------------------
1114 
1115 float LArHierarchyHelper::MCMatches::GetPurity(const RecoHierarchy::Node *pReco, const bool adcWeighted) const
1116 {
1117  auto iter{std::find(m_recoNodes.begin(), m_recoNodes.end(), pReco)};
1118  if (iter == m_recoNodes.end())
1119  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1120 
1121  const CaloHitList &recoHits{pReco->GetCaloHits()};
1122  const CaloHitList &mcHits{m_pMCParticle->GetCaloHits()};
1123  CaloHitVector intersection;
1124  std::set_intersection(mcHits.begin(), mcHits.end(), recoHits.begin(), recoHits.end(), std::back_inserter(intersection));
1125 
1126  return this->GetPurity(intersection, recoHits, adcWeighted);
1127 }
1128 
1129 //------------------------------------------------------------------------------------------------------------------------------------------
1130 
1131 float LArHierarchyHelper::MCMatches::GetPurity(const RecoHierarchy::Node *pReco, const HitType view, const bool adcWeighted) const
1132 {
1133  auto iter{std::find(m_recoNodes.begin(), m_recoNodes.end(), pReco)};
1134  if (iter == m_recoNodes.end())
1135  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1136 
1137  CaloHitList recoHits;
1138  for (const CaloHit *pCaloHit : pReco->GetCaloHits())
1139  if (pCaloHit->GetHitType() == view)
1140  recoHits.emplace_back(pCaloHit);
1141  CaloHitList mcHits;
1142  for (const CaloHit *pCaloHit : m_pMCParticle->GetCaloHits())
1143  if (pCaloHit->GetHitType() == view)
1144  mcHits.emplace_back(pCaloHit);
1145 
1146  CaloHitVector intersection;
1147  std::set_intersection(mcHits.begin(), mcHits.end(), recoHits.begin(), recoHits.end(), std::back_inserter(intersection));
1148 
1149  return this->GetPurity(intersection, recoHits, adcWeighted);
1150 }
1151 
1152 //------------------------------------------------------------------------------------------------------------------------------------------
1153 
1154 float LArHierarchyHelper::MCMatches::GetCompleteness(const RecoHierarchy::Node *pReco, const bool adcWeighted) const
1155 {
1156  auto iter{std::find(m_recoNodes.begin(), m_recoNodes.end(), pReco)};
1157  if (iter == m_recoNodes.end())
1158  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1159 
1160  const CaloHitList &recoHits{pReco->GetCaloHits()};
1161  const CaloHitList &mcHits{m_pMCParticle->GetCaloHits()};
1162  CaloHitVector intersection;
1163  std::set_intersection(mcHits.begin(), mcHits.end(), recoHits.begin(), recoHits.end(), std::back_inserter(intersection));
1164 
1165  return this->GetCompleteness(intersection, mcHits, adcWeighted);
1166 }
1167 
1168 //------------------------------------------------------------------------------------------------------------------------------------------
1169 
1170 float LArHierarchyHelper::MCMatches::GetCompleteness(const RecoHierarchy::Node *pReco, const HitType view, const bool adcWeighted) const
1171 {
1172  auto iter{std::find(m_recoNodes.begin(), m_recoNodes.end(), pReco)};
1173  if (iter == m_recoNodes.end())
1174  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1175 
1176  CaloHitList recoHits;
1177  for (const CaloHit *pCaloHit : pReco->GetCaloHits())
1178  if (pCaloHit->GetHitType() == view)
1179  recoHits.emplace_back(pCaloHit);
1180  CaloHitList mcHits;
1181  for (const CaloHit *pCaloHit : m_pMCParticle->GetCaloHits())
1182  if (pCaloHit->GetHitType() == view)
1183  mcHits.emplace_back(pCaloHit);
1184 
1185  CaloHitVector intersection;
1186  std::set_intersection(mcHits.begin(), mcHits.end(), recoHits.begin(), recoHits.end(), std::back_inserter(intersection));
1187 
1188  return this->GetCompleteness(intersection, mcHits, adcWeighted);
1189 }
1190 
1191 //------------------------------------------------------------------------------------------------------------------------------------------
1192 
1193 float LArHierarchyHelper::MCMatches::GetPurity(const CaloHitVector &intersection, const CaloHitList &recoHits, const bool adcWeighted) const
1194 {
1195  float purity{0.f};
1196  if (!intersection.empty())
1197  {
1198  if (adcWeighted)
1199  {
1200  float adcSum{0.f};
1201  for (const CaloHit *pCaloHit : recoHits)
1202  adcSum += pCaloHit->GetInputEnergy();
1203  if (adcSum > std::numeric_limits<float>::epsilon())
1204  {
1205  for (const CaloHit *pCaloHit : intersection)
1206  purity += pCaloHit->GetInputEnergy();
1207  purity /= adcSum;
1208  }
1209  }
1210  else
1211  {
1212  purity = intersection.size() / static_cast<float>(recoHits.size());
1213  }
1214  }
1215 
1216  return purity;
1217 }
1218 
1219 //------------------------------------------------------------------------------------------------------------------------------------------
1220 
1221 float LArHierarchyHelper::MCMatches::GetCompleteness(const CaloHitVector &intersection, const CaloHitList &mcHits, const bool adcWeighted) const
1222 {
1223  float completeness{0.f};
1224  if (!intersection.empty())
1225  {
1226  if (adcWeighted)
1227  {
1228  float adcSum{0.f};
1229  for (const CaloHit *pCaloHit : mcHits)
1230  adcSum += pCaloHit->GetInputEnergy();
1231  if (adcSum > std::numeric_limits<float>::epsilon())
1232  {
1233  for (const CaloHit *pCaloHit : intersection)
1234  completeness += pCaloHit->GetInputEnergy();
1235  completeness /= adcSum;
1236  }
1237  }
1238  else
1239  {
1240  completeness = intersection.size() / static_cast<float>(mcHits.size());
1241  }
1242  }
1243 
1244  return completeness;
1245 }
1246 
1247 //------------------------------------------------------------------------------------------------------------------------------------------
1248 
1250 {
1251  if (m_recoNodes.empty())
1252  return false;
1253 
1254  int nAboveThreshold{0};
1255  if (m_recoNodes.size() != 1)
1256  {
1257  for (const LArHierarchyHelper::RecoHierarchy::Node *const pNode : m_recoNodes)
1258  if (this->GetCompleteness(pNode) > 0.1f)
1259  ++nAboveThreshold;
1260  if (nAboveThreshold != 1)
1261  return false;
1262  }
1263 
1264  if (this->GetPurity(m_recoNodes.front()) < qualityCuts.m_minPurity)
1265  return false;
1266 
1267  if (this->GetCompleteness(m_recoNodes.front()) < qualityCuts.m_minCompleteness)
1268  return false;
1269 
1270  return true;
1271 }
1272 
1273 //------------------------------------------------------------------------------------------------------------------------------------------
1274 //------------------------------------------------------------------------------------------------------------------------------------------
1275 
1276 LArHierarchyHelper::MatchInfo::MatchInfo(const MCHierarchy &mcHierarchy, const RecoHierarchy &recoHierarchy) :
1277  MatchInfo(mcHierarchy, recoHierarchy, QualityCuts())
1278 {
1279 }
1280 
1281 //------------------------------------------------------------------------------------------------------------------------------------------
1282 
1283 LArHierarchyHelper::MatchInfo::MatchInfo(const MCHierarchy &mcHierarchy, const RecoHierarchy &recoHierarchy, const QualityCuts &qualityCuts) :
1284  m_mcHierarchy{mcHierarchy},
1285  m_recoHierarchy{recoHierarchy},
1286  m_qualityCuts{qualityCuts}
1287 {
1288 }
1289 
1290 //------------------------------------------------------------------------------------------------------------------------------------------
1291 
1293 {
1294  MCParticleList rootMCParticles;
1295  m_mcHierarchy.GetRootMCParticles(rootMCParticles);
1296  PfoList rootPfos;
1297  m_recoHierarchy.GetRootPfos(rootPfos);
1298  std::map<const MCHierarchy::Node *, MCMatches> mcToMatchMap;
1299 
1300  for (const MCParticle *const pRootMC : rootMCParticles)
1301  {
1302  MCHierarchy::NodeVector mcNodes;
1303  m_mcHierarchy.GetFlattenedNodes(pRootMC, mcNodes);
1304  for (const ParticleFlowObject *const pRootPfo : rootPfos)
1305  {
1306  RecoHierarchy::NodeVector recoNodes;
1307  m_recoHierarchy.GetFlattenedNodes(pRootPfo, recoNodes);
1308 
1309  std::sort(mcNodes.begin(), mcNodes.end(),
1310  [](const MCHierarchy::Node *lhs, const MCHierarchy::Node *rhs)
1311  { return lhs->GetCaloHits().size() > rhs->GetCaloHits().size(); });
1312  std::sort(recoNodes.begin(), recoNodes.end(),
1313  [](const RecoHierarchy::Node *lhs, const RecoHierarchy::Node *rhs)
1314  { return lhs->GetCaloHits().size() > rhs->GetCaloHits().size(); });
1315 
1316  for (const RecoHierarchy::Node *pRecoNode : recoNodes)
1317  {
1318  const CaloHitList &recoHits{pRecoNode->GetCaloHits()};
1319  const MCHierarchy::Node *pBestNode{nullptr};
1320  size_t bestSharedHits{0};
1321  for (const MCHierarchy::Node *pMCNode : mcNodes)
1322  {
1323  if (!pMCNode->IsReconstructable())
1324  continue;
1325  const CaloHitList &mcHits{pMCNode->GetCaloHits()};
1326  CaloHitVector intersection;
1327  std::set_intersection(mcHits.begin(), mcHits.end(), recoHits.begin(), recoHits.end(), std::back_inserter(intersection));
1328 
1329  if (!intersection.empty())
1330  {
1331  const size_t sharedHits{intersection.size()};
1332  if (sharedHits > bestSharedHits)
1333  {
1334  bestSharedHits = sharedHits;
1335  pBestNode = pMCNode;
1336  }
1337  }
1338  }
1339  if (pBestNode)
1340  {
1341  auto iter{mcToMatchMap.find(pBestNode)};
1342  if (iter != mcToMatchMap.end())
1343  {
1344  MCMatches &match(iter->second);
1345  match.AddRecoMatch(pRecoNode, static_cast<int>(bestSharedHits));
1346  }
1347  else
1348  {
1349  MCMatches match(pBestNode);
1350  match.AddRecoMatch(pRecoNode, static_cast<int>(bestSharedHits));
1351  mcToMatchMap.insert(std::make_pair(pBestNode, match));
1352  }
1353  }
1354  else
1355  {
1356  m_unmatchedReco.emplace_back(pRecoNode);
1357  }
1358  }
1359  }
1360  }
1361 
1362  for (auto [pMCNode, matches] : mcToMatchMap)
1363  {
1364  // We need to figure out which MC interaction hierarchy the matches belongs to
1365  for (const MCParticle *const pRootMC : rootMCParticles)
1366  {
1367  MCHierarchy::NodeVector mcNodes;
1368  m_mcHierarchy.GetFlattenedNodes(pRootMC, mcNodes);
1369  if (std::find(mcNodes.begin(), mcNodes.end(), pMCNode) != mcNodes.end())
1370  {
1371  m_matches[pRootMC].emplace_back(matches);
1372  break;
1373  }
1374  }
1375  }
1376 
1377  const auto predicate = [](const MCMatches &lhs, const MCMatches &rhs)
1378  { return lhs.GetMC()->GetCaloHits().size() > rhs.GetMC()->GetCaloHits().size(); };
1379 
1380  for (const MCParticle *const pRootMC : rootMCParticles)
1381  {
1382  std::sort(m_matches[pRootMC].begin(), m_matches[pRootMC].end(), predicate);
1383 
1384  MCHierarchy::NodeVector mcNodes;
1385  m_mcHierarchy.GetFlattenedNodes(pRootMC, mcNodes);
1386 
1387  for (const MCHierarchy::Node *pMCNode : mcNodes)
1388  {
1389  if (pMCNode->IsReconstructable() && mcToMatchMap.find(pMCNode) == mcToMatchMap.end())
1390  {
1391  MCMatches match(pMCNode);
1392  m_matches[pRootMC].emplace_back(match);
1393  }
1394  }
1395  }
1396 }
1397 
1398 //------------------------------------------------------------------------------------------------------------------------------------------
1399 
1400 unsigned int LArHierarchyHelper::MatchInfo::GetNMCNodes(const MCParticle *const pRoot) const
1401 {
1402  if (m_matches.find(pRoot) == m_matches.end())
1403  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1404 
1405  return static_cast<unsigned int>(m_matches.at(pRoot).size());
1406 }
1407 
1408 //------------------------------------------------------------------------------------------------------------------------------------------
1409 
1410 unsigned int LArHierarchyHelper::MatchInfo::GetNNeutrinoMCNodes(const MCParticle *const pRoot) const
1411 {
1412  if (m_matches.find(pRoot) == m_matches.end())
1413  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1414 
1415  unsigned int nNodes{0};
1416  for (const MCMatches &match : m_matches.at(pRoot))
1417  {
1418  const MCHierarchy::Node *pNode{match.GetMC()};
1419  if (!(pNode->IsCosmicRay() || pNode->IsTestBeamParticle()))
1420  ++nNodes;
1421  }
1422 
1423  return nNodes;
1424 }
1425 
1426 //------------------------------------------------------------------------------------------------------------------------------------------
1427 
1428 unsigned int LArHierarchyHelper::MatchInfo::GetNCosmicRayMCNodes(const MCParticle *const pRoot) const
1429 {
1430  if (m_matches.find(pRoot) == m_matches.end())
1431  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1432 
1433  unsigned int nNodes{0};
1434  for (const MCMatches &match : m_matches.at(pRoot))
1435  {
1436  const MCHierarchy::Node *pNode{match.GetMC()};
1437  if (pNode->IsCosmicRay())
1438  ++nNodes;
1439  }
1440 
1441  return nNodes;
1442 }
1443 
1444 //------------------------------------------------------------------------------------------------------------------------------------------
1445 
1446 unsigned int LArHierarchyHelper::MatchInfo::GetNTestBeamMCNodes(const MCParticle *const pRoot) const
1447 {
1448  if (m_matches.find(pRoot) == m_matches.end())
1449  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1450 
1451  unsigned int nNodes{0};
1452  for (const MCMatches &match : m_matches.at(pRoot))
1453  {
1454  const MCHierarchy::Node *pNode{match.GetMC()};
1455  if (pNode->IsTestBeamParticle())
1456  ++nNodes;
1457  }
1458 
1459  return nNodes;
1460 }
1461 
1462 //------------------------------------------------------------------------------------------------------------------------------------------
1463 
1465 {
1466  (void)mcHierarchy;
1467  MCParticleList rootMCParticles;
1468  mcHierarchy.GetRootMCParticles(rootMCParticles);
1469 
1470  for (const MCParticle *const pRootMC : rootMCParticles)
1471  {
1472  const LArHierarchyHelper::MCMatchesVector &matches{this->GetMatches(pRootMC)};
1473 
1474  MCParticleList primaries;
1475  for (const LArHierarchyHelper::MCMatches &match : matches)
1476  {
1477  const LArHierarchyHelper::MCHierarchy::Node *pMCNode{match.GetMC()};
1478  if (pMCNode->GetHierarchyTier() == 1)
1479  {
1480  const MCParticle *const pLeadingMC{pMCNode->GetLeadingMCParticle()};
1481  primaries.emplace_back(pLeadingMC);
1482  }
1483  }
1484  primaries.sort(LArMCParticleHelper::SortByMomentum);
1486 
1487  const LArMCParticle *const pLArRoot{dynamic_cast<const LArMCParticle *const>(pRootMC)};
1488  if (pLArRoot)
1489  std::cout << "=== MC Interaction : PDG " << std::to_string(pLArRoot->GetParticleId())
1490  << " Energy: " << std::to_string(pLArRoot->GetEnergy()) << " Type: " << descriptor.ToString() << std::endl;
1491  else
1492  std::cout << "=== MC Interaction : PDG " << std::to_string(pRootMC->GetParticleId())
1493  << " Energy: " << std::to_string(pRootMC->GetEnergy()) << " Type: " << descriptor.ToString() << std::endl;
1494 
1495  unsigned int nNeutrinoMCParticles{this->GetNNeutrinoMCNodes(pRootMC)}, nNeutrinoRecoParticles{0};
1496  unsigned int nCosmicMCParticles{this->GetNCosmicRayMCNodes(pRootMC)}, nCosmicRecoParticles{0};
1497  unsigned int nTestBeamMCParticles{this->GetNTestBeamMCNodes(pRootMC)}, nTestBeamRecoParticles{0};
1498  std::cout << " === Matches ===" << std::endl;
1499  std::cout << std::fixed << std::setprecision(2);
1500  for (const MCMatches &match : m_matches.at(pRootMC))
1501  {
1502  const MCHierarchy::Node *pMCNode{match.GetMC()};
1503  const int pdg{pMCNode->GetParticleId()};
1504  const size_t mcHits{pMCNode->GetCaloHits().size()};
1505  const std::string tag{pMCNode->IsTestBeamParticle() ? "(Beam) " : pMCNode->IsCosmicRay() ? "(Cosmic) " : ""};
1506  std::cout << " MC " << tag << pdg << " hits " << mcHits << std::endl;
1507  const RecoHierarchy::NodeVector &nodeVector{match.GetRecoMatches()};
1508 
1509  for (const RecoHierarchy::Node *pRecoNode : nodeVector)
1510  {
1511  const unsigned int recoHits{static_cast<unsigned int>(pRecoNode->GetCaloHits().size())};
1512  const unsigned int sharedHits{match.GetSharedHits(pRecoNode)};
1513  const float purity{match.GetPurity(pRecoNode)};
1514  const float completeness{match.GetCompleteness(pRecoNode)};
1515  if (completeness > 0.1f)
1516  std::cout << " Matched " << sharedHits << " out of " << recoHits << " with purity " << purity << " and completeness "
1517  << completeness << std::endl;
1518  else
1519  std::cout << " (Below threshold) " << sharedHits << " out of " << recoHits << " with purity " << purity
1520  << " and completeness " << completeness << std::endl;
1521  }
1522  if (nodeVector.empty())
1523  {
1524  std::cout << " Unmatched" << std::endl;
1525  }
1526  else if (match.IsQuality(this->GetQualityCuts()))
1527  {
1528  if (pMCNode->IsTestBeamParticle())
1529  ++nTestBeamRecoParticles;
1530  else if (pMCNode->IsCosmicRay())
1531  ++nCosmicRecoParticles;
1532  else
1533  ++nNeutrinoRecoParticles;
1534  }
1535  }
1536 
1537  if (LArMCParticleHelper::IsNeutrino(pRootMC))
1538  {
1539  std::cout << " Neutrino Interaction Summary:" << std::endl;
1540  if (nNeutrinoMCParticles)
1541  {
1542  std::cout << " Good final state particles: " << nNeutrinoRecoParticles << " of " << nNeutrinoMCParticles << " : "
1543  << (100 * nNeutrinoRecoParticles / static_cast<float>(nNeutrinoMCParticles)) << "%" << std::endl;
1544  }
1545  }
1546  else if (LArMCParticleHelper::IsCosmicRay(pRootMC))
1547  {
1548  std::cout << " Cosmic Ray Interaction Summary:" << std::endl;
1549  std::cout << std::fixed << std::setprecision(1);
1550  if (nCosmicMCParticles)
1551  {
1552  std::cout << " Good cosmics: " << nCosmicRecoParticles << " of " << nCosmicMCParticles << " : "
1553  << (100 * nCosmicRecoParticles / static_cast<float>(nCosmicMCParticles)) << "%" << std::endl;
1554  }
1555  }
1556  else if (LArMCParticleHelper::IsBeamParticle(pRootMC))
1557  {
1558  std::cout << " Test Beam Interaction Summary:" << std::endl;
1559  std::cout << std::fixed << std::setprecision(1);
1560  if (nTestBeamMCParticles)
1561  {
1562  std::cout << " Good test beam particles: " << nTestBeamRecoParticles << " of " << nTestBeamMCParticles << " : "
1563  << (100 * nTestBeamRecoParticles / static_cast<float>(nTestBeamMCParticles)) << "%" << std::endl;
1564  }
1565  if (nCosmicMCParticles)
1566  {
1567  std::cout << " Matched cosmics: " << nCosmicRecoParticles << " of " << nCosmicMCParticles << " : "
1568  << (100 * nCosmicRecoParticles / static_cast<float>(nCosmicMCParticles)) << "%" << std::endl;
1569  }
1570  }
1571  if (!this->GetUnmatchedReco().empty())
1572  std::cout << " Unmatched reco: " << this->GetUnmatchedReco().size() << std::endl;
1573  }
1574 }
1575 
1576 //------------------------------------------------------------------------------------------------------------------------------------------
1577 
1578 void LArHierarchyHelper::MatchInfo::GetRootMCParticles(MCParticleList &rootMCParticles) const
1579 {
1580  for (auto iter = m_matches.begin(); iter != m_matches.end(); ++iter)
1581  rootMCParticles.emplace_back(iter->first);
1582 }
1583 
1584 //------------------------------------------------------------------------------------------------------------------------------------------
1585 //------------------------------------------------------------------------------------------------------------------------------------------
1586 
1588  const MCParticleList &mcParticleList, const CaloHitList &caloHitList, const FoldingParameters &foldParameters, MCHierarchy &hierarchy)
1589 {
1590  hierarchy.FillHierarchy(mcParticleList, caloHitList, foldParameters);
1591 }
1592 
1593 //------------------------------------------------------------------------------------------------------------------------------------------
1594 
1595 void LArHierarchyHelper::FillRecoHierarchy(const PfoList &pfoList, const FoldingParameters &foldParameters, RecoHierarchy &hierarchy)
1596 {
1597  hierarchy.FillHierarchy(pfoList, foldParameters);
1598 }
1599 
1600 //------------------------------------------------------------------------------------------------------------------------------------------
1601 
1603 {
1604  matchInfo.Match();
1605 }
1606 
1607 //------------------------------------------------------------------------------------------------------------------------------------------
1608 // private
1609 
1610 void LArHierarchyHelper::GetMCPrimaries(const MCParticle *pRoot, MCParticleSet &primaries)
1611 {
1612  try
1613  {
1614  MCParticleList visible;
1616  // Check if we still need to use sets here
1617  for (const MCParticle *const pMCParticle : visible)
1618  primaries.insert(pMCParticle);
1619  }
1620  catch (const StatusCodeException &)
1621  {
1622  if (pRoot->GetParticleId() != 111 && pRoot->GetParticleId() < 1e9)
1623  std::cout << "LArHierarchyHelper::MCHierarchy::FillHierarchy: MC particle with PDG code " << pRoot->GetParticleId()
1624  << " at address " << pRoot << " has no associated primary particle" << std::endl;
1625  }
1626 }
1627 
1628 void LArHierarchyHelper::GetRecoPrimaries(const ParticleFlowObject *pRoot, PfoSet &primaries)
1629 {
1630  if (LArPfoHelper::IsNeutrino(pRoot))
1631  {
1632  const PfoList &children{pRoot->GetDaughterPfoList()};
1633  // Check if we still need to use sets here
1634  for (const ParticleFlowObject *const pPfo : children)
1635  primaries.insert(pPfo);
1636  }
1637  else
1638  {
1639  // Might want different handling here for test beam and cosmic ray particles, but for now, just treat
1640  // a non-neutrino root node as something to be stored in the primaries set directly
1641  primaries.insert(pRoot);
1642  }
1643 }
1644 
1645 } // namespace lar_content
static void GetRecoPrimaries(const pandora::ParticleFlowObject *pRoot, PfoSet &primaries)
Retrieves the primary PFOs from a list and returns the root (neutrino) for hierarchy, if it exists.
const pandora::CaloHitList & GetCaloHits() const
Retrieve the CaloHits associated with this node.
static void GetFirstVisibleMCParticles(const pandora::MCParticle *const pRoot, pandora::MCParticleList &visibleParticleList)
Get the first visible MC particles given a root particle. For example, given a neutrino this would re...
const unsigned int m_minHitsForGoodView
the minimum number of Hits for a good view
std::set< const pandora::MCParticle * > MCParticleSet
void GetRootMCParticles(pandora::MCParticleList &rootMCParticles) const
Retrieve the root MC particles of the interaction hierarchies.
static bool SortByNHits(const pandora::ParticleFlowObject *const pLhs, const pandora::ParticleFlowObject *const pRhs)
Sort pfos by number of constituent hits.
const MCMatchesVector & GetMatches(const pandora::MCParticle *const pRoot) const
Retrieve the vector of matches (this will include null matches - i.e. MC nodes with no corresponding ...
void GetFlattenedNodes(const pandora::MCParticle *const pRoot, NodeVector &nodeVector) const
Retrieve a flat vector of the ndoes in the hierarchy.
const MCHierarchy::Node * m_pMCParticle
MC node associated with any matches.
bool m_foldToTier
Whether or not to apply folding based on particle tier.
int m_tier
If folding to a tier, the tier to be combined with its child particles.
Header file for the interaction type helper class.
RecoHierarchy::NodeVector m_recoNodes
Matched reco nodes.
const NodeVector & GetInteractions(const pandora::ParticleFlowObject *pRoot) const
Retrieve the root nodes in the hierarchy for a given interaction.
void CollectContinuations(const pandora::MCParticle *pRoot, pandora::MCParticleList &continuingParticles, pandora::MCParticleList &childParticles, const float cosAngleTolerance) const
Identify downstream particles that represent continuations of the parent particle from a reconstructi...
const std::string ToString() const
Produce a string representation of the hierarchy.
pandora::IntVector m_sharedHits
Number of shared hits for each match.
void GetRootMCParticles(pandora::MCParticleList &rootMCParticles) const
Retrieve the root MC particles of the interaction hierarchies.
const std::string ToString(const std::string &prefix) const
Produce a string representation of the hierarchy.
std::map< const pandora::MCParticle *, pandora::CaloHitList > m_mcToHitsMap
The map between MC particles and calo hits.
constexpr auto abs(T v)
Returns the absolute value of the argument.
static int GetHierarchyTier(const pandora::MCParticle *const pMCParticle)
Determine the position in the hierarchy for the MCParticle.
static void GetMCPrimaries(const pandora::MCParticle *pRoot, MCParticleSet &primaries)
Retrieves the primary MC particles from a list and returns the root (neutrino) for hierarchy...
const RecoHierarchy & m_hierarchy
The parent reco hierarchy.
int m_pdg
The particle ID (track = muon, shower = electron)
void FillHierarchy(const pandora::MCParticleList &mcParticleList, const pandora::CaloHitList &caloHitList, const FoldingParameters &foldParameters)
Creates an MC hierarchy representation. Without folding this will be a mirror image of the standard M...
QualityCuts m_qualityCuts
The quality cuts to be applied to matches.
const pandora::CaloHitList & GetCaloHits() const
Retrieve the CaloHits associated with this node.
MatchInfo(const MCHierarchy &mcHierarchy, const RecoHierarchy &recoHierarchy)
Default constructor.
static void MatchHierarchies(MatchInfo &matchInfo)
Finds the matches between reconstructed and MC hierarchies.
Node(const RecoHierarchy &hierarchy, const pandora::ParticleFlowObject *pPfo)
Create a node with a primary PFO.
static void GetAllCaloHits(const pandora::ParticleFlowObject *pPfo, pandora::CaloHitList &caloHitList)
Get a list of all calo hits (including isolated) of all types from a given pfo.
Definition: LArPfoHelper.cc:76
bool IsQuality(const QualityCuts &qualityCuts) const
Get whether this match passes quality cuts.
void Match()
Match the nodes in the MC and reco hierarchies.
const MCHierarchy::Node * GetMC() const
Retrieve the MC node.
TFile f
Definition: plotHisto.C:6
MCMatches(const MCHierarchy::Node *pMCParticle)
Constructor.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
unsigned int GetNTestBeamMCNodes(const pandora::MCParticle *const pRoot) const
Retrieve the number of test beam derived MC nodes available to match.
static bool SortByMomentum(const pandora::MCParticle *const pLhs, const pandora::MCParticle *const pRhs)
Sort mc particles by their momentum.
void FillHierarchy(const pandora::PfoList &pfoList, const FoldingParameters &foldParameters)
Creates a reconstructed hierarchy representation. Without folding this will be a mirror image of the ...
LAr mc particle class.
Definition: LArMCParticle.h:94
static bool IsCosmicRay(const pandora::MCParticle *const pMCParticle)
Return true if passed a primary cosmic ray MCParticle.
static bool IsInelasticScatter(const pandora::MCParticle *const pMCParticle)
Check whether or not an MC particle came from an inelastic scattering process.
static void FillRecoHierarchy(const pandora::PfoList &pfoList, const FoldingParameters &foldParameters, RecoHierarchy &hierarchy)
Fill a reconstructed hierarchy based on the specified folding criteria (see RecoHierarchy::FillHierar...
const pandora::MCParticle * GetLeadingMCParticle() const
Retrieve the leading MC particle associated with this node.
static int GetHierarchyTier(const pandora::ParticleFlowObject *const pPfo)
Determine the position in the hierarchy for the MCParticle.
decltype(auto) constexpr to_string(T &&obj)
ADL-aware version of std::to_string.
void InterpretHierarchy(const pandora::MCParticle *const pRoot, pandora::MCParticleList &leadingParticles, pandora::MCParticleList &childParticles, const float cosAngleTolerance) const
Interpret the hierarchy below a particular particle to determine if and how it should be folded...
const float m_minPurity
The minimum purity for a match to be considered good.
static bool IsBeamParticle(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary beam MCParticle.
bool IsReconstructable() const
Return whether or not this node should be considered reconstructable.
bool IsCosmicRay() const
Check if this is a cosmic ray particle.
void FillHierarchy(const pandora::ParticleFlowObject *pRoot, const FoldingParameters &foldParameters)
Recursively fill the hierarchy based on the criteria established for this RecoHierarchy.
int GetParticleId() const
Retrieve the PDG code for the leading particle in this node.
unsigned int GetSharedHits(const RecoHierarchy::Node *pReco) const
Retrieve the number of shared hits in the match.
float GetPurity(const RecoHierarchy::Node *pReco, const bool adcWeighted=false) const
Retrieve the purity of the match.
pandora::PfoList m_pfos
The list of PFOs of which this node is composed.
const std::string ToString(const std::string &prefix) const
Produce a string representation of the hierarchy.
float m_cosAngleTolerance
Cosine of the maximum angle at which topologies can be considered continuous.
static bool IsNeutrino(const pandora::ParticleFlowObject *const pPfo)
Whether a pfo is a neutrino or (antineutrino)
NodeVector m_children
The child nodes of this nodea.
bool m_foldDynamic
Whether or not to use process and topological information to make folding decisions.
const RecoHierarchy::NodeVector & GetUnmatchedReco() const
Retrieve the vector of unmatched reco nodes.
const pandora::PfoList & GetRecoParticles() const
Retrieve the PFOs associated with this node.
static bool AreTopologicallyContinuous(const pandora::MCParticle *const pMCParent, const pandora::MCParticle *const pMCChild, const float cosAngleTolerance)
Determine if two MC particles are topologically continuous within a given tolerance. If the parent does not travel any distance, a travelling parent is sought and the comparison made between this and the child. If no travelling parent can be found, the particles are treated as continuous.
ReconstructabilityCriteria m_recoCriteria
The criteria used to determine if the node is reconstructable.
void Print(const MCHierarchy &mcHierarchy) const
Prints information about which reco nodes are matched to the MC nodes, information about hit sharing...
bool m_foldToLeadingShowers
Whether or not to fold shower children to the leading shower particle.
const bool m_removeNeutrons
whether to remove neutrons and their downstream particles
Header file for the lar hierarchy helper class.
const NodeVector & GetInteractions(const pandora::MCParticle *pRoot) const
Retrieve the root nodes in this hierarchy.
const std::string ToString() const
Produce a string representation of the hierarchy.
const pandora::ParticleFlowObject * m_mainPfo
The leading particle flow object for this node.
std::vector< MCMatches > MCMatchesVector
pandora::CaloHitList m_caloHits
The list of calo hits of which this node is composed.
HitType
Definition: HitType.h:12
unsigned int GetNNeutrinoMCNodes(const pandora::MCParticle *const pRoot) const
Retrieve the number of neutrino interaction derived MC nodes available to match.
const MCHierarchy & m_mcHierarchy
The MC hierarchy for the matching procedure.
bool IsReconstructable(const pandora::MCParticle *pMCParticle) const
Checks if an individual particle meets reconstructability criteria.
void SetLeadingLepton()
Tags the particle as the leading lepton.
Node(MCHierarchy &hierarchy, const pandora::MCParticle *pMCParticle, const int tier=1)
Create a node with a primary MC particle.
const RecoHierarchy & m_recoHierarchy
The Reco hierarchy for the matching procedure.
static void GetAllDownstreamPfos(const pandora::PfoList &inputPfoList, pandora::PfoList &outputPfoList)
Get a flat list of all pfos, recursively, of all daughters associated with those pfos in an input lis...
void GetRootPfos(pandora::PfoList &rootPfos) const
Retrieve the root particle flow objects of the interaction hierarchies.
static void FillMCHierarchy(const pandora::MCParticleList &mcParticleList, const pandora::CaloHitList &caloHitList, const FoldingParameters &foldParameters, MCHierarchy &hierarchy)
Fill an MC hierarchy based on the specified folding criteria (see MCHierarchy::FillHierarchy for deta...
unsigned int GetNMCNodes(const pandora::MCParticle *const pRoot) const
Retrieve the number of MC nodes available to match.
int GetParticleId() const
Retrieve the PDG code for the leading particle in this node Note, for reco objects the PDG codes repr...
void GetFlattenedNodes(const pandora::ParticleFlowObject *const pRoot, NodeVector &nodeVector) const
Retrieve a flat vector of the nodes in the hierarchy.
void FillFlat(const pandora::MCParticle *pRoot)
Fill this node by folding all descendent particles to this node.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
int GetId() const
Retrieve the unique ID of this node.
int m_nextNodeId
The ID to use for the next node.
const unsigned int m_minHits
the minimum number of primary good Hits
RecoHierarchy::NodeVector m_unmatchedReco
The vector of unmatched reco nodes.
const unsigned int m_minGoodViews
the minimum number of primary good views
static InteractionDescriptor GetInteractionDescriptor(const pandora::MCParticleList &mcPrimaryList)
Get the interaction descriptor of an event.
void RegisterNode(const Node *pNode)
Register a node with the hierarchy.
static void GetAllDescendentMCParticles(const pandora::MCParticle *const pMCParticle, pandora::MCParticleList &descendentMCParticleList)
Get all descendent mc particles.
MCNodeVectorMap m_interactions
Map from incident particles (e.g. neutrino) to primaries.
float GetCompleteness(const RecoHierarchy::Node *pReco, const bool adcWeighted=false) const
Retrieve the completeness of the match.
void FillHierarchy(const pandora::MCParticle *pRoot, const FoldingParameters &foldParameters)
Recursively fill the hierarchy based on the criteria established for this MCHierarchy.
bool IsTestBeamParticle() const
Check if this is a test beam particle.
std::set< const pandora::ParticleFlowObject * > PfoSet
const float m_minCompleteness
The minimum completeness for a match to be considered good.
static bool IsElasticScatter(const pandora::MCParticle *const pMCParticle)
Check whether or not an MC particle came from an elastic scattering process.
static bool IsNeutrino(const pandora::MCParticle *const pMCParticle)
Whether a mc particle is a neutrino or antineutrino.
void FillFlat(const pandora::ParticleFlowObject *pRoot)
Fill this node by folding all descendent particles to this node.
std::map< const Node *, int > m_nodeToIdMap
A map from nodes to unique ids.
InteractionInfo m_matches
The map between an interaction and the vector of good matches from MC to reco.
unsigned int GetNCosmicRayMCNodes(const pandora::MCParticle *const pRoot) const
Retrieve the number of cosmic ray derived MC nodes available to match.
void AddRecoMatch(const RecoHierarchy::Node *pReco, const int nSharedHits)
Add a reconstructed node as a match for this MC node.