9 #include "Pandora/PdgTable.h" 10 #include "Pandora/StatusCodes.h" 23 m_foldToLeadingShowers{
false},
53 std::cout <<
"LArHierarchyHelper: Error - attempting to fold to non-positive tier" << std::endl;
54 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
98 for (
const Node *pNode : nodeVector)
101 m_interactions.clear();
108 const auto predicate = [](
const MCParticle *pMCParticle) {
return std::abs(pMCParticle->GetParticleId()) == NEUTRON; };
110 for (
const CaloHit *pCaloHit : caloHitList)
114 const MCParticle *pMCParticle{MCParticleHelper::GetMainMCParticle(pCaloHit)};
117 catch (
const StatusCodeException &)
122 MCParticleList rootNodes;
123 for (
const MCParticle *pMCParticle : mcParticleList)
125 const MCParticleList &parentList{pMCParticle->GetParentList()};
126 if (parentList.empty())
128 rootNodes.emplace_back(pMCParticle);
132 for (
const MCParticle *pRoot : rootNodes)
136 MCParticleList primaries(primarySet.begin(), primarySet.end());
139 primaries.erase(std::remove_if(primaries.begin(), primaries.end(), predicate), primaries.end());
142 for (
const MCParticle *pPrimary : primaries)
144 MCParticleList allParticles{pPrimary};
152 MCParticleList dummy;
156 for (
const MCParticle *pMCParticle : allParticles)
162 allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
170 for (
const MCParticle *pPrimary : primaries)
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};
179 for (
const MCParticle *pMCParticle : allParticles)
185 allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
188 Node *pNode{
new Node(*
this, allParticles, allHits)};
190 if (!(isShower || isNeutron))
193 const MCParticleList &children{pPrimary->GetDaughterList()};
194 for (
const MCParticle *pChild : children)
195 pNode->FillHierarchy(pChild, foldParameters);
201 for (
const MCParticle *pPrimary : primaries)
203 MCParticleList leadingParticles, childParticles;
206 for (
const MCParticle *pMCParticle : leadingParticles)
212 allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
216 Node *pNode{
new Node(*
this, leadingParticles, allHits)};
218 for (
const MCParticle *pChild : childParticles)
219 pNode->FillHierarchy(pChild, foldParameters);
225 for (
const MCParticle *pPrimary : primaries)
227 MCParticleList allParticles{pPrimary};
229 for (
const MCParticle *pMCParticle : allParticles)
235 allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
238 Node *pNode{
new Node(*
this, allParticles, allHits)};
241 const MCParticleList &children{pPrimary->GetDaughterList()};
242 for (
const MCParticle *pChild : children)
243 pNode->FillHierarchy(pChild, foldParameters);
247 Node *pLeadingLepton{
nullptr};
248 float leadingLeptonEnergy{-std::numeric_limits<float>::max()};
254 const int pdg{
std::abs(pMC->GetParticleId())};
255 if ((pdg == MU_MINUS || pdg == E_MINUS || pdg == TAU_MINUS) && pMC->GetEnergy() > leadingLeptonEnergy)
257 pLeadingLepton =
const_cast<Node *
>(pNode);
258 leadingLeptonEnergy = pMC->GetEnergy();
272 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
282 rootMCParticles.emplace_back(iter->first);
288 const MCParticle *
const pRoot, MCParticleList &leadingParticles, MCParticleList &childParticles,
const float cosAngleTolerance)
const 290 leadingParticles.emplace_back(pRoot);
291 MCParticleList foldCandidates, childCandidates;
292 const MCParticleList &children{pRoot->GetDaughterList()};
293 for (
const MCParticle *pMCParticle : children)
302 if (pMCParticle->GetParticleId() == pRoot->GetParticleId())
303 foldCandidates.emplace_back(pMCParticle);
305 childCandidates.emplace_back(pMCParticle);
310 childCandidates.emplace_back(pMCParticle);
313 const MCParticle *pBestFoldCandidate{
nullptr};
314 float bestDp{std::numeric_limits<float>::max()};
315 for (
const MCParticle *pMCParticle : foldCandidates)
317 if (foldCandidates.size() == 1)
322 pBestFoldCandidate = pMCParticle;
324 childCandidates.emplace_back(pMCParticle);
331 const float dp{pRoot->GetMomentum().GetMagnitude() - pMCParticle->GetMomentum().GetMagnitude()};
334 pBestFoldCandidate = pMCParticle;
340 childCandidates.emplace_back(pMCParticle);
344 if (pBestFoldCandidate)
346 leadingParticles.emplace_back(pBestFoldCandidate);
349 this->
CollectContinuations(pBestFoldCandidate, leadingParticles, childCandidates, cosAngleTolerance);
351 for (
const MCParticle *pMCParticle : childCandidates)
355 childParticles.emplace_back(pMCParticle);
358 MCParticleList localHierarchy{pMCParticle};
359 CaloHitList localHits;
361 for (
const MCParticle *pLocalMCParticle : localHierarchy)
365 const CaloHitList &caloHits(
m_mcToHitsMap.at(pLocalMCParticle));
366 localHits.insert(localHits.begin(), caloHits.begin(), caloHits.end());
370 childParticles.emplace_back(pMCParticle);
378 const MCParticle *pRoot, MCParticleList &continuingParticles, MCParticleList &childParticles,
const float cosAngleTolerance)
const 380 const MCParticleList &children{pRoot->GetDaughterList()};
381 MCParticleList foldCandidates;
382 for (
const MCParticle *pMCParticle : children)
390 if (pMCParticle->GetParticleId() == pRoot->GetParticleId())
391 foldCandidates.emplace_back(pMCParticle);
396 childParticles.emplace_back(pMCParticle);
399 const MCParticle *pBestFoldCandidate{
nullptr};
400 float bestDp{std::numeric_limits<float>::max()};
401 for (
const MCParticle *pMCParticle : foldCandidates)
403 if (foldCandidates.size() == 1)
408 pBestFoldCandidate = pMCParticle;
415 const float dp{pRoot->GetMomentum().GetMagnitude() - pMCParticle->GetMomentum().GetMagnitude()};
418 pBestFoldCandidate = pMCParticle;
424 if (pBestFoldCandidate)
426 continuingParticles.emplace_back(pBestFoldCandidate);
427 const MCParticleList &newLeadingParticles{pBestFoldCandidate->GetDaughterList()};
429 childParticles.insert(childParticles.begin(), newLeadingParticles.begin(), newLeadingParticles.end());
431 const auto iter{std::find(childParticles.begin(), childParticles.end(), pBestFoldCandidate)};
432 if (iter != childParticles.end())
433 childParticles.erase(iter);
447 nodeVector.emplace_back(pNode);
448 queue.emplace_back(pNode);
450 while (!queue.empty())
452 const NodeVector &children{queue.front()->GetChildren()};
454 for (
const Node *pChild : children)
456 nodeVector.emplace_back(pChild);
457 queue.emplace_back(pChild);
479 str +=
"=== MC Interaction : PDG " +
std::to_string(pLArRoot->GetParticleId()) +
483 for (
const Node *pNode : nodeVector)
484 str +=
" " + pNode->ToString(
"") +
"\n";
496 unsigned int nHitsU{0}, nHitsV{0}, nHitsW{0};
497 for (
const CaloHit *pCaloHit :
m_mcToHitsMap.at(pMCParticle))
499 const HitType view{pCaloHit->GetHitType()};
500 if (view == TPC_VIEW_U)
502 else if (view == TPC_VIEW_V)
504 else if (view == TPC_VIEW_W)
507 const unsigned int nHits{nHitsU + nHitsV + nHitsW};
508 unsigned int nGoodViews{0};
523 unsigned int nHitsU{0}, nHitsV{0}, nHitsW{0};
524 for (
const CaloHit *pCaloHit : caloHits)
526 const HitType view{pCaloHit->GetHitType()};
527 if (view == TPC_VIEW_U)
529 else if (view == TPC_VIEW_V)
531 else if (view == TPC_VIEW_W)
534 const unsigned int nHits{nHitsU + nHitsV + nHitsW};
535 unsigned int nGoodViews{0};
547 m_hierarchy(hierarchy),
548 m_mainParticle(pMCParticle),
551 m_isLeadingLepton{
false}
555 m_pdg = pMCParticle->GetParticleId();
556 m_mcParticles.emplace_back(pMCParticle);
558 m_hierarchy.RegisterNode(
this);
564 m_hierarchy(hierarchy),
565 m_mcParticles(mcParticleList),
566 m_caloHits(caloHitList),
567 m_mainParticle(
nullptr),
570 m_isLeadingLepton{
false}
572 if (!mcParticleList.empty())
574 m_mainParticle = mcParticleList.front();
575 m_pdg = m_mainParticle->GetParticleId();
579 m_hierarchy.RegisterNode(
this);
586 m_mcParticles.clear();
588 for (
const Node *node : m_children)
599 MCParticleList leadingParticles, childParticles;
600 m_hierarchy.InterpretHierarchy(pRoot, leadingParticles, childParticles, foldParameters.
m_cosAngleTolerance);
602 for (
const MCParticle *pMCParticle : leadingParticles)
605 if (m_hierarchy.m_mcToHitsMap.find(pMCParticle) != m_hierarchy.m_mcToHitsMap.end())
607 const CaloHitList &caloHits(m_hierarchy.m_mcToHitsMap.at(pMCParticle));
608 allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
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)
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};
626 else if (foldParameters.
m_foldToLeadingShowers && (isShower || (isNeutron && !m_hierarchy.m_recoCriteria.m_removeNeutrons)))
628 else if (m_hierarchy.m_recoCriteria.m_removeNeutrons && isNeutron)
632 for (
const MCParticle *pMCParticle : allParticles)
635 if (m_hierarchy.m_mcToHitsMap.find(pMCParticle) != m_hierarchy.m_mcToHitsMap.end())
637 const CaloHitList &caloHits(m_hierarchy.m_mcToHitsMap.at(pMCParticle));
638 allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
642 if (!allParticles.empty())
648 if (hasChildren || (!hasChildren && !allHits.empty()))
650 Node *pNode{
new Node(m_hierarchy, allParticles, allHits, this->m_tier + 1)};
651 m_children.emplace_back(pNode);
655 const MCParticleList &children{pRoot->GetDaughterList()};
656 for (
const MCParticle *pChild : children)
668 MCParticleList allParticles{pRoot};
669 if (!m_hierarchy.m_recoCriteria.m_removeNeutrons)
679 for (
const MCParticle *pMCParticle : allParticles)
682 if (m_hierarchy.m_mcToHitsMap.find(pMCParticle) != m_hierarchy.m_mcToHitsMap.end())
684 const CaloHitList &caloHits(m_hierarchy.m_mcToHitsMap.at(pMCParticle));
685 allHits.insert(allHits.begin(), caloHits.begin(), caloHits.end());
688 if (!allParticles.empty())
690 Node *pNode{
new Node(m_hierarchy, allParticles, allHits, this->m_tier + 1)};
691 m_children.emplace_back(pNode);
699 return m_hierarchy.m_nodeToIdMap.at(
this);
706 const bool enoughHits{m_caloHits.size() >= m_hierarchy.m_recoCriteria.m_minHits};
709 bool enoughGoodViews{
false};
710 unsigned int nHitsU{0}, nHitsV{0}, nHitsW{0};
711 for (
const CaloHit *pCaloHit : m_caloHits)
713 switch (pCaloHit->GetHitType())
727 unsigned int nGoodViews{0};
728 if (nHitsU >= m_hierarchy.m_recoCriteria.m_minHitsForGoodView)
730 if (nHitsV >= m_hierarchy.m_recoCriteria.m_minHitsForGoodView)
732 if (nHitsW >= m_hierarchy.m_recoCriteria.m_minHitsForGoodView)
734 if (nGoodViews >= m_hierarchy.m_recoCriteria.m_minGoodViews)
736 enoughGoodViews =
true;
741 return enoughGoodViews;
768 std::string str(prefix +
"PDG: " +
std::to_string(m_pdg) +
" Energy: " +
std::to_string(m_mainParticle ? m_mainParticle->GetEnergy() : 0) +
770 for (
const Node *pChild : m_children)
771 str += pChild->ToString(prefix +
" ");
800 const unsigned int minHits,
const unsigned int minHitsForGoodView,
const unsigned int minGoodViews,
const bool removeNeutrons) :
821 for (
const Node *pNode : nodeVector)
824 m_interactions.clear();
832 for (
const ParticleFlowObject *pPfo : pfoList)
834 const PfoList &parentList{pPfo->GetParentPfoList()};
835 if (parentList.empty())
837 rootNodes.emplace_back(pPfo);
841 for (
const ParticleFlowObject *
const pRoot : rootNodes)
845 PfoList primaries(primarySet.begin(), primarySet.end());
849 for (
const ParticleFlowObject *pPrimary : primaries)
851 PfoList allParticles;
855 for (
const ParticleFlowObject *pPfo : allParticles)
862 for (
const ParticleFlowObject *pPrimary : primaries)
864 PfoList allParticles;
865 int pdg{
std::abs(pPrimary->GetParticleId())};
866 const bool isShower{pdg == E_MINUS};
871 allParticles.emplace_back(pPrimary);
874 for (
const ParticleFlowObject *pPfo : allParticles)
876 Node *pNode{
new Node(*
this, allParticles, allHits)};
881 const PfoList &children{pPrimary->GetDaughterPfoList()};
882 for (
const ParticleFlowObject *pChild : children)
883 pNode->FillHierarchy(pChild, foldParameters);
890 for (
const ParticleFlowObject *pPrimary : primaries)
892 PfoList allParticles{pPrimary};
894 for (
const ParticleFlowObject *pPfo : allParticles)
896 Node *pNode{
new Node(*
this, allParticles, allHits)};
899 const PfoList &children{pPrimary->GetDaughterPfoList()};
900 for (
const ParticleFlowObject *pChild : children)
912 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
922 rootPfos.emplace_back(iter->first);
932 nodeVector.emplace_back(pNode);
933 queue.emplace_back(pNode);
935 while (!queue.empty())
937 const NodeVector &children{queue.front()->GetChildren()};
939 for (
const Node *pChild : children)
941 nodeVector.emplace_back(pChild);
942 queue.emplace_back(pChild);
954 str +=
"=== Reco Interaction : PDG " +
std::to_string(pRoot->GetParticleId()) +
"\n";
955 for (
const Node *pNode : nodeVector)
956 str +=
" " + pNode->ToString(
"") +
"\n";
966 m_hierarchy{hierarchy},
972 m_pdg = pPfo->GetParticleId();
973 m_pfos.emplace_back(pPfo);
984 if (!pfoList.empty())
987 m_pdg = pfoList.front()->GetParticleId();
1010 PfoList allParticles;
1011 int pdg{
std::abs(pRoot->GetParticleId())};
1012 const bool isShower{pdg == E_MINUS};
1018 allParticles.emplace_back(pRoot);
1020 CaloHitList allHits;
1021 for (
const ParticleFlowObject *pPfo : allParticles)
1026 if (hasChildren || (!hasChildren && !allHits.empty()))
1033 const PfoList &children{pRoot->GetDaughterPfoList()};
1034 for (
const ParticleFlowObject *pChild : children)
1035 pNode->FillHierarchy(pChild, foldParameters);
1044 PfoList allParticles;
1046 CaloHitList allHits;
1047 for (
const ParticleFlowObject *pPfo : allParticles)
1080 str += pChild->ToString(prefix +
" ");
1089 m_pMCParticle{pMCParticle}
1107 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1119 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1121 const CaloHitList &recoHits{pReco->
GetCaloHits()};
1123 CaloHitVector intersection;
1124 std::set_intersection(mcHits.begin(), mcHits.end(), recoHits.begin(), recoHits.end(), std::back_inserter(intersection));
1126 return this->
GetPurity(intersection, recoHits, adcWeighted);
1135 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1137 CaloHitList recoHits;
1138 for (
const CaloHit *pCaloHit : pReco->
GetCaloHits())
1139 if (pCaloHit->GetHitType() == view)
1140 recoHits.emplace_back(pCaloHit);
1143 if (pCaloHit->GetHitType() == view)
1144 mcHits.emplace_back(pCaloHit);
1146 CaloHitVector intersection;
1147 std::set_intersection(mcHits.begin(), mcHits.end(), recoHits.begin(), recoHits.end(), std::back_inserter(intersection));
1149 return this->
GetPurity(intersection, recoHits, adcWeighted);
1158 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1160 const CaloHitList &recoHits{pReco->
GetCaloHits()};
1162 CaloHitVector intersection;
1163 std::set_intersection(mcHits.begin(), mcHits.end(), recoHits.begin(), recoHits.end(), std::back_inserter(intersection));
1174 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1176 CaloHitList recoHits;
1177 for (
const CaloHit *pCaloHit : pReco->
GetCaloHits())
1178 if (pCaloHit->GetHitType() == view)
1179 recoHits.emplace_back(pCaloHit);
1182 if (pCaloHit->GetHitType() == view)
1183 mcHits.emplace_back(pCaloHit);
1185 CaloHitVector intersection;
1186 std::set_intersection(mcHits.begin(), mcHits.end(), recoHits.begin(), recoHits.end(), std::back_inserter(intersection));
1196 if (!intersection.empty())
1201 for (
const CaloHit *pCaloHit : recoHits)
1202 adcSum += pCaloHit->GetInputEnergy();
1203 if (adcSum > std::numeric_limits<float>::epsilon())
1205 for (
const CaloHit *pCaloHit : intersection)
1206 purity += pCaloHit->GetInputEnergy();
1212 purity = intersection.size() /
static_cast<float>(recoHits.size());
1223 float completeness{0.f};
1224 if (!intersection.empty())
1229 for (
const CaloHit *pCaloHit : mcHits)
1230 adcSum += pCaloHit->GetInputEnergy();
1231 if (adcSum > std::numeric_limits<float>::epsilon())
1233 for (
const CaloHit *pCaloHit : intersection)
1234 completeness += pCaloHit->GetInputEnergy();
1235 completeness /= adcSum;
1240 completeness = intersection.size() /
static_cast<float>(mcHits.size());
1244 return completeness;
1254 int nAboveThreshold{0};
1260 if (nAboveThreshold != 1)
1294 MCParticleList rootMCParticles;
1298 std::map<const MCHierarchy::Node *, MCMatches> mcToMatchMap;
1300 for (
const MCParticle *
const pRootMC : rootMCParticles)
1304 for (
const ParticleFlowObject *
const pRootPfo : rootPfos)
1309 std::sort(mcNodes.begin(), mcNodes.end(),
1311 {
return lhs->
GetCaloHits().size() > rhs->GetCaloHits().size(); });
1312 std::sort(recoNodes.begin(), recoNodes.end(),
1314 {
return lhs->
GetCaloHits().size() > rhs->GetCaloHits().size(); });
1318 const CaloHitList &recoHits{pRecoNode->GetCaloHits()};
1320 size_t bestSharedHits{0};
1323 if (!pMCNode->IsReconstructable())
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));
1329 if (!intersection.empty())
1331 const size_t sharedHits{intersection.size()};
1332 if (sharedHits > bestSharedHits)
1334 bestSharedHits = sharedHits;
1335 pBestNode = pMCNode;
1341 auto iter{mcToMatchMap.find(pBestNode)};
1342 if (iter != mcToMatchMap.end())
1345 match.
AddRecoMatch(pRecoNode, static_cast<int>(bestSharedHits));
1350 match.
AddRecoMatch(pRecoNode, static_cast<int>(bestSharedHits));
1351 mcToMatchMap.insert(std::make_pair(pBestNode, match));
1362 for (
auto [pMCNode, matches] : mcToMatchMap)
1365 for (
const MCParticle *
const pRootMC : rootMCParticles)
1369 if (std::find(mcNodes.begin(), mcNodes.end(), pMCNode) != mcNodes.end())
1371 m_matches[pRootMC].emplace_back(matches);
1378 {
return lhs.
GetMC()->
GetCaloHits().size() > rhs.GetMC()->GetCaloHits().size(); };
1380 for (
const MCParticle *
const pRootMC : rootMCParticles)
1389 if (pMCNode->IsReconstructable() && mcToMatchMap.find(pMCNode) == mcToMatchMap.end())
1403 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1405 return static_cast<unsigned int>(
m_matches.at(pRoot).size());
1413 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1415 unsigned int nNodes{0};
1419 if (!(pNode->IsCosmicRay() || pNode->IsTestBeamParticle()))
1431 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1433 unsigned int nNodes{0};
1437 if (pNode->IsCosmicRay())
1449 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
1451 unsigned int nNodes{0};
1455 if (pNode->IsTestBeamParticle())
1467 MCParticleList rootMCParticles;
1470 for (
const MCParticle *
const pRootMC : rootMCParticles)
1474 MCParticleList primaries;
1478 if (pMCNode->GetHierarchyTier() == 1)
1481 primaries.emplace_back(pLeadingMC);
1489 std::cout <<
"=== MC Interaction : PDG " <<
std::to_string(pLArRoot->GetParticleId())
1490 <<
" Energy: " <<
std::to_string(pLArRoot->GetEnergy()) <<
" Type: " << descriptor.ToString() << std::endl;
1492 std::cout <<
"=== MC Interaction : PDG " <<
std::to_string(pRootMC->GetParticleId())
1493 <<
" Energy: " <<
std::to_string(pRootMC->GetEnergy()) <<
" Type: " << descriptor.ToString() << std::endl;
1495 unsigned int nNeutrinoMCParticles{this->
GetNNeutrinoMCNodes(pRootMC)}, nNeutrinoRecoParticles{0};
1497 unsigned int nTestBeamMCParticles{this->
GetNTestBeamMCNodes(pRootMC)}, nTestBeamRecoParticles{0};
1498 std::cout <<
" === Matches ===" << std::endl;
1499 std::cout << std::fixed << std::setprecision(2);
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;
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.1
f)
1516 std::cout <<
" Matched " << sharedHits <<
" out of " << recoHits <<
" with purity " << purity <<
" and completeness " 1517 << completeness << std::endl;
1519 std::cout <<
" (Below threshold) " << sharedHits <<
" out of " << recoHits <<
" with purity " << purity
1520 <<
" and completeness " << completeness << std::endl;
1522 if (nodeVector.empty())
1524 std::cout <<
" Unmatched" << std::endl;
1526 else if (match.IsQuality(this->GetQualityCuts()))
1528 if (pMCNode->IsTestBeamParticle())
1529 ++nTestBeamRecoParticles;
1530 else if (pMCNode->IsCosmicRay())
1531 ++nCosmicRecoParticles;
1533 ++nNeutrinoRecoParticles;
1539 std::cout <<
" Neutrino Interaction Summary:" << std::endl;
1540 if (nNeutrinoMCParticles)
1542 std::cout <<
" Good final state particles: " << nNeutrinoRecoParticles <<
" of " << nNeutrinoMCParticles <<
" : " 1543 << (100 * nNeutrinoRecoParticles /
static_cast<float>(nNeutrinoMCParticles)) <<
"%" << std::endl;
1548 std::cout <<
" Cosmic Ray Interaction Summary:" << std::endl;
1549 std::cout << std::fixed << std::setprecision(1);
1550 if (nCosmicMCParticles)
1552 std::cout <<
" Good cosmics: " << nCosmicRecoParticles <<
" of " << nCosmicMCParticles <<
" : " 1553 << (100 * nCosmicRecoParticles /
static_cast<float>(nCosmicMCParticles)) <<
"%" << std::endl;
1558 std::cout <<
" Test Beam Interaction Summary:" << std::endl;
1559 std::cout << std::fixed << std::setprecision(1);
1560 if (nTestBeamMCParticles)
1562 std::cout <<
" Good test beam particles: " << nTestBeamRecoParticles <<
" of " << nTestBeamMCParticles <<
" : " 1563 << (100 * nTestBeamRecoParticles /
static_cast<float>(nTestBeamMCParticles)) <<
"%" << std::endl;
1565 if (nCosmicMCParticles)
1567 std::cout <<
" Matched cosmics: " << nCosmicRecoParticles <<
" of " << nCosmicMCParticles <<
" : " 1568 << (100 * nCosmicRecoParticles /
static_cast<float>(nCosmicMCParticles)) <<
"%" << std::endl;
1572 std::cout <<
" Unmatched reco: " << this->
GetUnmatchedReco().size() << std::endl;
1581 rootMCParticles.emplace_back(iter->first);
1590 hierarchy.
FillHierarchy(mcParticleList, caloHitList, foldParameters);
1614 MCParticleList visible;
1617 for (
const MCParticle *
const pMCParticle : visible)
1618 primaries.insert(pMCParticle);
1620 catch (
const StatusCodeException &)
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;
1632 const PfoList &children{pRoot->GetDaughterPfoList()};
1634 for (
const ParticleFlowObject *
const pPfo : children)
1635 primaries.insert(pPfo);
1641 primaries.insert(pRoot);
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
FoldingParameters()
Default constructor.
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.
std::vector< const Node * > NodeVector
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.
MCHierarchy()
Default constructor.
constexpr auto abs(T v)
Returns the absolute value of the argument.
ReconstructabilityCriteria()
Default constructor.
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.
virtual ~RecoHierarchy()
Destructor.
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.
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.
MCMatches(const MCHierarchy::Node *pMCParticle)
Constructor.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
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 ...
virtual ~Node()
Destructor.
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.
std::vector< const Node * > NodeVector
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.
ReconstructabilityCriteria class.
bool IsCosmicRay() const
Check if this is a cosmic ray particle.
std::list< const Node * > NodeList
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.
virtual ~MCHierarchy()
Destructor.
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.
virtual ~Node()
Destructor.
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.
QualityCuts()
Default constructor.
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.
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.
std::list< const Node * > NodeList
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.
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.
RecoHierarchy()
Default constructor.
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.