Public Member Functions | Private Member Functions | Private Attributes | List of all members
lar_content::CosmicRayExtensionAlgorithm Class Reference

CosmicRayExtensionAlgorithm class. More...

#include <CosmicRayExtensionAlgorithm.h>

Inheritance diagram for lar_content::CosmicRayExtensionAlgorithm:
lar_content::ClusterExtensionAlgorithm lar_content::ClusterMergingAlgorithm

Public Member Functions

 CosmicRayExtensionAlgorithm ()
 Default constructor. More...
 

Private Member Functions

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. More...
 
void FillClusterAssociationMatrix (const pandora::ClusterVector &clusterVector, ClusterAssociationMatrix &clusterAssociationMatrix) const
 Fill the cluster association matrix. More...
 
void FillClusterMergeMap (const ClusterAssociationMatrix &clusterAssociationMatrix, ClusterMergeMap &clusterMergeMap) const
 Fill the cluster merge map. More...
 
void FillClusterAssociationMatrix (const LArPointingCluster &clusterI, const LArPointingCluster &clusterJ, ClusterAssociationMatrix &clusterAssociationMatrix) const
 Form association between two pointing clusters. More...
 
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. More...
 
pandora::StatusCode ReadSettings (const pandora::TiXmlHandle xmlHandle)
 

Private Attributes

float m_minClusterLength
 
float m_minSeedClusterLength
 
float m_maxLongitudinalDisplacement
 
float m_maxTransverseDisplacement
 
float m_minCosRelativeAngle
 
float m_maxAverageRms
 

Additional Inherited Members

- Protected Types inherited from lar_content::ClusterExtensionAlgorithm
typedef std::unordered_map< const pandora::Cluster *, ClusterAssociationClusterAssociationMap
 
typedef std::unordered_map< const pandora::Cluster *, ClusterAssociationMapClusterAssociationMatrix
 
- Protected Types inherited from lar_content::ClusterMergingAlgorithm
typedef std::unordered_map< const pandora::Cluster *, pandora::ClusterList > ClusterMergeMap
 
- Protected Member Functions inherited from lar_content::ClusterExtensionAlgorithm
void PopulateClusterMergeMap (const pandora::ClusterVector &clusterVector, ClusterMergeMap &clusterMergeMatrix) const
 Form associations between pointing clusters. More...
 
- Protected Member Functions inherited from lar_content::ClusterMergingAlgorithm
virtual pandora::StatusCode Run ()
 
void MergeClusters (pandora::ClusterVector &clusterVector, ClusterMergeMap &clusterMergeMap) const
 Merge associated clusters. More...
 
void CollectAssociatedClusters (const pandora::Cluster *const pSeedCluster, const ClusterMergeMap &clusterMergeMap, pandora::ClusterList &associatedClusterList) const
 Collect up all clusters associations related to a given seed cluster. More...
 
void CollectAssociatedClusters (const pandora::Cluster *const pSeedCluster, const pandora::Cluster *const pCurrentCluster, const ClusterMergeMap &clusterMergeMap, const pandora::ClusterSet &clusterVetoList, pandora::ClusterList &associatedClusterList) const
 Collect up all clusters associations related to a given seed cluster. More...
 
void GetSortedListOfCleanClusters (const pandora::ClusterVector &inputClusters, pandora::ClusterVector &outputClusters) const
 Sort the selected clusters, so that they have a well-defined ordering. More...
 
- Protected Attributes inherited from lar_content::ClusterMergingAlgorithm
std::string m_inputClusterListName
 The name of the input cluster list. If not specified, will access current list. More...
 

Detailed Description

CosmicRayExtensionAlgorithm class.

Definition at line 23 of file CosmicRayExtensionAlgorithm.h.

Constructor & Destructor Documentation

lar_content::CosmicRayExtensionAlgorithm::CosmicRayExtensionAlgorithm ( )

Member Function Documentation

float lar_content::CosmicRayExtensionAlgorithm::CalculateRms ( const pandora::Cluster *const  pCluster,
const pandora::CartesianVector &  position,
const pandora::CartesianVector &  direction 
) const
private

Calculate RMS deviation of a cluster with respect to the reference line.

Parameters
pClusterthe input cluster
positionthe intercept of the reference line
directionthe direction of the reference line

Definition at line 346 of file CosmicRayExtensionAlgorithm.cc.

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 }
intermediate_table::const_iterator const_iterator
void lar_content::CosmicRayExtensionAlgorithm::FillClusterAssociationMatrix ( const pandora::ClusterVector &  clusterVector,
ClusterAssociationMatrix clusterAssociationMatrix 
) const
privatevirtual

Fill the cluster association matrix.

Parameters
clusterVectorthe input vector of clusters
clusterAssociationMatrixthe matrix of associations

Implements lar_content::ClusterExtensionAlgorithm.

void lar_content::CosmicRayExtensionAlgorithm::FillClusterAssociationMatrix ( const LArPointingCluster clusterI,
const LArPointingCluster clusterJ,
ClusterAssociationMatrix clusterAssociationMatrix 
) const
private

Form association between two pointing clusters.

Parameters
clusterIthe first pointing cluster
clusterJthe second pointing cluster
clusterAssociationMatrixthe matrix of cluster associations

Definition at line 85 of file CosmicRayExtensionAlgorithm.cc.

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
161  const ClusterAssociation::VertexType vertexTypeI(clusterVertexI.IsInnerVertex() ? ClusterAssociation::INNER : ClusterAssociation::OUTER);
162  const ClusterAssociation::VertexType vertexTypeJ(clusterVertexJ.IsInnerVertex() ? ClusterAssociation::INNER : ClusterAssociation::OUTER);
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 }
double rms(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:40
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.
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.
static int max(int a, int b)
#define M_PI
Definition: includeROOT.h:54
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
static float GetLengthSquared(const LArPointingCluster &pointingCluster)
Calculate distance squared between inner and outer vertices of pointing cluster.
void lar_content::CosmicRayExtensionAlgorithm::FillClusterMergeMap ( const ClusterAssociationMatrix clusterAssociationMatrix,
ClusterMergeMap clusterMergeMap 
) const
privatevirtual

Fill the cluster merge map.

Parameters
clusterAssociationMatrixthe matrix of cluster associations
clusterMergeMapthe map of cluster merges

Implements lar_content::ClusterExtensionAlgorithm.

Definition at line 174 of file CosmicRayExtensionAlgorithm.cc.

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);
257  ClusterAssociation bestAssociationInner(ClusterAssociation::UNDEFINED, ClusterAssociation::UNDEFINED, ClusterAssociation::NONE, 0.f);
258 
259  const Cluster *pBestClusterOuter(nullptr);
260  ClusterAssociation bestAssociationOuter(ClusterAssociation::UNDEFINED, ClusterAssociation::UNDEFINED, ClusterAssociation::NONE, 0.f);
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 }
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 *, ClusterAssociation > ClusterAssociationMap
std::unordered_map< const pandora::Cluster *, ClusterAssociationMap > ClusterAssociationMatrix
std::vector< art::Ptr< recob::Cluster > > ClusterVector
void lar_content::CosmicRayExtensionAlgorithm::GetListOfCleanClusters ( const pandora::ClusterList *const  pClusterList,
pandora::ClusterVector &  clusterVector 
) const
privatevirtual

Populate cluster vector with subset of cluster list, containing clusters judged to be clean.

Parameters
pClusterListaddress of the cluster list
clusterVectorto receive the populated cluster vector

Implements lar_content::ClusterMergingAlgorithm.

Definition at line 33 of file CosmicRayExtensionAlgorithm.cc.

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 }
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
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
StatusCode lar_content::CosmicRayExtensionAlgorithm::ReadSettings ( const pandora::TiXmlHandle  xmlHandle)
privatevirtual

Reimplemented from lar_content::ClusterMergingAlgorithm.

Definition at line 373 of file CosmicRayExtensionAlgorithm.cc.

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 }
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)

Member Data Documentation

float lar_content::CosmicRayExtensionAlgorithm::m_maxAverageRms
private

Definition at line 62 of file CosmicRayExtensionAlgorithm.h.

float lar_content::CosmicRayExtensionAlgorithm::m_maxLongitudinalDisplacement
private

Definition at line 59 of file CosmicRayExtensionAlgorithm.h.

float lar_content::CosmicRayExtensionAlgorithm::m_maxTransverseDisplacement
private

Definition at line 60 of file CosmicRayExtensionAlgorithm.h.

float lar_content::CosmicRayExtensionAlgorithm::m_minClusterLength
private

Definition at line 57 of file CosmicRayExtensionAlgorithm.h.

float lar_content::CosmicRayExtensionAlgorithm::m_minCosRelativeAngle
private

Definition at line 61 of file CosmicRayExtensionAlgorithm.h.

float lar_content::CosmicRayExtensionAlgorithm::m_minSeedClusterLength
private

Definition at line 58 of file CosmicRayExtensionAlgorithm.h.


The documentation for this class was generated from the following files: