LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DeltaRayExtensionAlgorithm.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 DeltaRayExtensionAlgorithm::DeltaRayExtensionAlgorithm() :
21  m_minClusterLength(1.f),
22  m_maxClusterLength(10.f),
23  m_maxLongitudinalDisplacement(2.5f),
24  m_maxTransverseDisplacement(1.5f)
25 {
26 }
27 
28 //------------------------------------------------------------------------------------------------------------------------------------------
29 
30 void DeltaRayExtensionAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const
31 {
32  for (ClusterList::const_iterator iter = pClusterList->begin(), iterEnd = pClusterList->end(); iter != iterEnd; ++iter)
33  {
34  const Cluster *const pCluster = *iter;
35 
37  continue;
38 
39  clusterVector.push_back(pCluster);
40  }
41 
42  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
43 }
44 
45 //------------------------------------------------------------------------------------------------------------------------------------------
46 
47 void DeltaRayExtensionAlgorithm::FillClusterAssociationMatrix(const ClusterVector &clusterVector, ClusterAssociationMatrix &clusterAssociationMatrix) const
48 {
49  ClusterToCoordinateMap innerCoordinateMap, outerCoordinateMap;
50 
51  for (const Cluster *const pParentCluster : clusterVector)
52  {
53  for (const Cluster *const pDaughterCluster : clusterVector)
54  {
55  if (pParentCluster == pDaughterCluster)
56  continue;
57 
58  this->FillClusterAssociationMatrix(pParentCluster, pDaughterCluster, innerCoordinateMap, outerCoordinateMap, clusterAssociationMatrix);
59  }
60  }
61 }
62 
63 //------------------------------------------------------------------------------------------------------------------------------------------
64 
65 void DeltaRayExtensionAlgorithm::GetExtremalCoordinatesFromCache(const Cluster *const pCluster, ClusterToCoordinateMap &innerCoordinateMap,
66  ClusterToCoordinateMap &outerCoordinateMap, CartesianVector &innerCoordinate, CartesianVector &outerCoordinate) const
67 {
68  ClusterToCoordinateMap::const_iterator innerIter = innerCoordinateMap.find(pCluster);
69  ClusterToCoordinateMap::const_iterator outerIter = outerCoordinateMap.find(pCluster);
70 
71  if ((innerCoordinateMap.end() == innerIter) || (outerCoordinateMap.end() == outerIter))
72  {
73  LArClusterHelper::GetExtremalCoordinates(pCluster, innerCoordinate, outerCoordinate);
74  (void)innerCoordinateMap.insert(ClusterToCoordinateMap::value_type(pCluster, innerCoordinate));
75  (void)outerCoordinateMap.insert(ClusterToCoordinateMap::value_type(pCluster, outerCoordinate));
76  }
77  else
78  {
79  innerCoordinate = innerIter->second;
80  outerCoordinate = outerIter->second;
81  }
82 }
83 
84 //------------------------------------------------------------------------------------------------------------------------------------------
85 
86 void DeltaRayExtensionAlgorithm::FillClusterAssociationMatrix(const Cluster *const pParentCluster, const Cluster *const pDaughterCluster,
87  ClusterToCoordinateMap &innerCoordinateMap, ClusterToCoordinateMap &outerCoordinateMap, ClusterAssociationMatrix &clusterAssociationMatrix) const
88 {
89  // Daughter cluster must be available for any association to proceed
90  if (!PandoraContentApi::IsAvailable(*this, pDaughterCluster))
91  return;
92 
93  // Weak association: between parent cosmic-ray muon and daughter delta ray
94  // Strong association: between parent and daughter fragments of delta ray
95  // Figure of merit: distance between parent and daughter clusters
96  CartesianVector innerCoordinateP(0.f, 0.f, 0.f), outerCoordinateP(0.f, 0.f, 0.f);
97  this->GetExtremalCoordinatesFromCache(pParentCluster, innerCoordinateMap, outerCoordinateMap, innerCoordinateP, outerCoordinateP);
98 
99  CartesianVector innerCoordinateD(0.f, 0.f, 0.f), outerCoordinateD(0.f, 0.f, 0.f);
100  this->GetExtremalCoordinatesFromCache(pDaughterCluster, innerCoordinateMap, outerCoordinateMap, innerCoordinateD, outerCoordinateD);
101 
102  for (unsigned int useInnerD = 0; useInnerD < 2; ++useInnerD)
103  {
104  const CartesianVector daughterVertex(useInnerD == 1 ? innerCoordinateD : outerCoordinateD);
105  const CartesianVector daughterEnd(useInnerD == 1 ? outerCoordinateD : innerCoordinateD);
106 
107  const float daughterLengthSquared((daughterEnd - daughterVertex).GetMagnitudeSquared());
108 
109  // Daughter cluster must be available and below a length cut for any association
110  if (daughterLengthSquared > m_maxClusterLength * m_maxClusterLength)
111  continue;
112 
113  const CartesianVector projectedVertex(LArClusterHelper::GetClosestPosition(daughterVertex, pParentCluster));
114 
115  const float daughterVertexDistanceSquared((projectedVertex - daughterVertex).GetMagnitudeSquared());
116  const float daughterEndDistanceSquared((projectedVertex - daughterEnd).GetMagnitudeSquared());
117 
118  // Cut on proximity of daughter cluster to parent cluster
119  if (daughterVertexDistanceSquared > std::min(m_maxLongitudinalDisplacement * m_maxLongitudinalDisplacement,
120  std::min(daughterEndDistanceSquared, daughterLengthSquared)))
121  continue;
122 
123  const ClusterAssociation::VertexType daughterVertexType(useInnerD == 1 ? ClusterAssociation::INNER : ClusterAssociation::OUTER);
124  const float figureOfMerit(daughterVertexDistanceSquared);
125 
128 
129  for (unsigned int useInnerP = 0; useInnerP < 2; ++useInnerP)
130  {
131  const CartesianVector parentVertex(useInnerP == 1 ? innerCoordinateP : outerCoordinateP);
132  const CartesianVector parentEnd(useInnerP == 1 ? outerCoordinateP : innerCoordinateP);
133 
134  const float parentVertexDistanceSquared((projectedVertex - parentVertex).GetMagnitudeSquared());
135  const float parentEndDistanceSquared((projectedVertex - parentEnd).GetMagnitudeSquared());
136  const float parentLengthSquared((parentEnd - parentVertex).GetMagnitudeSquared());
137 
138  if (parentVertexDistanceSquared < parentEndDistanceSquared)
139  parentVertexType = (useInnerP == 1 ? ClusterAssociation::INNER : ClusterAssociation::OUTER);
140 
141  // Parent cluster must be available and below a length cut for a strong association
142  if (!PandoraContentApi::IsAvailable(*this, pParentCluster) || parentLengthSquared > m_maxClusterLength * m_maxClusterLength)
143  continue;
144 
145  // Require an end-to-end join between parent and daughter cluster
146  if (parentVertexDistanceSquared > std::min(m_maxTransverseDisplacement * m_maxTransverseDisplacement,
147  std::min(parentEndDistanceSquared, daughterEndDistanceSquared)))
148  continue;
149 
150  // Cut on pointing information
151  const CartesianVector daughterDirection((daughterEnd - daughterVertex).GetUnitVector());
152  const CartesianVector parentDirection((parentEnd - projectedVertex).GetUnitVector());
153 
154  const float forwardDistance(daughterDirection.GetDotProduct((daughterVertex - projectedVertex)));
155  const float sidewaysDistance(daughterDirection.GetCrossProduct((daughterVertex - projectedVertex)).GetMagnitude());
156 
157  if (forwardDistance < 0.f || forwardDistance > m_maxLongitudinalDisplacement || sidewaysDistance > m_maxTransverseDisplacement)
158  continue;
159 
160  if (-parentDirection.GetDotProduct(daughterDirection) < 0.25f)
161  continue;
162 
163  associationType = ClusterAssociation::STRONG;
164  break;
165  }
166 
167  if (parentVertexType > ClusterAssociation::UNDEFINED)
168  {
169  (void)clusterAssociationMatrix[pParentCluster].insert(ClusterAssociationMap::value_type(
170  pDaughterCluster, ClusterAssociation(parentVertexType, daughterVertexType, associationType, figureOfMerit)));
171  }
172  }
173 }
174 
175 //------------------------------------------------------------------------------------------------------------------------------------------
176 
177 void DeltaRayExtensionAlgorithm::FillClusterMergeMap(const ClusterAssociationMatrix &parentToDaughterMatrix, ClusterMergeMap &clusterMergeMap) const
178 {
179  // Merge parent and daughter clusters if they are strongly associated
180  // and the associations have the best figures of merit
181  // (i.e. the P --> D association is the best P --> X association,
182  // and the P <-- D association is the best X <-- D association).
183  ClusterAssociationMatrix daughterToParentMatrix;
184 
185  ClusterVector sortedParentClusters;
186  for (const auto &mapEntry : parentToDaughterMatrix)
187  sortedParentClusters.push_back(mapEntry.first);
188  std::sort(sortedParentClusters.begin(), sortedParentClusters.end(), LArClusterHelper::SortByNHits);
189 
190  for (const Cluster *const pParentCluster : sortedParentClusters)
191  {
192  const ClusterAssociationMap &daughterToAssociationMap(parentToDaughterMatrix.at(pParentCluster));
193 
194  ClusterVector sortedLocalDaughterClusters;
195  for (const auto &mapEntry : daughterToAssociationMap)
196  sortedLocalDaughterClusters.push_back(mapEntry.first);
197  std::sort(sortedLocalDaughterClusters.begin(), sortedLocalDaughterClusters.end(), LArClusterHelper::SortByNHits);
198 
199  for (const Cluster *const pDaughterCluster : sortedLocalDaughterClusters)
200  {
201  const ClusterAssociation &clusterAssociation(daughterToAssociationMap.at(pDaughterCluster));
202  (void)daughterToParentMatrix[pDaughterCluster].insert(ClusterAssociationMap::value_type(pParentCluster, clusterAssociation));
203  }
204  }
205 
206  ClusterAssociationMatrix reducedParentToDaughterMatrix;
207 
208  ClusterVector sortedDaughterClusters;
209  for (const auto &mapEntry : daughterToParentMatrix)
210  sortedDaughterClusters.push_back(mapEntry.first);
211  std::sort(sortedDaughterClusters.begin(), sortedDaughterClusters.end(), LArClusterHelper::SortByNHits);
212 
213  // Loop over parent clusters and select nearby daughter clusters that are closer than another parent cluster
214  for (const Cluster *const pDaughterCluster : sortedDaughterClusters)
215  {
216  const ClusterAssociationMap &parentToAssociationMap(daughterToParentMatrix.at(pDaughterCluster));
217 
218  const Cluster *pBestInner(NULL);
219  const Cluster *pBestOuter(NULL);
220 
221  float bestFomInner(std::numeric_limits<float>::max());
222  float bestFomOuter(std::numeric_limits<float>::max());
223 
224  ClusterVector sortedLocalParentClusters;
225  for (const auto &mapEntry : parentToAssociationMap)
226  sortedLocalParentClusters.push_back(mapEntry.first);
227  std::sort(sortedLocalParentClusters.begin(), sortedLocalParentClusters.end(), LArClusterHelper::SortByNHits);
228 
229  for (const Cluster *const pParentCluster : sortedLocalParentClusters)
230  {
231  const ClusterAssociation &clusterAssociation(parentToAssociationMap.at(pParentCluster));
232 
233  if (clusterAssociation.GetParent() == ClusterAssociation::INNER)
234  {
235  if (clusterAssociation.GetFigureOfMerit() < bestFomInner)
236  {
237  bestFomInner = clusterAssociation.GetFigureOfMerit();
238 
239  if (clusterAssociation.GetAssociation() == ClusterAssociation::STRONG)
240  {
241  pBestInner = pParentCluster;
242  }
243  else
244  {
245  pBestInner = NULL;
246  }
247  }
248  }
249 
250  if (clusterAssociation.GetParent() == ClusterAssociation::OUTER)
251  {
252  if (clusterAssociation.GetFigureOfMerit() < bestFomOuter)
253  {
254  bestFomOuter = clusterAssociation.GetFigureOfMerit();
255 
256  if (clusterAssociation.GetAssociation() == ClusterAssociation::STRONG)
257  {
258  pBestOuter = pParentCluster;
259  }
260  else
261  {
262  pBestOuter = NULL;
263  }
264  }
265  }
266  }
267 
268  if (pBestInner)
269  {
270  ClusterAssociationMatrix::const_iterator iter3A = parentToDaughterMatrix.find(pBestInner);
271 
272  if (parentToDaughterMatrix.end() == iter3A)
273  throw pandora::StatusCodeException(STATUS_CODE_FAILURE);
274 
275  const ClusterAssociationMap &parentToDaughterMap(iter3A->second);
276  ClusterAssociationMap::const_iterator iter3B = parentToDaughterMap.find(pDaughterCluster);
277 
278  if (parentToDaughterMap.end() == iter3B)
279  throw pandora::StatusCodeException(STATUS_CODE_FAILURE);
280 
281  const ClusterAssociation &bestAssociationInner(iter3B->second);
282  (void)reducedParentToDaughterMatrix[pBestInner].insert(ClusterAssociationMap::value_type(pDaughterCluster, bestAssociationInner));
283  }
284 
285  if (pBestOuter)
286  {
287  ClusterAssociationMatrix::const_iterator iter3A = parentToDaughterMatrix.find(pBestOuter);
288 
289  if (parentToDaughterMatrix.end() == iter3A)
290  throw pandora::StatusCodeException(STATUS_CODE_FAILURE);
291 
292  const ClusterAssociationMap &parentToDaughterMap(iter3A->second);
293  ClusterAssociationMap::const_iterator iter3B = parentToDaughterMap.find(pDaughterCluster);
294 
295  if (parentToDaughterMap.end() == iter3B)
296  throw pandora::StatusCodeException(STATUS_CODE_FAILURE);
297 
298  const ClusterAssociation &bestAssociationOuter(iter3B->second);
299  (void)reducedParentToDaughterMatrix[pBestOuter].insert(ClusterAssociationMap::value_type(pDaughterCluster, bestAssociationOuter));
300  }
301  }
302 
303  ClusterVector sortedReducedParentClusters;
304  for (const auto &mapEntry : reducedParentToDaughterMatrix)
305  sortedReducedParentClusters.push_back(mapEntry.first);
306  std::sort(sortedReducedParentClusters.begin(), sortedReducedParentClusters.end(), LArClusterHelper::SortByNHits);
307 
308  for (const Cluster *const pParentCluster : sortedReducedParentClusters)
309  {
310  const ClusterAssociationMap &daughterToAssociationMap(reducedParentToDaughterMatrix.at(pParentCluster));
311 
312  const Cluster *pBestInner(NULL);
313  const Cluster *pBestOuter(NULL);
314 
315  float bestFomInner(std::numeric_limits<float>::max());
316  float bestFomOuter(std::numeric_limits<float>::max());
317 
318  ClusterVector sortedLocalDaughterClusters;
319  for (const auto &mapEntry : daughterToAssociationMap)
320  sortedLocalDaughterClusters.push_back(mapEntry.first);
321  std::sort(sortedLocalDaughterClusters.begin(), sortedLocalDaughterClusters.end(), LArClusterHelper::SortByNHits);
322 
323  for (const Cluster *const pDaughterCluster : sortedLocalDaughterClusters)
324  {
325  const ClusterAssociation &clusterAssociation(daughterToAssociationMap.at(pDaughterCluster));
326 
327  if (clusterAssociation.GetParent() == ClusterAssociation::INNER)
328  {
329  if (clusterAssociation.GetFigureOfMerit() < bestFomInner)
330  {
331  bestFomInner = clusterAssociation.GetFigureOfMerit();
332 
333  if (clusterAssociation.GetAssociation() == ClusterAssociation::STRONG)
334  {
335  pBestInner = pDaughterCluster;
336  }
337  else
338  {
339  pBestInner = NULL;
340  }
341  }
342  }
343 
344  if (clusterAssociation.GetParent() == ClusterAssociation::OUTER)
345  {
346  if (clusterAssociation.GetFigureOfMerit() < bestFomOuter)
347  {
348  bestFomOuter = clusterAssociation.GetFigureOfMerit();
349 
350  if (clusterAssociation.GetAssociation() == ClusterAssociation::STRONG)
351  {
352  pBestOuter = pDaughterCluster;
353  }
354  else
355  {
356  pBestOuter = NULL;
357  }
358  }
359  }
360  }
361 
362  if (pBestInner)
363  {
364  ClusterList &parentList(clusterMergeMap[pParentCluster]);
365 
366  if (parentList.end() == std::find(parentList.begin(), parentList.end(), pBestInner))
367  parentList.push_back(pBestInner);
368 
369  ClusterList &bestInnerList(clusterMergeMap[pBestInner]);
370 
371  if (bestInnerList.end() == std::find(bestInnerList.begin(), bestInnerList.end(), pParentCluster))
372  bestInnerList.push_back(pParentCluster);
373  }
374 
375  if (pBestOuter)
376  {
377  ClusterList &parentList(clusterMergeMap[pParentCluster]);
378 
379  if (parentList.end() == std::find(parentList.begin(), parentList.end(), pBestOuter))
380  parentList.push_back(pBestOuter);
381 
382  ClusterList &bestOuterList(clusterMergeMap[pBestOuter]);
383 
384  if (bestOuterList.end() == std::find(bestOuterList.begin(), bestOuterList.end(), pParentCluster))
385  bestOuterList.push_back(pParentCluster);
386  }
387  }
388 }
389 
390 //------------------------------------------------------------------------------------------------------------------------------------------
391 
392 StatusCode DeltaRayExtensionAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
393 {
394  PANDORA_RETURN_RESULT_IF_AND_IF(
395  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterLength", m_minClusterLength));
396 
397  PANDORA_RETURN_RESULT_IF_AND_IF(
398  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxClusterLength", m_maxClusterLength));
399 
400  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
401  XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
402 
403  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
404  XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
405 
406  return ClusterExtensionAlgorithm::ReadSettings(xmlHandle);
407 }
408 
409 } // 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.
intermediate_table::const_iterator const_iterator
std::unordered_map< const pandora::Cluster *, pandora::CartesianVector > ClusterToCoordinateMap
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
TFile f
Definition: plotHisto.C:6
Header file for the delta ray extension algorithm class.
Header file for the cluster helper class.
std::unordered_map< const pandora::Cluster *, ClusterAssociation > ClusterAssociationMap
std::unordered_map< const pandora::Cluster *, ClusterAssociationMap > ClusterAssociationMatrix
void FillClusterMergeMap(const ClusterAssociationMatrix &clusterAssociationMatrix, ClusterMergeMap &clusterMergeMap) const
Fill the cluster merge map.
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 GetExtremalCoordinates(const pandora::ClusterList &clusterList, pandora::CartesianVector &innerCoordinate, pandora::CartesianVector &outerCoordinate)
Get positions of the two most distant calo hits in a list of cluster (ordered by Z) ...
void GetExtremalCoordinatesFromCache(const pandora::Cluster *const pCluster, ClusterToCoordinateMap &innerCoordinateMap, ClusterToCoordinateMap &outerCoordinateMap, pandora::CartesianVector &innerCoordinate, pandora::CartesianVector &outerCoordinate) const
Reduce number of extremal coordinates calculations by caching results when they are first obtained...
void FillClusterAssociationMatrix(const pandora::ClusterVector &clusterVector, ClusterAssociationMatrix &clusterAssociationMatrix) const
Fill the cluster association matrix.
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
static pandora::CartesianVector GetClosestPosition(const pandora::CartesianVector &position, const pandora::ClusterList &clusterList)
Get closest position in a list of clusters to a specified input position vector.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)