LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
CrossGapsAssociationAlgorithm.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 CrossGapsAssociationAlgorithm::CrossGapsAssociationAlgorithm() :
22  m_minClusterHits(10),
23  m_minClusterLayers(6),
24  m_slidingFitWindow(20),
25  m_maxSamplingPoints(1000),
26  m_sampleStepSize(0.5f),
27  m_maxUnmatchedSampleRun(8),
28  m_maxOnClusterDistance(1.5f),
29  m_minMatchedSamplingPoints(10),
30  m_minMatchedSamplingFraction(0.5f),
31  m_gapTolerance(0.f)
32 {
33 }
34 
35 //------------------------------------------------------------------------------------------------------------------------------------------
36 
37 void CrossGapsAssociationAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const
38 {
39  // ATTN May want to opt-out completely if no gap information available
40  // if (PandoraContentApi::GetGeometry(*this)->GetDetectorGapList().empty())
41  // return;
42 
43  for (const Cluster *const pCluster : *pClusterList)
44  {
45  if (pCluster->GetNCaloHits() < m_minClusterHits)
46  continue;
47 
48  if (1 + pCluster->GetOuterPseudoLayer() - pCluster->GetInnerPseudoLayer() < m_minClusterLayers)
49  continue;
50 
51  clusterVector.push_back(pCluster);
52  }
53 
54  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByInnerLayer);
55 }
56 
57 //------------------------------------------------------------------------------------------------------------------------------------------
58 
60 {
61  TwoDSlidingFitResultMap slidingFitResultMap;
62  const float slidingFitPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));
63 
64  for (const Cluster *const pCluster : clusterVector)
65  {
66  try {(void) slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(pCluster, TwoDSlidingFitResult(pCluster, m_slidingFitWindow, slidingFitPitch)));}
67  catch (StatusCodeException &) {}
68  }
69 
70  // ATTN This method assumes that clusters have been sorted by layer
71  for (ClusterVector::const_iterator iterI = clusterVector.begin(), iterIEnd = clusterVector.end(); iterI != iterIEnd; ++iterI)
72  {
73  const Cluster *const pInnerCluster = *iterI;
74  TwoDSlidingFitResultMap::const_iterator fitIterI = slidingFitResultMap.find(pInnerCluster);
75 
76  if (slidingFitResultMap.end() == fitIterI)
77  continue;
78 
79  for (ClusterVector::const_iterator iterJ = iterI, iterJEnd = clusterVector.end(); iterJ != iterJEnd; ++iterJ)
80  {
81  const Cluster *const pOuterCluster = *iterJ;
82 
83  if (pInnerCluster == pOuterCluster)
84  continue;
85 
86  TwoDSlidingFitResultMap::const_iterator fitIterJ = slidingFitResultMap.find(pOuterCluster);
87 
88  if (slidingFitResultMap.end() == fitIterJ)
89  continue;
90 
91  if (!this->AreClustersAssociated(fitIterI->second, fitIterJ->second))
92  continue;
93 
94  clusterAssociationMap[pInnerCluster].m_forwardAssociations.insert(pOuterCluster);
95  clusterAssociationMap[pOuterCluster].m_backwardAssociations.insert(pInnerCluster);
96  }
97  }
98 }
99 
100 //------------------------------------------------------------------------------------------------------------------------------------------
101 
102 bool CrossGapsAssociationAlgorithm::IsExtremalCluster(const bool isForward, const Cluster *const pCurrentCluster, const Cluster *const pTestCluster) const
103 {
104  const unsigned int currentLayer(isForward ? pCurrentCluster->GetOuterPseudoLayer() : pCurrentCluster->GetInnerPseudoLayer());
105  const unsigned int testLayer(isForward ? pTestCluster->GetOuterPseudoLayer() : pTestCluster->GetInnerPseudoLayer());
106 
107  if (isForward && ((testLayer > currentLayer) || ((testLayer == currentLayer) && LArClusterHelper::SortByNHits(pTestCluster, pCurrentCluster))))
108  return true;
109 
110  if (!isForward && ((testLayer < currentLayer) || ((testLayer == currentLayer) && LArClusterHelper::SortByNHits(pTestCluster, pCurrentCluster))))
111  return true;
112 
113  return false;
114 }
115 
116 //------------------------------------------------------------------------------------------------------------------------------------------
117 
119 {
120  if (outerFitResult.GetCluster()->GetInnerPseudoLayer() < innerFitResult.GetCluster()->GetInnerPseudoLayer())
121  throw pandora::StatusCodeException(STATUS_CODE_NOT_ALLOWED);
122 
123  if (outerFitResult.GetCluster()->GetInnerPseudoLayer() < innerFitResult.GetCluster()->GetOuterPseudoLayer())
124  return false;
125 
126  return (this->IsAssociated(innerFitResult.GetGlobalMaxLayerPosition(), innerFitResult.GetGlobalMaxLayerDirection(), outerFitResult) &&
127  this->IsAssociated(outerFitResult.GetGlobalMinLayerPosition(), outerFitResult.GetGlobalMinLayerDirection() * -1.f, innerFitResult));
128 }
129 
130 //------------------------------------------------------------------------------------------------------------------------------------------
131 
132 bool CrossGapsAssociationAlgorithm::IsAssociated(const CartesianVector &startPosition, const CartesianVector &startDirection,
133  const TwoDSlidingFitResult &targetFitResult) const
134 {
135  const HitType hitType(LArClusterHelper::GetClusterHitType(targetFitResult.GetCluster()));
136  unsigned int nSamplingPoints(0), nGapSamplingPoints(0), nMatchedSamplingPoints(0), nUnmatchedSampleRun(0);
137 
138  for (unsigned int iSample = 0; iSample < m_maxSamplingPoints; ++iSample)
139  {
140  ++nSamplingPoints;
141  const CartesianVector samplingPoint(startPosition + startDirection * static_cast<float>(iSample) * m_sampleStepSize);
142 
143  if (LArGeometryHelper::IsInGap(this->GetPandora(), samplingPoint, hitType, m_gapTolerance))
144  {
145  ++nGapSamplingPoints;
146  nUnmatchedSampleRun = 0; // ATTN Choose to also reset run when entering gap region
147  continue;
148  }
149 
150  if (this->IsNearCluster(samplingPoint, targetFitResult))
151  {
152  ++nMatchedSamplingPoints;
153  nUnmatchedSampleRun = 0;
154  }
155  else if (++nUnmatchedSampleRun > m_maxUnmatchedSampleRun)
156  {
157  break;
158  }
159  }
160 
161  const float expectation((targetFitResult.GetGlobalMaxLayerPosition() - targetFitResult.GetGlobalMinLayerPosition()).GetMagnitude() / m_sampleStepSize);
162  const float matchedSamplingFraction(expectation > 0.f ? static_cast<float>(nMatchedSamplingPoints) / expectation : 0.f);
163 
164  if ((nMatchedSamplingPoints > m_minMatchedSamplingPoints) || (matchedSamplingFraction > m_minMatchedSamplingFraction))
165  return true;
166 
167  return false;
168 }
169 
170 //------------------------------------------------------------------------------------------------------------------------------------------
171 
172 bool CrossGapsAssociationAlgorithm::IsNearCluster(const CartesianVector &samplingPoint, const TwoDSlidingFitResult &targetFitResult) const
173 {
175  targetFitResult.GetLocalPosition(samplingPoint, rL, rT);
176 
177  CartesianVector fitPosition(0.f, 0.f, 0.f);
178 
179  if (STATUS_CODE_SUCCESS == targetFitResult.GetGlobalFitPosition(rL, fitPosition))
180  {
181  if ((fitPosition - samplingPoint).GetMagnitudeSquared() < m_maxOnClusterDistance * m_maxOnClusterDistance)
182  return true;
183  }
184 
185  CartesianVector fitPositionAtX(0.f, 0.f, 0.f);
186 
187  if (STATUS_CODE_SUCCESS == targetFitResult.GetGlobalFitPositionAtX(samplingPoint.GetX(), fitPositionAtX))
188  {
189  if ((fitPositionAtX - samplingPoint).GetMagnitudeSquared() < m_maxOnClusterDistance * m_maxOnClusterDistance)
190  return true;
191  }
192 
193  return false;
194 }
195 
196 //------------------------------------------------------------------------------------------------------------------------------------------
197 
198 StatusCode CrossGapsAssociationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
199 {
200  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
201  "MinClusterHits", m_minClusterHits));
202 
203  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
204  "MinClusterLayers", m_minClusterLayers));
205 
206  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
207  "SlidingFitWindow", m_slidingFitWindow));
208 
209  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
210  "MaxSamplingPoints", m_maxSamplingPoints));
211 
212  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
213  "SampleStepSize", m_sampleStepSize));
214 
215  if (m_sampleStepSize < std::numeric_limits<float>::epsilon())
216  {
217  std::cout << "CrossGapsAssociationAlgorithm: Invalid value for SampleStepSize " << m_sampleStepSize << std::endl;
218  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
219  }
220 
221  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
222  "MaxUnmatchedSampleRun", m_maxUnmatchedSampleRun));
223 
224  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
225  "MaxOnClusterDistance", m_maxOnClusterDistance));
226 
227  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
228  "MinMatchedSamplingPoints", m_minMatchedSamplingPoints));
229 
230  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
231  "MinMatchedSamplingFraction", m_minMatchedSamplingFraction));
232 
233  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
234  "GapTolerance", m_gapTolerance));
235 
237 }
238 
239 } // namespace lar_content
float m_minMatchedSamplingFraction
Minimum ratio between matched sampling points and expectation to declare association.
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 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.
void PopulateClusterAssociationMap(const pandora::ClusterVector &clusterVector, ClusterAssociationMap &clusterAssociationMap) const
Populate the cluster association map.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
Header file for the cross gaps association algorithm class.
unsigned int m_minClusterLayers
The minimum allowed number of layers for a clean cluster.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
pandora::CartesianVector GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
unsigned int m_minClusterHits
The minimum allowed number of hits in a clean cluster.
unsigned int m_maxSamplingPoints
The maximum number of extension sampling points considered per association check. ...
bool IsAssociated(const pandora::CartesianVector &startPosition, const pandora::CartesianVector &startDirection, const TwoDSlidingFitResult &targetFitResult) const
Sample points along the extrapolation from a starting position to a target fit result to declare clus...
unsigned int m_slidingFitWindow
The layer window for the sliding linear fits.
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.
Int_t max
Definition: plot.C:27
float m_gapTolerance
The tolerance to use when querying whether a sampling point is in a gap, units cm.
intermediate_table::const_iterator const_iterator
Header file for the cluster helper class.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
unsigned int m_maxUnmatchedSampleRun
The maximum run of unmatched (and non-gap) samples to consider before stopping.
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
float m_maxOnClusterDistance
The maximum distance between a sampling point and sliding fit to target cluster.
const pandora::Cluster * GetCluster() const
Get the address of the cluster, if originally provided.
pandora::StatusCode GetGlobalFitPositionAtX(const float x, pandora::CartesianVector &position) const
Get global fit position for a given input x coordinate.
unsigned int m_minMatchedSamplingPoints
Minimum number of matched sampling points to declare association.
pandora::CartesianVector GetGlobalMaxLayerDirection() const
Get global direction corresponding to the fit result in maximum fit layer.
bool AreClustersAssociated(const TwoDSlidingFitResult &innerFitResult, const TwoDSlidingFitResult &outerFitResult) const
Determine whether two clusters are associated.
std::unordered_map< const pandora::Cluster *, ClusterAssociation > ClusterAssociationMap
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
bool IsNearCluster(const pandora::CartesianVector &samplingPoint, const TwoDSlidingFitResult &targetFitResult) const
Whether a sampling point lies near a target 2d sliding fit result.
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
bool IsExtremalCluster(const bool isForward, const pandora::Cluster *const pCurrentCluster, const pandora::Cluster *const pTestCluster) const
Determine which of two clusters is extremal.
float m_sampleStepSize
The sampling step size used in association checks, units cm.
static bool SortByInnerLayer(const pandora::Cluster *const pLhs, const pandora::Cluster *const pRhs)
Sort clusters by inner layer, then position, then pulse-height.