LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LongitudinalAssociationAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
13 
15 
16 using namespace pandora;
17 
18 namespace lar_content
19 {
20 
21 LongitudinalAssociationAlgorithm::LongitudinalAssociationAlgorithm() :
22  m_minClusterLayers(4),
23  m_maxGapLayers(7),
24  m_fitLayers(30),
25  m_maxGapDistanceSquared(10.f),
26  m_minCosRelativeAngle(0.985f),
27  m_maxTransverseDisplacement(2.f),
28  m_maxLongitudinalDisplacement(2.f),
29  m_hitSizeZ(0.3f),
30  m_hitSizeX(0.5f),
31  m_view(TPC_3D)
32 {
33 }
34 
35 //------------------------------------------------------------------------------------------------------------------------------------------
36 
37 void LongitudinalAssociationAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const
38 {
39  if (!pClusterList->empty())
40  m_view = LArClusterHelper::GetClusterHitType(pClusterList->front());
41 
42  for (ClusterList::const_iterator iter = pClusterList->begin(), iterEnd = pClusterList->end(); iter != iterEnd; ++iter)
43  {
44  const Cluster *const pCluster = *iter;
45 
46  if (1 + pCluster->GetOuterPseudoLayer() - pCluster->GetInnerPseudoLayer() < m_minClusterLayers)
47  continue;
48 
49  clusterVector.push_back(pCluster);
50  }
51 
52  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByInnerLayer);
53 }
54 
55 //------------------------------------------------------------------------------------------------------------------------------------------
56 
58 {
59  // ATTN This method assumes that clusters have been sorted by layer
60  for (ClusterVector::const_iterator iterI = clusterVector.begin(), iterIEnd = clusterVector.end(); iterI != iterIEnd; ++iterI)
61  {
62  const Cluster *const pInnerCluster = *iterI;
63 
64  for (ClusterVector::const_iterator iterJ = iterI, iterJEnd = clusterVector.end(); iterJ != iterJEnd; ++iterJ)
65  {
66  const Cluster *const pOuterCluster = *iterJ;
67 
68  if (pInnerCluster == pOuterCluster)
69  continue;
70 
71  if (!this->AreClustersAssociated(pInnerCluster, pOuterCluster))
72  continue;
73 
74  clusterAssociationMap[pInnerCluster].m_forwardAssociations.insert(pOuterCluster);
75  clusterAssociationMap[pOuterCluster].m_backwardAssociations.insert(pInnerCluster);
76  }
77  }
78 }
79 
80 //------------------------------------------------------------------------------------------------------------------------------------------
81 
82 bool LongitudinalAssociationAlgorithm::IsExtremalCluster(const bool isForward, const Cluster *const pCurrentCluster, const Cluster *const pTestCluster) const
83 {
84  const unsigned int currentLayer(isForward ? pCurrentCluster->GetOuterPseudoLayer() : pCurrentCluster->GetInnerPseudoLayer());
85  const unsigned int testLayer(isForward ? pTestCluster->GetOuterPseudoLayer() : pTestCluster->GetInnerPseudoLayer());
86 
87  if (isForward && ((testLayer > currentLayer) || ((testLayer == currentLayer) && LArClusterHelper::SortByNHits(pTestCluster, pCurrentCluster))))
88  return true;
89 
90  if (!isForward && ((testLayer < currentLayer) || ((testLayer == currentLayer) && LArClusterHelper::SortByNHits(pTestCluster, pCurrentCluster))))
91  return true;
92 
93  return false;
94 }
95 
96 //------------------------------------------------------------------------------------------------------------------------------------------
97 
98 bool LongitudinalAssociationAlgorithm::AreClustersAssociated(const Cluster *const pInnerCluster, const Cluster *const pOuterCluster) const
99 {
100  if (pOuterCluster->GetInnerPseudoLayer() < pInnerCluster->GetInnerPseudoLayer())
101  throw pandora::StatusCodeException(STATUS_CODE_NOT_ALLOWED);
102 
103  // TODO Remove hardcoded numbers
104  if ((pOuterCluster->GetInnerPseudoLayer() < pInnerCluster->GetInnerPseudoLayer() + 3) ||
105  (pInnerCluster->GetOuterPseudoLayer() + 3 > pOuterCluster->GetOuterPseudoLayer()))
106  {
107  return false;
108  }
109 
110  if ((pInnerCluster->GetOuterPseudoLayer() > pOuterCluster->GetInnerPseudoLayer() + 1) ||
111  (pOuterCluster->GetInnerPseudoLayer() > pInnerCluster->GetOuterPseudoLayer() + m_maxGapLayers))
112  {
113  return false;
114  }
115 
116  if ((2 * pInnerCluster->GetOuterPseudoLayer() < pOuterCluster->GetInnerPseudoLayer() + pInnerCluster->GetInnerPseudoLayer()) ||
117  (pInnerCluster->GetOuterPseudoLayer() + pOuterCluster->GetOuterPseudoLayer() < 2 * pOuterCluster->GetInnerPseudoLayer()))
118  {
119  return false;
120  }
121 
122  const CartesianVector innerEndCentroid(pInnerCluster->GetCentroid(pInnerCluster->GetOuterPseudoLayer()));
123  const CartesianVector outerStartCentroid(pOuterCluster->GetCentroid(pOuterCluster->GetInnerPseudoLayer()));
124 
125  const float ratio{LArGeometryHelper::GetWirePitchRatio(this->GetPandora(), m_view)};
126  const float maxGapDistanceSquaredAdjusted{ratio * ratio * m_maxGapDistanceSquared};
127 
128  if ((innerEndCentroid - outerStartCentroid).GetMagnitudeSquared() > maxGapDistanceSquaredAdjusted)
129  return false;
130 
131  ClusterFitResult innerEndFit, outerStartFit;
132  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, ClusterFitHelper::FitEnd(pInnerCluster, m_fitLayers, innerEndFit));
133  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, ClusterFitHelper::FitStart(pOuterCluster, m_fitLayers, outerStartFit));
134 
135  if (this->AreClustersAssociated(innerEndCentroid, outerStartCentroid, innerEndFit, outerStartFit))
136  return true;
137 
138  return false;
139 }
140 
141 //------------------------------------------------------------------------------------------------------------------------------------------
142 
143 bool LongitudinalAssociationAlgorithm::AreClustersAssociated(const CartesianVector &innerClusterEnd,
144  const CartesianVector &outerClusterStart, const ClusterFitResult &innerFit, const ClusterFitResult &outerFit) const
145 {
146  if (!innerFit.IsFitSuccessful() || !outerFit.IsFitSuccessful())
147  return false;
148 
149  if (innerFit.GetDirection().GetCosOpeningAngle(outerFit.GetDirection()) < m_minCosRelativeAngle)
150  return false;
151 
152  const float ratio{LArGeometryHelper::GetWirePitchRatio(this->GetPandora(), m_view)};
153  const float maxTransverseDisplacementAdjusted{ratio * m_maxTransverseDisplacement};
154  const float maxLongitudinalDisplacementAdjusted{ratio * m_maxLongitudinalDisplacement};
155 
156  const CartesianVector innerEndFit1(
157  innerFit.GetIntercept() + innerFit.GetDirection() * (innerFit.GetDirection().GetDotProduct(innerClusterEnd - innerFit.GetIntercept())));
158  const CartesianVector innerEndFit2(
159  outerFit.GetIntercept() + outerFit.GetDirection() * (outerFit.GetDirection().GetDotProduct(innerClusterEnd - outerFit.GetIntercept())));
160 
161  const CartesianVector outerStartFit1(outerFit.GetIntercept() +
162  outerFit.GetDirection() * (outerFit.GetDirection().GetDotProduct(outerClusterStart - outerFit.GetIntercept())));
163  const CartesianVector outerStartFit2(innerFit.GetIntercept() +
164  innerFit.GetDirection() * (innerFit.GetDirection().GetDotProduct(outerClusterStart - innerFit.GetIntercept())));
165 
166  const CartesianVector clusterSeparation(outerClusterStart - innerClusterEnd);
167 
168  if ((std::fabs(clusterSeparation.GetX()) < m_hitSizeX * maxTransverseDisplacementAdjusted) &&
169  (std::fabs(clusterSeparation.GetZ()) < m_hitSizeZ * maxLongitudinalDisplacementAdjusted))
170  return true;
171 
172  const CartesianVector fittedSeparation(outerStartFit1 - innerEndFit1);
173 
174  if ((std::fabs(fittedSeparation.GetX()) < m_hitSizeX * maxTransverseDisplacementAdjusted) &&
175  (std::fabs(fittedSeparation.GetZ()) < m_hitSizeZ * maxLongitudinalDisplacementAdjusted))
176  return true;
177 
178  const CartesianVector fittedInnerSeparation(innerEndFit2 - innerEndFit1);
179 
180  if ((std::fabs(fittedInnerSeparation.GetX()) < m_hitSizeX * maxTransverseDisplacementAdjusted) &&
181  (std::fabs(fittedInnerSeparation.GetZ()) < m_hitSizeZ * maxLongitudinalDisplacementAdjusted))
182  return true;
183 
184  const CartesianVector fittedOuterSeparation(outerStartFit2 - outerStartFit1);
185 
186  if ((std::fabs(fittedOuterSeparation.GetX()) < m_hitSizeX * maxTransverseDisplacementAdjusted) &&
187  (std::fabs(fittedOuterSeparation.GetZ()) < m_hitSizeZ * maxLongitudinalDisplacementAdjusted))
188  return true;
189 
190  return false;
191 }
192 
193 //------------------------------------------------------------------------------------------------------------------------------------------
194 
195 StatusCode LongitudinalAssociationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
196 {
197  PANDORA_RETURN_RESULT_IF_AND_IF(
198  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterLayers", m_minClusterLayers));
199 
200  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxGapLayers", m_maxGapLayers));
201 
202  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "FitLayers", m_fitLayers));
203 
204  float maxGapDistance = std::sqrt(m_maxGapDistanceSquared);
205  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxGapDistance", maxGapDistance));
206  m_maxGapDistanceSquared = maxGapDistance * maxGapDistance;
207 
208  PANDORA_RETURN_RESULT_IF_AND_IF(
209  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCosRelativeAngle", m_minCosRelativeAngle));
210 
211  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
212  XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
213 
214  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
215  XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
216 
217  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "HitSizeZ", m_hitSizeZ));
218 
219  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "HitSizeX", m_hitSizeX));
220 
222 }
223 
224 } // namespace lar_content
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.
Header file for the longitudinal association algorithm class.
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 PopulateClusterAssociationMap(const pandora::ClusterVector &clusterVector, ClusterAssociationMap &clusterAssociationMap) const
Populate the cluster association map.
static float GetWirePitchRatio(const pandora::Pandora &pandora, const pandora::HitType view)
Return the ratio of the wire pitch of the specified view to the minimum wire pitch for the detector...
float m_maxLongitudinalDisplacement
maximum allowed longitudinal displacement after extrapolation (normalised to cell size) ...
float m_minCosRelativeAngle
maximum allowed relative angle between associated clusters
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
std::unordered_map< const pandora::Cluster *, ClusterAssociation > ClusterAssociationMap
pandora::HitType m_view
The view to which the hits under consideration belong.
intermediate_table::const_iterator const_iterator
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
unsigned int m_minClusterLayers
minimum allowed number of layers for a clean cluster
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
Header file for the cluster helper class.
float m_hitSizeX
estimated hit size in x (drift time) dimension, units cm
float m_maxGapDistanceSquared
maximum allowed distance (squared) between associated clusters
bool IsExtremalCluster(const bool isForward, const pandora::Cluster *const pCurrentCluster, const pandora::Cluster *const pTestCluster) const
Determine which of two clusters is extremal.
unsigned int m_maxGapLayers
maximum allowed number of layers between associated clusters
float m_hitSizeZ
estimated hit size in z (wire number) dimension, units cm
bool AreClustersAssociated(const pandora::Cluster *const pInnerCluster, const pandora::Cluster *const pOuterCluster) const
Determine whether two clusters are associated.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
unsigned int m_fitLayers
number of layers to fit at start and end of cluster
float m_maxTransverseDisplacement
maximum allowed transverse displacement after extrapolation (normalised to cell size) ...
static bool SortByInnerLayer(const pandora::Cluster *const pLhs, const pandora::Cluster *const pRhs)
Sort clusters by inner layer, then position, then pulse-height.