LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
OvershootSplittingAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
13 
15 
17 
18 using namespace pandora;
19 
20 namespace lar_content
21 {
22 
23 OvershootSplittingAlgorithm::OvershootSplittingAlgorithm() :
25  m_minClusterLength(5.f),
26  m_maxClusterSeparation(5.f),
27  m_minVertexDisplacement(2.f),
28  m_maxIntersectDisplacement(1.5f),
29  m_minSplitDisplacement(10.f)
30 {
31 }
32 
33 //------------------------------------------------------------------------------------------------------------------------------------------
34 
35 void OvershootSplittingAlgorithm::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 
42  continue;
43 
44  clusterVector.push_back(pCluster);
45  }
46 
47  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
48 }
49 
50 //------------------------------------------------------------------------------------------------------------------------------------------
51 
52 void OvershootSplittingAlgorithm::FindBestSplitPositions(const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterPositionMap &clusterSplittingMap) const
53 {
54  // Use sliding fit results to build a list of intersection points
55  ClusterPositionMap clusterIntersectionMap;
56  this->BuildIntersectionMap(slidingFitResultMap, clusterIntersectionMap);
57 
58  // Sort intersection points according to their position along the sliding fit
59  ClusterPositionMap sortedIntersectionMap;
60  this->BuildSortedIntersectionMap(slidingFitResultMap, clusterIntersectionMap, sortedIntersectionMap);
61 
62  // Use intersection points to decide where/if to split cluster
63  this->PopulateSplitPositionMap(sortedIntersectionMap, clusterSplittingMap);
64 }
65 
66 //------------------------------------------------------------------------------------------------------------------------------------------
67 
68 void OvershootSplittingAlgorithm::BuildIntersectionMap(const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterPositionMap &clusterIntersectionMap) const
69 {
70  ClusterList clusterList;
71  for (const auto &mapEntry : slidingFitResultMap)
72  clusterList.push_back(mapEntry.first);
73  clusterList.sort(LArClusterHelper::SortByNHits);
74 
75  for (const Cluster *const pCluster1 : clusterList)
76  {
77  const TwoDSlidingFitResult &slidingFitResult1(slidingFitResultMap.at(pCluster1));
78 
79  for (const Cluster *const pCluster2 : clusterList)
80  {
81  if (pCluster1 == pCluster2)
82  continue;
83 
84  const TwoDSlidingFitResult &slidingFitResult2(slidingFitResultMap.at(pCluster2));
85 
86  try
87  {
88  const LArPointingCluster pointingCluster(slidingFitResult2);
89 
90  // Project pointing cluster onto target cluster
91  const CartesianVector innerPosition(pointingCluster.GetInnerVertex().GetPosition());
92  const CartesianVector outerPosition(pointingCluster.GetOuterVertex().GetPosition());
93  const float innerDisplacement(LArClusterHelper::GetClosestDistance(innerPosition, pCluster1));
94  const float outerDisplacement(LArClusterHelper::GetClosestDistance(outerPosition, pCluster1));
95  const bool useInner((innerDisplacement < outerDisplacement) ? true : false);
96 
97  const LArPointingCluster::Vertex &clusterVertex = (useInner ? pointingCluster.GetInnerVertex() : pointingCluster.GetOuterVertex());
98 
99  float rL2(0.f), rT2(0.f);
100  CartesianVector intersectPosition2(0.f, 0.f, 0.f);
101 
102  try
103  {
104  LArPointingClusterHelper::GetIntersection(clusterVertex, pCluster1, intersectPosition2, rL2, rT2);
105  }
106  catch (const StatusCodeException &)
107  {
108  continue;
109  }
110 
111  if (rL2 < -m_maxIntersectDisplacement || rL2 > m_maxClusterSeparation)
112  continue;
113 
114  // Find projected position and direction on target cluster
115  float rL1(0.f), rT1(0.f);
116  CartesianVector projectedPosition1(0.f, 0.f, 0.f), projectedDirection1(0.f, 0.f, 0.f);
117  slidingFitResult1.GetLocalPosition(intersectPosition2, rL1, rT1);
118 
119  const StatusCode statusCodePosition(slidingFitResult1.GetGlobalFitPosition(rL1, projectedPosition1));
120  if (STATUS_CODE_SUCCESS != statusCodePosition)
121  throw pandora::StatusCodeException(statusCodePosition);
122 
123  const StatusCode statusCodeDirection(slidingFitResult1.GetGlobalFitDirection(rL1, projectedDirection1));
124  if (STATUS_CODE_SUCCESS != statusCodeDirection)
125  throw pandora::StatusCodeException(statusCodeDirection);
126 
127  const CartesianVector projectedPosition2(clusterVertex.GetPosition());
128  const CartesianVector projectedDirection2(clusterVertex.GetDirection());
129 
130  // Find intersection of pointing cluster and target cluster
131  float firstDisplacement(0.f), secondDisplacement(0.f);
132  CartesianVector intersectPosition1(0.f, 0.f, 0.f);
133 
134  try
135  {
136  LArPointingClusterHelper::GetIntersection(projectedPosition1, projectedDirection1, projectedPosition2,
137  projectedDirection2, intersectPosition1, firstDisplacement, secondDisplacement);
138  }
139  catch (const StatusCodeException &)
140  {
141  continue;
142  }
143 
144  // Store intersections if they're sufficiently far along the cluster trajectory
145  const float closestDisplacement1(LArClusterHelper::GetClosestDistance(intersectPosition1, pCluster1));
146  const float closestDisplacement2(LArClusterHelper::GetClosestDistance(intersectPosition1, pCluster2));
147 
148  if (std::max(closestDisplacement1, closestDisplacement2) > m_maxClusterSeparation)
149  continue;
150 
151  const CartesianVector minPosition(slidingFitResult1.GetGlobalMinLayerPosition());
152  const CartesianVector maxPosition(slidingFitResult1.GetGlobalMaxLayerPosition());
153  const float lengthSquared((maxPosition - minPosition).GetMagnitudeSquared());
154 
155  const float minDisplacementSquared((minPosition - intersectPosition1).GetMagnitudeSquared());
156  const float maxDisplacementSquared((maxPosition - intersectPosition1).GetMagnitudeSquared());
157 
158  if (std::min(minDisplacementSquared, maxDisplacementSquared) < (m_minVertexDisplacement * m_minVertexDisplacement) ||
159  std::max(minDisplacementSquared, maxDisplacementSquared) > lengthSquared)
160  continue;
161 
162  clusterIntersectionMap[pCluster1].push_back(intersectPosition1);
163  }
164  catch (StatusCodeException &statusCodeException)
165  {
166  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
167  throw statusCodeException;
168  }
169  }
170  }
171 }
172 
173 //------------------------------------------------------------------------------------------------------------------------------------------
174 
176  const ClusterPositionMap &clusterIntersectionMap, ClusterPositionMap &sortedIntersectionMap) const
177 {
178  ClusterList clusterList;
179  for (const auto &mapEntry : clusterIntersectionMap)
180  clusterList.push_back(mapEntry.first);
181  clusterList.sort(LArClusterHelper::SortByNHits);
182 
183  for (const Cluster *const pCluster : clusterList)
184  {
185  const CartesianPointVector &inputPositionVector(clusterIntersectionMap.at(pCluster));
186 
187  if (inputPositionVector.empty())
188  continue;
189 
190  TwoDSlidingFitResultMap::const_iterator sIter = slidingFitResultMap.find(pCluster);
191  if (slidingFitResultMap.end() == sIter)
192  throw StatusCodeException(STATUS_CODE_FAILURE);
193 
194  const TwoDSlidingFitResult &slidingFitResult = sIter->second;
195 
196  MyTrajectoryPointList trajectoryPointList;
197  for (CartesianPointVector::const_iterator pIter = inputPositionVector.begin(), pIterEnd = inputPositionVector.end(); pIter != pIterEnd; ++pIter)
198  {
199  const CartesianVector &position = *pIter;
200  float rL(0.f), rT(0.f);
201  slidingFitResult.GetLocalPosition(position, rL, rT);
202  trajectoryPointList.push_back(MyTrajectoryPoint(rL, position));
203  }
204 
205  std::sort(trajectoryPointList.begin(), trajectoryPointList.end(), OvershootSplittingAlgorithm::SortByHitProjection);
206 
207  if (trajectoryPointList.empty())
208  throw StatusCodeException(STATUS_CODE_FAILURE);
209 
210  for (MyTrajectoryPointList::const_iterator qIter = trajectoryPointList.begin(), qIterEnd = trajectoryPointList.end(); qIter != qIterEnd; ++qIter)
211  {
212  const CartesianVector &clusterPosition = qIter->second;
213  sortedIntersectionMap[pCluster].push_back(clusterPosition);
214  }
215  }
216 }
217 
218 //------------------------------------------------------------------------------------------------------------------------------------------
219 
220 void OvershootSplittingAlgorithm::PopulateSplitPositionMap(const ClusterPositionMap &clusterIntersectionMap, ClusterPositionMap &clusterSplittingMap) const
221 {
222  ClusterList clusterList;
223  for (const auto &mapEntry : clusterIntersectionMap)
224  clusterList.push_back(mapEntry.first);
225  clusterList.sort(LArClusterHelper::SortByNHits);
226 
227  for (const Cluster *const pCluster : clusterList)
228  {
229  const CartesianPointVector &inputPositionVector(clusterIntersectionMap.at(pCluster));
230 
231  if (inputPositionVector.empty())
232  continue;
233 
234  // Select pairs of positions within a given separation, and calculate their average position
235  MyTrajectoryPointList candidatePositionList;
236 
237  bool foundPrevPosition(false);
238  CartesianVector prevPosition(0.f, 0.f, 0.f);
239 
240  for (CartesianPointVector::const_iterator pIter = inputPositionVector.begin(), pIterEnd = inputPositionVector.end(); pIter != pIterEnd; ++pIter)
241  {
242  const CartesianVector &nextPosition = *pIter;
243 
244  if (foundPrevPosition)
245  {
246  const CartesianVector averagePosition((nextPosition + prevPosition) * 0.5f);
247  const float displacementSquared((nextPosition - prevPosition).GetMagnitudeSquared());
248 
249  if (displacementSquared < m_maxIntersectDisplacement * m_maxIntersectDisplacement)
250  candidatePositionList.push_back(MyTrajectoryPoint(displacementSquared, averagePosition));
251  }
252 
253  prevPosition = nextPosition;
254  foundPrevPosition = true;
255  }
256 
257  if (candidatePositionList.empty())
258  continue;
259 
260  std::sort(candidatePositionList.begin(), candidatePositionList.end(), OvershootSplittingAlgorithm::SortByHitProjection);
261 
262  // Use the average positions of the closest pairs of points as the split position
263  bool foundPrevCandidate(false);
264  CartesianVector prevCandidate(0.f, 0.f, 0.f);
265 
266  for (MyTrajectoryPointList::const_iterator pIter = candidatePositionList.begin(), pIterEnd = candidatePositionList.end();
267  pIter != pIterEnd; ++pIter)
268  {
269  const CartesianVector &nextCandidate = pIter->second;
270 
271  if (foundPrevCandidate)
272  {
273  if ((nextCandidate - prevCandidate).GetMagnitudeSquared() < m_minSplitDisplacement * m_minSplitDisplacement)
274  continue;
275  }
276 
277  clusterSplittingMap[pCluster].push_back(nextCandidate);
278  prevCandidate = nextCandidate;
279  foundPrevCandidate = true;
280  }
281  }
282 }
283 
284 //------------------------------------------------------------------------------------------------------------------------------------------
285 
287 {
288  if (lhs.first != rhs.first)
289  return (lhs.first < rhs.first);
290 
291  return (lhs.second.GetMagnitudeSquared() > rhs.second.GetMagnitudeSquared());
292 }
293 
294 //------------------------------------------------------------------------------------------------------------------------------------------
295 
296 StatusCode OvershootSplittingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
297 {
298  PANDORA_RETURN_RESULT_IF_AND_IF(
299  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterLength", m_minClusterLength));
300 
301  PANDORA_RETURN_RESULT_IF_AND_IF(
302  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxClusterSeparation", m_maxClusterSeparation));
303 
304  PANDORA_RETURN_RESULT_IF_AND_IF(
305  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinVertexDisplacement", m_minVertexDisplacement));
306 
307  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
308  XmlHelper::ReadValue(xmlHandle, "MaxIntersectDisplacement", m_maxIntersectDisplacement));
309 
310  PANDORA_RETURN_RESULT_IF_AND_IF(
311  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinSplitDisplacement", m_minSplitDisplacement));
312 
314 }
315 
316 } // 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 lar pointing cluster 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...
static void GetIntersection(const LArPointingCluster::Vertex &firstVertex, const LArPointingCluster::Vertex &secondVertex, pandora::CartesianVector &intersectPosition, float &firstDisplacement, float &secondDisplacement)
Get intersection of two vertices.
Header file for the overshoot splitting algorithm class.
intermediate_table::const_iterator const_iterator
LArPointingCluster class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
TFile f
Definition: plotHisto.C:6
std::unordered_map< const pandora::Cluster *, pandora::CartesianPointVector > ClusterPositionMap
Header file for the cluster helper class.
const Vertex & GetOuterVertex() const
Get the outer vertex.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
const Vertex & GetInnerVertex() const
Get the inner vertex.
void BuildIntersectionMap(const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterPositionMap &clusterIntersectionMap) const
Use sliding fit results to calculate intersections of clusters.
const pandora::CartesianVector & GetDirection() const
Get the vertex direction.
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
void FindBestSplitPositions(const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterPositionMap &clusterSplittingMap) const
Determine best split positions based on sliding fit result.
std::vector< MyTrajectoryPoint > MyTrajectoryPointList
std::pair< float, pandora::CartesianVector > MyTrajectoryPoint
static bool SortByHitProjection(const MyTrajectoryPoint &lhs, const MyTrajectoryPoint &rhs)
Sort pfos by number of constituent hits.
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
void PopulateSplitPositionMap(const ClusterPositionMap &sortedIntersectionMap, ClusterPositionMap &clusterSplittingMap) const
Select split positions from sorted list of candidate positions.
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
void BuildSortedIntersectionMap(const TwoDSlidingFitResultMap &slidingFitResultMap, const ClusterPositionMap &clusterIntersectionMap, ClusterPositionMap &sortedIntersectionMap) const
Use intersection points to decide on splitting points.
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.