LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LongitudinalExtensionAlgorithm.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 LongitudinalExtensionAlgorithm::LongitudinalExtensionAlgorithm() :
22  m_clusterMinLength(5.f),
23  m_clusterMinLayerOccupancy(0.75f),
24  m_nodeMaxDisplacement(1.5f),
25  m_nodeMaxCosRelativeAngle(0.906f),
26  m_emissionMaxLongitudinalDisplacement(15.f),
27  m_emissionMaxTransverseDisplacement(2.5f),
28  m_emissionMaxCosRelativeAngle(0.985f)
29 {
30 }
31 
32 //------------------------------------------------------------------------------------------------------------------------------------------
33 
34 void LongitudinalExtensionAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const
35 {
36  for (ClusterList::const_iterator iter = pClusterList->begin(), iterEnd = pClusterList->end(); iter != iterEnd; ++iter)
37  {
38  const Cluster *const pCluster = *iter;
39 
41  continue;
42 
44  continue;
45 
46  clusterVector.push_back(pCluster);
47  }
48 
49  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
50 }
51 
52 //------------------------------------------------------------------------------------------------------------------------------------------
53 
54 void LongitudinalExtensionAlgorithm::FillClusterAssociationMatrix(const ClusterVector &clusterVector, ClusterAssociationMatrix &clusterAssociationMatrix) const
55 {
56  // Convert each input cluster into a pointing cluster
57  LArPointingClusterList pointingClusterList;
58 
59  for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
60  {
61  try
62  {
63  pointingClusterList.push_back(LArPointingCluster(*iter));
64  }
65  catch (StatusCodeException &)
66  {
67  }
68  }
69 
70  // Form associations between pairs of pointing clusters
71  for (LArPointingClusterList::const_iterator iterI = pointingClusterList.begin(), iterEndI = pointingClusterList.end(); iterI != iterEndI; ++iterI)
72  {
73  const LArPointingCluster &clusterI = *iterI;
74 
75  for (LArPointingClusterList::const_iterator iterJ = iterI, iterEndJ = pointingClusterList.end(); iterJ != iterEndJ; ++iterJ)
76  {
77  const LArPointingCluster &clusterJ = *iterJ;
78 
79  if (clusterI.GetCluster() == clusterJ.GetCluster())
80  continue;
81 
82  this->FillClusterAssociationMatrix(clusterI, clusterJ, clusterAssociationMatrix);
83  }
84  }
85 }
86 
87 //------------------------------------------------------------------------------------------------------------------------------------------
88 
90  const LArPointingCluster &clusterI, const LArPointingCluster &clusterJ, ClusterAssociationMatrix &clusterAssociationMatrix) const
91 {
92  const Cluster *const pClusterI(clusterI.GetCluster());
93  const Cluster *const pClusterJ(clusterJ.GetCluster());
94 
95  if (pClusterI == pClusterJ)
96  return;
97 
98  // Check that new layer occupancy would be reasonable
100  return;
101 
102  // Identify closest pair of vertices
103  LArPointingCluster::Vertex targetVertexI, targetVertexJ;
104 
105  try
106  {
107  LArPointingClusterHelper::GetClosestVertices(clusterI, clusterJ, targetVertexI, targetVertexJ);
108  }
109  catch (StatusCodeException &)
110  {
111  return;
112  }
113 
114  // (Just in case...)
115  if (!(targetVertexI.IsInitialized() && targetVertexJ.IsInitialized()))
116  throw StatusCodeException(STATUS_CODE_FAILURE);
117 
118  const CartesianVector &vertexPositionI(targetVertexI.GetPosition());
119  const CartesianVector &vertexPositionJ(targetVertexJ.GetPosition());
120  const CartesianVector &vertexDirectionI(targetVertexI.GetDirection());
121  const CartesianVector &vertexDirectionJ(targetVertexJ.GetDirection());
122 
123  // Check for reasonable proximity between vertices
124  const float distanceSquared((vertexPositionI - vertexPositionJ).GetMagnitudeSquared());
125 
127  return;
128 
129  // Check that vertices have a reasonable linear fit
130  if (targetVertexI.GetRms() > 1.f || targetVertexJ.GetRms() > 1.f)
131  return;
132 
133  // Association type
135 
136  // Requirements for Nodes
137  if (distanceSquared < 2.f * m_nodeMaxDisplacement * m_nodeMaxDisplacement)
138  {
139  associationType = ClusterAssociation::WEAK;
140 
141  if (distanceSquared < m_nodeMaxDisplacement * m_nodeMaxDisplacement)
142  {
143  const float cosTheta(-vertexDirectionI.GetDotProduct(vertexDirectionJ));
144 
145  if (cosTheta > m_nodeMaxCosRelativeAngle)
146  {
147  associationType = ClusterAssociation::STRONG;
148  }
149  }
150  }
151 
152  // Requirements for Emissions
153  const float clusterLengthI(LArPointingClusterHelper::GetLength(clusterI));
154  const float clusterLengthJ(LArPointingClusterHelper::GetLength(clusterJ));
155 
156  if (associationType < ClusterAssociation::STRONG)
157  {
158  const float cosTheta(-vertexDirectionI.GetDotProduct(vertexDirectionJ));
159  const float cosThetaI((vertexPositionI - vertexPositionJ).GetUnitVector().GetDotProduct(vertexDirectionI));
160  const float cosThetaJ((vertexPositionJ - vertexPositionI).GetUnitVector().GetDotProduct(vertexDirectionJ));
161 
162  float rT1(0.f), rL1(0.f), rT2(0.f), rL2(0.f);
163  LArPointingClusterHelper::GetImpactParameters(vertexPositionI, vertexDirectionI, vertexPositionJ, rL1, rT1);
164  LArPointingClusterHelper::GetImpactParameters(vertexPositionJ, vertexDirectionJ, vertexPositionI, rL2, rT2);
165 
166  if ((rL1 > -2.5f && rL1 < std::min(0.66f * clusterLengthJ, m_emissionMaxLongitudinalDisplacement)) &&
167  (rL2 > -2.5f && rL2 < std::min(0.66f * clusterLengthI, m_emissionMaxLongitudinalDisplacement)) && (rT1 < m_emissionMaxTransverseDisplacement) &&
168  (rT2 < m_emissionMaxTransverseDisplacement) && (targetVertexI.GetRms() < 0.5f && targetVertexJ.GetRms() < 0.5f) &&
169  (cosTheta > m_emissionMaxCosRelativeAngle) && (std::fabs(cosThetaI) > 0.25f) && (std::fabs(cosThetaJ) > 0.25f))
170  {
171  associationType = ClusterAssociation::STRONG;
172  }
173  }
174 
175  // Record the association
176  if (ClusterAssociation::NONE != associationType)
177  {
180  (void)clusterAssociationMatrix[pClusterI].insert(
181  ClusterAssociationMap::value_type(pClusterJ, ClusterAssociation(vertexTypeI, vertexTypeJ, associationType, clusterLengthJ)));
182  (void)clusterAssociationMatrix[pClusterJ].insert(
183  ClusterAssociationMap::value_type(pClusterI, ClusterAssociation(vertexTypeJ, vertexTypeI, associationType, clusterLengthI)));
184  return;
185  }
186 }
187 
188 //------------------------------------------------------------------------------------------------------------------------------------------
189 
190 void LongitudinalExtensionAlgorithm::FillClusterMergeMap(const ClusterAssociationMatrix &inputAssociationMatrix, ClusterMergeMap &clusterMergeMap) const
191 {
192  // Decide which associations will become merges
193  // To make the merge A <-> B, both A -> B and B -> A must be strong associations
194  // with the largest figures of merit of all the A -> X and B -> Y associations
195 
196  // First step: remove double-counting from the map of associations
197  // i.e. if the map has A <-> B, B <-> C, A <-> C, then remove A <-> C
198  ClusterAssociationMatrix clusterAssociationMatrix;
199 
200  ClusterVector sortedInputClusters;
201  for (const auto &mapEntry : inputAssociationMatrix)
202  sortedInputClusters.push_back(mapEntry.first);
203  std::sort(sortedInputClusters.begin(), sortedInputClusters.end(), LArClusterHelper::SortByNHits);
204 
205  for (const Cluster *const pCluster1 : sortedInputClusters)
206  {
207  const ClusterAssociationMap &associationMap1(inputAssociationMatrix.at(pCluster1));
208 
209  for (const Cluster *const pCluster2 : sortedInputClusters)
210  {
211  if (pCluster1 == pCluster2)
212  continue;
213 
214  const ClusterAssociationMap &associationMap2(inputAssociationMatrix.at(pCluster2));
215 
216  ClusterAssociationMap::const_iterator iter12 = associationMap1.find(pCluster2);
217  if (associationMap1.end() == iter12)
218  continue;
219 
220  ClusterAssociationMap::const_iterator iter21 = associationMap2.find(pCluster1);
221  if (associationMap2.end() == iter21)
222  continue;
223 
224  const ClusterAssociation &association12(iter12->second);
225  const ClusterAssociation &association21(iter21->second);
226 
227  bool isAssociated(true);
228 
229  ClusterVector sortedAssociationClusters;
230  for (const auto &mapEntry : associationMap1)
231  sortedAssociationClusters.push_back(mapEntry.first);
232  std::sort(sortedAssociationClusters.begin(), sortedAssociationClusters.end(), LArClusterHelper::SortByNHits);
233 
234  for (const Cluster *const pCluster3 : sortedAssociationClusters)
235  {
236  const ClusterAssociation &association13(associationMap1.at(pCluster3));
237 
238  ClusterAssociationMap::const_iterator iter23 = associationMap2.find(pCluster3);
239  if (associationMap2.end() == iter23)
240  continue;
241 
242  const ClusterAssociation &association23(iter23->second);
243 
244  if (association12.GetParent() == association13.GetParent() && association23.GetParent() == association21.GetParent() &&
245  association13.GetDaughter() != association23.GetDaughter())
246  {
247  isAssociated = false;
248  break;
249  }
250  }
251 
252  if (isAssociated)
253  {
254  (void)clusterAssociationMatrix[pCluster1].insert(ClusterAssociationMap::value_type(pCluster2, association12));
255  (void)clusterAssociationMatrix[pCluster2].insert(ClusterAssociationMap::value_type(pCluster1, association21));
256  }
257  }
258  }
259 
260  // Second step: find the best associations A -> X and B -> Y
261  ClusterAssociationMatrix intermediateAssociationMatrix;
262 
263  ClusterVector sortedClusters;
264  for (const auto &mapEntry : clusterAssociationMatrix)
265  sortedClusters.push_back(mapEntry.first);
266  std::sort(sortedClusters.begin(), sortedClusters.end(), LArClusterHelper::SortByNHits);
267 
268  for (const Cluster *const pParentCluster : sortedClusters)
269  {
270  const ClusterAssociationMap &clusterAssociationMap(clusterAssociationMatrix.at(pParentCluster));
271 
272  const Cluster *pBestClusterInner(nullptr);
274 
275  const Cluster *pBestClusterOuter(nullptr);
277 
278  ClusterVector sortedAssociationClusters;
279  for (const auto &mapEntry : clusterAssociationMap)
280  sortedAssociationClusters.push_back(mapEntry.first);
281  std::sort(sortedAssociationClusters.begin(), sortedAssociationClusters.end(), LArClusterHelper::SortByNHits);
282 
283  for (const Cluster *const pDaughterCluster : sortedAssociationClusters)
284  {
285  const ClusterAssociation &clusterAssociation(clusterAssociationMap.at(pDaughterCluster));
286 
287  // Inner associations
288  if (clusterAssociation.GetParent() == ClusterAssociation::INNER)
289  {
290  if (clusterAssociation.GetFigureOfMerit() > bestAssociationInner.GetFigureOfMerit())
291  {
292  bestAssociationInner = clusterAssociation;
293 
294  if (clusterAssociation.GetAssociation() == ClusterAssociation::STRONG)
295  pBestClusterInner = pDaughterCluster;
296  else
297  pBestClusterInner = nullptr;
298  }
299  }
300 
301  // Outer associations
302  if (clusterAssociation.GetParent() == ClusterAssociation::OUTER)
303  {
304  if (clusterAssociation.GetFigureOfMerit() > bestAssociationOuter.GetFigureOfMerit())
305  {
306  bestAssociationOuter = clusterAssociation;
307 
308  if (clusterAssociation.GetAssociation() == ClusterAssociation::STRONG)
309  pBestClusterOuter = pDaughterCluster;
310  else
311  pBestClusterOuter = nullptr;
312  }
313  }
314  }
315 
316  if (pBestClusterInner)
317  (void)intermediateAssociationMatrix[pParentCluster].insert(ClusterAssociationMap::value_type(pBestClusterInner, bestAssociationInner));
318 
319  if (pBestClusterOuter)
320  (void)intermediateAssociationMatrix[pParentCluster].insert(ClusterAssociationMap::value_type(pBestClusterOuter, bestAssociationOuter));
321  }
322 
323  // Third step: make the merge if A -> X and B -> Y is in fact A -> B and B -> A
324  ClusterVector intermediateSortedClusters;
325  for (const auto &mapEntry : intermediateAssociationMatrix)
326  intermediateSortedClusters.push_back(mapEntry.first);
327  std::sort(intermediateSortedClusters.begin(), intermediateSortedClusters.end(), LArClusterHelper::SortByNHits);
328 
329  for (const Cluster *const pParentCluster : intermediateSortedClusters)
330  {
331  const ClusterAssociationMap &parentAssociationMap(intermediateAssociationMatrix.at(pParentCluster));
332 
333  ClusterVector sortedAssociationClusters;
334  for (const auto &mapEntry : parentAssociationMap)
335  sortedAssociationClusters.push_back(mapEntry.first);
336  std::sort(sortedAssociationClusters.begin(), sortedAssociationClusters.end(), LArClusterHelper::SortByNHits);
337 
338  for (const Cluster *const pDaughterCluster : sortedAssociationClusters)
339  {
340  const ClusterAssociation &parentToDaughterAssociation(parentAssociationMap.at(pDaughterCluster));
341 
342  ClusterAssociationMatrix::const_iterator iter5 = intermediateAssociationMatrix.find(pDaughterCluster);
343 
344  if (intermediateAssociationMatrix.end() == iter5)
345  continue;
346 
347  const ClusterAssociationMap &daughterAssociationMap(iter5->second);
348 
349  ClusterAssociationMap::const_iterator iter6 = daughterAssociationMap.find(pParentCluster);
350 
351  if (daughterAssociationMap.end() == iter6)
352  continue;
353 
354  const ClusterAssociation &daughterToParentAssociation(iter6->second);
355 
356  if (parentToDaughterAssociation.GetParent() == daughterToParentAssociation.GetDaughter() &&
357  parentToDaughterAssociation.GetDaughter() == daughterToParentAssociation.GetParent())
358  {
359  ClusterList &parentList(clusterMergeMap[pParentCluster]);
360 
361  if (parentList.end() == std::find(parentList.begin(), parentList.end(), pDaughterCluster))
362  parentList.push_back(pDaughterCluster);
363  }
364  }
365  }
366 }
367 
368 //------------------------------------------------------------------------------------------------------------------------------------------
369 
370 StatusCode LongitudinalExtensionAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
371 {
372  PANDORA_RETURN_RESULT_IF_AND_IF(
373  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ClusterMinLength", m_clusterMinLength));
374 
375  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
376  XmlHelper::ReadValue(xmlHandle, "ClusterMinLayerOccupancy", m_clusterMinLayerOccupancy));
377 
378  PANDORA_RETURN_RESULT_IF_AND_IF(
379  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "NodeMaxDisplacement", m_nodeMaxDisplacement));
380 
381  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
382  XmlHelper::ReadValue(xmlHandle, "NodeMaxCosRelativeAngle", m_nodeMaxCosRelativeAngle));
383 
384  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
385  XmlHelper::ReadValue(xmlHandle, "EmissionMaxLongitudinalDisplacement", m_emissionMaxLongitudinalDisplacement));
386 
387  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
388  XmlHelper::ReadValue(xmlHandle, "EmissionMaxTransverseDisplacement", m_emissionMaxTransverseDisplacement));
389 
390  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
391  XmlHelper::ReadValue(xmlHandle, "EmissionMaxCosRelativeAngle", m_emissionMaxCosRelativeAngle));
392 
393  return ClusterExtensionAlgorithm::ReadSettings(xmlHandle);
394 }
395 
396 } // 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.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
void FillClusterAssociationMatrix(const pandora::ClusterVector &clusterVector, ClusterAssociationMatrix &clusterAssociationMatrix) const
Fill the cluster association matrix.
std::vector< LArPointingCluster > LArPointingClusterList
static void GetImpactParameters(const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
Calculate impact parameters between a pair of pointing vertices.
static float GetLength(const LArPointingCluster &pointingCluster)
Calculate distance squared between inner and outer vertices of pointing cluster.
intermediate_table::const_iterator const_iterator
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.
TFile f
Definition: plotHisto.C:6
Header file for the cluster helper class.
bool IsInitialized() const
Whether the vertex has been initialized.
std::unordered_map< const pandora::Cluster *, ClusterAssociation > ClusterAssociationMap
std::unordered_map< const pandora::Cluster *, ClusterAssociationMap > ClusterAssociationMatrix
static float GetLayerOccupancy(const pandora::Cluster *const pCluster)
Fraction of occupied layers in cluster.
const pandora::CartesianVector & GetDirection() const
Get the vertex direction.
float GetRms() const
Get rms from vertex fit.
void FillClusterMergeMap(const ClusterAssociationMatrix &clusterAssociationMatrix, ClusterMergeMap &clusterMergeMap) const
Fill the cluster merge map.
Header file for the longitudinal extension algorithm class.
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
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...
std::unordered_map< const pandora::Cluster *, pandora::ClusterList > ClusterMergeMap
bool IsInnerVertex() const
Is this the inner vertex.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)