CosmicRaySplittingAlgorithm class. More...

#include <CosmicRaySplittingAlgorithm.h>

Inheritance diagram for lar_content::CosmicRaySplittingAlgorithm:

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

Private Member Functions
pandora::StatusCode	Run ()

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

void	GetListOfCleanClusters (const pandora::ClusterList *const pClusterList, pandora::ClusterVector &clusterVector) const
	Populate cluster vector with subset of cluster list, containing clusters judged to be clean. More...

void	BuildSlidingFitResultMap (const pandora::ClusterVector &clusterVector, TwoDSlidingFitResultMap &slidingFitResultMap) const
	Build the map of sliding fit results. More...

pandora::StatusCode	FindBestSplitPosition (const TwoDSlidingFitResult &slidingFitResult, pandora::CartesianVector &splitPosition, pandora::CartesianVector &splitDirection1, pandora::CartesianVector &splitDirection2) const
	Find the position of greatest scatter along a sliding linear fit. More...

pandora::StatusCode	ConfirmSplitPosition (const TwoDSlidingFitResult &branchSlidingFitResult, const TwoDSlidingFitResult &replacementSlidingFitResult, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2, float &lengthSquared1, float &lengthSquared2) const
	Find a second replacement cluster that aligns with the scatter of the first branch cluster. More...

pandora::StatusCode	PerformSingleSplit (const pandora::Cluster const pBranchCluster, const pandora::Cluster const pReplacementCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &forwardDirection, const pandora::CartesianVector &backwardDirection) const
	Split a branch cluster for case of one replacement cluster. More...

pandora::StatusCode	PerformDoubleSplit (const pandora::Cluster const pBranchCluster, const pandora::Cluster const pReplacementCluster1, const pandora::Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2) const
	Split a branch cluster for case of two replacement clusters. More...

void	GetCaloHitListToMove (const pandora::Cluster const pBranchCluster, const pandora::Cluster const pReplacementCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &forwardDirection, const pandora::CartesianVector &backwardDirection, pandora::CaloHitList &caloHitList) const
	Get list of calo hits to move in order to split a branch cluster into two segments for case of one replacement cluster. More...

void	GetCaloHitListsToMove (const pandora::Cluster *const pBranchCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2, pandora::CaloHitList &caloHitList1, pandora::CaloHitList &caloHitList2) const
	Get lists of calo hits to move in order to split a branch cluster into two segments for case of two replacement clusters. More...

bool	IdentifyCrossedTracks (const pandora::Cluster const pBranchCluster, const pandora::Cluster const pReplacementCluster1, const pandora::Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition) const
	Identify crossed tracks formed from the branch cluster and its replacement cluster. More...

pandora::StatusCode	GetCaloHitListToKeep (const pandora::Cluster *const pBranchCluster, const pandora::CaloHitList &caloHitListToMove, pandora::CaloHitList &caloHitListToKeep) const
	Split the branch cluster and add hits to the replacement cluster. More...

pandora::StatusCode	SplitCluster (const pandora::Cluster const pBranchCluster, const pandora::Cluster const pReplacementCluster, const pandora::CaloHitList &caloHitListToMove) const
	Split the branch cluster and add hits to the replacement cluster. More...

Private Attributes
float	m_clusterMinLength
	minimum length of clusters for this algorithm More...

unsigned int	m_halfWindowLayers
	number of layers to use for half-window of sliding fit More...

float	m_samplingPitch
	sampling pitch for walking along sliding linear fit More...

float	m_maxCosSplittingAngle
	smallest scatter angle allowed when splitting cluster More...

float	m_minCosMergingAngle
	largest relative angle allowed when merging clusters More...

float	m_maxTransverseDisplacement
	maximum transverse displacement of associated clusters More...

float	m_maxLongitudinalDisplacement
	maximum longitudinal displacement of associated clusters More...

float	m_maxLongitudinalDisplacementSquared
	longitudinal displacement squared More...

Detailed Description

CosmicRaySplittingAlgorithm class.

Definition at line 21 of file CosmicRaySplittingAlgorithm.h.

Constructor & Destructor Documentation

lar_content::CosmicRaySplittingAlgorithm::CosmicRaySplittingAlgorithm ( )

Default constructor.

Definition at line 22 of file CosmicRaySplittingAlgorithm.cc.

                                                          :
     m_clusterMinLength(10.f),
     m_halfWindowLayers(30),
     m_samplingPitch(1.f),
     m_maxCosSplittingAngle(0.9925f),
     m_minCosMergingAngle(0.94f),
     m_maxTransverseDisplacement(5.f),
     m_maxLongitudinalDisplacement(30.f),
     m_maxLongitudinalDisplacementSquared(m_maxLongitudinalDisplacement * m_maxLongitudinalDisplacement)
 {
 }

Member Function Documentation

void lar_content::CosmicRaySplittingAlgorithm::BuildSlidingFitResultMap	(	const pandora::ClusterVector &	clusterVector,
		TwoDSlidingFitResultMap &	slidingFitResultMap
	)		const

private

Build the map of sliding fit results.

Parameters

clusterVector	the input cluster vector
slidingFitResultMap	the output sliding fit result map

Definition at line 181 of file CosmicRaySplittingAlgorithm.cc.

 {
     const float slidingFitPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));
 
     for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
     {
         if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
         {
             try
             {
                 const TwoDSlidingFitResult slidingFitResult(*iter, m_halfWindowLayers, slidingFitPitch);
 
                 if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
                     throw StatusCodeException(STATUS_CODE_FAILURE);
             }
             catch (StatusCodeException &statusCodeException)
             {
                 if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
                     throw statusCodeException;
             }
         }
     }
 }

StatusCode lar_content::CosmicRaySplittingAlgorithm::ConfirmSplitPosition	(	const TwoDSlidingFitResult &	branchSlidingFitResult,
		const TwoDSlidingFitResult &	replacementSlidingFitResult,
		const pandora::CartesianVector &	splitPosition,
		const pandora::CartesianVector &	splitDirection1,
		const pandora::CartesianVector &	splitDirection2,
		float &	lengthSquared1,
		float &	lengthSquared2
	)		const

private

Find a second replacement cluster that aligns with the scatter of the first branch cluster.

Parameters

branchSlidingFitResult	the sliding fit result for the branch cluster
replacementSlidingFitResult	the sliding fit result for the replacement cluster
splitPosition	the candidate split position on the branch cluster
splitDirection1	the first track direction just above the split position
splitDirection2	the second track direction just below the split position
lengthSquared1	figure of merit for association with first track direction
lengthSquared2	figure of merit for association with second track direction

Definition at line 270 of file CosmicRaySplittingAlgorithm.cc.

 {
     // Check if the replacement cluster points to the split position on the branch cluster
     bestLengthSquared1 = std::numeric_limits<float>::max();
     bestLengthSquared2 = std::numeric_limits<float>::max();
 
     bool foundMatch(false);
 
     const CartesianVector minPosition(replacementSlidingFitResult.GetGlobalMinLayerPosition());
     const CartesianVector maxPosition(replacementSlidingFitResult.GetGlobalMaxLayerPosition());
     const CartesianVector minDirection(replacementSlidingFitResult.GetGlobalMinLayerDirection());
     const CartesianVector maxDirection(replacementSlidingFitResult.GetGlobalMaxLayerDirection());
 
     const CartesianVector branchVertex1(branchSlidingFitResult.GetGlobalMinLayerPosition());
     const CartesianVector branchVertex2(branchSlidingFitResult.GetGlobalMaxLayerPosition());
     const CartesianVector branchDirection1(splitDirection1 * -1.f);
     const CartesianVector branchDirection2(splitDirection2 * -1.f);
 
     const float cosSplittingAngle(-splitDirection1.GetDotProduct(splitDirection2));
     const float branchLengthSquared((branchVertex2 - branchVertex1).GetMagnitudeSquared());
     const float replacementLengthSquared((maxPosition - minPosition).GetMagnitudeSquared());
 
     // Loop over each end of the replacement cluster
     for (unsigned int iFwd = 0; iFwd < 2; ++iFwd)
     {
         const CartesianVector pAxis((0 == iFwd) ? (maxPosition - minPosition) : (minPosition - maxPosition));
         const CartesianVector vtxPosition((0 == iFwd) ? minPosition : maxPosition);
         const CartesianVector endPosition((0 == iFwd) ? maxPosition : minPosition);
         const CartesianVector vtxDirection((0 == iFwd) ? (pAxis.GetDotProduct(minDirection) > 0 ? minDirection : minDirection * -1.f)
                                                        : (pAxis.GetDotProduct(maxDirection) > 0 ? maxDirection : maxDirection * -1.f));
 
         // Choose the correct end of the replacement cluster and require proximity to the branch cluster
         const float vtxDisplacement(LArClusterHelper::GetClosestDistance(vtxPosition, branchSlidingFitResult.GetCluster()));
         const float endDisplacement(LArClusterHelper::GetClosestDistance(endPosition, branchSlidingFitResult.GetCluster()));
 
         const float lengthSquared((vtxPosition - splitPosition).GetMagnitudeSquared());
         const float lengthSquared1((vtxPosition - branchVertex1).GetMagnitudeSquared());
         const float lengthSquared2((vtxPosition - branchVertex2).GetMagnitudeSquared());
 
         if (vtxDisplacement > endDisplacement)
             continue;
 
         if (std::min(lengthSquared, std::min(lengthSquared1, lengthSquared2)) > std::min(branchLengthSquared, replacementLengthSquared))
             continue;
 
         // Require good pointing information between replacement cluster and candidate split position
         float impactL(0.f), impactT(0.f);
         LArPointingClusterHelper::GetImpactParameters(vtxPosition, vtxDirection, splitPosition, impactL, impactT);
 
         if (impactT > m_maxTransverseDisplacement || impactL > m_maxLongitudinalDisplacement || impactL < -1.f ||
             impactT > std::max(1.5f, 0.577f * (1.f + impactL)))
             continue;
 
         // Check the segment of the branch cluster above the split position
         if (vtxDirection.GetDotProduct(branchDirection1) > std::max(m_minCosMergingAngle, cosSplittingAngle))
         {
             if (lengthSquared < bestLengthSquared1)
             {
                 bestLengthSquared1 = lengthSquared;
                 foundMatch = true;
             }
         }
 
         // Check the segment of the branch cluster below the split position
         if (vtxDirection.GetDotProduct(branchDirection2) > std::max(m_minCosMergingAngle, cosSplittingAngle))
         {
             if (lengthSquared < bestLengthSquared2)
             {
                 bestLengthSquared2 = lengthSquared;
                 foundMatch = true;
             }
         }
     }
 
     if (!foundMatch)
         return STATUS_CODE_NOT_FOUND;
 
     return STATUS_CODE_SUCCESS;
 }

StatusCode lar_content::CosmicRaySplittingAlgorithm::FindBestSplitPosition	(	const TwoDSlidingFitResult &	slidingFitResult,
		pandora::CartesianVector &	splitPosition,
		pandora::CartesianVector &	splitDirection1,
		pandora::CartesianVector &	splitDirection2
	)		const

private

Find the position of greatest scatter along a sliding linear fit.

Parameters

slidingFitResult	the input sliding linear fit result
splitPosition	the position of greatest scatter
splitDirection1	the direction vector just above the scatter position
splitDirection2	the direction vector just below the scatter position

Definition at line 207 of file CosmicRaySplittingAlgorithm.cc.

 {
     // Find position of greatest scatter for this cluster
     float splitCosTheta(m_maxCosSplittingAngle);
     bool foundSplit(false);
 
     const CartesianVector &minPosition(branchSlidingFitResult.GetGlobalMinLayerPosition());
     const CartesianVector &maxPosition(branchSlidingFitResult.GetGlobalMaxLayerPosition());
     const float halfWindowLength(branchSlidingFitResult.GetLayerFitHalfWindowLength());
 
     float minL(0.f), maxL(0.f), minT(0.f), maxT(0.f);
     branchSlidingFitResult.GetLocalPosition(minPosition, minL, minT);
     branchSlidingFitResult.GetLocalPosition(maxPosition, maxL, maxT);
 
     const unsigned int nSamplingPoints = static_cast<unsigned int>((maxL - minL) / m_samplingPitch);
 
     for (unsigned int n = 0; n < nSamplingPoints; ++n)
     {
         const float alpha((0.5f + static_cast<float>(n)) / static_cast<float>(nSamplingPoints));
         const float rL(minL + (maxL - minL) * alpha);
 
         CartesianVector centralPosition(0.f, 0.f, 0.f), forwardDirection(0.f, 0.f, 0.f), backwardDirection(0.f, 0.f, 0.f);
 
         if ((STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitPosition(rL, centralPosition)) ||
             (STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitDirection(rL + halfWindowLength, forwardDirection)) ||
             (STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitDirection(rL - halfWindowLength, backwardDirection)))
         {
             continue;
         }
 
         const float cosTheta(forwardDirection.GetDotProduct(backwardDirection));
 
         if (cosTheta < splitCosTheta)
         {
             CartesianVector forwardPosition(0.f, 0.f, 0.f);
             CartesianVector backwardPosition(0.f, 0.f, 0.f);
 
             if ((STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitPosition(rL + halfWindowLength, forwardPosition)) ||
                 (STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitPosition(rL - halfWindowLength, backwardPosition)))
             {
                 continue;
             }
 
             CartesianVector forwardVectorOutwards(forwardPosition - centralPosition);
             CartesianVector backwardVectorOutwards(backwardPosition - centralPosition);
 
             splitPosition = centralPosition;
             splitDirection1 = (forwardDirection.GetDotProduct(forwardVectorOutwards) > 0.f) ? forwardDirection : forwardDirection * -1.f;
             splitDirection2 = (backwardDirection.GetDotProduct(backwardVectorOutwards) > 0.f) ? backwardDirection : backwardDirection * -1.f;
             splitCosTheta = cosTheta;
             foundSplit = true;
         }
     }
 
     if (!foundSplit)
         return STATUS_CODE_NOT_FOUND;
 
     return STATUS_CODE_SUCCESS;
 }

void lar_content::CosmicRaySplittingAlgorithm::GetCaloHitListsToMove	(	const pandora::Cluster *const	pBranchCluster,
		const pandora::CartesianVector &	splitPosition,
		const pandora::CartesianVector &	splitDirection1,
		const pandora::CartesianVector &	splitDirection2,
		pandora::CaloHitList &	caloHitList1,
		pandora::CaloHitList &	caloHitList2
	)		const

private

Get lists of calo hits to move in order to split a branch cluster into two segments for case of two replacement clusters.

Parameters

pBranchCluster	the branch cluster
splitPosition	the split position
splitDirection1	the direction of the branch cluster just above the split position
splitDirection2	the direction of the branch cluster just below the split position
caloHitList1	the first segment to be split from the branch cluster
caloHitList2	the second segment to be split from the branch cluster

Definition at line 449 of file CosmicRaySplittingAlgorithm.cc.

 {
     CaloHitList caloHitListToSort;
     pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitListToSort);
 
     for (CaloHitList::const_iterator iter = caloHitListToSort.begin(), iterEnd = caloHitListToSort.end(); iter != iterEnd; ++iter)
     {
         const CaloHit *const pCaloHit = *iter;
 
         const float displacement1(splitDirection1.GetDotProduct(pCaloHit->GetPositionVector() - splitPosition));
         const float displacement2(splitDirection2.GetDotProduct(pCaloHit->GetPositionVector() - splitPosition));
 
         if (displacement1 > displacement2)
         {
             caloHitListToMove1.push_back(pCaloHit);
         }
         else
         {
             caloHitListToMove2.push_back(pCaloHit);
         }
     }
 }

StatusCode lar_content::CosmicRaySplittingAlgorithm::GetCaloHitListToKeep	(	const pandora::Cluster *const	pBranchCluster,
		const pandora::CaloHitList &	caloHitListToMove,
		pandora::CaloHitList &	caloHitListToKeep
	)		const

private

Split the branch cluster and add hits to the replacement cluster.

Parameters

pBranchCluster	the branch cluster
caloHitListToMove	the list of hits to be removed from the branch cluster and added to the replacement cluster
caloHitListToKeep	to receive the list of calo hits to keep

Definition at line 491 of file CosmicRaySplittingAlgorithm.cc.

 {
     if (caloHitListToMove.empty())
         return STATUS_CODE_FAILURE;
 
     CaloHitList caloHitList;
     pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
 
     for (CaloHitList::const_iterator iter = caloHitList.begin(), iterEnd = caloHitList.end(); iter != iterEnd; ++iter)
     {
         const CaloHit *const pCaloHit = *iter;
 
         if (caloHitListToMove.end() == std::find(caloHitListToMove.begin(), caloHitListToMove.end(), pCaloHit))
             caloHitListToKeep.push_back(pCaloHit);
     }
 
     return STATUS_CODE_SUCCESS;
 }

void lar_content::CosmicRaySplittingAlgorithm::GetCaloHitListToMove	(	const pandora::Cluster *const	pBranchCluster,
		const pandora::Cluster *const	pReplacementCluster,
		const pandora::CartesianVector &	splitPosition,
		const pandora::CartesianVector &	forwardDirection,
		const pandora::CartesianVector &	backwardDirection,
		pandora::CaloHitList &	caloHitList
	)		const

private

Get list of calo hits to move in order to split a branch cluster into two segments for case of one replacement cluster.

Parameters

pBranchCluster	the branch cluster
pReplacementCluster	the replacement cluster
splitPosition	the split position
forwardDirection	the direction of the branch cluster just above the split position
backwardDirection	the direction of the branch cluster just below the split position
caloHitList	the output hits to be removed from the branch cluster

Definition at line 397 of file CosmicRaySplittingAlgorithm.cc.

 {
     const CartesianVector forwardPosition(LArClusterHelper::GetClosestPosition(splitPosition, pBranchCluster));
     const CartesianVector vtxPosition(LArClusterHelper::GetClosestPosition(splitPosition, pReplacementCluster));
     const CartesianVector vtxDirection((forwardPosition - vtxPosition).GetUnitVector());
 
     CaloHitList caloHitListToSort, caloHitListToCheck;
     pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitListToSort);
 
     for (CaloHitList::const_iterator iter = caloHitListToSort.begin(), iterEnd = caloHitListToSort.end(); iter != iterEnd; ++iter)
     {
         const CaloHit *const pCaloHit = *iter;
 
         if (forwardDirection.GetDotProduct(pCaloHit->GetPositionVector() - forwardPosition) > 0.f)
         {
             caloHitListToMove.push_back(pCaloHit);
         }
         else if (forwardDirection.GetDotProduct(pCaloHit->GetPositionVector() - vtxPosition) > -1.25f)
         {
             caloHitListToCheck.push_back(pCaloHit);
         }
     }
 
     float closestLengthSquared(std::numeric_limits<float>::max());
 
     for (CaloHitList::const_iterator iter = caloHitListToCheck.begin(), iterEnd = caloHitListToCheck.end(); iter != iterEnd; ++iter)
     {
         const CaloHit *const pCaloHit = *iter;
 
         if (vtxDirection.GetCrossProduct(pCaloHit->GetPositionVector() - forwardPosition).GetMagnitude() >
             backwardDirection.GetCrossProduct(pCaloHit->GetPositionVector() - forwardPosition).GetMagnitude())
             continue;
 
         const float lengthSquared((pCaloHit->GetPositionVector() - vtxPosition).GetMagnitudeSquared());
 
         if (lengthSquared < closestLengthSquared)
             closestLengthSquared = lengthSquared;
     }
 
     for (CaloHitList::const_iterator iter = caloHitListToCheck.begin(), iterEnd = caloHitListToCheck.end(); iter != iterEnd; ++iter)
     {
         const CaloHit *const pCaloHit = *iter;
 
         if ((pCaloHit->GetPositionVector() - vtxPosition).GetMagnitudeSquared() >= closestLengthSquared)
             caloHitListToMove.push_back(pCaloHit);
     }
 }

void lar_content::CosmicRaySplittingAlgorithm::GetListOfCleanClusters	(	const pandora::ClusterList *const	pClusterList,
		pandora::ClusterVector &	clusterVector
	)		const

private

Populate cluster vector with subset of cluster list, containing clusters judged to be clean.

Parameters

pClusterList	address of the cluster list
clusterVector	to receive the populated cluster vector

Definition at line 164 of file CosmicRaySplittingAlgorithm.cc.

 {
     for (ClusterList::const_iterator iter = pClusterList->begin(), iterEnd = pClusterList->end(); iter != iterEnd; ++iter)
     {
         const Cluster *const pCluster = *iter;
 
         if (LArClusterHelper::GetLengthSquared(pCluster) < m_clusterMinLength * m_clusterMinLength)
             continue;
 
         clusterVector.push_back(pCluster);
     }
 
     std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
 }

bool lar_content::CosmicRaySplittingAlgorithm::IdentifyCrossedTracks	(	const pandora::Cluster *const	pBranchCluster,
		const pandora::Cluster *const	pReplacementCluster1,
		const pandora::Cluster *const	pReplacementCluster2,
		const pandora::CartesianVector &	splitPosition
	)		const

private

Identify crossed tracks formed from the branch cluster and its replacement cluster.

Parameters

pBranchCluster	the branch cluster
pReplacementCluster1	the first replacement cluster
pReplacementCluster2	the second replacement cluster
splitPosition	the split position

Definition at line 474 of file CosmicRaySplittingAlgorithm.cc.

 {
     CartesianVector branchVertex1(0.f, 0.f, 0.f), branchVertex2(0.f, 0.f, 0.f);
     LArClusterHelper::GetExtremalCoordinates(pBranchCluster, branchVertex1, branchVertex2);
 
     const CartesianVector replacementVertex1(LArClusterHelper::GetClosestPosition(splitPosition, pReplacementCluster1));
     const CartesianVector replacementVertex2(LArClusterHelper::GetClosestPosition(splitPosition, pReplacementCluster2));
 
     if ((replacementVertex2 - replacementVertex1).GetMagnitudeSquared() > (branchVertex2 - branchVertex1).GetMagnitudeSquared())
         return false;
 
     return true;
 }

StatusCode lar_content::CosmicRaySplittingAlgorithm::PerformDoubleSplit	(	const pandora::Cluster *const	pBranchCluster,
		const pandora::Cluster *const	pReplacementCluster1,
		const pandora::Cluster *const	pReplacementCluster2,
		const pandora::CartesianVector &	splitPosition,
		const pandora::CartesianVector &	splitDirection1,
		const pandora::CartesianVector &	splitDirection2
	)		const

private

Split a branch cluster for case of two replacement clusters.

Parameters

pBranchCluster	the branch cluster
pReplacementCluster1	the first replacement cluster
pReplacementCluster2	the second replacement cluster
splitPosition	the split position
splitDirection1	the direction of the branch cluster just above the split position
splitDirection2	the direction of the branch cluster just below the split position

Definition at line 371 of file CosmicRaySplittingAlgorithm.cc.

 {
     CaloHitList caloHitListToMove1, caloHitListToMove2;
     this->GetCaloHitListsToMove(pBranchCluster, splitPosition, splitDirection1, splitDirection2, caloHitListToMove1, caloHitListToMove2);
 
     CaloHitList caloHitListToKeep1;
     this->GetCaloHitListToKeep(pBranchCluster, caloHitListToMove1, caloHitListToKeep1);
 
     if (caloHitListToKeep1.empty())
         return STATUS_CODE_FAILURE;
 
     PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->SplitCluster(pBranchCluster, pReplacementCluster1, caloHitListToMove1));
 
     CaloHitList caloHitListToKeep2;
     this->GetCaloHitListToKeep(pBranchCluster, caloHitListToMove2, caloHitListToKeep2);
 
     if (caloHitListToKeep2.empty())
         return PandoraContentApi::MergeAndDeleteClusters(*this, pReplacementCluster2, pBranchCluster);
 
     return this->SplitCluster(pBranchCluster, pReplacementCluster2, caloHitListToMove2);
 }

StatusCode lar_content::CosmicRaySplittingAlgorithm::PerformSingleSplit	(	const pandora::Cluster *const	pBranchCluster,
		const pandora::Cluster *const	pReplacementCluster,
		const pandora::CartesianVector &	splitPosition,
		const pandora::CartesianVector &	forwardDirection,
		const pandora::CartesianVector &	backwardDirection
	)		const

private

Split a branch cluster for case of one replacement cluster.

Parameters

pBranchCluster	the branch cluster
pReplacementCluster	the replacement cluster
splitPosition	the split position
forwardDirection	the direction of the branch cluster just above the split position
backwardDirection	the direction of the branch cluster just below the split position

Definition at line 354 of file CosmicRaySplittingAlgorithm.cc.

 {
     CaloHitList caloHitListToMove;
     this->GetCaloHitListToMove(pBranchCluster, pReplacementCluster, splitPosition, forwardDirection, backwardDirection, caloHitListToMove);
 
     CaloHitList caloHitListToKeep;
     this->GetCaloHitListToKeep(pBranchCluster, caloHitListToMove, caloHitListToKeep);
 
     if (caloHitListToKeep.empty())
         return PandoraContentApi::MergeAndDeleteClusters(*this, pReplacementCluster, pBranchCluster);
 
     return this->SplitCluster(pBranchCluster, pReplacementCluster, caloHitListToMove);
 }

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

private

Definition at line 531 of file CosmicRaySplittingAlgorithm.cc.

 {
     PANDORA_RETURN_RESULT_IF_AND_IF(
         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ClusterMinLength", m_clusterMinLength));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(
         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SlidingFitHalfWindow", m_halfWindowLayers));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SamplingPitch", m_samplingPitch));
 
     if (m_samplingPitch < std::numeric_limits<float>::epsilon())
         return STATUS_CODE_INVALID_PARAMETER;
 
     PANDORA_RETURN_RESULT_IF_AND_IF(
         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxCosSplittingAngle", m_maxCosSplittingAngle));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(
         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCosMergingAngle", m_minCosMergingAngle));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
         XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
         XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
     m_maxLongitudinalDisplacementSquared = m_maxLongitudinalDisplacement * m_maxLongitudinalDisplacement;
 
     return STATUS_CODE_SUCCESS;
 }

StatusCode lar_content::CosmicRaySplittingAlgorithm::Run ( )

private

Definition at line 36 of file CosmicRaySplittingAlgorithm.cc.

 {
     const ClusterList *pClusterList = NULL;
     PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pClusterList));
 
     // Get ordered list of clean clusters
     ClusterVector clusterVector;
     this->GetListOfCleanClusters(pClusterList, clusterVector);
 
     // Calculate sliding fit results for clean clusters
     TwoDSlidingFitResultMap slidingFitResultMap;
     this->BuildSlidingFitResultMap(clusterVector, slidingFitResultMap);
 
     // Loop over clusters, identify and perform splits
     ClusterSet splitClusters;
 
     for (ClusterVector::const_iterator bIter = clusterVector.begin(), bIterEnd1 = clusterVector.end(); bIter != bIterEnd1; ++bIter)
     {
         if (splitClusters.count(*bIter) > 0)
             continue;
 
         TwoDSlidingFitResultMap::const_iterator bFitIter = slidingFitResultMap.find(*bIter);
 
         if (slidingFitResultMap.end() == bFitIter)
             continue;
 
         const TwoDSlidingFitResult &branchSlidingFitResult(bFitIter->second);
 
         // Find best split position for candidate branch cluster
         CartesianVector splitPosition(0.f, 0.f, 0.f);
         CartesianVector splitDirection1(0.f, 0.f, 0.f);
         CartesianVector splitDirection2(0.f, 0.f, 0.f);
 
         if (STATUS_CODE_SUCCESS != this->FindBestSplitPosition(branchSlidingFitResult, splitPosition, splitDirection1, splitDirection2))
             continue;
 
         // Find candidate replacement clusters to merge into branch cluster at the split position
         TwoDSlidingFitResultMap::const_iterator bestReplacementIter1(slidingFitResultMap.end());
         TwoDSlidingFitResultMap::const_iterator bestReplacementIter2(slidingFitResultMap.end());
 
         float bestLengthSquared1(m_maxLongitudinalDisplacementSquared);
         float bestLengthSquared2(m_maxLongitudinalDisplacementSquared);
 
         for (ClusterVector::const_iterator rIter = clusterVector.begin(), rIterEnd = clusterVector.end(); rIter != rIterEnd; ++rIter)
         {
             if (splitClusters.count(*rIter) > 0)
                 continue;
 
             TwoDSlidingFitResultMap::const_iterator rFitIter = slidingFitResultMap.find(*rIter);
 
             if (slidingFitResultMap.end() == rFitIter)
                 continue;
 
             const TwoDSlidingFitResult &replacementSlidingFitResult(rFitIter->second);
 
             if (branchSlidingFitResult.GetCluster() == replacementSlidingFitResult.GetCluster())
                 continue;
 
             float lengthSquared1(std::numeric_limits<float>::max());
             float lengthSquared2(std::numeric_limits<float>::max());
 
             if (STATUS_CODE_SUCCESS != this->ConfirmSplitPosition(branchSlidingFitResult, replacementSlidingFitResult, splitPosition,
                                            splitDirection1, splitDirection2, lengthSquared1, lengthSquared2))
                 continue;
 
             if (lengthSquared1 < bestLengthSquared1)
             {
                 bestLengthSquared1 = lengthSquared1;
                 bestReplacementIter1 = rFitIter;
             }
 
             if (lengthSquared2 < bestLengthSquared2)
             {
                 bestLengthSquared2 = lengthSquared2;
                 bestReplacementIter2 = rFitIter;
             }
         }
 
         const Cluster *pReplacementCluster1 = NULL;
         const Cluster *pReplacementCluster2 = NULL;
 
         if (slidingFitResultMap.end() != bestReplacementIter1)
             pReplacementCluster1 = bestReplacementIter1->first;
 
         if (slidingFitResultMap.end() != bestReplacementIter2)
             pReplacementCluster2 = bestReplacementIter2->first;
 
         if (NULL == pReplacementCluster1 && NULL == pReplacementCluster2)
             continue;
 
         const Cluster *const pBranchCluster = branchSlidingFitResult.GetCluster();
 
         // Crossed tracks
         if (pReplacementCluster1 && pReplacementCluster2)
         {
             if (!(this->IdentifyCrossedTracks(pBranchCluster, pReplacementCluster1, pReplacementCluster2, splitPosition)))
                 continue;
 
             PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
                 this->PerformDoubleSplit(pBranchCluster, pReplacementCluster1, pReplacementCluster2, splitPosition, splitDirection1, splitDirection2));
         }
         // Single scatter
         else if (pReplacementCluster1)
         {
             PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
                 this->PerformSingleSplit(pBranchCluster, pReplacementCluster1, splitPosition, splitDirection1, splitDirection2));
         }
         else if (pReplacementCluster2)
         {
             PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
                 this->PerformSingleSplit(pBranchCluster, pReplacementCluster2, splitPosition, splitDirection2, splitDirection1));
         }
 
         // Choose not to re-use clusters (for now)
         if (pReplacementCluster1)
             splitClusters.insert(pReplacementCluster1);
 
         if (pReplacementCluster2)
             splitClusters.insert(pReplacementCluster2);
 
         splitClusters.insert(pBranchCluster);
     }
 
     return STATUS_CODE_SUCCESS;
 }

StatusCode lar_content::CosmicRaySplittingAlgorithm::SplitCluster	(	const pandora::Cluster *const	pBranchCluster,
		const pandora::Cluster *const	pReplacementCluster,
		const pandora::CaloHitList &	caloHitListToMove
	)		const

private

Split the branch cluster and add hits to the replacement cluster.

Parameters

pBranchCluster	the branch cluster
pReplacementCluster	the replacement cluster
caloHitListToMove	the list of hits to be removed from the branch cluster and added to the replacement cluster

Definition at line 513 of file CosmicRaySplittingAlgorithm.cc.

 {
     if (caloHitListToMove.empty())
         return STATUS_CODE_FAILURE;
 
     for (CaloHitList::const_iterator iter = caloHitListToMove.begin(), iterEnd = caloHitListToMove.end(); iter != iterEnd; ++iter)
     {
         const CaloHit *const pCaloHit = *iter;
         PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::RemoveFromCluster(*this, pBranchCluster, pCaloHit));
         PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToCluster(*this, pReplacementCluster, pCaloHit));
     }
 
     return STATUS_CODE_SUCCESS;
 }