LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CrossGapsExtensionAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
14 
16 
17 using namespace pandora;
18 
19 namespace lar_content
20 {
21 
22 CrossGapsExtensionAlgorithm::CrossGapsExtensionAlgorithm() :
23  m_minClusterLength(5.f),
24  m_minGapFraction(0.5f),
25  m_maxGapTolerance(2.f),
26  m_maxTransverseDisplacement(2.5f),
27  m_maxRelativeAngle(10.f)
28 {
29 }
30 
31 //------------------------------------------------------------------------------------------------------------------------------------------
32 
33 void CrossGapsExtensionAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const
34 {
35  // ATTN May want to opt-out completely if no gap information available
36  // if (PandoraContentApi::GetGeometry(*this)->GetDetectorGapList().empty())
37  // return;
38 
39  for (const Cluster *const pCluster : *pClusterList)
40  {
42  continue;
43 
44  clusterVector.push_back(pCluster);
45  }
46 
47  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
48 }
49 
50 //------------------------------------------------------------------------------------------------------------------------------------------
51 
53 {
54  // Build lists of pointing clusters in proximity to gaps
55  LArPointingClusterList innerPointingClusterList, outerPointingClusterList;
56  this->BuildPointingClusterList(clusterVector, innerPointingClusterList, outerPointingClusterList);
57 
58  // Form associations between pairs of pointing clusters
59  for (const LArPointingCluster &pointingClusterInner : innerPointingClusterList)
60  {
61  const LArPointingCluster::Vertex &pointingVertexInner(pointingClusterInner.GetInnerVertex());
62  const float zInner(pointingVertexInner.GetPosition().GetZ());
63 
64  for (const LArPointingCluster &pointingClusterOuter : outerPointingClusterList)
65  {
66  const LArPointingCluster::Vertex &pointingVertexOuter(pointingClusterOuter.GetOuterVertex());
67  const float zOuter(pointingVertexOuter.GetPosition().GetZ());
68 
69  if (!this->IsAcrossGap(zOuter, zInner, LArClusterHelper::GetClusterHitType(pointingClusterInner.GetCluster())))
70  continue;
71 
72  if (!this->IsAssociated(pointingVertexInner, pointingVertexOuter))
73  continue;
74 
75  const Cluster *const pClusterInner(pointingClusterInner.GetCluster());
76  const Cluster *const pClusterOuter(pointingClusterOuter.GetCluster());
77 
78  const float lengthSquaredInner(LArClusterHelper::GetLengthSquared(pClusterInner));
79  const float lengthSquaredOuter(LArClusterHelper::GetLengthSquared(pClusterOuter));
80 
81  (void)clusterAssociationMatrix[pClusterInner].insert(ClusterAssociationMap::value_type(pClusterOuter,
83  (void)clusterAssociationMatrix[pClusterOuter].insert(ClusterAssociationMap::value_type(pClusterInner,
85  }
86  }
87 }
88 
89 //------------------------------------------------------------------------------------------------------------------------------------------
90 
92  LArPointingClusterList &innerPointingClusterList, LArPointingClusterList &outerPointingClusterList) const
93 {
94  // Convert each input cluster into a pointing cluster
95  LArPointingClusterList pointingClusterList;
96 
97  for (const Cluster *const pCluster : clusterVector)
98  {
99  try
100  {
101  pointingClusterList.push_back(LArPointingCluster(pCluster));
102  }
103  catch (StatusCodeException &)
104  {
105  }
106  }
107 
108  // Identify clusters adjacent to detector gaps
109  this->BuildPointingClusterList(true, pointingClusterList, innerPointingClusterList);
110  this->BuildPointingClusterList(false, pointingClusterList, outerPointingClusterList);
111 }
112 
113 //------------------------------------------------------------------------------------------------------------------------------------------
114 
116  const bool useInner, const LArPointingClusterList &inputPointingClusterList, LArPointingClusterList &outputPointingClusterList) const
117 {
118  for (const LArPointingCluster &pointingCluster : inputPointingClusterList)
119  {
120  const LArPointingCluster::Vertex &pointingVertex(useInner ? pointingCluster.GetInnerVertex() : pointingCluster.GetOuterVertex());
121  const HitType hitType(LArClusterHelper::GetClusterHitType(pointingCluster.GetCluster()));
122 
123  if (LArGeometryHelper::IsInGap(this->GetPandora(), pointingVertex.GetPosition(), hitType, m_maxGapTolerance))
124  outputPointingClusterList.push_back(pointingCluster);
125  }
126 }
127 
128 //------------------------------------------------------------------------------------------------------------------------------------------
129 
131 {
132  const float ratio{LArGeometryHelper::GetWirePitchRatio(this->GetPandora(), LArClusterHelper::GetClusterHitType(pointingVertex1.GetCluster()))};
133  const float maxTransverseDisplacementAdjusted{ratio * m_maxTransverseDisplacement};
134  const float maxLongitudinalDisplacement((pointingVertex2.GetPosition() - pointingVertex1.GetPosition()).GetMagnitude());
135 
136  const bool isAssociated1(LArPointingClusterHelper::IsEmission(pointingVertex1.GetPosition(), pointingVertex2, -1.f,
137  maxLongitudinalDisplacement + 1.f, maxTransverseDisplacementAdjusted, m_maxRelativeAngle));
138  const bool isAssociated2(LArPointingClusterHelper::IsEmission(pointingVertex2.GetPosition(), pointingVertex1, -1.f,
139  maxLongitudinalDisplacement + 1.f, maxTransverseDisplacementAdjusted, m_maxRelativeAngle));
140 
141  return (isAssociated1 && isAssociated2);
142 }
143 
144 //------------------------------------------------------------------------------------------------------------------------------------------
145 
146 bool CrossGapsExtensionAlgorithm::IsAcrossGap(const float minZ, const float maxZ, const HitType hitType) const
147 {
148  if (maxZ - minZ < std::numeric_limits<float>::epsilon())
149  return false;
150 
151  const float gapDeltaZ(LArGeometryHelper::CalculateGapDeltaZ(this->GetPandora(), minZ, maxZ, hitType));
152 
153  if (gapDeltaZ / (maxZ - minZ) < m_minGapFraction)
154  return false;
155 
156  return true;
157 }
158 
159 //------------------------------------------------------------------------------------------------------------------------------------------
160 
161 void CrossGapsExtensionAlgorithm::FillClusterMergeMap(const ClusterAssociationMatrix &inputAssociationMatrix, ClusterMergeMap &clusterMergeMap) const
162 {
163  // Decide which associations will become merges
164  // To make the merge A <-> B, both A -> B and B -> A must be strong associations
165  // with the largest figures of merit of all the A -> X and B -> Y associations
166 
167  // First step: remove double-counting from the map of associations
168  // i.e. if the map has A <-> B, B <-> C, A <-> C, then remove A <-> C
169  ClusterAssociationMatrix clusterAssociationMatrix;
170 
171  ClusterVector sortedInputClusters;
172  for (const auto &mapEntry : inputAssociationMatrix)
173  sortedInputClusters.push_back(mapEntry.first);
174  std::sort(sortedInputClusters.begin(), sortedInputClusters.end(), LArClusterHelper::SortByNHits);
175 
176  for (const Cluster *const pCluster1 : sortedInputClusters)
177  {
178  const ClusterAssociationMap &associationMap1(inputAssociationMatrix.at(pCluster1));
179 
180  for (const Cluster *const pCluster2 : sortedInputClusters)
181  {
182  if (pCluster1 == pCluster2)
183  continue;
184 
185  const ClusterAssociationMap &associationMap2(inputAssociationMatrix.at(pCluster2));
186 
187  ClusterAssociationMap::const_iterator iter12 = associationMap1.find(pCluster2);
188  if (associationMap1.end() == iter12)
189  continue;
190 
191  ClusterAssociationMap::const_iterator iter21 = associationMap2.find(pCluster1);
192  if (associationMap2.end() == iter21)
193  continue;
194 
195  const ClusterAssociation &association12(iter12->second);
196  const ClusterAssociation &association21(iter21->second);
197 
198  bool isAssociated(true);
199 
200  ClusterVector sortedAssociationClusters;
201  for (const auto &mapEntry : associationMap1)
202  sortedAssociationClusters.push_back(mapEntry.first);
203  std::sort(sortedAssociationClusters.begin(), sortedAssociationClusters.end(), LArClusterHelper::SortByNHits);
204 
205  for (const Cluster *const pCluster3 : sortedAssociationClusters)
206  {
207  const ClusterAssociation &association13(associationMap1.at(pCluster3));
208 
209  ClusterAssociationMap::const_iterator iter23 = associationMap2.find(pCluster3);
210  if (associationMap2.end() == iter23)
211  continue;
212 
213  const ClusterAssociation &association23(iter23->second);
214 
215  if (association12.GetParent() == association13.GetParent() && association23.GetParent() == association21.GetParent() &&
216  association13.GetDaughter() != association23.GetDaughter())
217  {
218  isAssociated = false;
219  break;
220  }
221  }
222 
223  if (isAssociated)
224  {
225  (void)clusterAssociationMatrix[pCluster1].insert(ClusterAssociationMap::value_type(pCluster2, association12));
226  (void)clusterAssociationMatrix[pCluster2].insert(ClusterAssociationMap::value_type(pCluster1, association21));
227  }
228  }
229  }
230 
231  // Second step: find the best associations A -> X and B -> Y
232  ClusterAssociationMatrix intermediateAssociationMatrix;
233 
234  ClusterVector sortedClusters;
235  for (const auto &mapEntry : clusterAssociationMatrix)
236  sortedClusters.push_back(mapEntry.first);
237  std::sort(sortedClusters.begin(), sortedClusters.end(), LArClusterHelper::SortByNHits);
238 
239  for (const Cluster *const pParentCluster : sortedClusters)
240  {
241  const ClusterAssociationMap &clusterAssociationMap(clusterAssociationMatrix.at(pParentCluster));
242 
243  const Cluster *pBestClusterInner(nullptr);
245 
246  const Cluster *pBestClusterOuter(nullptr);
248 
249  ClusterVector sortedAssociationClusters;
250  for (const auto &mapEntry : clusterAssociationMap)
251  sortedAssociationClusters.push_back(mapEntry.first);
252  std::sort(sortedAssociationClusters.begin(), sortedAssociationClusters.end(), LArClusterHelper::SortByNHits);
253 
254  for (const Cluster *const pDaughterCluster : sortedAssociationClusters)
255  {
256  const ClusterAssociation &clusterAssociation(clusterAssociationMap.at(pDaughterCluster));
257 
258  // Inner associations
259  if (clusterAssociation.GetParent() == ClusterAssociation::INNER)
260  {
261  if (clusterAssociation.GetFigureOfMerit() > bestAssociationInner.GetFigureOfMerit())
262  {
263  bestAssociationInner = clusterAssociation;
264  pBestClusterInner = pDaughterCluster;
265  }
266  }
267 
268  // Outer associations
269  if (clusterAssociation.GetParent() == ClusterAssociation::OUTER)
270  {
271  if (clusterAssociation.GetFigureOfMerit() > bestAssociationOuter.GetFigureOfMerit())
272  {
273  bestAssociationOuter = clusterAssociation;
274  pBestClusterOuter = pDaughterCluster;
275  }
276  }
277  }
278 
279  if (pBestClusterInner)
280  (void)intermediateAssociationMatrix[pParentCluster].insert(ClusterAssociationMap::value_type(pBestClusterInner, bestAssociationInner));
281 
282  if (pBestClusterOuter)
283  (void)intermediateAssociationMatrix[pParentCluster].insert(ClusterAssociationMap::value_type(pBestClusterOuter, bestAssociationOuter));
284  }
285 
286  // Third step: make the merge if A -> X and B -> Y is in fact A -> B and B -> A
287  ClusterVector intermediateSortedClusters;
288  for (const auto &mapEntry : intermediateAssociationMatrix)
289  intermediateSortedClusters.push_back(mapEntry.first);
290  std::sort(intermediateSortedClusters.begin(), intermediateSortedClusters.end(), LArClusterHelper::SortByNHits);
291 
292  for (const Cluster *const pParentCluster : intermediateSortedClusters)
293  {
294  const ClusterAssociationMap &parentAssociationMap(intermediateAssociationMatrix.at(pParentCluster));
295 
296  ClusterVector sortedAssociationClusters;
297  for (const auto &mapEntry : parentAssociationMap)
298  sortedAssociationClusters.push_back(mapEntry.first);
299  std::sort(sortedAssociationClusters.begin(), sortedAssociationClusters.end(), LArClusterHelper::SortByNHits);
300 
301  for (const Cluster *const pDaughterCluster : sortedAssociationClusters)
302  {
303  const ClusterAssociation &parentToDaughterAssociation(parentAssociationMap.at(pDaughterCluster));
304 
305  ClusterAssociationMatrix::const_iterator iter5 = intermediateAssociationMatrix.find(pDaughterCluster);
306 
307  if (intermediateAssociationMatrix.end() == iter5)
308  continue;
309 
310  const ClusterAssociationMap &daughterAssociationMap(iter5->second);
311 
312  ClusterAssociationMap::const_iterator iter6 = daughterAssociationMap.find(pParentCluster);
313 
314  if (daughterAssociationMap.end() == iter6)
315  continue;
316 
317  const ClusterAssociation &daughterToParentAssociation(iter6->second);
318 
319  if (parentToDaughterAssociation.GetParent() == daughterToParentAssociation.GetDaughter() &&
320  parentToDaughterAssociation.GetDaughter() == daughterToParentAssociation.GetParent())
321  {
322  ClusterList &parentList(clusterMergeMap[pParentCluster]);
323 
324  if (parentList.end() == std::find(parentList.begin(), parentList.end(), pDaughterCluster))
325  parentList.push_back(pDaughterCluster);
326  }
327  }
328  }
329 }
330 
331 //------------------------------------------------------------------------------------------------------------------------------------------
332 
333 StatusCode CrossGapsExtensionAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
334 {
335  PANDORA_RETURN_RESULT_IF_AND_IF(
336  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterLength", m_minClusterLength));
337 
338  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinGapFraction", m_minGapFraction));
339 
340  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxGapTolerance", m_maxGapTolerance));
341 
342  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
343  XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
344 
345  float maxCosRelativeAngle(std::cos(m_maxRelativeAngle * M_PI / 180.0));
346  PANDORA_RETURN_RESULT_IF_AND_IF(
347  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxCosRelativeAngle", maxCosRelativeAngle));
348  m_maxRelativeAngle = (180.0 / M_PI) * std::acos(maxCosRelativeAngle);
349 
350  return ClusterExtensionAlgorithm::ReadSettings(xmlHandle);
351 }
352 
353 } // 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.
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...
static bool IsEmission(const pandora::CartesianVector &parentVertex, const LArPointingCluster::Vertex &daughterVertex, const float minLongitudinalDistance, const float maxLongitudinalDistance, const float maxTransverseDistance, const float angularAllowance)
Whether pointing vertex is emitted from a given position.
void BuildPointingClusterList(const pandora::ClusterVector &clusterVector, LArPointingClusterList &innerPointingClusterList, LArPointingClusterList &outerPointingClusterList) const
Build lists of pointing clusters that are adjacent to a detector gap.
std::vector< LArPointingCluster > LArPointingClusterList
static float CalculateGapDeltaZ(const pandora::Pandora &pandora, const float minZ, const float maxZ, const pandora::HitType hitType)
Calculate the total distance within a given 2D region that is composed of detector gaps...
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 FillClusterMergeMap(const ClusterAssociationMatrix &clusterAssociationMatrix, ClusterMergeMap &clusterMergeMap) const
Fill the cluster merge map.
intermediate_table::const_iterator const_iterator
LArPointingCluster class.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
static bool IsInGap(const pandora::Pandora &pandora, const pandora::CartesianVector &testPoint2D, const pandora::HitType hitType, const float gapTolerance=0.f)
Whether a 2D test point lies in a registered gap with the associated hit type.
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
const pandora::Cluster * GetCluster() const
Get the address of the cluster.
Header file for the cluster helper class.
void FillClusterAssociationMatrix(const pandora::ClusterVector &clusterVector, ClusterAssociationMatrix &clusterAssociationMatrix) const
Fill the cluster association matrix.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
std::unordered_map< const pandora::Cluster *, ClusterAssociation > ClusterAssociationMap
std::unordered_map< const pandora::Cluster *, ClusterAssociationMap > ClusterAssociationMatrix
bool IsAcrossGap(const float minZ, const float maxZ, const pandora::HitType hitType) const
Determine whether a start and end position sit either side of a gap.
HitType
Definition: HitType.h:12
bool IsAssociated(const LArPointingCluster::Vertex &pointingVertex1, const LArPointingCluster::Vertex &pointingVertex2) const
Use pointing information to determine whether two clusters are associated.
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
std::unordered_map< const pandora::Cluster *, pandora::ClusterList > ClusterMergeMap
std::vector< art::Ptr< recob::Cluster > > ClusterVector
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
Header file for the cross gaps extension algorithm class.