VertexRefinementAlgorithm class. More...

#include <VertexRefinementAlgorithm.h>

Inheritance diagram for lar_content::VertexRefinementAlgorithm:

Public Member Functions
	VertexRefinementAlgorithm ()
	Default constructor. More...

Private Member Functions
pandora::StatusCode	Run ()

void	GetClusterLists (const pandora::StringVector &inputClusterListNames, pandora::ClusterList &clusterListU, pandora::ClusterList &clusterListV, pandora::ClusterList &clusterListW) const
	Get the input cluster lists. More...

void	RefineVertices (const pandora::VertexList *const pVertexList, const pandora::ClusterList &clusterListU, const pandora::ClusterList &clusterListV, const pandora::ClusterList &clusterListW) const
	Perform the refinement proceduce on a list of vertices. More...

pandora::CartesianVector	RefineVertexTwoD (const pandora::ClusterList &clusterList, const pandora::CartesianVector &originalVtxPos) const
	Refine the position of a two dimensional projection of a vertex using the clusters in that view. More...

void	GetBestFitPoint (const pandora::CartesianPointVector &intercepts, const pandora::CartesianPointVector &directions, const pandora::FloatVector &weights, pandora::CartesianVector &bestFitPoint) const
	Calculate the best fit point of a set of lines using a matrix equation. More...

pandora::StatusCode	ReadSettings (const pandora::TiXmlHandle xmlHandle)

Private Attributes
pandora::StringVector	m_inputClusterListNames
	The list of input cluster list names. More...

std::string	m_inputVertexListName
	The initial vertex list. More...

std::string	m_outputVertexListName
	The refined vertex list to be outputted. More...

float	m_chiSquaredCut
	The maximum chi2 value a refined vertex can have to be kept. More...

float	m_distanceCut
	The maximum distance a refined vertex can be from the original position to be kept. More...

unsigned int	m_minimumHitsCut
	The minimum size of a cluster to be used in refinement. More...

float	m_twoDDistanceCut
	The maximum distance a cluster can be from the original position to be used in refinement. More...

Detailed Description

VertexRefinementAlgorithm class.

Definition at line 19 of file VertexRefinementAlgorithm.h.

Constructor & Destructor Documentation

lar_content::VertexRefinementAlgorithm::VertexRefinementAlgorithm ( )

Default constructor.

Definition at line 24 of file VertexRefinementAlgorithm.cc.

                                                      :
     m_chiSquaredCut(2.f),
     m_distanceCut(5.f),
     m_minimumHitsCut(5),
     m_twoDDistanceCut(10.f)
 {
 }

Member Function Documentation

void lar_content::VertexRefinementAlgorithm::GetBestFitPoint	(	const pandora::CartesianPointVector &	intercepts,
		const pandora::CartesianPointVector &	directions,
		const pandora::FloatVector &	weights,
		pandora::CartesianVector &	bestFitPoint
	)		const

private

Calculate the best fit point of a set of lines using a matrix equation.

Parameters

intercepts	the vector of the defining points of the lines
directions	the vector of line directions
weights	the vector of weights for each line
bestFitPoint	the resulting best fit point

Definition at line 195 of file VertexRefinementAlgorithm.cc.

 {
     const int n(intercepts.size());
 
     Eigen::VectorXd d = Eigen::VectorXd::Zero(3 * n);
     Eigen::MatrixXd G = Eigen::MatrixXd::Zero(3 * n, n + 3);
 
     for (int i = 0; i < n; ++i)
     {
         d(3 * i) = intercepts[i].GetX() * weights[i];
         d(3 * i + 1) = intercepts[i].GetY() * weights[i];
         d(3 * i + 2) = intercepts[i].GetZ() * weights[i];
 
         G(3 * i, 0) = weights[i];
         G(3 * i + 1, 1) = weights[i];
         G(3 * i + 2, 2) = weights[i];
 
         G(3 * i, i + 3) = -directions[i].GetX();
         G(3 * i + 1, i + 3) = -directions[i].GetY();
         G(3 * i + 2, i + 3) = -directions[i].GetZ();
     }
 
     if ((G.transpose() * G).determinant() < std::numeric_limits<float>::epsilon())
     {
         bestFitPoint = CartesianVector(std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max());
         return;
     }
 
     Eigen::VectorXd m = (G.transpose() * G).inverse() * G.transpose() * d;
 
     bestFitPoint = CartesianVector(m[0], m[1], m[2]);
 }

void lar_content::VertexRefinementAlgorithm::GetClusterLists	(	const pandora::StringVector &	inputClusterListNames,
		pandora::ClusterList &	clusterListU,
		pandora::ClusterList &	clusterListV,
		pandora::ClusterList &	clusterListW
	)		const

private

Get the input cluster lists.

Parameters

inputClusterListNames	the input cluster list names
clusterListU	the U-view cluster list to populate
clusterListV	the V-view cluster list to populate
clusterListW	the W-view cluster list to populate

Definition at line 67 of file VertexRefinementAlgorithm.cc.

 {
     for (const std::string &clusterListName : inputClusterListNames)
     {
         const ClusterList *pClusterList(nullptr);
         PANDORA_THROW_RESULT_IF_AND_IF(
             STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, clusterListName, pClusterList));
 
         if (!pClusterList || pClusterList->empty())
         {
             if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
                 std::cout << "VertexRefinementAlgorithm: unable to find cluster list " << clusterListName << std::endl;
 
             continue;
         }
 
         const HitType hitType(LArClusterHelper::GetClusterHitType(pClusterList->front()));
 
         if ((TPC_VIEW_U != hitType) && (TPC_VIEW_V != hitType) && (TPC_VIEW_W != hitType))
             throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
 
         ClusterList &outputClusterList((TPC_VIEW_U == hitType) ? clusterListU : (TPC_VIEW_V == hitType) ? clusterListV : clusterListW);
         outputClusterList.insert(outputClusterList.end(), pClusterList->begin(), pClusterList->end());
     }
 }

StatusCode lar_content::VertexRefinementAlgorithm::ReadSettings ( const pandora::TiXmlHandle xmlHandle )

private

Definition at line 231 of file VertexRefinementAlgorithm.cc.

 {
     PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "InputClusterListNames", m_inputClusterListNames));
 
     PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputVertexListName", m_inputVertexListName));
 
     PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputVertexListName", m_outputVertexListName));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ChiSquaredCut", m_chiSquaredCut));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DistanceCut", m_distanceCut));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinimumHitsCut", m_minimumHitsCut));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "TwoDDistanceCut", m_twoDDistanceCut));
 
     return STATUS_CODE_SUCCESS;
 }

CartesianVector lar_content::VertexRefinementAlgorithm::RefineVertexTwoD	(	const pandora::ClusterList &	clusterList,
		const pandora::CartesianVector &	originalVtxPos
	)		const

private

Refine the position of a two dimensional projection of a vertex using the clusters in that view.

Parameters

clusterList	the list of two dimensional clusters
originalVtxPos	the original vertex position projected into two dimensions

Returns: the new refined position

Definition at line 157 of file VertexRefinementAlgorithm.cc.

 {
     CartesianPointVector intercepts, directions;
     FloatVector weights;
 
     for (const Cluster *const pCluster : clusterList)
     {
         if (LArClusterHelper::GetClosestDistance(originalVtxPos, pCluster) > 10)
             continue;
 
         if (pCluster->GetNCaloHits() < m_minimumHitsCut)
             continue;
 
         CartesianVector centroid(0.f, 0.f, 0.f);
         LArPcaHelper::EigenValues eigenValues(0.f, 0.f, 0.f);
         LArPcaHelper::EigenVectors eigenVectors;
         CartesianPointVector pointVector;
 
         LArClusterHelper::GetCoordinateVector(pCluster, pointVector);
         LArPcaHelper::RunPca(pointVector, centroid, eigenValues, eigenVectors);
 
         intercepts.push_back(LArClusterHelper::GetClosestPosition(originalVtxPos, pCluster));
         directions.push_back(eigenVectors.at(0).GetUnitVector());
         weights.push_back(1.f / ((LArClusterHelper::GetClosestPosition(originalVtxPos, pCluster) - originalVtxPos).GetMagnitudeSquared() + 1));
     }
 
     CartesianVector newVtxPos(originalVtxPos);
     if (intercepts.size() > 1)
         GetBestFitPoint(intercepts, directions, weights, newVtxPos);
 
     if ((newVtxPos - originalVtxPos).GetMagnitude() > m_twoDDistanceCut)
         return originalVtxPos;
 
     return newVtxPos;
 }

void lar_content::VertexRefinementAlgorithm::RefineVertices	(	const pandora::VertexList *const	pVertexList,
		const pandora::ClusterList &	clusterListU,
		const pandora::ClusterList &	clusterListV,
		const pandora::ClusterList &	clusterListW
	)		const

private

Perform the refinement proceduce on a list of vertices.

Parameters

pVertexList	address of the vertex list
clusterListU	the list of U-view clusters
clusterListV	the list of V-view clusters
clusterListW	the list of W-view clusters

Definition at line 96 of file VertexRefinementAlgorithm.cc.

 {
     for (const Vertex *const pVertex : *pVertexList)
     {
         const CartesianVector originalPosition(pVertex->GetPosition());
 
         const CartesianVector vtxU(
             this->RefineVertexTwoD(clusterListU, LArGeometryHelper::ProjectPosition(this->GetPandora(), originalPosition, TPC_VIEW_U)));
         const CartesianVector vtxV(
             this->RefineVertexTwoD(clusterListV, LArGeometryHelper::ProjectPosition(this->GetPandora(), originalPosition, TPC_VIEW_V)));
         const CartesianVector vtxW(
             this->RefineVertexTwoD(clusterListW, LArGeometryHelper::ProjectPosition(this->GetPandora(), originalPosition, TPC_VIEW_W)));
 
         CartesianVector vtxUV(0.f, 0.f, 0.f), vtxUW(0.f, 0.f, 0.f), vtxVW(0.f, 0.f, 0.f), vtx3D(0.f, 0.f, 0.f), position3D(0.f, 0.f, 0.f);
         float chi2UV(0.f), chi2UW(0.f), chi2VW(0.f), chi23D(0.f), chi2(0.f);
 
         LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), TPC_VIEW_U, TPC_VIEW_V, vtxU, vtxV, vtxUV, chi2UV);
         LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), TPC_VIEW_U, TPC_VIEW_W, vtxU, vtxW, vtxUW, chi2UW);
         LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), TPC_VIEW_V, TPC_VIEW_W, vtxV, vtxW, vtxVW, chi2VW);
         LArGeometryHelper::MergeThreePositions3D(this->GetPandora(), TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W, vtxU, vtxV, vtxW, vtx3D, chi23D);
 
         if (chi2UV < chi2UW && chi2UV < chi2VW && chi2UV < chi23D)
         {
             position3D = vtxUV;
             chi2 = chi2UV;
         }
         else if (chi2UW < chi2VW && chi2UW < chi23D)
         {
             position3D = vtxUW;
             chi2 = chi2UW;
         }
         else if (chi2VW < chi23D)
         {
             position3D = vtxVW;
             chi2 = chi2VW;
         }
         else
         {
             position3D = vtx3D;
             chi2 = chi23D;
         }
 
         if (chi2 > m_chiSquaredCut)
             position3D = originalPosition;
 
         if ((position3D - originalPosition).GetMagnitude() > m_distanceCut)
             position3D = originalPosition;
 
         PandoraContentApi::Vertex::Parameters parameters;
         parameters.m_position = position3D;
         parameters.m_vertexLabel = VERTEX_INTERACTION;
         parameters.m_vertexType = VERTEX_3D;
 
         const Vertex *pNewVertex(NULL);
         PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Vertex::Create(*this, parameters, pNewVertex));
     }
 }

StatusCode lar_content::VertexRefinementAlgorithm::Run ( )

private

Definition at line 34 of file VertexRefinementAlgorithm.cc.

 {
     const VertexList *pInputVertexList(nullptr);
     PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pInputVertexList));
 
     if (!pInputVertexList || pInputVertexList->empty())
     {
         if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
             std::cout << "VertexRefinementAlgorithm: unable to find current vertex list " << std::endl;
 
         return STATUS_CODE_SUCCESS;
     }
 
     const VertexList *pOutputVertexList(NULL);
     std::string temporaryListName;
     PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pOutputVertexList, temporaryListName));
 
     ClusterList clusterListU, clusterListV, clusterListW;
     this->GetClusterLists(m_inputClusterListNames, clusterListU, clusterListV, clusterListW);
 
     this->RefineVertices(pInputVertexList, clusterListU, clusterListV, clusterListW);
 
     if (!pOutputVertexList->empty())
     {
         PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Vertex>(*this, m_outputVertexListName));
         PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<Vertex>(*this, m_outputVertexListName));
     }
 
     return STATUS_CODE_SUCCESS;
 }