9 #include "Pandora/AlgorithmHeaders.h"
26 VertexBasedPfoMopUpAlgorithm::VertexBasedPfoMopUpAlgorithm() :
27 m_minVertexLongitudinalDistance(-2.5f),
28 m_maxVertexTransverseDistance(1.5f),
29 m_minVertexAssociatedHitTypes(2),
30 m_coneAngleCentile(0.8f),
31 m_maxConeCosHalfAngle(0.95f),
32 m_maxConeLengthMultiplier(3.f),
33 m_directionTanAngle(1.732f),
34 m_directionApexShift(0.333f),
35 m_meanBoundedFractionCut(0.6f),
36 m_maxBoundedFractionCut(0.7f),
37 m_minBoundedFractionCut(0.3f),
38 m_minConsistentDirections(2),
39 m_minConsistentDirectionsTrack(3)
48 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetCurrentList(*
this, pVertexList));
50 const Vertex *
const pSelectedVertex(
51 (pVertexList && (pVertexList->size() == 1) && (VERTEX_3D == (*(pVertexList->begin()))->GetVertexType())) ? *(pVertexList->begin()) :
nullptr);
55 if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
56 std::cout <<
"VertexBasedPfoMopUp: unable to find vertex in current list " << std::endl;
58 return STATUS_CODE_SUCCESS;
63 PfoList vertexPfos, nonVertexPfos;
64 this->
GetInputPfos(pSelectedVertex, vertexPfos, nonVertexPfos);
67 this->
GetPfoAssociations(pSelectedVertex, vertexPfos, nonVertexPfos, pfoAssociationList);
69 std::sort(pfoAssociationList.begin(), pfoAssociationList.end());
76 return STATUS_CODE_SUCCESS;
90 const Pfo *
const pVertexPfo,
const Pfo *
const pDaughterPfo, HitTypeToAssociationMap &hitTypeToAssociationMap)
const
92 if ((pVertexPfo->GetClusterList().size() != pDaughterPfo->GetClusterList().size()) || (3 != pVertexPfo->GetClusterList().size()))
93 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
95 return PfoAssociation(pVertexPfo, pDaughterPfo, hitTypeToAssociationMap.at(TPC_VIEW_U), hitTypeToAssociationMap.at(TPC_VIEW_V),
96 hitTypeToAssociationMap.at(TPC_VIEW_W));
103 StringVector listNames;
107 for (
const std::string &listName : listNames)
109 const PfoList *pPfoList(
nullptr);
111 if (STATUS_CODE_SUCCESS != PandoraContentApi::GetList(*
this, listName, pPfoList))
114 for (
const Pfo *
const pPfo : *pPfoList)
116 PfoList &pfoTargetList(this->
IsVertexAssociated(pPfo, pVertex) ? vertexPfos : nonVertexPfos);
117 pfoTargetList.push_back(pPfo);
126 if (VERTEX_3D != pVertex->GetVertexType())
127 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
131 for (
const Cluster *
const pCluster : pPfo->GetClusterList())
135 if ((TPC_VIEW_U != hitType) && (TPC_VIEW_V != hitType) && (TPC_VIEW_W != hitType))
145 hitTypeSet.insert(hitType);
147 catch (StatusCodeException &)
152 const unsigned int nVertexAssociatedHitTypes(hitTypeSet.size());
159 const Vertex *
const pVertex,
const PfoList &vertexPfos,
const PfoList &nonVertexPfos,
PfoAssociationList &pfoAssociationList)
const
161 for (
const Pfo *
const pVertexPfo : vertexPfos)
163 for (
const Pfo *
const pDaughterPfo : nonVertexPfos)
168 pfoAssociationList.push_back(pfoAssociation);
170 catch (StatusCodeException &)
180 const Vertex *
const pVertex,
const Pfo *
const pVertexPfo,
const Pfo *
const pDaughterPfo)
const
182 if (pVertexPfo->GetClusterList().empty() || pDaughterPfo->GetClusterList().empty())
183 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
187 for (
const Cluster *
const pVertexCluster : pVertexPfo->GetClusterList())
191 for (
const Cluster *
const pDaughterCluster : pDaughterPfo->GetClusterList())
195 if (vertexHitType != daughterHitType)
199 hitTypeToAssociationMap[vertexHitType] = clusterAssociation;
203 return this->
GetPfoAssociation(pVertexPfo, pDaughterPfo, hitTypeToAssociationMap);
209 const Vertex *
const pVertex,
const Cluster *
const pVertexCluster,
const Cluster *
const pDaughterCluster)
const
221 const bool isConsistentDirection(vertexClusterDirection == daughterClusterDirection);
223 return ClusterAssociation(pVertexCluster, pDaughterCluster, boundedFraction, isConsistentDirection);
230 const PfoList *pTrackPfoList(
nullptr);
244 if ((pTrackPfoList->end() != std::find(pTrackPfoList->begin(), pTrackPfoList->end(), pfoAssociation.GetVertexPfo())) &&
245 (pTrackPfoList->end() != std::find(pTrackPfoList->begin(), pTrackPfoList->end(), pfoAssociation.GetDaughterPfo())))
250 if (((pTrackPfoList->end() != std::find(pTrackPfoList->begin(), pTrackPfoList->end(), pfoAssociation.GetVertexPfo())) ||
251 (pTrackPfoList->end() != std::find(pTrackPfoList->begin(), pTrackPfoList->end(), pfoAssociation.GetDaughterPfo()))) &&
269 const PfoList *pTrackPfoList(
nullptr), *pShowerPfoList(
nullptr);
273 if (!pTrackPfoList && !pShowerPfoList)
274 throw StatusCodeException(STATUS_CODE_FAILURE);
276 const Pfo *
const pVertexPfo(pfoAssociation.
GetVertexPfo());
277 const bool isvertexTrack(pTrackPfoList && (pTrackPfoList->end() != std::find(pTrackPfoList->begin(), pTrackPfoList->end(), pVertexPfo)));
279 const bool isDaughterShower(pShowerPfoList && (pShowerPfoList->end() != std::find(pShowerPfoList->begin(), pShowerPfoList->end(), pDaughterPfo)));
283 if (isvertexTrack && isDaughterShower)
285 const PfoList vertexPfoList(1, pVertexPfo);
287 PandoraContentApi::ParticleFlowObject::Metadata metadata;
288 metadata.m_particleId = E_MINUS;
289 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::AlterMetadata(*
this, pVertexPfo, metadata));
299 m_pVertexCluster(nullptr),
300 m_pDaughterCluster(nullptr),
301 m_boundedFraction(0.f),
302 m_isConsistentDirection(
false)
309 const Cluster *
const pDaughterCluster,
const float boundedFraction,
const bool isConsistentDirection) :
310 m_pVertexCluster(pVertexCluster),
311 m_pDaughterCluster(pDaughterCluster),
312 m_boundedFraction(boundedFraction),
313 m_isConsistentDirection(isConsistentDirection)
322 m_pVertexPfo(pVertexPfo),
323 m_pDaughterPfo(pDaughterPfo),
324 m_clusterAssociationU(clusterAssociationU),
325 m_clusterAssociationV(clusterAssociationV),
326 m_clusterAssociationW(clusterAssociationW)
334 return ((this->GetClusterAssociationU().GetBoundedFraction() + this->GetClusterAssociationV().GetBoundedFraction() +
335 this->GetClusterAssociationW().GetBoundedFraction()) /
343 return std::max(this->GetClusterAssociationU().GetBoundedFraction(),
344 std::max(this->GetClusterAssociationV().GetBoundedFraction(), this->GetClusterAssociationW().GetBoundedFraction()));
351 return std::min(this->GetClusterAssociationU().GetBoundedFraction(),
352 std::min(this->GetClusterAssociationV().GetBoundedFraction(), this->GetClusterAssociationW().GetBoundedFraction()));
359 unsigned int nConsistentDirections(0);
361 if (this->GetClusterAssociationU().IsConsistentDirection())
362 ++nConsistentDirections;
364 if (this->GetClusterAssociationV().IsConsistentDirection())
365 ++nConsistentDirections;
367 if (this->GetClusterAssociationW().IsConsistentDirection())
368 ++nConsistentDirections;
370 return nConsistentDirections;
377 if (std::fabs(this->GetMeanBoundedFraction() - rhs.
GetMeanBoundedFraction()) > std::numeric_limits<float>::epsilon())
390 const Cluster *
const pCluster,
const CartesianVector &vertexPosition2D,
const float coneAngleCentile,
const float maxCosHalfAngle) :
391 m_pCluster(pCluster),
392 m_apex(vertexPosition2D),
393 m_direction(0.f, 0.f, 0.f),
395 m_coneCosHalfAngle(0.f)
401 if (
m_coneLength < std::numeric_limits<float>::epsilon())
414 unsigned int nMatchedHits(0);
415 const OrderedCaloHitList &orderedCaloHitList(pDaughterCluster->GetOrderedCaloHitList());
417 for (OrderedCaloHitList::const_iterator iter = orderedCaloHitList.begin(), iterEnd = orderedCaloHitList.end(); iter != iterEnd; ++iter)
419 for (CaloHitList::const_iterator hIter = iter->second->begin(), hIterEnd = iter->second->end(); hIter != hIterEnd; ++hIter)
421 const CartesianVector &positionVector((*hIter)->GetPositionVector());
423 if (m_direction.GetCosOpeningAngle(positionVector - m_apex) < m_coneCosHalfAngle)
426 if (m_direction.GetDotProduct(positionVector - m_apex) > coneLengthMultiplier * m_coneLength)
433 if (0 == pDaughterCluster->GetNCaloHits())
434 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
436 return (static_cast<float>(nMatchedHits) /
static_cast<float>(pDaughterCluster->GetNCaloHits()));
443 const OrderedCaloHitList &orderedCaloHitList(m_pCluster->GetOrderedCaloHitList());
444 float sumDxDz(0.f), sumDxDx(0.f);
446 for (OrderedCaloHitList::const_iterator iter = orderedCaloHitList.begin(), iterEnd = orderedCaloHitList.end(); iter != iterEnd; ++iter)
448 for (CaloHitList::const_iterator hitIter = iter->second->begin(), hitIterEnd = iter->second->end(); hitIter != hitIterEnd; ++hitIter)
450 const CartesianVector apexDisplacement((*hitIter)->GetPositionVector() - m_apex);
451 sumDxDz += apexDisplacement.GetX() * apexDisplacement.GetZ();
452 sumDxDx += apexDisplacement.GetX() * apexDisplacement.GetX();
456 if (sumDxDx < std::numeric_limits<float>::epsilon())
457 return CartesianVector(0.f, 0.f, 1.f);
459 return CartesianVector(1.f, 0.f, sumDxDz / sumDxDx).GetUnitVector();
466 float maxProjectedLength(0.f);
467 const OrderedCaloHitList &orderedCaloHitList(m_pCluster->GetOrderedCaloHitList());
469 for (OrderedCaloHitList::const_iterator iter = orderedCaloHitList.begin(), iterEnd = orderedCaloHitList.end(); iter != iterEnd; ++iter)
471 for (CaloHitList::const_iterator hitIter = iter->second->begin(), hitIterEnd = iter->second->end(); hitIter != hitIterEnd; ++hitIter)
473 const float projectedLength(m_direction.GetDotProduct((*hitIter)->GetPositionVector() - m_apex));
475 if (std::fabs(projectedLength) > std::fabs(maxProjectedLength))
476 maxProjectedLength = projectedLength;
480 return maxProjectedLength;
487 FloatVector halfAngleValues;
488 const OrderedCaloHitList &orderedCaloHitList(m_pCluster->GetOrderedCaloHitList());
490 for (OrderedCaloHitList::const_iterator iter = orderedCaloHitList.begin(), iterEnd = orderedCaloHitList.end(); iter != iterEnd; ++iter)
492 for (CaloHitList::const_iterator hitIter = iter->second->begin(), hitIterEnd = iter->second->end(); hitIter != hitIterEnd; ++hitIter)
493 halfAngleValues.push_back(m_direction.GetOpeningAngle((*hitIter)->GetPositionVector() - m_apex));
496 std::sort(halfAngleValues.begin(), halfAngleValues.end());
498 if (halfAngleValues.empty())
499 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
501 const unsigned int halfAngleBin(coneAngleCentile * halfAngleValues.size());
502 return std::cos(halfAngleValues.at(halfAngleBin));
510 PANDORA_RETURN_RESULT_IF_AND_IF(
511 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"TrackPfoListName",
m_trackPfoListName));
513 PANDORA_RETURN_RESULT_IF_AND_IF(
514 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"ShowerPfoListName",
m_showerPfoListName));
516 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
519 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
522 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
525 PANDORA_RETURN_RESULT_IF_AND_IF(
526 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"ConeAngleCentile",
m_coneAngleCentile));
528 PANDORA_RETURN_RESULT_IF_AND_IF(
529 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxConeCosHalfAngle",
m_maxConeCosHalfAngle));
531 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
534 PANDORA_RETURN_RESULT_IF_AND_IF(
535 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"DirectionTanAngle",
m_directionTanAngle));
537 PANDORA_RETURN_RESULT_IF_AND_IF(
538 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"DirectionApexShift",
m_directionApexShift));
540 PANDORA_RETURN_RESULT_IF_AND_IF(
541 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MeanBoundedFractionCut",
m_meanBoundedFractionCut));
543 PANDORA_RETURN_RESULT_IF_AND_IF(
544 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxBoundedFractionCut",
m_maxBoundedFractionCut));
546 PANDORA_RETURN_RESULT_IF_AND_IF(
547 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinBoundedFractionCut",
m_minBoundedFractionCut));
549 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
552 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
std::set< pandora::HitType > HitTypeSet
float GetBoundedFraction(const pandora::Cluster *const pDaughterCluster, const float coneLengthMultiplier) const
Get the fraction of hits in a candidate daughter cluster bounded by the cone.
void GetPfoAssociations(const pandora::Vertex *const pVertex, const pandora::PfoList &vertexPfos, const pandora::PfoList &nonVertexPfos, PfoAssociationList &pfoAssociationList) const
Get the list of associations between vertex-associated pfos and non-vertex-associated pfos...
static bool SortByNHits(const pandora::ParticleFlowObject *const pLhs, const pandora::ParticleFlowObject *const pRhs)
Sort pfos by number of constituent hits.
Header file for the pfo helper class.
unsigned int m_minVertexAssociatedHitTypes
The min number of vertex associated hit types for a vertex associated pfo.
float m_maxConeLengthMultiplier
Consider hits as bound if inside cone, with projected distance less than N times cone length...
ClusterAssociation class.
PfoAssociation(const pandora::Pfo *const pVertexPfo, const pandora::Pfo *const pDaughterPfo, const ClusterAssociation &clusterAssociationU, const ClusterAssociation &clusterAssociationV, const ClusterAssociation &clusterAssociationW)
Constructor.
ClusterAssociation GetClusterAssociation(const pandora::Vertex *const pVertex, const pandora::Cluster *const pVertexCluster, const pandora::Cluster *const pDaughterCluster) const
Get cluster association details between a vertex-associated cluster and a non-vertex associated daugh...
ConeParameters(const pandora::Cluster *const pCluster, const pandora::CartesianVector &vertexPosition2D, const float coneAngleCentile, const float maxConeCosHalfAngle)
Constructor.
Header file for the lar pointing cluster class.
unsigned int GetNConsistentDirections() const
Get the number of views for which the vertex and daughter cluster directions are consistent.
float m_coneAngleCentile
Cluster cone angle is defined using specified centile of distribution of hit half angles...
float GetCosHalfAngleEstimate(const float coneAngleCentile) const
Get the cone cos half angle estimate.
unsigned int m_minConsistentDirections
The minimum number of consistent cluster directions to allow a pfo merge.
pandora::StringVector m_daughterListNames
The list of potential daughter object list names.
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
pandora::CartesianVector m_direction
The cone direction.
void MergePfos(const PfoAssociation &pfoAssociation) const
Merge the vertex and daughter pfos (deleting daughter pfo, merging clusters, etc.) described in the s...
std::vector< PfoAssociation > PfoAssociationList
LArPointingCluster class.
float GetMaxBoundedFraction() const
Get the maximum bounded fraction from the u, v and w views.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
float m_maxVertexTransverseDistance
Vertex association check: max transverse distance cut.
float m_minVertexLongitudinalDistance
Vertex association check: min longitudinal distance cut.
float m_maxBoundedFractionCut
Cut on association info (max bounded fraction) for determining pfo merges.
float GetSignedConeLength() const
Get the cone length (signed, by projections of hits onto initial direction estimate) ...
bool ProcessPfoAssociations(const PfoAssociationList &pfoAssociationList) const
Process the list of pfo associations, merging the best-matching pfo.
ClusterAssociation()
Default constructor.
Header file for the geometry helper class.
std::map< pandora::HitType, ClusterAssociation > HitTypeToAssociationMap
Header file for the cluster helper class.
float m_directionApexShift
Direction determination, look for vertex inside triangle with apex shifted along the cluster length...
const Vertex & GetOuterVertex() const
Get the outer vertex.
float m_meanBoundedFractionCut
Cut on association info (mean bounded fraction) for determining pfo merges.
pandora::StatusCode Run()
float m_coneLength
The cone length.
const Vertex & GetInnerVertex() const
Get the inner vertex.
virtual PfoAssociation GetPfoAssociation(const pandora::Pfo *const pVertexPfo, const pandora::Pfo *const pDaughterPfo, HitTypeToAssociationMap &hitTypeToAssociationMap) const
Get pfo association details between a vertex-associated pfo and a non-vertex associated daughter cand...
float m_coneCosHalfAngle
The cone cos half angle.
pandora::CartesianVector GetDirectionEstimate() const
Get the cone direction estimate, with apex fixed at the 2d vertex position.
const pandora::Pfo * m_pDaughterPfo
The address of the non-vertex-associated candidate daughter pfo.
float m_directionTanAngle
Direction determination, look for vertex inside triangle with apex shifted along the cluster length...
void GetInputPfos(const pandora::Vertex *const pVertex, pandora::PfoList &vertexPfos, pandora::PfoList &nonVertexPfos) const
Get the list of input pfos and divide them into vertex-associated and non-vertex-associated lists...
Header file for the vertex helper class.
virtual bool IsVertexAssociated(const pandora::CartesianVector &vertex2D, const LArPointingCluster &pointingCluster) const
Whether a specified pfo is associated with a specified vertex.
unsigned int m_minConsistentDirectionsTrack
The minimum number of consistent cluster directions to allow a merge involving a track pfo...
required by fuzzyCluster table::sbnd_g4_services gaushitTruthMatch pandora
float m_maxConeCosHalfAngle
Maximum value for cosine of cone half angle.
const pandora::Pfo * GetDaughterPfo() const
Get the address of the non-vertex-associated candidate daughter pfo.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
std::string m_trackPfoListName
The input track pfo list name.
virtual void MergeAndDeletePfos(const pandora::ParticleFlowObject *const pPfoToEnlarge, const pandora::ParticleFlowObject *const pPfoToDelete) const
Merge and delete a pair of pfos, with a specific set of conventions for cluster merging, vertex use, etc.
ClusterAssociation class.
const pandora::Pfo * GetVertexPfo() const
Get the address of the vertex-associated pfo.
float GetMinBoundedFraction() const
Get the minimum bounded fraction from the u, v and w views.
const pandora::Pfo * m_pVertexPfo
The address of the vertex-associated pfo.
float GetMeanBoundedFraction() const
Get the mean bounded fraction, averaging over the u, v and w views.
bool operator<(const PfoAssociation &rhs) const
operator<
std::list< Vertex > VertexList
BEGIN_PROLOG could also be cout
Header file for the vertex based pfo mop up algorithm class.
static bool IsNode(const pandora::CartesianVector &parentVertex, const LArPointingCluster::Vertex &daughterVertex, const float minLongitudinalDistance, const float maxTransverseDistance)
Whether pointing vertex is adjacent to a given position.
static ClusterDirection GetClusterDirectionInZ(const pandora::Pandora &pandora, const pandora::Vertex *const pVertex, const pandora::Cluster *const pCluster, const float tanAngle, const float apexShift)
Get the direction of the cluster in z, using a projection of the provided vertex. ...
std::string m_showerPfoListName
The input shower pfo list name.
float m_minBoundedFractionCut
Cut on association info (min bounded fraction) for determining pfo merges.