LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArMCParticleHelper.cc
Go to the documentation of this file.
1 
9 #include "Helpers/MCParticleHelper.h"
10 
11 #include "Objects/CaloHit.h"
12 #include "Objects/Cluster.h"
13 #include "Objects/MCParticle.h"
14 
15 #include "Pandora/PdgTable.h"
16 #include "Pandora/StatusCodes.h"
17 
22 
23 #include <algorithm>
24 #include <cstdlib>
25 
26 namespace lar_content
27 {
28 
29 using namespace pandora;
30 
32  m_minPrimaryGoodHits(15),
33  m_minHitsForGoodView(5),
34  m_minPrimaryGoodViews(2),
35  m_selectInputHits(true),
36  m_maxPhotonPropagation(2.5f),
37  m_minHitSharingFraction(0.9f),
38  m_foldBackHierarchy(true)
39 {
40 }
41 
42 //------------------------------------------------------------------------------------------------------------------------------------------
43 //------------------------------------------------------------------------------------------------------------------------------------------
44 
45 bool LArMCParticleHelper::DoesPrimaryMeetCriteria(const MCParticle *const pMCParticle, std::function<bool(const MCParticle *const)> fCriteria)
46 {
47  try
48  {
49  const MCParticle *const pPrimaryMCParticle = LArMCParticleHelper::GetPrimaryMCParticle(pMCParticle);
50  return fCriteria(pPrimaryMCParticle);
51  }
52  catch (const StatusCodeException &)
53  {
54  }
55 
56  return false;
57 }
58 
59 //------------------------------------------------------------------------------------------------------------------------------------------
60 
61 bool LArMCParticleHelper::DoesLeadingMeetCriteria(const MCParticle *const pMCParticle, std::function<bool(const MCParticle *const)> fCriteria)
62 {
63  try
64  {
65  const MCParticle *const pLeadingMCParticle = LArMCParticleHelper::GetLeadingMCParticle(pMCParticle);
66  return fCriteria(pLeadingMCParticle);
67  }
68  catch (const StatusCodeException &)
69  {
70  }
71 
72  return false;
73 }
74 
75 //------------------------------------------------------------------------------------------------------------------------------------------
76 
77 bool LArMCParticleHelper::IsBeamNeutrinoFinalState(const MCParticle *const pMCParticle)
78 {
79  const MCParticle *const pParentMCParticle(LArMCParticleHelper::GetParentMCParticle(pMCParticle));
80  return (LArMCParticleHelper::IsPrimary(pMCParticle) && LArMCParticleHelper::IsNeutrino(pParentMCParticle));
81 }
82 
83 //------------------------------------------------------------------------------------------------------------------------------------------
84 
85 bool LArMCParticleHelper::IsTriggeredBeamParticle(const MCParticle *const pMCParticle)
86 {
87  const int nuance(LArMCParticleHelper::GetNuanceCode(pMCParticle));
88  return (LArMCParticleHelper::IsPrimary(pMCParticle) && (nuance == 2001));
89 }
90 
91 //------------------------------------------------------------------------------------------------------------------------------------------
92 
93 bool LArMCParticleHelper::IsBeamParticle(const MCParticle *const pMCParticle)
94 {
95  const int nuance(LArMCParticleHelper::GetNuanceCode(pMCParticle));
96  return (LArMCParticleHelper::IsPrimary(pMCParticle) && ((nuance == 2000) || (nuance == 2001)));
97 }
98 
99 //------------------------------------------------------------------------------------------------------------------------------------------
100 
101 bool LArMCParticleHelper::IsLeadingBeamParticle(const MCParticle *const pMCParticle)
102 {
103  // ATTN: Only the parent triggered beam particle has nuance code 2001
105  return (LArMCParticleHelper::IsLeading(pMCParticle) && (parentNuance == 2001 || parentNuance == 2000));
106 }
107 
108 //------------------------------------------------------------------------------------------------------------------------------------------
109 
110 bool LArMCParticleHelper::IsCosmicRay(const MCParticle *const pMCParticle)
111 {
112  const int nuance(LArMCParticleHelper::GetNuanceCode(pMCParticle));
113  return (LArMCParticleHelper::IsPrimary(pMCParticle) &&
114  ((nuance == 3000) || ((nuance == 0) && !LArMCParticleHelper::IsBeamNeutrinoFinalState(pMCParticle))));
115 }
116 
117 //------------------------------------------------------------------------------------------------------------------------------------------
118 
119 unsigned int LArMCParticleHelper::GetNuanceCode(const MCParticle *const pMCParticle)
120 {
121  const LArMCParticle *const pLArMCParticle(dynamic_cast<const LArMCParticle *>(pMCParticle));
122  if (pLArMCParticle)
123  return pLArMCParticle->GetNuanceCode();
124 
125  std::cout << "LArMCParticleHelper::GetNuanceCode - Error: Can't cast to LArMCParticle" << std::endl;
126  throw StatusCodeException(STATUS_CODE_NOT_ALLOWED);
127 }
128 
129 //------------------------------------------------------------------------------------------------------------------------------------------
130 
131 bool LArMCParticleHelper::IsNeutrino(const MCParticle *const pMCParticle)
132 {
133  const int nuance(LArMCParticleHelper::GetNuanceCode(pMCParticle));
134  if ((nuance == 0) || (nuance == 2000) || (nuance == 2001) || (nuance == 3000))
135  return false;
136 
137  const int absoluteParticleId(std::abs(pMCParticle->GetParticleId()));
138  return ((NU_E == absoluteParticleId) || (NU_MU == absoluteParticleId) || (NU_TAU == absoluteParticleId));
139 }
140 
141 //------------------------------------------------------------------------------------------------------------------------------------------
142 
143 bool LArMCParticleHelper::IsPrimary(const pandora::MCParticle *const pMCParticle)
144 {
145  try
146  {
147  return (LArMCParticleHelper::GetPrimaryMCParticle(pMCParticle) == pMCParticle);
148  }
149  catch (const StatusCodeException &)
150  {
151  }
152 
153  return false;
154 }
155 
156 //------------------------------------------------------------------------------------------------------------------------------------------
157 
158 bool LArMCParticleHelper::IsLeading(const pandora::MCParticle *const pMCParticle)
159 {
160  try
161  {
162  return (LArMCParticleHelper::GetLeadingMCParticle(pMCParticle) == pMCParticle);
163  }
164  catch (const StatusCodeException &)
165  {
166  }
167 
168  return false;
169 }
170 
171 //------------------------------------------------------------------------------------------------------------------------------------------
172 
173 int LArMCParticleHelper::GetHierarchyTier(const pandora::MCParticle *const pMCParticle)
174 {
175  const MCParticle *pParentMCParticle = pMCParticle;
176  int tier(0);
177 
178  while (pParentMCParticle->GetParentList().empty() == false)
179  {
180  if (1 != pParentMCParticle->GetParentList().size())
181  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
182 
183  pParentMCParticle = *(pParentMCParticle->GetParentList().begin());
184  ++tier;
185  }
186 
187  return tier;
188 }
189 
190 //------------------------------------------------------------------------------------------------------------------------------------------
191 
192 bool LArMCParticleHelper::IsVisible(const MCParticle *const pMCParticle)
193 {
194  const int absoluteParticleId(std::abs(pMCParticle->GetParticleId()));
195 
196  if ((E_MINUS == absoluteParticleId) || (MU_MINUS == absoluteParticleId) || (PI_PLUS == absoluteParticleId) || (K_PLUS == absoluteParticleId) ||
197  (SIGMA_MINUS == absoluteParticleId) || (SIGMA_PLUS == absoluteParticleId) || (HYPERON_MINUS == absoluteParticleId) ||
198  (PROTON == absoluteParticleId) || (PHOTON == absoluteParticleId) || (NEUTRON == absoluteParticleId))
199  return true;
200 
201  return false;
202 }
203 
204 //------------------------------------------------------------------------------------------------------------------------------------------
205 
206 void LArMCParticleHelper::GetTrueNeutrinos(const MCParticleList *const pMCParticleList, MCParticleVector &trueNeutrinos)
207 {
208  for (const MCParticle *const pMCParticle : *pMCParticleList)
209  {
210  if (LArMCParticleHelper::IsNeutrino(pMCParticle))
211  trueNeutrinos.push_back(pMCParticle);
212  }
213 
214  std::sort(trueNeutrinos.begin(), trueNeutrinos.end(), LArMCParticleHelper::SortByMomentum);
215 }
216 
217 //------------------------------------------------------------------------------------------------------------------------------------------
218 
219 void LArMCParticleHelper::GetTrueTestBeamParticles(const MCParticleList *const pMCParticleList, MCParticleVector &trueTestBeamParticles)
220 {
221  for (const MCParticle *const pMCParticle : *pMCParticleList)
222  {
224  trueTestBeamParticles.push_back(pMCParticle);
225  }
226 
227  std::sort(trueTestBeamParticles.begin(), trueTestBeamParticles.end(), LArMCParticleHelper::SortByMomentum);
228 }
229 
230 //------------------------------------------------------------------------------------------------------------------------------------------
231 
232 const MCParticle *LArMCParticleHelper::GetParentMCParticle(const MCParticle *const pMCParticle)
233 {
234  const MCParticle *pParentMCParticle = pMCParticle;
235 
236  while (pParentMCParticle->GetParentList().empty() == false)
237  {
238  if (1 != pParentMCParticle->GetParentList().size())
239  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
240 
241  pParentMCParticle = *(pParentMCParticle->GetParentList().begin());
242  }
243 
244  return pParentMCParticle;
245 }
246 
247 //------------------------------------------------------------------------------------------------------------------------------------------
248 
249 void LArMCParticleHelper::GetAllDescendentMCParticles(const pandora::MCParticle *const pMCParticle, pandora::MCParticleList &descendentMCParticleList)
250 {
251  for (const MCParticle *pDaughterMCParticle : pMCParticle->GetDaughterList())
252  {
253  if (std::find(descendentMCParticleList.begin(), descendentMCParticleList.end(), pDaughterMCParticle) == descendentMCParticleList.end())
254  {
255  descendentMCParticleList.emplace_back(pDaughterMCParticle);
256  LArMCParticleHelper::GetAllDescendentMCParticles(pDaughterMCParticle, descendentMCParticleList);
257  }
258  }
259 }
260 
261 //------------------------------------------------------------------------------------------------------------------------------------------
262 
263 void LArMCParticleHelper::GetAllDescendentMCParticles(const MCParticle *const pMCParticle, MCParticleList &descendentTrackParticles,
264  MCParticleList &leadingShowerParticles, MCParticleList &leadingNeutrons)
265 {
266  for (const MCParticle *pDaughterMCParticle : pMCParticle->GetDaughterList())
267  {
268  if (std::find(descendentTrackParticles.begin(), descendentTrackParticles.end(), pDaughterMCParticle) == descendentTrackParticles.end())
269  {
270  const int pdg{std::abs(pDaughterMCParticle->GetParticleId())};
271  if (pdg == E_MINUS || pdg == PHOTON)
272  {
273  leadingShowerParticles.emplace_back(pDaughterMCParticle);
274  }
275  else if (pdg == NEUTRON)
276  {
277  leadingNeutrons.emplace_back(pDaughterMCParticle);
278  }
279  else
280  {
281  descendentTrackParticles.emplace_back(pDaughterMCParticle);
282  LArMCParticleHelper::GetAllDescendentMCParticles(pDaughterMCParticle, descendentTrackParticles, leadingShowerParticles, leadingNeutrons);
283  }
284  }
285  }
286 }
287 
288 //------------------------------------------------------------------------------------------------------------------------------------------
289 
290 void LArMCParticleHelper::GetAllAncestorMCParticles(const pandora::MCParticle *const pMCParticle, pandora::MCParticleList &ancestorMCParticleList)
291 {
292  const MCParticleList &parentMCParticleList = pMCParticle->GetParentList();
293  if (parentMCParticleList.empty())
294  return;
295  if (parentMCParticleList.size() != 1)
296  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
297 
298  const MCParticle *pParentMCParticle = *parentMCParticleList.begin();
299  if (std::find(ancestorMCParticleList.begin(), ancestorMCParticleList.end(), pParentMCParticle) == ancestorMCParticleList.end())
300  {
301  ancestorMCParticleList.push_back(pParentMCParticle);
302  LArMCParticleHelper::GetAllAncestorMCParticles(pParentMCParticle, ancestorMCParticleList);
303  }
304 }
305 
306 //------------------------------------------------------------------------------------------------------------------------------------------
307 
308 void LArMCParticleHelper::GetPrimaryMCParticleList(const MCParticleList *const pMCParticleList, MCParticleVector &mcPrimaryVector)
309 {
310  for (const MCParticle *const pMCParticle : *pMCParticleList)
311  {
312  if (LArMCParticleHelper::IsPrimary(pMCParticle))
313  mcPrimaryVector.push_back(pMCParticle);
314  }
315 
316  std::sort(mcPrimaryVector.begin(), mcPrimaryVector.end(), LArMCParticleHelper::SortByMomentum);
317 }
318 
319 //------------------------------------------------------------------------------------------------------------------------------------------
320 
321 void LArMCParticleHelper::GetLeadingMCParticleList(const MCParticleList *const pMCParticleList, MCParticleVector &mcLeadingVector)
322 {
323  for (const MCParticle *const pMCParticle : *pMCParticleList)
324  {
325  const bool isBeamParticle(LArMCParticleHelper::IsBeamParticle(LArMCParticleHelper::GetParentMCParticle(pMCParticle)));
326 
327  if ((isBeamParticle && LArMCParticleHelper::IsLeading(pMCParticle)) || (!isBeamParticle && LArMCParticleHelper::IsPrimary(pMCParticle)))
328  {
329  mcLeadingVector.push_back(pMCParticle);
330  }
331  }
332 
333  std::sort(mcLeadingVector.begin(), mcLeadingVector.end(), LArMCParticleHelper::SortByMomentum);
334 }
335 
336 //------------------------------------------------------------------------------------------------------------------------------------------
337 
338 void LArMCParticleHelper::GetFirstVisibleMCParticles(const MCParticle *const pRoot, MCParticleList &visibleParticleList)
339 {
341  {
342  visibleParticleList.emplace_back(pRoot);
343  }
344  else
345  {
346  for (const MCParticle *const pMCParticle : pRoot->GetDaughterList())
347  LArMCParticleHelper::GetFirstVisibleMCParticles(pMCParticle, visibleParticleList);
348  }
349 }
350 
351 //------------------------------------------------------------------------------------------------------------------------------------------
352 
353 const MCParticle *LArMCParticleHelper::GetPrimaryMCParticle(const MCParticle *const pMCParticle)
354 {
355  // Navigate upward through MC daughter/parent links - collect this particle and all its parents
356  MCParticleVector mcVector;
357 
358  const MCParticle *pParentMCParticle = pMCParticle;
359  mcVector.push_back(pParentMCParticle);
360 
361  while (!pParentMCParticle->GetParentList().empty())
362  {
363  if (1 != pParentMCParticle->GetParentList().size())
364  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
365 
366  pParentMCParticle = *(pParentMCParticle->GetParentList().begin());
367  mcVector.push_back(pParentMCParticle);
368  }
369 
370  // Navigate downward through MC parent/daughter links - return the first long-lived charged particle
371  for (MCParticleVector::const_reverse_iterator iter = mcVector.rbegin(), iterEnd = mcVector.rend(); iter != iterEnd; ++iter)
372  {
373  const MCParticle *const pNextParticle = *iter;
374 
375  if (LArMCParticleHelper::IsVisible(pNextParticle))
376  return pNextParticle;
377  }
378 
379  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
380 }
381 
382 //------------------------------------------------------------------------------------------------------------------------------------------
383 
384 const MCParticle *LArMCParticleHelper::GetLeadingMCParticle(const MCParticle *const pMCParticle, const int hierarchyTierLimit)
385 {
386  // ATTN: If not beam particle return primary particle
387  const bool isBeamParticle(LArMCParticleHelper::IsBeamParticle(LArMCParticleHelper::GetParentMCParticle(pMCParticle)));
388 
389  if (!isBeamParticle)
390  return LArMCParticleHelper::GetPrimaryMCParticle(pMCParticle);
391 
392  // Navigate upward through MC daughter/parent links - collect this particle and all its parents
393  MCParticleVector mcVector;
394 
395  const MCParticle *pParentMCParticle = pMCParticle;
396  mcVector.push_back(pParentMCParticle);
397 
398  while (!pParentMCParticle->GetParentList().empty())
399  {
400  if (1 != pParentMCParticle->GetParentList().size())
401  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
402 
403  pParentMCParticle = *(pParentMCParticle->GetParentList().begin());
404  mcVector.push_back(pParentMCParticle);
405  }
406 
407  int hierarchyTier(0);
408  const MCParticle *pLeadingMCParticle(nullptr);
409 
410  // Navigate downward through MC parent/daughter links - return the first long-lived charged particle
411  for (MCParticleVector::const_reverse_iterator iter = mcVector.rbegin(), iterEnd = mcVector.rend(); iter != iterEnd; ++iter)
412  {
413  const MCParticle *const pNextParticle = *iter;
414 
415  // ATTN: Squash any invisible particles (e.g. pizero)
416  if (!LArMCParticleHelper::IsVisible(pNextParticle))
417  continue;
418 
419  if (hierarchyTier <= hierarchyTierLimit)
420  pLeadingMCParticle = pNextParticle;
421 
422  hierarchyTier++;
423  }
424 
425  if (!pLeadingMCParticle)
426  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
427 
428  return pLeadingMCParticle;
429 }
430 
431 //------------------------------------------------------------------------------------------------------------------------------------------
432 
433 void LArMCParticleHelper::GetMCPrimaryMap(const MCParticleList *const pMCParticleList, MCRelationMap &mcPrimaryMap)
434 {
435  for (const MCParticle *const pMCParticle : *pMCParticleList)
436  {
437  try
438  {
439  const MCParticle *const pPrimaryMCParticle = LArMCParticleHelper::GetPrimaryMCParticle(pMCParticle);
440  mcPrimaryMap[pMCParticle] = pPrimaryMCParticle;
441  }
442  catch (const StatusCodeException &)
443  {
444  }
445  }
446 }
447 
448 //------------------------------------------------------------------------------------------------------------------------------------------
449 
450 void LArMCParticleHelper::GetMCLeadingMap(const MCParticleList *const pMCParticleList, MCRelationMap &mcLeadingMap)
451 {
452  for (const MCParticle *const pMCParticle : *pMCParticleList)
453  {
454  try
455  {
456  const MCParticle *const pLeadingMCParticle = LArMCParticleHelper::GetLeadingMCParticle(pMCParticle);
457  mcLeadingMap[pMCParticle] = pLeadingMCParticle;
458  }
459  catch (const StatusCodeException &)
460  {
461  }
462  }
463 }
464 
465 //------------------------------------------------------------------------------------------------------------------------------------------
466 
467 void LArMCParticleHelper::GetMCToSelfMap(const MCParticleList *const pMCParticleList, MCRelationMap &mcToSelfMap)
468 {
469  for (const MCParticle *const pMCParticle : *pMCParticleList)
470  {
471  mcToSelfMap[pMCParticle] = pMCParticle;
472  }
473 }
474 
475 //------------------------------------------------------------------------------------------------------------------------------------------
476 
477 const MCParticle *LArMCParticleHelper::GetMainMCParticle(const ParticleFlowObject *const pPfo)
478 {
479  ClusterList clusterList;
480  LArPfoHelper::GetTwoDClusterList(pPfo, clusterList);
481  return MCParticleHelper::GetMainMCParticle(&clusterList);
482 }
483 
484 //------------------------------------------------------------------------------------------------------------------------------------------
485 
486 bool LArMCParticleHelper::SortByMomentum(const MCParticle *const pLhs, const MCParticle *const pRhs)
487 {
488  // Sort by momentum (prefer higher momentum)
489  const float momentumLhs(pLhs->GetMomentum().GetMagnitudeSquared());
490  const float momentumRhs(pRhs->GetMomentum().GetMagnitudeSquared());
491 
492  if (std::fabs(momentumLhs - momentumRhs) > std::numeric_limits<float>::epsilon())
493  return (momentumLhs > momentumRhs);
494 
495  // Sort by energy (prefer lighter particles)
496  if (std::fabs(pLhs->GetEnergy() - pRhs->GetEnergy()) > std::numeric_limits<float>::epsilon())
497  return (pLhs->GetEnergy() < pRhs->GetEnergy());
498 
499  // Sort by PDG code (prefer smaller numbers)
500  if (pLhs->GetParticleId() != pRhs->GetParticleId())
501  return (pLhs->GetParticleId() < pRhs->GetParticleId());
502 
503  // Sort by vertex position (tie-breaker)
504  const float positionLhs(pLhs->GetVertex().GetMagnitudeSquared());
505  const float positionRhs(pRhs->GetVertex().GetMagnitudeSquared());
506 
507  return (positionLhs < positionRhs);
508 }
509 
510 //------------------------------------------------------------------------------------------------------------------------------------------
511 
512 void LArMCParticleHelper::GetMCParticleToCaloHitMatches(const CaloHitList *const pCaloHitList, const MCRelationMap &mcToTargetMCMap,
513  CaloHitToMCMap &hitToMCMap, MCContributionMap &mcToTrueHitListMap)
514 {
515  for (const CaloHit *const pCaloHit : *pCaloHitList)
516  {
517  try
518  {
519  const MCParticle *const pHitParticle(MCParticleHelper::GetMainMCParticle(pCaloHit));
520  const MCParticle *pTargetParticle(pHitParticle);
521 
522  // ATTN Do not map back to target if mc to primary mc map or mc to self map not provided
523  if (!mcToTargetMCMap.empty())
524  {
525  MCRelationMap::const_iterator mcIter = mcToTargetMCMap.find(pHitParticle);
526 
527  if (mcToTargetMCMap.end() == mcIter)
528  continue;
529 
530  pTargetParticle = mcIter->second;
531  }
532 
533  mcToTrueHitListMap[pTargetParticle].push_back(pCaloHit);
534  hitToMCMap[pCaloHit] = pTargetParticle;
535  }
536  catch (StatusCodeException &statusCodeException)
537  {
538  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
539  throw statusCodeException;
540  }
541  }
542 }
543 
544 //------------------------------------------------------------------------------------------------------------------------------------------
545 
546 void LArMCParticleHelper::SelectReconstructableMCParticles(const MCParticleList *pMCParticleList, const CaloHitList *pCaloHitList,
547  const PrimaryParameters &parameters, std::function<bool(const MCParticle *const)> fCriteria, MCContributionMap &selectedMCParticlesToHitsMap)
548 {
549  // Obtain map: [mc particle -> target mc particle]
550  LArMCParticleHelper::MCRelationMap mcToTargetMCMap;
551  parameters.m_foldBackHierarchy ? LArMCParticleHelper::GetMCPrimaryMap(pMCParticleList, mcToTargetMCMap)
552  : LArMCParticleHelper::GetMCToSelfMap(pMCParticleList, mcToTargetMCMap);
553 
554  // Remove non-reconstructable hits, e.g. those downstream of a neutron
555  // Unless selectInputHits == false
556  CaloHitList selectedCaloHitList;
557  LArMCParticleHelper::SelectCaloHits(pCaloHitList, mcToTargetMCMap, selectedCaloHitList, parameters.m_selectInputHits, parameters.m_maxPhotonPropagation);
558 
559  // Obtain maps: [hit -> target mc particle], [target mc particle -> list of hits]
560  CaloHitToMCMap trueHitToTargetMCMap;
561  MCContributionMap targetMCToTrueHitListMap;
562  LArMCParticleHelper::GetMCParticleToCaloHitMatches(&selectedCaloHitList, mcToTargetMCMap, trueHitToTargetMCMap, targetMCToTrueHitListMap);
563 
564  // Obtain vector: target mc particles
565  MCParticleVector targetMCVector;
566  if (parameters.m_foldBackHierarchy)
567  {
568  LArMCParticleHelper::GetPrimaryMCParticleList(pMCParticleList, targetMCVector);
569  }
570  else
571  {
572  std::copy(pMCParticleList->begin(), pMCParticleList->end(), std::back_inserter(targetMCVector));
573  }
574 
575  // Select MCParticles matching criteria
576  MCParticleVector candidateTargets;
577  LArMCParticleHelper::SelectParticlesMatchingCriteria(targetMCVector, fCriteria, candidateTargets, parameters, false);
578 
579  // Ensure the MCParticles have enough "good" hits to be reconstructed
580  LArMCParticleHelper::SelectParticlesByHitCount(candidateTargets, targetMCToTrueHitListMap, mcToTargetMCMap, parameters, selectedMCParticlesToHitsMap);
581 }
582 
583 //------------------------------------------------------------------------------------------------------------------------------------------
584 
585 void LArMCParticleHelper::SelectReconstructableTestBeamHierarchyMCParticles(const MCParticleList *pMCParticleList, const CaloHitList *pCaloHitList,
586  const PrimaryParameters &parameters, std::function<bool(const MCParticle *const)> fCriteria, MCContributionMap &selectedMCParticlesToHitsMap)
587 {
588  // Obtain map: [mc particle -> target mc particle]
589  LArMCParticleHelper::MCRelationMap mcToTargetMCMap;
590  parameters.m_foldBackHierarchy ? LArMCParticleHelper::GetMCLeadingMap(pMCParticleList, mcToTargetMCMap)
591  : LArMCParticleHelper::GetMCToSelfMap(pMCParticleList, mcToTargetMCMap);
592 
593  // Remove non-reconstructable hits, e.g. those downstream of a neutron
594  // Unless selectInputHits == false
595  CaloHitList selectedCaloHitList;
596  LArMCParticleHelper::SelectCaloHits(pCaloHitList, mcToTargetMCMap, selectedCaloHitList, parameters.m_selectInputHits, parameters.m_maxPhotonPropagation);
597 
598  // Obtain maps: [hit -> target mc particle], [target mc particle -> list of hits]
599  CaloHitToMCMap trueHitToTargetMCMap;
600  MCContributionMap targetMCToTrueHitListMap;
601  LArMCParticleHelper::GetMCParticleToCaloHitMatches(&selectedCaloHitList, mcToTargetMCMap, trueHitToTargetMCMap, targetMCToTrueHitListMap);
602 
603  // Obtain vector: target mc particles
604  MCParticleVector targetMCVector;
605  if (parameters.m_foldBackHierarchy)
606  {
607  LArMCParticleHelper::GetLeadingMCParticleList(pMCParticleList, targetMCVector);
608  }
609  else
610  {
611  std::copy(pMCParticleList->begin(), pMCParticleList->end(), std::back_inserter(targetMCVector));
612  }
613 
614  // Select MCParticles matching criteria
615  MCParticleVector candidateTargets;
616  LArMCParticleHelper::SelectParticlesMatchingCriteria(targetMCVector, fCriteria, candidateTargets, parameters, true);
617 
618  // Ensure the MCParticles have enough "good" hits to be reconstructed
619  LArMCParticleHelper::SelectParticlesByHitCount(candidateTargets, targetMCToTrueHitListMap, mcToTargetMCMap, parameters, selectedMCParticlesToHitsMap);
620 
621  // ATTN: The parent particle must be in the hierarchy map, event if not reconstructable
622  for (const MCParticle *const pMCParticle : candidateTargets)
623  {
624  if (!LArMCParticleHelper::IsBeamParticle(pMCParticle))
625  continue;
626 
627  const MCParticle *const pParentMCParticle(LArMCParticleHelper::GetParentMCParticle(pMCParticle));
628  if (selectedMCParticlesToHitsMap.find(pParentMCParticle) == selectedMCParticlesToHitsMap.end())
629  {
630  CaloHitList caloHitList;
631  selectedMCParticlesToHitsMap.insert(MCContributionMap::value_type(pParentMCParticle, caloHitList));
632  }
633  }
634 }
635 
636 //------------------------------------------------------------------------------------------------------------------------------------------
637 
638 void LArMCParticleHelper::GetPfoToReconstructable2DHitsMap(const PfoList &pfoList, const MCContributionMap &selectedMCParticleToHitsMap,
639  PfoContributionMap &pfoToReconstructable2DHitsMap, const bool foldBackHierarchy)
640 {
642  pfoList, MCContributionMapVector({selectedMCParticleToHitsMap}), pfoToReconstructable2DHitsMap, foldBackHierarchy);
643 }
644 
645 //------------------------------------------------------------------------------------------------------------------------------------------
646 
648  const MCContributionMap &selectedMCParticleToHitsMap, PfoContributionMap &pfoToReconstructable2DHitsMap, const bool foldBackHierarchy)
649 {
651  pfoList, MCContributionMapVector({selectedMCParticleToHitsMap}), pfoToReconstructable2DHitsMap, foldBackHierarchy);
652 }
653 
654 //------------------------------------------------------------------------------------------------------------------------------------------
655 
656 void LArMCParticleHelper::GetPfoToReconstructable2DHitsMap(const PfoList &pfoList, const MCContributionMapVector &selectedMCParticleToHitsMaps,
657  PfoContributionMap &pfoToReconstructable2DHitsMap, const bool foldBackHierarchy)
658 {
659  for (const ParticleFlowObject *const pPfo : pfoList)
660  {
661  CaloHitList pfoHitList;
662  LArMCParticleHelper::CollectReconstructable2DHits(pPfo, selectedMCParticleToHitsMaps, pfoHitList, foldBackHierarchy);
663 
664  if (!pfoToReconstructable2DHitsMap.insert(PfoContributionMap::value_type(pPfo, pfoHitList)).second)
665  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
666  }
667 }
668 
669 //------------------------------------------------------------------------------------------------------------------------------------------
670 
672  const MCContributionMapVector &selectedMCParticleToHitsMaps, PfoContributionMap &pfoToReconstructable2DHitsMap, const bool foldBackHierarchy)
673 {
674  for (const ParticleFlowObject *const pPfo : pfoList)
675  {
676  CaloHitList pfoHitList;
677  LArMCParticleHelper::CollectReconstructableTestBeamHierarchy2DHits(pPfo, selectedMCParticleToHitsMaps, pfoHitList, foldBackHierarchy);
678 
679  if (!pfoToReconstructable2DHitsMap.insert(PfoContributionMap::value_type(pPfo, pfoHitList)).second)
680  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
681  }
682 }
683 
684 //------------------------------------------------------------------------------------------------------------------------------------------
685 
687  const MCContributionMapVector &selectedMCParticleToHitsMaps, PfoToMCParticleHitSharingMap &pfoToMCParticleHitSharingMap,
688  MCParticleToPfoHitSharingMap &mcParticleToPfoHitSharingMap)
689 {
690  PfoVector sortedPfos;
691  for (const auto &mapEntry : pfoToReconstructable2DHitsMap)
692  sortedPfos.push_back(mapEntry.first);
693  std::sort(sortedPfos.begin(), sortedPfos.end(), LArPfoHelper::SortByNHits);
694 
695  for (const ParticleFlowObject *const pPfo : sortedPfos)
696  {
697  for (const MCContributionMap &mcParticleToHitsMap : selectedMCParticleToHitsMaps)
698  {
699  MCParticleVector sortedMCParticles;
700  for (const auto &mapEntry : mcParticleToHitsMap)
701  sortedMCParticles.push_back(mapEntry.first);
702  std::sort(sortedMCParticles.begin(), sortedMCParticles.end(), PointerLessThan<MCParticle>());
703 
704  for (const MCParticle *const pMCParticle : sortedMCParticles)
705  {
706  // Add map entries for this Pfo & MCParticle if required
707  if (pfoToMCParticleHitSharingMap.find(pPfo) == pfoToMCParticleHitSharingMap.end())
708  if (!pfoToMCParticleHitSharingMap.insert(PfoToMCParticleHitSharingMap::value_type(pPfo, MCParticleToSharedHitsVector())).second)
709  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT); // ATTN maybe overkill
710 
711  if (mcParticleToPfoHitSharingMap.find(pMCParticle) == mcParticleToPfoHitSharingMap.end())
712  if (!mcParticleToPfoHitSharingMap.insert(MCParticleToPfoHitSharingMap::value_type(pMCParticle, PfoToSharedHitsVector())).second)
713  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
714 
715  // Check this Pfo & MCParticle pairing hasn't already been checked
716  MCParticleToSharedHitsVector &mcHitPairs(pfoToMCParticleHitSharingMap.at(pPfo));
717  PfoToSharedHitsVector &pfoHitPairs(mcParticleToPfoHitSharingMap.at(pMCParticle));
718 
719  if (std::any_of(mcHitPairs.begin(), mcHitPairs.end(),
720  [&](const MCParticleCaloHitListPair &pair) { return (pair.first == pMCParticle); }))
721  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
722 
723  if (std::any_of(pfoHitPairs.begin(), pfoHitPairs.end(), [&](const PfoCaloHitListPair &pair) { return (pair.first == pPfo); }))
724  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
725 
726  // Add records to maps if there are any shared hits
727  const CaloHitList sharedHits(
728  LArMCParticleHelper::GetSharedHits(pfoToReconstructable2DHitsMap.at(pPfo), mcParticleToHitsMap.at(pMCParticle)));
729 
730  if (!sharedHits.empty())
731  {
732  mcHitPairs.push_back(MCParticleCaloHitListPair(pMCParticle, sharedHits));
733  pfoHitPairs.push_back(PfoCaloHitListPair(pPfo, sharedHits));
734 
735  std::sort(mcHitPairs.begin(), mcHitPairs.end(),
736  [](const MCParticleCaloHitListPair &a, const MCParticleCaloHitListPair &b) -> bool
737  {
738  return ((a.second.size() != b.second.size()) ? a.second.size() > b.second.size()
739  : LArMCParticleHelper::SortByMomentum(a.first, b.first));
740  });
741 
742  std::sort(pfoHitPairs.begin(), pfoHitPairs.end(),
743  [](const PfoCaloHitListPair &a, const PfoCaloHitListPair &b) -> bool {
744  return ((a.second.size() != b.second.size()) ? a.second.size() > b.second.size()
745  : LArPfoHelper::SortByNHits(a.first, b.first));
746  });
747  }
748  }
749  }
750  }
751 }
752 
753 //------------------------------------------------------------------------------------------------------------------------------------------
754 
755 void LArMCParticleHelper::GetClusterToReconstructable2DHitsMap(const pandora::ClusterList &clusterList,
756  const MCContributionMap &selectedMCToHitsMap, ClusterContributionMap &clusterToReconstructable2DHitsMap)
757 {
758  LArMCParticleHelper::GetClusterToReconstructable2DHitsMap(clusterList, MCContributionMapVector({selectedMCToHitsMap}), clusterToReconstructable2DHitsMap);
759 }
760 
761 //------------------------------------------------------------------------------------------------------------------------------------------
762 
763 void LArMCParticleHelper::GetClusterToReconstructable2DHitsMap(const pandora::ClusterList &clusterList,
764  const MCContributionMapVector &selectedMCToHitsMaps, ClusterContributionMap &clusterToReconstructable2DHitsMap)
765 {
766  for (const Cluster *const pCluster : clusterList)
767  {
768  CaloHitList caloHitList;
769  LArMCParticleHelper::CollectReconstructable2DHits(pCluster, selectedMCToHitsMaps, caloHitList);
770 
771  if (!clusterToReconstructable2DHitsMap.insert(ClusterContributionMap::value_type(pCluster, caloHitList)).second)
772  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
773  }
774 }
775 
776 //------------------------------------------------------------------------------------------------------------------------------------------
777 
778 bool LArMCParticleHelper::IsBremsstrahlung(const MCParticle *const pMCParticle)
779 {
780  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
781  if (!pLArMCParticle)
782  return false;
783 
784  switch (pLArMCParticle->GetProcess())
785  {
786  case MC_PROC_E_BREM:
787  case MC_PROC_MU_BREM:
788  case MC_PROC_HAD_BREM:
789  return true;
790  default:
791  return false;
792  }
793 
794  return false;
795 }
796 
797 //------------------------------------------------------------------------------------------------------------------------------------------
798 
799 bool LArMCParticleHelper::IsCapture(const MCParticle *const pMCParticle)
800 {
801  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
802  if (!pLArMCParticle)
803  return false;
804 
805  switch (pLArMCParticle->GetProcess())
806  {
808  case MC_PROC_N_CAPTURE:
812  return true;
813  default:
814  return false;
815  }
816 
817  return false;
818 }
819 
820 //------------------------------------------------------------------------------------------------------------------------------------------
821 
822 bool LArMCParticleHelper::IsDecay(const MCParticle *const pMCParticle)
823 {
824  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
825  if (!pLArMCParticle)
826  return false;
827 
828  switch (pLArMCParticle->GetProcess())
829  {
830  case MC_PROC_DECAY:
831  return true;
832  default:
833  return false;
834  }
835 
836  return false;
837 }
838 
839 //------------------------------------------------------------------------------------------------------------------------------------------
840 
841 bool LArMCParticleHelper::IsElasticScatter(const MCParticle *const pMCParticle)
842 {
843  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
844  if (!pLArMCParticle)
845  return false;
846 
847  switch (pLArMCParticle->GetProcess())
848  {
851  case MC_PROC_HAD_ELASTIC:
852  case MC_PROC_RAYLEIGH:
853  return true;
854  default:
855  return false;
856  }
857 
858  return false;
859 }
860 
861 //------------------------------------------------------------------------------------------------------------------------------------------
862 
863 bool LArMCParticleHelper::IsInelasticScatter(const MCParticle *const pMCParticle)
864 {
865  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
866  if (!pLArMCParticle)
867  return false;
868 
869  switch (pLArMCParticle->GetProcess())
870  {
871  case MC_PROC_COMPT:
890  return true;
891  default:
892  return false;
893  }
894 
895  return false;
896 }
897 
898 //------------------------------------------------------------------------------------------------------------------------------------------
899 
900 bool LArMCParticleHelper::IsIonisation(const MCParticle *const pMCParticle)
901 {
902  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
903  if (!pLArMCParticle)
904  return false;
905 
906  switch (pLArMCParticle->GetProcess())
907  {
908  case MC_PROC_E_IONI:
909  case MC_PROC_MU_IONI:
910  case MC_PROC_HAD_IONI:
911  case MC_PROC_ION_IONI:
912  return true;
913  default:
914  return false;
915  }
916 
917  return false;
918 }
919 
920 //------------------------------------------------------------------------------------------------------------------------------------------
921 
922 bool LArMCParticleHelper::IsNuclear(const MCParticle *const pMCParticle)
923 {
924  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
925  if (!pLArMCParticle)
926  return false;
927 
928  switch (pLArMCParticle->GetProcess())
929  {
932  case MC_PROC_MU_NUCLEAR:
933  return true;
934  default:
935  return false;
936  }
937 
938  return false;
939 }
940 
941 //------------------------------------------------------------------------------------------------------------------------------------------
942 
943 bool LArMCParticleHelper::IsPairProduction(const MCParticle *const pMCParticle)
944 {
945  const LArMCParticle *pLArMCParticle{dynamic_cast<const LArMCParticle *>(pMCParticle)};
946  if (!pLArMCParticle)
947  return false;
948 
949  switch (pLArMCParticle->GetProcess())
950  {
953  return true;
954  default:
955  return false;
956  }
957 
958  return false;
959 }
960 
961 //------------------------------------------------------------------------------------------------------------------------------------------
962 
963 CaloHitList LArMCParticleHelper::GetSharedHits(const CaloHitList &hitListA, const CaloHitList &hitListB)
964 {
965  CaloHitList sharedHits;
966 
967  for (const CaloHit *const pCaloHit : hitListA)
968  {
969  if (std::find(hitListB.begin(), hitListB.end(), pCaloHit) != hitListB.end())
970  sharedHits.push_back(pCaloHit);
971  }
972 
973  return sharedHits;
974 }
975 
976 //------------------------------------------------------------------------------------------------------------------------------------------
977 
978 bool LArMCParticleHelper::AreTopologicallyContinuous(const MCParticle *const pMCParent, const MCParticle *const pMCChild, const float cosAngleTolerance)
979 {
980  CartesianVector childDirection{pMCChild->GetEndpoint() - pMCChild->GetVertex()};
981  if (childDirection.GetMagnitude() < std::numeric_limits<float>::epsilon())
982  return true;
983  childDirection = childDirection.GetUnitVector();
984 
985  const MCParticle *pMCUpstream{pMCParent};
986  while (true)
987  {
988  CartesianVector parentDirection{pMCUpstream->GetEndpoint() - pMCUpstream->GetVertex()};
989  if (parentDirection.GetMagnitude() > std::numeric_limits<float>::epsilon())
990  {
991  parentDirection = parentDirection.GetUnitVector();
992  return parentDirection.GetDotProduct(childDirection) >= cosAngleTolerance;
993  }
994  else
995  {
996  const MCParticleList &parentList{pMCUpstream->GetParentList()};
997  const size_t size{parentList.size()};
998  if (size == 1)
999  pMCUpstream = parentList.front();
1000  else if (size == 0)
1001  return true;
1002  else
1003  return false;
1004  }
1005  }
1006 
1007  return false;
1008 }
1009 
1010 // private
1011 //------------------------------------------------------------------------------------------------------------------------------------------
1012 
1013 void LArMCParticleHelper::CollectReconstructable2DHits(const ParticleFlowObject *const pPfo,
1014  const MCContributionMapVector &selectedMCParticleToHitsMaps, CaloHitList &reconstructableCaloHitList2D, const bool foldBackHierarchy)
1015 {
1016 
1017  PfoList pfoList;
1018 
1019  // If foldBackHierarchy collect all 2D calo hits in pfo hierarchy
1020  // else collect hits directly belonging to pfo
1021  if (foldBackHierarchy)
1022  {
1023  LArPfoHelper::GetAllDownstreamPfos(pPfo, pfoList);
1024  }
1025  else
1026  {
1027  pfoList.push_back(pPfo);
1028  }
1029 
1030  LArMCParticleHelper::CollectReconstructable2DHits(pfoList, selectedMCParticleToHitsMaps, reconstructableCaloHitList2D);
1031 }
1032 
1033 //------------------------------------------------------------------------------------------------------------------------------------------
1034 
1036  const MCContributionMapVector &selectedMCParticleToHitsMaps, CaloHitList &reconstructableCaloHitList2D, const bool foldBackHierarchy)
1037 {
1038 
1039  PfoList pfoList;
1040  // If foldBackHierarchy collect all 2D calo hits in pfo hierarchy
1041  // else collect hits directly belonging to pfo
1042  if (foldBackHierarchy)
1043  {
1044  // ATTN: Only collect downstream pfos for daughter test beam particles & cosmics
1045  if (pPfo->GetParentPfoList().empty() && LArPfoHelper::IsTestBeam(pPfo))
1046  {
1047  pfoList.push_back(pPfo);
1048  }
1049  else
1050  {
1051  LArPfoHelper::GetAllDownstreamPfos(pPfo, pfoList);
1052  }
1053  }
1054  else
1055  {
1056  pfoList.push_back(pPfo);
1057  }
1058 
1059  LArMCParticleHelper::CollectReconstructable2DHits(pfoList, selectedMCParticleToHitsMaps, reconstructableCaloHitList2D);
1060 }
1061 
1062 //------------------------------------------------------------------------------------------------------------------------------------------
1063 
1065  const PfoList &pfoList, const MCContributionMapVector &selectedMCParticleToHitsMaps, CaloHitList &reconstructableCaloHitList2D)
1066 {
1067  CaloHitList caloHitList2D;
1068  LArPfoHelper::GetCaloHits(pfoList, TPC_VIEW_U, caloHitList2D);
1069  LArPfoHelper::GetCaloHits(pfoList, TPC_VIEW_V, caloHitList2D);
1070  LArPfoHelper::GetCaloHits(pfoList, TPC_VIEW_W, caloHitList2D);
1071  LArPfoHelper::GetIsolatedCaloHits(pfoList, TPC_VIEW_U, caloHitList2D); // TODO check isolated usage throughout
1072  LArPfoHelper::GetIsolatedCaloHits(pfoList, TPC_VIEW_V, caloHitList2D);
1073  LArPfoHelper::GetIsolatedCaloHits(pfoList, TPC_VIEW_W, caloHitList2D);
1074 
1075  // Filter for only reconstructable hits
1076  for (const CaloHit *const pCaloHit : caloHitList2D)
1077  {
1078  bool isTargetHit(false);
1079  for (const MCContributionMap &mcParticleToHitsMap : selectedMCParticleToHitsMaps)
1080  {
1081  // ATTN This map is unordered, but this does not impact search for specific target hit
1082  for (const MCContributionMap::value_type &mapEntry : mcParticleToHitsMap)
1083  {
1084  if (std::find(mapEntry.second.begin(), mapEntry.second.end(), pCaloHit) != mapEntry.second.end())
1085  {
1086  isTargetHit = true;
1087  break;
1088  }
1089  }
1090  if (isTargetHit)
1091  break;
1092  }
1093 
1094  if (isTargetHit)
1095  reconstructableCaloHitList2D.push_back(pCaloHit);
1096  }
1097 }
1098 
1099 //------------------------------------------------------------------------------------------------------------------------------------------
1100 
1101 void LArMCParticleHelper::CollectReconstructable2DHits(const pandora::Cluster *const pCluster,
1102  const MCContributionMapVector &selectedMCToHitsMaps, pandora::CaloHitList &reconstructableCaloHitList2D)
1103 {
1104  const CaloHitList &isolatedCaloHitList{pCluster->GetIsolatedCaloHitList()};
1105  CaloHitList caloHitList;
1106  pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
1107  for (const CaloHit *pCaloHit : isolatedCaloHitList)
1108  caloHitList.push_back(pCaloHit);
1109 
1110  // Filter for only reconstructable hits
1111  for (const CaloHit *const pCaloHit : caloHitList)
1112  {
1113  bool isTargetHit{false};
1114  for (const MCContributionMap &mcParticleToHitsMap : selectedMCToHitsMaps)
1115  { // ATTN This map is unordered, but this does not impact search for specific target hit
1116  for (const MCContributionMap::value_type &mapEntry : mcParticleToHitsMap)
1117  {
1118  if (std::find(mapEntry.second.begin(), mapEntry.second.end(), pCaloHit) != mapEntry.second.end())
1119  {
1120  isTargetHit = true;
1121  break;
1122  }
1123  }
1124  if (isTargetHit)
1125  break;
1126  }
1127 
1128  if (isTargetHit)
1129  reconstructableCaloHitList2D.push_back(pCaloHit);
1130  }
1131 }
1132 
1133 //------------------------------------------------------------------------------------------------------------------------------------------
1134 
1135 void LArMCParticleHelper::SelectCaloHits(const CaloHitList *const pCaloHitList, const LArMCParticleHelper::MCRelationMap &mcToTargetMCMap,
1136  CaloHitList &selectedCaloHitList, const bool selectInputHits, const float maxPhotonPropagation)
1137 {
1138  if (!selectInputHits)
1139  {
1140  selectedCaloHitList.insert(selectedCaloHitList.end(), pCaloHitList->begin(), pCaloHitList->end());
1141  return;
1142  }
1143 
1144  for (const CaloHit *const pCaloHit : *pCaloHitList)
1145  {
1146  try
1147  {
1148  const MCParticle *const pHitParticle(MCParticleHelper::GetMainMCParticle(pCaloHit));
1149 
1150  LArMCParticleHelper::MCRelationMap::const_iterator mcIter = mcToTargetMCMap.find(pHitParticle);
1151 
1152  if (mcToTargetMCMap.end() == mcIter)
1153  continue;
1154 
1155  // ATTN With folding on or off, still require primary particle to review hierarchy details
1156  const MCParticle *const pPrimaryParticle = LArMCParticleHelper::GetPrimaryMCParticle(pHitParticle);
1157 
1158  if (PassMCParticleChecks(pPrimaryParticle, pPrimaryParticle, pHitParticle, maxPhotonPropagation))
1159  selectedCaloHitList.push_back(pCaloHit);
1160  }
1161  catch (const StatusCodeException &)
1162  {
1163  }
1164  }
1165 }
1166 
1167 //------------------------------------------------------------------------------------------------------------------------------------------
1168 
1169 bool LArMCParticleHelper::IsDescendentOf(const MCParticle *const pMCParticle, const int pdg, const bool isChargeSensitive)
1170 {
1171  const MCParticle *pCurrentParticle = pMCParticle;
1172  while (!pCurrentParticle->GetParentList().empty())
1173  {
1174  if (pCurrentParticle->GetParentList().size() > 1)
1175  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
1176 
1177  const MCParticle *pParent = *(pCurrentParticle->GetParentList().begin());
1178  const bool found{isChargeSensitive ? pParent->GetParticleId() == pdg : std::abs(pParent->GetParticleId()) == std::abs(pdg)};
1179  if (found)
1180  return true;
1181  pCurrentParticle = pParent;
1182  }
1183 
1184  return false;
1185 }
1186 
1187 //------------------------------------------------------------------------------------------------------------------------------------------
1188 
1189 void LArMCParticleHelper::GetBreadthFirstHierarchyRepresentation(const MCParticle *const pMCParticle, MCParticleList &mcParticleList)
1190 {
1191  const MCParticle *const pRoot{LArMCParticleHelper::GetParentMCParticle(pMCParticle)};
1192  MCParticleList queue;
1193  mcParticleList.emplace_back(pRoot);
1194  queue.emplace_back(pRoot);
1195 
1196  while (!queue.empty())
1197  {
1198  const MCParticleList &daughters{queue.front()->GetDaughterList()};
1199  queue.pop_front();
1200  for (const MCParticle *pDaughter : daughters)
1201  {
1202  mcParticleList.emplace_back(pDaughter);
1203  queue.emplace_back(pDaughter);
1204  }
1205  }
1206 }
1207 
1208 //------------------------------------------------------------------------------------------------------------------------------------------
1209 
1210 void LArMCParticleHelper::SelectParticlesByHitCount(const MCParticleVector &candidateTargets, const MCContributionMap &mcToTrueHitListMap,
1211  const MCRelationMap &mcToTargetMCMap, const PrimaryParameters &parameters, MCContributionMap &selectedMCParticlesToHitsMap)
1212 {
1213  // Apply restrictions on the number of good hits associated with the MCParticles
1214  for (const MCParticle *const pMCTarget : candidateTargets)
1215  {
1216  MCContributionMap::const_iterator trueHitsIter = mcToTrueHitListMap.find(pMCTarget);
1217  if (mcToTrueHitListMap.end() == trueHitsIter)
1218  continue;
1219 
1220  const CaloHitList &caloHitList(trueHitsIter->second);
1221 
1222  // Remove shared hits where target particle deposits below threshold energy fraction
1223  CaloHitList goodCaloHitList;
1225  &caloHitList, mcToTargetMCMap, goodCaloHitList, parameters.m_selectInputHits, parameters.m_minHitSharingFraction);
1226 
1227  if (goodCaloHitList.size() < parameters.m_minPrimaryGoodHits)
1228  continue;
1229 
1230  unsigned int nGoodViews(0);
1231  if (LArMonitoringHelper::CountHitsByType(TPC_VIEW_U, goodCaloHitList) >= parameters.m_minHitsForGoodView)
1232  ++nGoodViews;
1233 
1234  if (LArMonitoringHelper::CountHitsByType(TPC_VIEW_V, goodCaloHitList) >= parameters.m_minHitsForGoodView)
1235  ++nGoodViews;
1236 
1237  if (LArMonitoringHelper::CountHitsByType(TPC_VIEW_W, goodCaloHitList) >= parameters.m_minHitsForGoodView)
1238  ++nGoodViews;
1239 
1240  if (nGoodViews < parameters.m_minPrimaryGoodViews)
1241  continue;
1242 
1243  if (!selectedMCParticlesToHitsMap.insert(MCContributionMap::value_type(pMCTarget, caloHitList)).second)
1244  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
1245  }
1246 }
1247 
1248 //------------------------------------------------------------------------------------------------------------------------------------------
1249 
1250 void LArMCParticleHelper::SelectGoodCaloHits(const CaloHitList *const pSelectedCaloHitList, const LArMCParticleHelper::MCRelationMap &mcToTargetMCMap,
1251  CaloHitList &selectedGoodCaloHitList, const bool selectInputHits, const float minHitSharingFraction)
1252 {
1253  if (!selectInputHits)
1254  {
1255  selectedGoodCaloHitList.insert(selectedGoodCaloHitList.end(), pSelectedCaloHitList->begin(), pSelectedCaloHitList->end());
1256  return;
1257  }
1258 
1259  for (const CaloHit *const pCaloHit : *pSelectedCaloHitList)
1260  {
1261  MCParticleVector mcParticleVector;
1262  for (const auto &mapEntry : pCaloHit->GetMCParticleWeightMap())
1263  mcParticleVector.push_back(mapEntry.first);
1264  std::sort(mcParticleVector.begin(), mcParticleVector.end(), PointerLessThan<MCParticle>());
1265 
1266  MCParticleWeightMap targetWeightMap;
1267 
1268  for (const MCParticle *const pMCParticle : mcParticleVector)
1269  {
1270  const float weight(pCaloHit->GetMCParticleWeightMap().at(pMCParticle));
1271  LArMCParticleHelper::MCRelationMap::const_iterator mcIter = mcToTargetMCMap.find(pMCParticle);
1272 
1273  if (mcToTargetMCMap.end() != mcIter)
1274  targetWeightMap[mcIter->second] += weight;
1275  }
1276 
1277  MCParticleVector mcTargetVector;
1278  for (const auto &mapEntry : targetWeightMap)
1279  mcTargetVector.push_back(mapEntry.first);
1280  std::sort(mcTargetVector.begin(), mcTargetVector.end(), PointerLessThan<MCParticle>());
1281 
1282  const MCParticle *pBestTargetParticle(nullptr);
1283  float bestTargetWeight(0.f), targetWeightSum(0.f);
1284 
1285  for (const MCParticle *const pTargetMCParticle : mcTargetVector)
1286  {
1287  const float targetWeight(targetWeightMap.at(pTargetMCParticle));
1288  targetWeightSum += targetWeight;
1289 
1290  if (targetWeight > bestTargetWeight)
1291  {
1292  bestTargetWeight = targetWeight;
1293  pBestTargetParticle = pTargetMCParticle;
1294  }
1295  }
1296 
1297  if (!pBestTargetParticle || (targetWeightSum < std::numeric_limits<float>::epsilon()) || ((bestTargetWeight / targetWeightSum) < minHitSharingFraction))
1298  continue;
1299 
1300  selectedGoodCaloHitList.push_back(pCaloHit);
1301  }
1302 }
1303 
1304 //------------------------------------------------------------------------------------------------------------------------------------------
1305 
1307  std::function<bool(const MCParticle *const)> fCriteria, MCParticleVector &selectedParticles, const PrimaryParameters &parameters, const bool isTestBeam)
1308 {
1309  for (const MCParticle *const pMCParticle : inputMCParticles)
1310  {
1311  if (parameters.m_foldBackHierarchy)
1312  {
1313  if (!fCriteria(pMCParticle))
1314  continue;
1315  }
1316  else
1317  {
1318  if (isTestBeam)
1319  {
1320  if (!LArMCParticleHelper::DoesLeadingMeetCriteria(pMCParticle, fCriteria))
1321  continue;
1322  }
1323  else
1324  {
1325  if (!LArMCParticleHelper::DoesPrimaryMeetCriteria(pMCParticle, fCriteria))
1326  continue;
1327  }
1328  }
1329  selectedParticles.push_back(pMCParticle);
1330  }
1331 }
1332 
1333 //------------------------------------------------------------------------------------------------------------------------------------------
1334 
1335 bool LArMCParticleHelper::PassMCParticleChecks(const MCParticle *const pOriginalPrimary, const MCParticle *const pThisMCParticle,
1336  const MCParticle *const pHitMCParticle, const float maxPhotonPropagation)
1337 {
1338  if (NEUTRON == std::abs(pThisMCParticle->GetParticleId()))
1339  return false;
1340 
1341  if ((PHOTON == pThisMCParticle->GetParticleId()) && (PHOTON != pOriginalPrimary->GetParticleId()) &&
1342  (E_MINUS != std::abs(pOriginalPrimary->GetParticleId())))
1343  {
1344  if ((pThisMCParticle->GetEndpoint() - pThisMCParticle->GetVertex()).GetMagnitude() > maxPhotonPropagation)
1345  return false;
1346  }
1347 
1348  if (pThisMCParticle == pHitMCParticle)
1349  return true;
1350 
1351  for (const MCParticle *const pDaughterMCParticle : pThisMCParticle->GetDaughterList())
1352  {
1353  if (PassMCParticleChecks(pOriginalPrimary, pDaughterMCParticle, pHitMCParticle, maxPhotonPropagation))
1354  return true;
1355  }
1356 
1357  return false;
1358 }
1359 
1360 } // namespace lar_content
static bool IsVisible(const pandora::MCParticle *const pMCParticle)
Whether a mc particle is visible (i.e. long-lived charged particle)
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...
static bool PassMCParticleChecks(const pandora::MCParticle *const pOriginalPrimary, const pandora::MCParticle *const pThisMCParticle, const pandora::MCParticle *const pHitMCParticle, const float maxPhotonPropagation)
Whether it is possible to navigate from a primary mc particle to a downstream mc particle without "pa...
unsigned int m_minPrimaryGoodViews
the minimum number of primary good views
static bool SortByNHits(const pandora::ParticleFlowObject *const pLhs, const pandora::ParticleFlowObject *const pRhs)
Sort pfos by number of constituent hits.
static bool IsNuclear(const pandora::MCParticle *const pMCParticle)
Check whether or not an MC particle comes from a nuclear interaction process.
Header file for the pfo helper class.
std::unordered_map< const pandora::MCParticle *, pandora::CaloHitList > MCContributionMap
bool m_selectInputHits
whether to select input hits
std::pair< const pandora::MCParticle *, pandora::CaloHitList > MCParticleCaloHitListPair
static void GetTestBeamHierarchyPfoToReconstructable2DHitsMap(const pandora::PfoList &pfoList, const MCContributionMap &selectedMCParticleToHitsMap, PfoContributionMap &pfoToReconstructable2DHitsMap, const bool foldBackHierarchy)
Get mapping from Pfo in reconstructed test beam hierarchy to reconstructable 2D hits (=good hits belo...
static bool IsPrimary(const pandora::MCParticle *const pMCParticle)
Whether a provided mc particle matches the implemented definition of being primary.
static void GetTwoDClusterList(const pandora::ParticleFlowObject *const pPfo, pandora::ClusterList &clusterList)
Get the list of 2D clusters from an input pfo.
static void GetLeadingMCParticleList(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &mcLeadingVector)
Get vector of leading MC particles from an input list of MC particles.
static bool IsLeadingBeamParticle(const pandora::MCParticle *const pMCParticle)
Returns true if passed a leading beam MCParticle.
unsigned int m_minPrimaryGoodHits
the minimum number of primary good Hits
int GetNuanceCode() const
Get the nuance code.
static void GetPfoToReconstructable2DHitsMap(const pandora::PfoList &pfoList, const MCContributionMap &selectedMCParticleToHitsMap, PfoContributionMap &pfoToReconstructable2DHitsMap, const bool foldBackHierarchy)
Get mapping from Pfo to reconstructable 2D hits (=good hits belonging to a selected reconstructable M...
std::unordered_map< const pandora::CaloHit *, const pandora::MCParticle * > CaloHitToMCMap
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 GetPfoMCParticleHitSharingMaps(const PfoContributionMap &pfoToReconstructable2DHitsMap, const MCContributionMapVector &selectedMCParticleToHitsMaps, PfoToMCParticleHitSharingMap &pfoToMCParticleHitSharingMap, MCParticleToPfoHitSharingMap &mcParticleToPfoHitSharingMap)
Get the mappings from Pfo -> pair (reconstructable MCparticles, number of reconstructable 2D hits sha...
static const pandora::MCParticle * GetPrimaryMCParticle(const pandora::MCParticle *const pMCParticle)
Get the primary parent mc particle.
std::map< const pandora::MCParticle *, PfoToSharedHitsVector > MCParticleToPfoHitSharingMap
std::vector< PfoCaloHitListPair > PfoToSharedHitsVector
intermediate_table::const_iterator const_iterator
static void SelectGoodCaloHits(const pandora::CaloHitList *const pSelectedCaloHitList, const MCRelationMap &mcToTargetMCMap, pandora::CaloHitList &selectedGoodCaloHitList, const bool selectInputHits, const float minHitSharingFraction)
Apply further selection criteria to end up with a collection of "good" calo hits that can be use to d...
bool m_foldBackHierarchy
whether to fold the hierarchy back to the primary (neutrino) or leading particles (test beam) ...
static void SelectReconstructableTestBeamHierarchyMCParticles(const pandora::MCParticleList *pMCParticleList, const pandora::CaloHitList *pCaloHitList, const PrimaryParameters &parameters, std::function< bool(const pandora::MCParticle *const)> fCriteria, MCContributionMap &selectedMCParticlesToHitsMap)
Select target, reconstructable mc particles in the relevant hierarchy that match given criteria...
static void GetMCPrimaryMap(const pandora::MCParticleList *const pMCParticleList, MCRelationMap &mcPrimaryMap)
Get mapping from individual mc particles (in a provided list) and their primary parent mc particles...
static bool IsPairProduction(const pandora::MCParticle *const pMCParticle)
Check whether or not an MC particle comes from a pair production process.
Header file for the lar monitoring helper helper class.
static const pandora::MCParticle * GetMainMCParticle(const pandora::ParticleFlowObject *const pPfo)
Find the mc particle making the largest contribution to 2D clusters in a specified pfo...
static unsigned int GetNuanceCode(const pandora::MCParticle *const pMCParticle)
Get the nuance code of an MCParticle.
unsigned int m_minHitsForGoodView
the minimum number of Hits for a good view
static const pandora::MCParticle * GetLeadingMCParticle(const pandora::MCParticle *const pMCParticle, const int hierarchyTierLimit=1)
Get the leading particle in the hierarchy, for use at ProtoDUNE.
float m_maxPhotonPropagation
the maximum photon propagation length
TFile f
Definition: plotHisto.C:6
static void SelectParticlesByHitCount(const pandora::MCParticleVector &candidateTargets, const MCContributionMap &mcToTrueHitListMap, const MCRelationMap &mcToTargetMCMap, const PrimaryParameters &parameters, MCContributionMap &selectedMCParticlesToHitsMap)
Filter an input vector of MCParticles to ensure they have sufficient good hits to be reconstructable...
static void CollectReconstructableTestBeamHierarchy2DHits(const pandora::ParticleFlowObject *const pPfo, const MCContributionMapVector &selectedMCParticleToHitsMaps, pandora::CaloHitList &reconstructableCaloHitList2D, const bool foldBackHierarchy)
For a given Pfo, collect the hits which are reconstructable (=good hits belonging to a selected recon...
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
static bool SortByMomentum(const pandora::MCParticle *const pLhs, const pandora::MCParticle *const pRhs)
Sort mc particles by their momentum.
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 bool IsLeading(const pandora::MCParticle *const pMCParticle)
Whether a provided mc particle matches the implemented definition of being leading.
static void SelectReconstructableMCParticles(const pandora::MCParticleList *pMCParticleList, const pandora::CaloHitList *pCaloHitList, const PrimaryParameters &parameters, std::function< bool(const pandora::MCParticle *const)> fCriteria, MCContributionMap &selectedMCParticlesToHitsMap)
Select target, reconstructable mc particles that match given criteria.
static bool IsTestBeam(const pandora::ParticleFlowObject *const pPfo)
Whether a pfo is a test beam particle.
static void GetMCToSelfMap(const pandora::MCParticleList *const pMCParticleList, MCRelationMap &mcToSelfMap)
Get mapping from individual mc particles (in a provided list) to themselves (to be used when not fold...
Header file for the lar monte carlo particle helper helper class.
Header file for the cluster helper class.
static void GetMCParticleToCaloHitMatches(const pandora::CaloHitList *const pCaloHitList, const MCRelationMap &mcToTargetMCMap, CaloHitToMCMap &hitToMCMap, MCContributionMap &mcToTrueHitListMap)
Match calo hits to their parent particles.
static bool IsBeamParticle(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary beam MCParticle.
std::unordered_map< const pandora::Cluster *, pandora::CaloHitList > ClusterContributionMap
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
static const pandora::MCParticle * GetParentMCParticle(const pandora::MCParticle *const pMCParticle)
Get the parent mc particle.
static bool IsTriggeredBeamParticle(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary triggered beam MCParticle.
float m_minHitSharingFraction
the minimum Hit sharing fraction
static void GetMCLeadingMap(const pandora::MCParticleList *const pMCParticleList, MCRelationMap &mcLeadingMap)
Get mapping from individual mc particles (in a provided list) and their leading parent mc particles...
static bool IsDecay(const pandora::MCParticle *const pMCParticle)
Check whether or not an MC particle comes from a decay process.
std::vector< MCContributionMap > MCContributionMapVector
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.
static bool IsBremsstrahlung(const pandora::MCParticle *const pMCParticle)
static void GetPrimaryMCParticleList(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &mcPrimaryVector)
Get vector of primary MC particles from an input list of MC particles.
double weight
Definition: plottest35.C:25
static bool DoesPrimaryMeetCriteria(const pandora::MCParticle *const pMCParticle, std::function< bool(const pandora::MCParticle *const)> fCriteria)
Returns true if passed particle whose primary meets the passed criteria.
static void GetAllAncestorMCParticles(const pandora::MCParticle *const pMCParticle, pandora::MCParticleList &ancestorMCParticleList)
Get all ancestor mc particles.
static void CollectReconstructable2DHits(const pandora::ParticleFlowObject *const pPfo, const MCContributionMapVector &selectedMCParticleToHitsMaps, pandora::CaloHitList &reconstructableCaloHitList2D, const bool foldBackHierarchy)
For a given Pfo, collect the hits which are reconstructable (=good hits belonging to a selected recon...
static void GetTrueNeutrinos(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &trueNeutrinos)
Get neutrino MC particles from an input MC particle list.
std::pair< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoCaloHitListPair
static bool IsCapture(const pandora::MCParticle *const pMCParticle)
Check whether or not an MC particle comes from a capture process.
static void GetIsolatedCaloHits(const pandora::PfoList &pfoList, const pandora::HitType &hitType, pandora::CaloHitList &caloHitList)
Get a list of isolated calo hits of a particular hit type from a list of pfos.
static void GetClusterToReconstructable2DHitsMap(const pandora::ClusterList &clusterList, const MCContributionMap &selectedMCToHitsMap, ClusterContributionMap &clusterToReconstructable2DHitsMap)
Get mapping from cluster to reconstructable 2D hits (=good hits belonging to a selected reconstructab...
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...
static bool IsIonisation(const pandora::MCParticle *const pMCParticle)
Check whether or not an MC particle comes from an ionisation process.
static bool DoesLeadingMeetCriteria(const pandora::MCParticle *const pMCParticle, std::function< bool(const pandora::MCParticle *const)> fCriteria)
Returns true if passed particle whose leading meets the passed criteria.
std::vector< MCParticleCaloHitListPair > MCParticleToSharedHitsVector
static void SelectParticlesMatchingCriteria(const pandora::MCParticleVector &inputMCParticles, std::function< bool(const pandora::MCParticle *const)> fCriteria, pandora::MCParticleVector &selectedParticles, const PrimaryParameters &parameters, const bool isTestBeam)
Select mc particles matching given criteria from an input list.
static void GetTrueTestBeamParticles(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &trueTestBeamParticles)
Get triggered test beam MC particles from an input MC particle list.
static bool IsDescendentOf(const pandora::MCParticle *const pMCParticle, const int pdg, const bool isChargeSensitive=false)
Determine if the MC particle is a descendent of a particle with the given PDG code.
std::unordered_map< const pandora::ParticleFlowObject *, pandora::CaloHitList > PfoContributionMap
static void GetAllDescendentMCParticles(const pandora::MCParticle *const pMCParticle, pandora::MCParticleList &descendentMCParticleList)
Get all descendent mc particles.
static void GetCaloHits(const pandora::PfoList &pfoList, const pandora::HitType &hitType, pandora::CaloHitList &caloHitList)
Get a list of calo hits of a particular hit type from a list of pfos.
static pandora::CaloHitList GetSharedHits(const pandora::CaloHitList &hitListA, const pandora::CaloHitList &hitListB)
Get the hits in the intersection of two hit lists.
static bool IsBeamNeutrinoFinalState(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary neutrino final state MCParticle.
static unsigned int CountHitsByType(const pandora::HitType hitType, const pandora::CaloHitList &caloHitList)
Count the number of calo hits, in a provided list, of a specified type.
std::unordered_map< const pandora::MCParticle *, const pandora::MCParticle * > MCRelationMap
static void GetBreadthFirstHierarchyRepresentation(const pandora::MCParticle *const pMCParticle, pandora::MCParticleList &mcParticleList)
Retrieve a linearised representation of the MC particle hierarchy in breadth first order...
static bool IsElasticScatter(const pandora::MCParticle *const pMCParticle)
Check whether or not an MC particle came from an elastic scattering process.
std::map< const pandora::ParticleFlowObject *, MCParticleToSharedHitsVector > PfoToMCParticleHitSharingMap
static bool IsNeutrino(const pandora::MCParticle *const pMCParticle)
Whether a mc particle is a neutrino or antineutrino.
static void SelectCaloHits(const pandora::CaloHitList *const pCaloHitList, const MCRelationMap &mcToTargetMCMap, pandora::CaloHitList &selectedCaloHitList, const bool selectInputHits, const float maxPhotonPropagation)
Select a subset of calo hits representing those that represent "reconstructable" regions of the event...