LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CosmicRayExtensionAlgorithm.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 CosmicRayExtensionAlgorithm::CosmicRayExtensionAlgorithm() :
22  m_minClusterLength(3.f),
23  m_minSeedClusterLength(6.f),
24  m_maxLongitudinalDisplacement(30.f),
25  m_maxTransverseDisplacement(2.f),
26  m_minCosRelativeAngle(0.966f),
27  m_maxAverageRms(1.f)
28 {
29 }
30 
31 //------------------------------------------------------------------------------------------------------------------------------------------
32 
33 void CosmicRayExtensionAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const
34 {
35  for (ClusterList::const_iterator iter = pClusterList->begin(), iterEnd = pClusterList->end(); iter != iterEnd; ++iter)
36  {
37  const Cluster *const pCluster = *iter;
38 
40  continue;
41 
42  clusterVector.push_back(pCluster);
43  }
44 
45  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
46 }
47 
48 //------------------------------------------------------------------------------------------------------------------------------------------
49 
50 void CosmicRayExtensionAlgorithm::FillClusterAssociationMatrix(const ClusterVector &clusterVector, ClusterAssociationMatrix &clusterAssociationMatrix) const
51 {
52  // Convert each input cluster into a pointing cluster
53  LArPointingClusterList pointingClusterList;
54 
55  for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
56  {
57  try
58  {
59  pointingClusterList.push_back(LArPointingCluster(*iter));
60  }
61  catch (StatusCodeException &)
62  {
63  }
64  }
65 
66  // Form associations between pairs of pointing clusters
67  for (LArPointingClusterList::const_iterator iterI = pointingClusterList.begin(), iterEndI = pointingClusterList.end(); iterI != iterEndI; ++iterI)
68  {
69  const LArPointingCluster &clusterI = *iterI;
70 
71  for (LArPointingClusterList::const_iterator iterJ = iterI, iterEndJ = pointingClusterList.end(); iterJ != iterEndJ; ++iterJ)
72  {
73  const LArPointingCluster &clusterJ = *iterJ;
74 
75  if (clusterI.GetCluster() == clusterJ.GetCluster())
76  continue;
77 
78  this->FillClusterAssociationMatrix(clusterI, clusterJ, clusterAssociationMatrix);
79  }
80  }
81 }
82 
83 //------------------------------------------------------------------------------------------------------------------------------------------
84 
86  const LArPointingCluster &clusterI, const LArPointingCluster &clusterJ, ClusterAssociationMatrix &clusterAssociationMatrix) const
87 {
88  const Cluster *const pClusterI(clusterI.GetCluster());
89  const Cluster *const pClusterJ(clusterJ.GetCluster());
90 
91  if (pClusterI == pClusterJ)
92  return;
93 
94  // Get closest pair of vertices
95  LArPointingCluster::Vertex clusterVertexI, clusterVertexJ;
96 
97  try
98  {
99  LArPointingClusterHelper::GetClosestVertices(clusterI, clusterJ, clusterVertexI, clusterVertexJ);
100  }
101  catch (StatusCodeException &)
102  {
103  return;
104  }
105 
106  // (Just in case...)
107  if (!(clusterVertexI.IsInitialized() && clusterVertexJ.IsInitialized()))
108  throw StatusCodeException(STATUS_CODE_FAILURE);
109 
110  const CartesianVector vertexI(clusterVertexI.GetPosition());
111  const CartesianVector vertexJ(clusterVertexJ.GetPosition());
112 
113  const CartesianVector endI(clusterVertexI.IsInnerVertex() ? clusterI.GetOuterVertex().GetPosition() : clusterI.GetInnerVertex().GetPosition());
114  const CartesianVector endJ(clusterVertexJ.IsInnerVertex() ? clusterJ.GetOuterVertex().GetPosition() : clusterJ.GetInnerVertex().GetPosition());
115 
116  // Requirements on length
117  const float lengthSquaredI(LArPointingClusterHelper::GetLengthSquared(clusterI));
118  const float lengthSquaredJ(LArPointingClusterHelper::GetLengthSquared(clusterJ));
119 
120  if (std::max(lengthSquaredI, lengthSquaredJ) < m_minSeedClusterLength * m_minSeedClusterLength)
121  return;
122 
123  // Requirements on proximity
124  const float distanceSquaredIJ((vertexI - vertexJ).GetMagnitudeSquared());
125 
127  return;
128 
129  // Requirements on pointing information
130  const CartesianVector directionI((endI - vertexI).GetUnitVector());
131  const CartesianVector directionJ((endJ - vertexJ).GetUnitVector());
132 
133  const float cosTheta(-directionI.GetDotProduct(directionJ));
134  const float cosThetaCut(-1.f + 2.f * m_minCosRelativeAngle);
135 
136  if (cosTheta < cosThetaCut)
137  return;
138 
139  // Requirements on overlap between clusters
140  const CartesianVector directionIJ((endJ - endI).GetUnitVector());
141  const CartesianVector directionJI((endI - endJ).GetUnitVector());
142 
143  const float overlapL(directionIJ.GetDotProduct(vertexJ - vertexI));
144  const float overlapT(directionIJ.GetCrossProduct(vertexJ - vertexI).GetMagnitude());
145 
146  if (overlapL < -1.f || overlapL * overlapL > 2.f * std::min(lengthSquaredI, lengthSquaredJ) ||
147  overlapT > m_maxTransverseDisplacement + std::fabs(overlapL) * std::tan(5.f * M_PI / 180.f))
148  return;
149 
150  // Calculate RMS deviations on composite cluster
151  const float rms1(this->CalculateRms(pClusterI, endI, directionIJ));
152  const float rms2(this->CalculateRms(pClusterJ, endJ, directionJI));
153 
154  const float rms(0.5f * (rms1 + rms2));
155  const float rmsCut(2.f * m_maxAverageRms * (cosTheta - cosThetaCut) / (1.0 - cosThetaCut));
156 
157  if (rms > rmsCut)
158  return;
159 
160  // Record the association
164 
165  (void)clusterAssociationMatrix[pClusterI].insert(
166  ClusterAssociationMap::value_type(pClusterJ, ClusterAssociation(vertexTypeI, vertexTypeJ, associationType, lengthSquaredJ)));
167 
168  (void)clusterAssociationMatrix[pClusterJ].insert(
169  ClusterAssociationMap::value_type(pClusterI, ClusterAssociation(vertexTypeJ, vertexTypeI, associationType, lengthSquaredI)));
170 }
171 
172 //------------------------------------------------------------------------------------------------------------------------------------------
173 
174 void CosmicRayExtensionAlgorithm::FillClusterMergeMap(const ClusterAssociationMatrix &inputAssociationMatrix, ClusterMergeMap &clusterMergeMap) const
175 {
176  // Decide which associations will become merges
177  // To make the merge A <-> B, both A -> B and B -> A must be strong associations
178  // with the largest figures of merit of all the A -> X and B -> Y associations
179 
180  // First step: remove double-counting from the map of associations
181  // i.e. if the map has A <-> B, B <-> C, A <-> C, then remove A <-> C
182  ClusterAssociationMatrix clusterAssociationMatrix;
183 
184  ClusterVector sortedInputClusters;
185  for (const auto &mapEntry : inputAssociationMatrix)
186  sortedInputClusters.push_back(mapEntry.first);
187  std::sort(sortedInputClusters.begin(), sortedInputClusters.end(), LArClusterHelper::SortByNHits);
188 
189  for (const Cluster *const pCluster1 : sortedInputClusters)
190  {
191  const ClusterAssociationMap &associationMap1(inputAssociationMatrix.at(pCluster1));
192 
193  for (const Cluster *const pCluster2 : sortedInputClusters)
194  {
195  if (pCluster1 == pCluster2)
196  continue;
197 
198  const ClusterAssociationMap &associationMap2(inputAssociationMatrix.at(pCluster2));
199 
200  ClusterAssociationMap::const_iterator iter12 = associationMap1.find(pCluster2);
201  if (associationMap1.end() == iter12)
202  continue;
203 
204  ClusterAssociationMap::const_iterator iter21 = associationMap2.find(pCluster1);
205  if (associationMap2.end() == iter21)
206  continue;
207 
208  const ClusterAssociation &association12(iter12->second);
209  const ClusterAssociation &association21(iter21->second);
210 
211  bool isAssociated(true);
212 
213  ClusterVector sortedAssociationClusters;
214  for (const auto &mapEntry : associationMap1)
215  sortedAssociationClusters.push_back(mapEntry.first);
216  std::sort(sortedAssociationClusters.begin(), sortedAssociationClusters.end(), LArClusterHelper::SortByNHits);
217 
218  for (const Cluster *const pCluster3 : sortedAssociationClusters)
219  {
220  const ClusterAssociation &association13(associationMap1.at(pCluster3));
221 
222  ClusterAssociationMap::const_iterator iter23 = associationMap2.find(pCluster3);
223  if (associationMap2.end() == iter23)
224  continue;
225 
226  const ClusterAssociation &association23(iter23->second);
227 
228  if (association12.GetParent() == association13.GetParent() && association23.GetParent() == association21.GetParent() &&
229  association13.GetDaughter() != association23.GetDaughter())
230  {
231  isAssociated = false;
232  break;
233  }
234  }
235 
236  if (isAssociated)
237  {
238  (void)clusterAssociationMatrix[pCluster1].insert(ClusterAssociationMap::value_type(pCluster2, association12));
239  (void)clusterAssociationMatrix[pCluster2].insert(ClusterAssociationMap::value_type(pCluster1, association21));
240  }
241  }
242  }
243 
244  // Second step: find the best associations A -> X and B -> Y
245  ClusterAssociationMatrix intermediateAssociationMatrix;
246 
247  ClusterVector sortedClusters;
248  for (const auto &mapEntry : clusterAssociationMatrix)
249  sortedClusters.push_back(mapEntry.first);
250  std::sort(sortedClusters.begin(), sortedClusters.end(), LArClusterHelper::SortByNHits);
251 
252  for (const Cluster *const pParentCluster : sortedClusters)
253  {
254  const ClusterAssociationMap &clusterAssociationMap(clusterAssociationMatrix.at(pParentCluster));
255 
256  const Cluster *pBestClusterInner(nullptr);
258 
259  const Cluster *pBestClusterOuter(nullptr);
261 
262  ClusterVector sortedAssociationClusters;
263  for (const auto &mapEntry : clusterAssociationMap)
264  sortedAssociationClusters.push_back(mapEntry.first);
265  std::sort(sortedAssociationClusters.begin(), sortedAssociationClusters.end(), LArClusterHelper::SortByNHits);
266 
267  for (const Cluster *const pDaughterCluster : sortedAssociationClusters)
268  {
269  const ClusterAssociation &clusterAssociation(clusterAssociationMap.at(pDaughterCluster));
270 
271  // Inner associations
272  if (clusterAssociation.GetParent() == ClusterAssociation::INNER)
273  {
274  if (clusterAssociation.GetFigureOfMerit() > bestAssociationInner.GetFigureOfMerit())
275  {
276  bestAssociationInner = clusterAssociation;
277  pBestClusterInner = pDaughterCluster;
278  }
279  }
280 
281  // Outer associations
282  if (clusterAssociation.GetParent() == ClusterAssociation::OUTER)
283  {
284  if (clusterAssociation.GetFigureOfMerit() > bestAssociationOuter.GetFigureOfMerit())
285  {
286  bestAssociationOuter = clusterAssociation;
287  pBestClusterOuter = pDaughterCluster;
288  }
289  }
290  }
291 
292  if (pBestClusterInner)
293  (void)intermediateAssociationMatrix[pParentCluster].insert(ClusterAssociationMap::value_type(pBestClusterInner, bestAssociationInner));
294 
295  if (pBestClusterOuter)
296  (void)intermediateAssociationMatrix[pParentCluster].insert(ClusterAssociationMap::value_type(pBestClusterOuter, bestAssociationOuter));
297  }
298 
299  // Third step: make the merge if A -> X and B -> Y is in fact A -> B and B -> A
300  ClusterVector intermediateSortedClusters;
301  for (const auto &mapEntry : intermediateAssociationMatrix)
302  intermediateSortedClusters.push_back(mapEntry.first);
303  std::sort(intermediateSortedClusters.begin(), intermediateSortedClusters.end(), LArClusterHelper::SortByNHits);
304 
305  for (const Cluster *const pParentCluster : intermediateSortedClusters)
306  {
307  const ClusterAssociationMap &parentAssociationMap(intermediateAssociationMatrix.at(pParentCluster));
308 
309  ClusterVector sortedAssociationClusters;
310  for (const auto &mapEntry : parentAssociationMap)
311  sortedAssociationClusters.push_back(mapEntry.first);
312  std::sort(sortedAssociationClusters.begin(), sortedAssociationClusters.end(), LArClusterHelper::SortByNHits);
313 
314  for (const Cluster *const pDaughterCluster : sortedAssociationClusters)
315  {
316  const ClusterAssociation &parentToDaughterAssociation(parentAssociationMap.at(pDaughterCluster));
317 
318  ClusterAssociationMatrix::const_iterator iter5 = intermediateAssociationMatrix.find(pDaughterCluster);
319 
320  if (intermediateAssociationMatrix.end() == iter5)
321  continue;
322 
323  const ClusterAssociationMap &daughterAssociationMap(iter5->second);
324 
325  ClusterAssociationMap::const_iterator iter6 = daughterAssociationMap.find(pParentCluster);
326 
327  if (daughterAssociationMap.end() == iter6)
328  continue;
329 
330  const ClusterAssociation &daughterToParentAssociation(iter6->second);
331 
332  if (parentToDaughterAssociation.GetParent() == daughterToParentAssociation.GetDaughter() &&
333  parentToDaughterAssociation.GetDaughter() == daughterToParentAssociation.GetParent())
334  {
335  ClusterList &parentList(clusterMergeMap[pParentCluster]);
336 
337  if (parentList.end() == std::find(parentList.begin(), parentList.end(), pDaughterCluster))
338  parentList.push_back(pDaughterCluster);
339  }
340  }
341  }
342 }
343 
344 //------------------------------------------------------------------------------------------------------------------------------------------
345 
346 float CosmicRayExtensionAlgorithm::CalculateRms(const Cluster *const pCluster, const CartesianVector &position, const CartesianVector &direction) const
347 {
348  float totalChi2(0.f);
349  float totalHits(0.f);
350 
351  CaloHitList caloHitList;
352  pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
353 
354  for (CaloHitList::const_iterator iter = caloHitList.begin(), iterEnd = caloHitList.end(); iter != iterEnd; ++iter)
355  {
356  const CaloHit *const pCaloHit = *iter;
357 
358  const CartesianVector hitPosition(pCaloHit->GetPositionVector());
359  const CartesianVector predictedPosition(position + direction * direction.GetDotProduct(hitPosition - position));
360 
361  totalChi2 += (predictedPosition - hitPosition).GetMagnitudeSquared();
362  totalHits += 1.f;
363  }
364 
365  if (totalHits > 0.f)
366  return std::sqrt(totalChi2 / totalHits);
367 
368  return 0.f;
369 }
370 
371 //------------------------------------------------------------------------------------------------------------------------------------------
372 
373 StatusCode CosmicRayExtensionAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
374 {
375  PANDORA_RETURN_RESULT_IF_AND_IF(
376  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterLength", m_minClusterLength));
377 
378  PANDORA_RETURN_RESULT_IF_AND_IF(
379  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinSeedClusterLength", m_minSeedClusterLength));
380 
381  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
382  XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
383 
384  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
385  XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
386 
387  PANDORA_RETURN_RESULT_IF_AND_IF(
388  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCosRelativeAngle", m_minCosRelativeAngle));
389 
390  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxAverageRms", m_maxAverageRms));
391 
392  return ClusterExtensionAlgorithm::ReadSettings(xmlHandle);
393 }
394 
395 } // 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.
std::vector< LArPointingCluster > LArPointingClusterList
Header file for the cosmic-ray extension algorithm class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
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.
const Vertex & GetOuterVertex() const
Get the outer vertex.
const Vertex & GetInnerVertex() const
Get the inner vertex.
std::unordered_map< const pandora::Cluster *, ClusterAssociation > ClusterAssociationMap
std::unordered_map< const pandora::Cluster *, ClusterAssociationMap > ClusterAssociationMatrix
float CalculateRms(const pandora::Cluster *const pCluster, const pandora::CartesianVector &position, const pandora::CartesianVector &direction) const
Calculate RMS deviation of a cluster with respect to the reference line.
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 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
bool IsInnerVertex() const
Is this the inner vertex.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
static float GetLengthSquared(const LArPointingCluster &pointingCluster)
Calculate distance squared between inner and outer vertices of pointing cluster.
void FillClusterMergeMap(const ClusterAssociationMatrix &clusterAssociationMatrix, ClusterMergeMap &clusterMergeMap) const
Fill the cluster merge map.
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)