9 #include "Pandora/AlgorithmHeaders.h" 21 CosmicRayExtensionAlgorithm::CosmicRayExtensionAlgorithm() :
22 m_minClusterLength(3.
f),
23 m_minSeedClusterLength(6.
f),
24 m_maxLongitudinalDisplacement(30.
f),
25 m_maxTransverseDisplacement(2.
f),
26 m_minCosRelativeAngle(0.966
f),
37 const Cluster *
const pCluster = *iter;
42 clusterVector.push_back(pCluster);
61 catch (StatusCodeException &)
88 const Cluster *
const pClusterI(clusterI.
GetCluster());
89 const Cluster *
const pClusterJ(clusterJ.
GetCluster());
91 if (pClusterI == pClusterJ)
101 catch (StatusCodeException &)
108 throw StatusCodeException(STATUS_CODE_FAILURE);
110 const CartesianVector vertexI(clusterVertexI.
GetPosition());
111 const CartesianVector vertexJ(clusterVertexJ.
GetPosition());
124 const float distanceSquaredIJ((vertexI - vertexJ).GetMagnitudeSquared());
130 const CartesianVector directionI((endI - vertexI).GetUnitVector());
131 const CartesianVector directionJ((endJ - vertexJ).GetUnitVector());
133 const float cosTheta(-directionI.GetDotProduct(directionJ));
136 if (cosTheta < cosThetaCut)
140 const CartesianVector directionIJ((endJ - endI).GetUnitVector());
141 const CartesianVector directionJI((endI - endJ).GetUnitVector());
143 const float overlapL(directionIJ.GetDotProduct(vertexJ - vertexI));
144 const float overlapT(directionIJ.GetCrossProduct(vertexJ - vertexI).GetMagnitude());
146 if (overlapL < -1.f || overlapL * overlapL > 2.
f * std::min(lengthSquaredI, lengthSquaredJ) ||
151 const float rms1(this->
CalculateRms(pClusterI, endI, directionIJ));
152 const float rms2(this->
CalculateRms(pClusterJ, endJ, directionJI));
154 const float rms(0.5
f * (rms1 + rms2));
155 const float rmsCut(2.
f *
m_maxAverageRms * (cosTheta - cosThetaCut) / (1.0 - cosThetaCut));
165 (void)clusterAssociationMatrix[pClusterI].insert(
166 ClusterAssociationMap::value_type(pClusterJ,
ClusterAssociation(vertexTypeI, vertexTypeJ, associationType, lengthSquaredJ)));
168 (void)clusterAssociationMatrix[pClusterJ].insert(
169 ClusterAssociationMap::value_type(pClusterI,
ClusterAssociation(vertexTypeJ, vertexTypeI, associationType, lengthSquaredI)));
185 for (
const auto &mapEntry : inputAssociationMatrix)
186 sortedInputClusters.push_back(mapEntry.first);
189 for (
const Cluster *
const pCluster1 : sortedInputClusters)
193 for (
const Cluster *
const pCluster2 : sortedInputClusters)
195 if (pCluster1 == pCluster2)
201 if (associationMap1.end() == iter12)
205 if (associationMap2.end() == iter21)
211 bool isAssociated(
true);
214 for (
const auto &mapEntry : associationMap1)
215 sortedAssociationClusters.push_back(mapEntry.first);
218 for (
const Cluster *
const pCluster3 : sortedAssociationClusters)
223 if (associationMap2.end() == iter23)
229 association13.GetDaughter() != association23.
GetDaughter())
231 isAssociated =
false;
238 (void)clusterAssociationMatrix[pCluster1].insert(ClusterAssociationMap::value_type(pCluster2, association12));
239 (void)clusterAssociationMatrix[pCluster2].insert(ClusterAssociationMap::value_type(pCluster1, association21));
248 for (
const auto &mapEntry : clusterAssociationMatrix)
249 sortedClusters.push_back(mapEntry.first);
252 for (
const Cluster *
const pParentCluster : sortedClusters)
256 const Cluster *pBestClusterInner(
nullptr);
259 const Cluster *pBestClusterOuter(
nullptr);
263 for (
const auto &mapEntry : clusterAssociationMap)
264 sortedAssociationClusters.push_back(mapEntry.first);
267 for (
const Cluster *
const pDaughterCluster : sortedAssociationClusters)
269 const ClusterAssociation &clusterAssociation(clusterAssociationMap.at(pDaughterCluster));
274 if (clusterAssociation.GetFigureOfMerit() > bestAssociationInner.
GetFigureOfMerit())
276 bestAssociationInner = clusterAssociation;
277 pBestClusterInner = pDaughterCluster;
284 if (clusterAssociation.GetFigureOfMerit() > bestAssociationOuter.
GetFigureOfMerit())
286 bestAssociationOuter = clusterAssociation;
287 pBestClusterOuter = pDaughterCluster;
292 if (pBestClusterInner)
293 (void)intermediateAssociationMatrix[pParentCluster].insert(ClusterAssociationMap::value_type(pBestClusterInner, bestAssociationInner));
295 if (pBestClusterOuter)
296 (void)intermediateAssociationMatrix[pParentCluster].insert(ClusterAssociationMap::value_type(pBestClusterOuter, bestAssociationOuter));
301 for (
const auto &mapEntry : intermediateAssociationMatrix)
302 intermediateSortedClusters.push_back(mapEntry.first);
305 for (
const Cluster *
const pParentCluster : intermediateSortedClusters)
310 for (
const auto &mapEntry : parentAssociationMap)
311 sortedAssociationClusters.push_back(mapEntry.first);
314 for (
const Cluster *
const pDaughterCluster : sortedAssociationClusters)
316 const ClusterAssociation &parentToDaughterAssociation(parentAssociationMap.at(pDaughterCluster));
320 if (intermediateAssociationMatrix.end() == iter5)
327 if (daughterAssociationMap.end() == iter6)
332 if (parentToDaughterAssociation.GetParent() == daughterToParentAssociation.
GetDaughter() &&
333 parentToDaughterAssociation.GetDaughter() == daughterToParentAssociation.
GetParent())
335 ClusterList &parentList(clusterMergeMap[pParentCluster]);
337 if (parentList.end() == std::find(parentList.begin(), parentList.end(), pDaughterCluster))
338 parentList.push_back(pDaughterCluster);
348 float totalChi2(0.
f);
349 float totalHits(0.
f);
351 CaloHitList caloHitList;
352 pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
356 const CaloHit *
const pCaloHit = *iter;
358 const CartesianVector hitPosition(pCaloHit->GetPositionVector());
359 const CartesianVector predictedPosition(position + direction * direction.GetDotProduct(hitPosition - position));
361 totalChi2 += (predictedPosition - hitPosition).GetMagnitudeSquared();
366 return std::sqrt(totalChi2 / totalHits);
375 PANDORA_RETURN_RESULT_IF_AND_IF(
376 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinClusterLength",
m_minClusterLength));
378 PANDORA_RETURN_RESULT_IF_AND_IF(
379 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinSeedClusterLength",
m_minSeedClusterLength));
381 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
384 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
387 PANDORA_RETURN_RESULT_IF_AND_IF(
388 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinCosRelativeAngle",
m_minCosRelativeAngle));
390 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxAverageRms",
m_maxAverageRms));
static bool SortByNHits(const pandora::Cluster *const pLhs, const pandora::Cluster *const pRhs)
Sort clusters by number of hits, then layer span, then inner layer, then position, then pulse-height.
std::vector< LArPointingCluster > LArPointingClusterList
Header file for the cosmic-ray extension algorithm class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
ClusterAssociation class.
float m_minSeedClusterLength
LArPointingCluster class.
const pandora::Cluster * GetCluster() const
Get the address of the cluster.
static void GetClosestVertices(const bool useX, const bool useY, const bool useZ, const LArPointingCluster &pointingClusterI, const LArPointingCluster &pointingClusterJ, LArPointingCluster::Vertex &closestVertexI, LArPointingCluster::Vertex &closestVertexJ)
Given a pair of pointing clusters, receive the closest or farthest pair of vertices.
Header file for the cluster helper class.
bool IsInitialized() const
Whether the vertex has been initialized.
const Vertex & GetOuterVertex() const
Get the outer vertex.
const Vertex & GetInnerVertex() const
Get the inner vertex.
AssociationType
Association enumeration.
std::unordered_map< const pandora::Cluster *, ClusterAssociation > ClusterAssociationMap
std::unordered_map< const pandora::Cluster *, ClusterAssociationMap > ClusterAssociationMatrix
float CalculateRms(const pandora::Cluster *const pCluster, const pandora::CartesianVector &position, const pandora::CartesianVector &direction) const
Calculate RMS deviation of a cluster with respect to the reference line.
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...
void FillClusterAssociationMatrix(const pandora::ClusterVector &clusterVector, ClusterAssociationMatrix &clusterAssociationMatrix) const
Fill the cluster association matrix.
VertexType GetParent() const
Get parent.
float m_minCosRelativeAngle
float m_maxTransverseDisplacement
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
VertexType GetDaughter() const
Get daughter.
std::unordered_map< const pandora::Cluster *, pandora::ClusterList > ClusterMergeMap
bool IsInnerVertex() const
Is this the inner vertex.
float m_maxLongitudinalDisplacement
std::vector< art::Ptr< recob::Cluster > > ClusterVector
static float GetLengthSquared(const LArPointingCluster &pointingCluster)
Calculate distance squared between inner and outer vertices of pointing cluster.
void FillClusterMergeMap(const ClusterAssociationMatrix &clusterAssociationMatrix, ClusterMergeMap &clusterMergeMap) const
Fill the cluster merge map.
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
VertexType
Vertex enumeration.
float GetFigureOfMerit() const
Get figure of merit.