LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
LongitudinalAssociationAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
14 
15 using namespace pandora;
16 
17 namespace lar_content
18 {
19 
20 LongitudinalAssociationAlgorithm::LongitudinalAssociationAlgorithm() :
21  m_minClusterLayers(4),
22  m_maxGapLayers(7),
23  m_fitLayers(30),
24  m_maxGapDistanceSquared(10.f),
25  m_minCosRelativeAngle(0.985f),
26  m_maxTransverseDisplacement(2.f),
27  m_maxLongitudinalDisplacement(2.f),
28  m_hitSizeZ(0.3f),
29  m_hitSizeX(0.5f)
30 {
31 }
32 
33 //------------------------------------------------------------------------------------------------------------------------------------------
34 
35 void LongitudinalAssociationAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const
36 {
37  for (ClusterList::const_iterator iter = pClusterList->begin(), iterEnd = pClusterList->end(); iter != iterEnd; ++iter)
38  {
39  const Cluster *const pCluster = *iter;
40 
41  if (1 + pCluster->GetOuterPseudoLayer() - pCluster->GetInnerPseudoLayer() < m_minClusterLayers)
42  continue;
43 
44  clusterVector.push_back(pCluster);
45  }
46 
47  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByInnerLayer);
48 }
49 
50 //------------------------------------------------------------------------------------------------------------------------------------------
51 
53 {
54  // ATTN This method assumes that clusters have been sorted by layer
55  for (ClusterVector::const_iterator iterI = clusterVector.begin(), iterIEnd = clusterVector.end(); iterI != iterIEnd; ++iterI)
56  {
57  const Cluster *const pInnerCluster = *iterI;
58 
59  for (ClusterVector::const_iterator iterJ = iterI, iterJEnd = clusterVector.end(); iterJ != iterJEnd; ++iterJ)
60  {
61  const Cluster *const pOuterCluster = *iterJ;
62 
63  if (pInnerCluster == pOuterCluster)
64  continue;
65 
66  if (!this->AreClustersAssociated(pInnerCluster, pOuterCluster))
67  continue;
68 
69  clusterAssociationMap[pInnerCluster].m_forwardAssociations.insert(pOuterCluster);
70  clusterAssociationMap[pOuterCluster].m_backwardAssociations.insert(pInnerCluster);
71  }
72  }
73 }
74 
75 //------------------------------------------------------------------------------------------------------------------------------------------
76 
77 bool LongitudinalAssociationAlgorithm::IsExtremalCluster(const bool isForward, const Cluster *const pCurrentCluster, const Cluster *const pTestCluster) const
78 {
79  const unsigned int currentLayer(isForward ? pCurrentCluster->GetOuterPseudoLayer() : pCurrentCluster->GetInnerPseudoLayer());
80  const unsigned int testLayer(isForward ? pTestCluster->GetOuterPseudoLayer() : pTestCluster->GetInnerPseudoLayer());
81 
82  if (isForward && ((testLayer > currentLayer) || ((testLayer == currentLayer) && LArClusterHelper::SortByNHits(pTestCluster, pCurrentCluster))))
83  return true;
84 
85  if (!isForward && ((testLayer < currentLayer) || ((testLayer == currentLayer) && LArClusterHelper::SortByNHits(pTestCluster, pCurrentCluster))))
86  return true;
87 
88  return false;
89 }
90 
91 //------------------------------------------------------------------------------------------------------------------------------------------
92 
93 bool LongitudinalAssociationAlgorithm::AreClustersAssociated(const Cluster *const pInnerCluster, const Cluster *const pOuterCluster) const
94 {
95  if (pOuterCluster->GetInnerPseudoLayer() < pInnerCluster->GetInnerPseudoLayer())
96  throw pandora::StatusCodeException(STATUS_CODE_NOT_ALLOWED);
97 
98  // TODO Remove hardcoded numbers
99  if ((pOuterCluster->GetInnerPseudoLayer() < pInnerCluster->GetInnerPseudoLayer() + 3) ||
100  (pInnerCluster->GetOuterPseudoLayer() + 3 > pOuterCluster->GetOuterPseudoLayer()))
101  {
102  return false;
103  }
104 
105  if ((pInnerCluster->GetOuterPseudoLayer() > pOuterCluster->GetInnerPseudoLayer() + 1) ||
106  (pOuterCluster->GetInnerPseudoLayer() > pInnerCluster->GetOuterPseudoLayer() + m_maxGapLayers))
107  {
108  return false;
109  }
110 
111  if ((2 * pInnerCluster->GetOuterPseudoLayer() < pOuterCluster->GetInnerPseudoLayer() + pInnerCluster->GetInnerPseudoLayer()) ||
112  (pInnerCluster->GetOuterPseudoLayer() + pOuterCluster->GetOuterPseudoLayer() < 2 * pOuterCluster->GetInnerPseudoLayer()))
113  {
114  return false;
115  }
116 
117  const CartesianVector innerEndCentroid(pInnerCluster->GetCentroid(pInnerCluster->GetOuterPseudoLayer()));
118  const CartesianVector outerStartCentroid(pOuterCluster->GetCentroid(pOuterCluster->GetInnerPseudoLayer()));
119 
120  if ((innerEndCentroid - outerStartCentroid).GetMagnitudeSquared() > m_maxGapDistanceSquared)
121  return false;
122 
123  ClusterFitResult innerEndFit, outerStartFit;
124  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, ClusterFitHelper::FitEnd(pInnerCluster, m_fitLayers, innerEndFit));
125  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, ClusterFitHelper::FitStart(pOuterCluster, m_fitLayers, outerStartFit));
126 
127  if (this->AreClustersAssociated(innerEndCentroid, outerStartCentroid, innerEndFit, outerStartFit))
128  return true;
129 
130  return false;
131 }
132 
133 //------------------------------------------------------------------------------------------------------------------------------------------
134 
135 bool LongitudinalAssociationAlgorithm::AreClustersAssociated(const CartesianVector &innerClusterEnd, const CartesianVector &outerClusterStart,
136  const ClusterFitResult &innerFit, const ClusterFitResult &outerFit) const
137 {
138  if (!innerFit.IsFitSuccessful() || !outerFit.IsFitSuccessful())
139  return false;
140 
141  if (innerFit.GetDirection().GetCosOpeningAngle(outerFit.GetDirection()) < m_minCosRelativeAngle)
142  return false;
143 
144  const CartesianVector innerEndFit1(innerFit.GetIntercept() + innerFit.GetDirection() * (innerFit.GetDirection().GetDotProduct(innerClusterEnd - innerFit.GetIntercept())));
145  const CartesianVector innerEndFit2(outerFit.GetIntercept() + outerFit.GetDirection() * (outerFit.GetDirection().GetDotProduct(innerClusterEnd - outerFit.GetIntercept())));
146 
147  const CartesianVector outerStartFit1(outerFit.GetIntercept() + outerFit.GetDirection() * (outerFit.GetDirection().GetDotProduct(outerClusterStart - outerFit.GetIntercept())));
148  const CartesianVector outerStartFit2(innerFit.GetIntercept() + innerFit.GetDirection() * (innerFit.GetDirection().GetDotProduct(outerClusterStart - innerFit.GetIntercept())));
149 
150  const CartesianVector clusterSeparation(outerClusterStart - innerClusterEnd);
151 
152  if ((std::fabs(clusterSeparation.GetX()) < m_hitSizeX * m_maxTransverseDisplacement) &&
153  (std::fabs(clusterSeparation.GetZ()) < m_hitSizeZ * m_maxLongitudinalDisplacement))
154  return true;
155 
156  const CartesianVector fittedSeparation(outerStartFit1 - innerEndFit1);
157 
158  if ((std::fabs(fittedSeparation.GetX()) < m_hitSizeX * m_maxTransverseDisplacement) &&
159  (std::fabs(fittedSeparation.GetZ()) < m_hitSizeZ * m_maxLongitudinalDisplacement))
160  return true;
161 
162  const CartesianVector fittedInnerSeparation(innerEndFit2 - innerEndFit1);
163 
164  if ((std::fabs(fittedInnerSeparation.GetX()) < m_hitSizeX * m_maxTransverseDisplacement) &&
165  (std::fabs(fittedInnerSeparation.GetZ()) < m_hitSizeZ * m_maxLongitudinalDisplacement))
166  return true;
167 
168  const CartesianVector fittedOuterSeparation(outerStartFit2 - outerStartFit1);
169 
170  if ((std::fabs(fittedOuterSeparation.GetX()) < m_hitSizeX * m_maxTransverseDisplacement) &&
171  (std::fabs(fittedOuterSeparation.GetZ()) < m_hitSizeZ * m_maxLongitudinalDisplacement))
172  return true;
173 
174  return false;
175 }
176 
177 //------------------------------------------------------------------------------------------------------------------------------------------
178 
179 StatusCode LongitudinalAssociationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
180 {
181  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
182  "MinClusterLayers", m_minClusterLayers));
183 
184  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
185  "MaxGapLayers", m_maxGapLayers));
186 
187  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
188  "FitLayers", m_fitLayers));
189 
190  float maxGapDistance = std::sqrt(m_maxGapDistanceSquared);
191  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
192  "MaxGapDistance", maxGapDistance));
193  m_maxGapDistanceSquared = maxGapDistance * maxGapDistance;
194 
195  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
196  "MinCosRelativeAngle", m_minCosRelativeAngle));
197 
198  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
199  "MaxTransverseDisplacement", m_maxTransverseDisplacement));
200 
201  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
202  "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
203 
204  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
205  "HitSizeZ", m_hitSizeZ));
206 
207  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
208  "HitSizeX", m_hitSizeX));
209 
211 }
212 
213 } // 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.
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)
unsigned int m_minClusterLayers
minimum allowed number of layers for a clean cluster
TFile f
Definition: plotHisto.C:6
intermediate_table::const_iterator const_iterator
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
std::vector< art::Ptr< recob::Cluster > > ClusterVector
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::unordered_map< const pandora::Cluster *, ClusterAssociation > ClusterAssociationMap
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.