9 #include "Pandora/PdgTable.h" 27 void LArPfoHelper::GetCoordinateVector(
const ParticleFlowObject *
const pPfo,
const HitType &hitType, CartesianPointVector &coordinateVector)
29 ClusterList clusterList;
30 LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
32 for (
const Cluster *
const pCluster : clusterList)
33 LArClusterHelper::GetCoordinateVector(pCluster, coordinateVector);
38 void LArPfoHelper::GetCaloHits(
const PfoList &pfoList,
const HitType &hitType, CaloHitList &caloHitList)
40 for (
const ParticleFlowObject *
const pPfo : pfoList)
41 LArPfoHelper::GetCaloHits(pPfo, hitType, caloHitList);
46 void LArPfoHelper::GetCaloHits(
const ParticleFlowObject *
const pPfo,
const HitType &hitType, CaloHitList &caloHitList)
48 ClusterList clusterList;
49 LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
51 for (
const Cluster *
const pCluster : clusterList)
52 pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
57 void LArPfoHelper::GetIsolatedCaloHits(
const PfoList &pfoList,
const HitType &hitType, CaloHitList &caloHitList)
59 for (
const ParticleFlowObject *
const pPfo : pfoList)
60 LArPfoHelper::GetIsolatedCaloHits(pPfo, hitType, caloHitList);
65 void LArPfoHelper::GetIsolatedCaloHits(
const ParticleFlowObject *
const pPfo,
const HitType &hitType, CaloHitList &caloHitList)
67 ClusterList clusterList;
68 LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
70 for (
const Cluster *
const pCluster : clusterList)
71 caloHitList.insert(caloHitList.end(), pCluster->GetIsolatedCaloHitList().begin(), pCluster->GetIsolatedCaloHitList().end());
76 void LArPfoHelper::GetAllCaloHits(
const ParticleFlowObject *
const pPfo, CaloHitList &caloHitList)
78 LArPfoHelper::GetCaloHits(pPfo, TPC_VIEW_U, caloHitList);
79 LArPfoHelper::GetCaloHits(pPfo, TPC_VIEW_V, caloHitList);
80 LArPfoHelper::GetCaloHits(pPfo, TPC_VIEW_W, caloHitList);
81 LArPfoHelper::GetIsolatedCaloHits(pPfo, TPC_VIEW_U, caloHitList);
82 LArPfoHelper::GetIsolatedCaloHits(pPfo, TPC_VIEW_V, caloHitList);
83 LArPfoHelper::GetIsolatedCaloHits(pPfo, TPC_VIEW_W, caloHitList);
88 void LArPfoHelper::GetClusters(
const PfoList &pfoList,
const HitType &hitType, ClusterList &clusterList)
90 for (
const ParticleFlowObject *
const pPfo : pfoList)
91 LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
96 void LArPfoHelper::GetClusters(
const ParticleFlowObject *
const pPfo,
const HitType &hitType, ClusterList &clusterList)
98 for (
const Cluster *
const pCluster : pPfo->GetClusterList())
100 if (hitType != LArClusterHelper::GetClusterHitType(pCluster))
103 clusterList.push_back(pCluster);
109 unsigned int LArPfoHelper::GetNumberOfTwoDHits(
const ParticleFlowObject *
const pPfo)
111 unsigned int totalHits(0);
113 ClusterList clusterList;
114 LArPfoHelper::GetTwoDClusterList(pPfo, clusterList);
115 for (
const Cluster *
const pCluster : clusterList)
116 totalHits += pCluster->GetNCaloHits();
123 void LArPfoHelper::GetTwoDClusterList(
const ParticleFlowObject *
const pPfo, ClusterList &clusterList)
125 for (
const Cluster *
const pCluster : pPfo->GetClusterList())
127 if (TPC_3D == LArClusterHelper::GetClusterHitType(pCluster))
130 clusterList.push_back(pCluster);
136 void LArPfoHelper::GetThreeDClusterList(
const ParticleFlowObject *
const pPfo, ClusterList &clusterList)
138 for (
const Cluster *
const pCluster : pPfo->GetClusterList())
140 if (TPC_3D != LArClusterHelper::GetClusterHitType(pCluster))
143 clusterList.push_back(pCluster);
149 void LArPfoHelper::GetAllConnectedPfos(
const PfoList &inputPfoList, PfoList &outputPfoList)
151 for (
const ParticleFlowObject *
const pPfo : inputPfoList)
152 LArPfoHelper::GetAllConnectedPfos(pPfo, outputPfoList);
157 void LArPfoHelper::GetAllConnectedPfos(
const ParticleFlowObject *
const pPfo, PfoList &outputPfoList)
159 if (outputPfoList.end() != std::find(outputPfoList.begin(), outputPfoList.end(), pPfo))
162 outputPfoList.push_back(pPfo);
163 LArPfoHelper::GetAllConnectedPfos(pPfo->GetParentPfoList(), outputPfoList);
164 LArPfoHelper::GetAllConnectedPfos(pPfo->GetDaughterPfoList(), outputPfoList);
169 void LArPfoHelper::GetAllDownstreamPfos(
const PfoList &inputPfoList, PfoList &outputPfoList)
171 for (
const ParticleFlowObject *
const pPfo : inputPfoList)
172 LArPfoHelper::GetAllDownstreamPfos(pPfo, outputPfoList);
177 void LArPfoHelper::GetAllDownstreamPfos(
const ParticleFlowObject *
const pPfo, PfoList &outputPfoList)
179 if (outputPfoList.end() != std::find(outputPfoList.begin(), outputPfoList.end(), pPfo))
182 outputPfoList.push_back(pPfo);
183 LArPfoHelper::GetAllDownstreamPfos(pPfo->GetDaughterPfoList(), outputPfoList);
188 void LArPfoHelper::GetAllDownstreamPfos(
189 const pandora::ParticleFlowObject *
const pPfo, pandora::PfoList &outputTrackPfoList, pandora::PfoList &outputLeadingShowerPfoList)
191 if (LArPfoHelper::IsTrack(pPfo))
193 outputTrackPfoList.emplace_back(pPfo);
194 for (
const ParticleFlowObject *pChild : pPfo->GetDaughterPfoList())
196 if (std::find(outputTrackPfoList.begin(), outputTrackPfoList.end(), pChild) == outputTrackPfoList.end())
198 const int pdg{
std::abs(pChild->GetParticleId())};
201 outputLeadingShowerPfoList.emplace_back(pChild);
205 outputTrackPfoList.emplace_back(pChild);
206 LArPfoHelper::GetAllDownstreamPfos(pChild, outputTrackPfoList, outputLeadingShowerPfoList);
213 outputLeadingShowerPfoList.emplace_back(pPfo);
219 int LArPfoHelper::GetHierarchyTier(
const ParticleFlowObject *
const pPfo)
221 const ParticleFlowObject *pParentPfo = pPfo;
224 while (pParentPfo->GetParentPfoList().empty() ==
false)
226 if (1 != pParentPfo->GetParentPfoList().size())
227 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
229 pParentPfo = *(pParentPfo->GetParentPfoList().begin());
238 float LArPfoHelper::GetTwoDLengthSquared(
const ParticleFlowObject *
const pPfo)
240 if (!LArPfoHelper::IsTwoD(pPfo))
241 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
243 float lengthSquared(0.
f);
245 for (
const Cluster *
const pCluster : pPfo->GetClusterList())
247 if (TPC_3D == LArClusterHelper::GetClusterHitType(pCluster))
250 lengthSquared += LArClusterHelper::GetLengthSquared(pCluster);
253 return lengthSquared;
258 float LArPfoHelper::GetThreeDLengthSquared(
const ParticleFlowObject *
const pPfo)
260 if (!LArPfoHelper::IsThreeD(pPfo))
261 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
263 float lengthSquared(0.
f);
265 for (
const Cluster *
const pCluster : pPfo->GetClusterList())
267 if (TPC_3D != LArClusterHelper::GetClusterHitType(pCluster))
270 lengthSquared += LArClusterHelper::GetLengthSquared(pCluster);
273 return lengthSquared;
278 float LArPfoHelper::GetClosestDistance(
const ParticleFlowObject *
const pPfo,
const Cluster *
const pCluster)
280 const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));
282 ClusterList clusterList;
283 LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
285 if (clusterList.empty())
286 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
288 float bestDistance(std::numeric_limits<float>::max());
290 for (
const Cluster *
const pPfoCluster : clusterList)
292 const float thisDistance(LArClusterHelper::GetClosestDistance(pCluster, pPfoCluster));
294 if (thisDistance < bestDistance)
295 bestDistance = thisDistance;
303 float LArPfoHelper::GetThreeDSeparation(
const ParticleFlowObject *
const pPfo1,
const ParticleFlowObject *
const pPfo2)
305 ClusterList clusterList1, clusterList2;
307 LArPfoHelper::GetClusters(pPfo1, TPC_3D, clusterList1);
308 LArPfoHelper::GetClusters(pPfo2, TPC_3D, clusterList2);
310 if (clusterList1.empty() || clusterList2.empty())
311 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
313 return LArClusterHelper::GetClosestDistance(clusterList1, clusterList2);
318 bool LArPfoHelper::IsTwoD(
const ParticleFlowObject *
const pPfo)
320 for (
const Cluster *
const pCluster : pPfo->GetClusterList())
322 if (TPC_3D != LArClusterHelper::GetClusterHitType(pCluster))
331 bool LArPfoHelper::IsThreeD(
const ParticleFlowObject *
const pPfo)
333 for (
const Cluster *
const pCluster : pPfo->GetClusterList())
335 if (TPC_3D == LArClusterHelper::GetClusterHitType(pCluster))
344 bool LArPfoHelper::IsTrack(
const ParticleFlowObject *
const pPfo)
346 const int pdg(pPfo->GetParticleId());
354 bool LArPfoHelper::IsShower(
const ParticleFlowObject *
const pPfo)
356 const int pdg(pPfo->GetParticleId());
364 int LArPfoHelper::GetPrimaryNeutrino(
const ParticleFlowObject *
const pPfo)
368 const ParticleFlowObject *
const pParentPfo = LArPfoHelper::GetParentNeutrino(pPfo);
369 return pParentPfo->GetParticleId();
371 catch (
const StatusCodeException &)
379 bool LArPfoHelper::IsFinalState(
const ParticleFlowObject *
const pPfo)
381 if (pPfo->GetParentPfoList().empty() && !LArPfoHelper::IsNeutrino(pPfo))
384 if (LArPfoHelper::IsNeutrinoFinalState(pPfo))
392 bool LArPfoHelper::IsNeutrinoFinalState(
const ParticleFlowObject *
const pPfo)
394 return ((pPfo->GetParentPfoList().size() == 1) && (LArPfoHelper::IsNeutrino(*(pPfo->GetParentPfoList().begin()))));
399 bool LArPfoHelper::IsNeutrino(
const ParticleFlowObject *
const pPfo)
401 const int absoluteParticleId(
std::abs(pPfo->GetParticleId()));
403 if ((NU_E == absoluteParticleId) || (NU_MU == absoluteParticleId) || (NU_TAU == absoluteParticleId))
411 bool LArPfoHelper::IsTestBeamFinalState(
const ParticleFlowObject *
const pPfo)
413 return LArPfoHelper::IsTestBeam(LArPfoHelper::GetParentPfo(pPfo));
418 bool LArPfoHelper::IsTestBeam(
const ParticleFlowObject *
const pPfo)
420 const PropertiesMap &properties(pPfo->GetPropertiesMap());
423 if (iter != properties.end())
424 return ((iter->second > 0.f) ?
true :
false);
431 void LArPfoHelper::GetRecoNeutrinos(
const PfoList *
const pPfoList, PfoList &recoNeutrinos)
436 for (
const ParticleFlowObject *
const pPfo : *pPfoList)
438 if (LArPfoHelper::IsNeutrino(pPfo))
439 recoNeutrinos.push_back(pPfo);
445 const ParticleFlowObject *LArPfoHelper::GetParentPfo(
const ParticleFlowObject *
const pPfo)
447 const ParticleFlowObject *pParentPfo = pPfo;
449 while (pParentPfo->GetParentPfoList().empty() ==
false)
451 pParentPfo = *(pParentPfo->GetParentPfoList().begin());
459 const ParticleFlowObject *LArPfoHelper::GetParentNeutrino(
const ParticleFlowObject *
const pPfo)
461 const ParticleFlowObject *
const pParentPfo = LArPfoHelper::GetParentPfo(pPfo);
463 if (!LArPfoHelper::IsNeutrino(pParentPfo))
464 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
471 const Vertex *LArPfoHelper::GetVertex(
const ParticleFlowObject *
const pPfo)
473 if (pPfo->GetVertexList().empty())
474 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
476 const Vertex *pVertex(
nullptr);
479 if (LArPfoHelper::IsTestBeam(pPfo) && pPfo->GetParentPfoList().empty())
481 pVertex = LArPfoHelper::GetVertexWithLabel(pPfo->GetVertexList(), VERTEX_START);
485 if (pPfo->GetVertexList().size() != 1)
486 throw StatusCodeException(STATUS_CODE_FAILURE);
488 pVertex = *(pPfo->GetVertexList().begin());
491 if (VERTEX_3D != pVertex->GetVertexType())
492 throw StatusCodeException(STATUS_CODE_FAILURE);
499 const Vertex *LArPfoHelper::GetTestBeamInteractionVertex(
const ParticleFlowObject *
const pPfo)
501 if (pPfo->GetVertexList().empty() || !pPfo->GetParentPfoList().empty() || !LArPfoHelper::IsTestBeam(pPfo))
502 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
504 const Vertex *pInteractionVertex(LArPfoHelper::GetVertexWithLabel(pPfo->GetVertexList(), VERTEX_INTERACTION));
506 if (VERTEX_3D != pInteractionVertex->GetVertexType())
507 throw StatusCodeException(STATUS_CODE_FAILURE);
509 return pInteractionVertex;
514 const Vertex *LArPfoHelper::GetVertexWithLabel(
const VertexList &vertexList,
const VertexLabel vertexLabel)
516 const Vertex *pTargetVertex(
nullptr);
518 for (
const Vertex *pCandidateVertex : vertexList)
520 if (pCandidateVertex->GetVertexLabel() == vertexLabel)
524 pTargetVertex = pCandidateVertex;
528 throw StatusCodeException(STATUS_CODE_FAILURE);
534 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
536 return pTargetVertex;
541 void LArPfoHelper::GetSlidingFitTrajectory(
const CartesianPointVector &pointVector,
const CartesianVector &vertexPosition,
544 LArPfoHelper::SlidingFitTrajectoryImpl(&pointVector, vertexPosition, layerWindow, layerPitch, trackStateVector, pIndexVector);
549 void LArPfoHelper::GetSlidingFitTrajectory(
const ParticleFlowObject *
const pPfo,
const Vertex *
const pVertex,
550 const unsigned int layerWindow,
const float layerPitch,
LArTrackStateVector &trackStateVector)
552 CaloHitList caloHitList;
553 LArPfoHelper::GetCaloHits(pPfo, TPC_3D, caloHitList);
554 LArPfoHelper::SlidingFitTrajectoryImpl(&caloHitList, pVertex->GetPosition(), layerWindow, layerPitch, trackStateVector);
559 LArShowerPCA LArPfoHelper::GetPrincipalComponents(
const CartesianPointVector &pointVector,
const CartesianVector &vertexPosition)
562 CartesianVector centroid(0.
f, 0.
f, 0.
f);
565 LArPcaHelper::RunPca(pointVector, centroid, eigenValues, eigenVecs);
568 if (eigenValues.GetX() < std::numeric_limits<float>::epsilon())
569 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
572 const float testProjection(eigenVecs.at(0).GetDotProduct(vertexPosition - centroid));
573 const float directionScaleFactor((testProjection > std::numeric_limits<float>::epsilon()) ? -1.
f : 1.
f);
575 const CartesianVector primaryAxis(eigenVecs.at(0) * directionScaleFactor);
576 const CartesianVector secondaryAxis(eigenVecs.at(1) * directionScaleFactor);
577 const CartesianVector tertiaryAxis(eigenVecs.at(2) * directionScaleFactor);
579 return LArShowerPCA(centroid, primaryAxis, secondaryAxis, tertiaryAxis, eigenValues);
584 LArShowerPCA LArPfoHelper::GetPrincipalComponents(
const ParticleFlowObject *
const pPfo,
const Vertex *
const pVertex)
586 CartesianPointVector pointVector;
587 LArPfoHelper::GetCoordinateVector(pPfo, TPC_3D, pointVector);
588 return LArPfoHelper::GetPrincipalComponents(pointVector, pVertex->GetPosition());
595 if (lhs.first != rhs.first)
596 return (lhs.first < rhs.first);
602 const float dx(lhs.second.GetPosition().GetX() - rhs.second.GetPosition().GetX());
603 const float dy(lhs.second.GetPosition().GetY() - rhs.second.GetPosition().GetY());
604 const float dz(lhs.second.GetPosition().GetZ() - rhs.second.GetPosition().GetZ());
605 return (dx + dy + dz > 0.
f);
610 bool LArPfoHelper::SortByNHits(
const ParticleFlowObject *
const pLhs,
const ParticleFlowObject *
const pRhs)
612 unsigned int nTwoDHitsLhs(0), nThreeDHitsLhs(0);
613 float energyLhs(0.
f);
614 for (
ClusterList::const_iterator iter = pLhs->GetClusterList().begin(), iterEnd = pLhs->GetClusterList().end(); iter != iterEnd; ++iter)
616 const Cluster *
const pClusterLhs = *iter;
618 if (TPC_3D != LArClusterHelper::GetClusterHitType(pClusterLhs))
619 nTwoDHitsLhs += pClusterLhs->GetNCaloHits();
621 nThreeDHitsLhs += pClusterLhs->GetNCaloHits();
623 energyLhs += pClusterLhs->GetHadronicEnergy();
626 unsigned int nTwoDHitsRhs(0), nThreeDHitsRhs(0);
627 float energyRhs(0.
f);
628 for (
ClusterList::const_iterator iter = pRhs->GetClusterList().begin(), iterEnd = pRhs->GetClusterList().end(); iter != iterEnd; ++iter)
630 const Cluster *
const pClusterRhs = *iter;
632 if (TPC_3D != LArClusterHelper::GetClusterHitType(pClusterRhs))
633 nTwoDHitsRhs += pClusterRhs->GetNCaloHits();
635 nThreeDHitsRhs += pClusterRhs->GetNCaloHits();
637 energyRhs += pClusterRhs->GetHadronicEnergy();
640 if (nTwoDHitsLhs != nTwoDHitsRhs)
641 return (nTwoDHitsLhs > nTwoDHitsRhs);
643 if (nThreeDHitsLhs != nThreeDHitsRhs)
644 return (nThreeDHitsLhs > nThreeDHitsRhs);
647 return (energyLhs > energyRhs);
652 void LArPfoHelper::GetBreadthFirstHierarchyRepresentation(
const pandora::ParticleFlowObject *
const pPfo, pandora::PfoList &pfoList)
654 const ParticleFlowObject *pRoot{pPfo};
655 PfoList parents{pRoot->GetParentPfoList()};
656 while (!parents.empty())
658 if (parents.size() > 1)
659 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
660 pRoot = parents.front();
661 parents = pRoot->GetParentPfoList();
664 pfoList.emplace_back(pRoot);
665 queue.emplace_back(pRoot);
667 while (!queue.empty())
669 const PfoList &daughters{queue.front()->GetDaughterPfoList()};
671 for (
const ParticleFlowObject *pDaughter : daughters)
673 pfoList.emplace_back(pDaughter);
674 queue.emplace_back(pDaughter);
681 template <
typename T>
682 void LArPfoHelper::SlidingFitTrajectoryImpl(
const T *
const pT,
const CartesianVector &vertexPosition,
const unsigned int layerWindow,
685 CartesianPointVector pointVector;
687 for (
const auto &nextPoint : *pT)
688 pointVector.push_back(LArObjectHelper::TypeAdaptor::GetPosition(nextPoint));
690 if (pointVector.empty())
691 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
693 std::sort(pointVector.begin(), pointVector.end(), LArClusterHelper::SortCoordinatesByPosition);
698 pIndexVector->clear();
706 if ((maxPosition - minPosition).GetMagnitudeSquared() < std::numeric_limits<float>::epsilon())
707 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
709 const CartesianVector seedPosition((maxPosition + minPosition) * 0.5
f);
710 const CartesianVector seedDirection((maxPosition - minPosition).GetUnitVector());
712 const float scaleFactor((seedDirection.GetDotProduct(seedPosition - vertexPosition) > 0.f) ? +1.f : -1.f);
715 for (
const auto &nextPoint : *pT)
723 CartesianVector position(0.f, 0.f, 0.f);
726 if (positionStatusCode != STATUS_CODE_SUCCESS)
727 throw StatusCodeException(positionStatusCode);
729 CartesianVector direction(0.f, 0.f, 0.f);
732 if (directionStatusCode != STATUS_CODE_SUCCESS)
733 throw StatusCodeException(directionStatusCode);
735 const float projection(seedDirection.GetDotProduct(position - seedPosition));
738 LArTrackState(position, direction * scaleFactor, LArObjectHelper::TypeAdaptor::GetCaloHit(nextPoint)), index));
740 catch (
const StatusCodeException &statusCodeException1)
742 indicesWithoutSpacePoints.push_back(index);
744 if (STATUS_CODE_FAILURE == statusCodeException1.GetStatusCode())
745 throw statusCodeException1;
749 catch (
const StatusCodeException &statusCodeException2)
751 if (STATUS_CODE_FAILURE == statusCodeException2.GetStatusCode())
752 throw statusCodeException2;
756 std::sort(trackTrajectory.begin(), trackTrajectory.end(), LArPfoHelper::SortByHitProjection);
760 trackStateVector.push_back(larTrackTrajectoryPoint.second);
762 pIndexVector->push_back(larTrackTrajectoryPoint.GetIndex());
768 for (
const int index : indicesWithoutSpacePoints)
769 pIndexVector->push_back(index);
775 template void LArPfoHelper::SlidingFitTrajectoryImpl(
777 template void LArPfoHelper::SlidingFitTrajectoryImpl(
Header file for the pfo helper class.
pandora::CartesianVector EigenValues
constexpr auto abs(T v)
Returns the absolute value of the argument.
Header file for the principal curve analysis helper class.
std::vector< int > IntVector
float GetLongitudinalDisplacement(const pandora::CartesianVector &position) const
Get longitudinal projection onto primary axis.
std::vector< LArTrackTrajectoryPoint > LArTrackTrajectory
Header file for the object helper class.
Header file for the cluster helper class.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
const pandora::CartesianVector & GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
Header file for the lar three dimensional sliding fit result class.
std::vector< pandora::CartesianVector > EigenVectors
ThreeDSlidingFitResult class.
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
std::vector< LArTrackState > LArTrackStateVector
LArTrackTrajectoryPoint class.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
std::list< Vertex > VertexList
const pandora::CartesianVector & GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.