LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
SlidingConeClusterMopUpAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
14 
16 
18 
19 using namespace pandora;
20 
21 namespace lar_content
22 {
23 
24 SlidingConeClusterMopUpAlgorithm::SlidingConeClusterMopUpAlgorithm() :
25  m_useVertex(true),
26  m_maxIterations(1000),
27  m_maxHitsToConsider3DTrack(100),
28  m_minHitsToConsider3DShower(20),
29  m_maxHitsToConsider2DCluster(50),
30  m_halfWindowLayers(20),
31  m_nConeFitLayers(40),
32  m_nConeFits(5),
33  m_coneLengthMultiplier(3.f),
34  m_maxConeLength(126.f),
35  m_coneTanHalfAngle(0.2f),
36  m_coneBoundedFraction(0.5f)
37 {
38 }
39 
40 //------------------------------------------------------------------------------------------------------------------------------------------
41 
43 {
44  const Vertex *pVertex(nullptr);
45  this->GetInteractionVertex(pVertex);
46 
47  if (m_useVertex && !pVertex)
48  {
49  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
50  std::cout << "SlidingConeClusterMopUpAlgorithm - interaction vertex not available for use." << std::endl;
51  return STATUS_CODE_SUCCESS;
52  }
53 
54  ClusterVector clusters3D;
55  ClusterToPfoMap clusterToPfoMap;
56  this->GetThreeDClusters(clusters3D, clusterToPfoMap);
57 
58  ClusterVector availableClusters2D;
59  this->GetAvailableTwoDClusters(availableClusters2D);
60 
61  ClusterMergeMap clusterMergeMap;
62  this->GetClusterMergeMap(pVertex, clusters3D, availableClusters2D, clusterMergeMap);
63 
64  this->MakeClusterMerges(clusterToPfoMap, clusterMergeMap);
65 
66  return STATUS_CODE_SUCCESS;
67 }
68 
69 //------------------------------------------------------------------------------------------------------------------------------------------
70 
72 {
73  if (!m_useVertex)
74  return;
75 
76  const VertexList *pVertexList = nullptr;
77  PANDORA_THROW_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetCurrentList(*this, pVertexList));
78 
79  pVertex = ((pVertexList && (pVertexList->size() == 1) && (VERTEX_3D == (*(pVertexList->begin()))->GetVertexType())) ? *(pVertexList->begin()) : nullptr);
80 }
81 
82 //------------------------------------------------------------------------------------------------------------------------------------------
83 
85 {
86  for (const std::string &pfoListName : m_inputPfoListNames)
87  {
88  const PfoList *pPfoList(nullptr);
89 
90  if (STATUS_CODE_SUCCESS != PandoraContentApi::GetList(*this, pfoListName, pPfoList))
91  continue;
92 
93  for (const Pfo *const pPfo : *pPfoList)
94  {
95  ClusterList pfoClusters3D;
96  LArPfoHelper::GetThreeDClusterList(pPfo, pfoClusters3D);
97 
98  for (const Cluster *const pCluster3D : pfoClusters3D)
99  {
100  if (LArPfoHelper::IsTrack(pPfo) && (pCluster3D->GetNCaloHits() > m_maxHitsToConsider3DTrack))
101  continue;
102 
103  if (LArPfoHelper::IsShower(pPfo) && (pCluster3D->GetNCaloHits() < m_minHitsToConsider3DShower))
104  continue;
105 
106  if (!clusterToPfoMap.insert(ClusterToPfoMap::value_type(pCluster3D, pPfo)).second)
107  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
108 
109  clusters3D.push_back(pCluster3D);
110  }
111  }
112  }
113 
114  std::sort(clusters3D.begin(), clusters3D.end(), LArClusterHelper::SortByNHits);
115 }
116 
117 //------------------------------------------------------------------------------------------------------------------------------------------
118 
120 {
121  for (const std::string &clusterListName : m_daughterListNames)
122  {
123  const ClusterList *pClusterList(nullptr);
124 
125  if (STATUS_CODE_SUCCESS != PandoraContentApi::GetList(*this, clusterListName, pClusterList))
126  continue;
127 
128  for (const Cluster *const pCluster : *pClusterList)
129  {
130  if (!PandoraContentApi::IsAvailable(*this, pCluster))
131  continue;
132 
133  if (pCluster->GetNCaloHits() > m_maxHitsToConsider2DCluster)
134  continue;
135 
136  availableClusters2D.push_back(pCluster);
137  }
138  }
139 
140  std::sort(availableClusters2D.begin(), availableClusters2D.end(), LArClusterHelper::SortByNHits);
141 }
142 
143 //------------------------------------------------------------------------------------------------------------------------------------------
144 
145 void SlidingConeClusterMopUpAlgorithm::GetClusterMergeMap(const Vertex *const pVertex, const ClusterVector &clusters3D,
146  const ClusterVector &availableClusters2D, ClusterMergeMap &clusterMergeMap) const
147 {
148  const float layerPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));
149 
150  for (const Cluster *const pShowerCluster : clusters3D)
151  {
152  float coneLength3D(0.f);
153  SimpleConeList simpleConeList3D;
154 
155  try
156  {
157  const ThreeDSlidingConeFitResult slidingConeFitResult3D(pShowerCluster, m_halfWindowLayers, layerPitch);
158 
159  const CartesianVector &minLayerPosition(slidingConeFitResult3D.GetSlidingFitResult().GetGlobalMinLayerPosition());
160  const CartesianVector &maxLayerPosition(slidingConeFitResult3D.GetSlidingFitResult().GetGlobalMaxLayerPosition());
161  coneLength3D = std::min(m_coneLengthMultiplier * (maxLayerPosition - minLayerPosition).GetMagnitude(), m_maxConeLength);
162 
163  const float vertexToMinLayer(!pVertex ? 0.f : (pVertex->GetPosition() - minLayerPosition).GetMagnitude());
164  const float vertexToMaxLayer(!pVertex ? 0.f : (pVertex->GetPosition() - maxLayerPosition).GetMagnitude());
165  const ConeSelection coneSelection(!pVertex ? CONE_BOTH_DIRECTIONS : (vertexToMaxLayer > vertexToMinLayer) ? CONE_FORWARD_ONLY : CONE_BACKWARD_ONLY);
166 
167  slidingConeFitResult3D.GetSimpleConeList(m_nConeFitLayers, m_nConeFits, coneSelection, simpleConeList3D);
168  }
169  catch (const StatusCodeException &) {continue;}
170 
171  for (const Cluster *const pNearbyCluster2D : availableClusters2D)
172  {
173  ClusterMerge bestClusterMerge(nullptr, 0.f, 0.f);
174  const HitType hitType(LArClusterHelper::GetClusterHitType(pNearbyCluster2D));
175 
176  for (const SimpleCone &simpleCone3D : simpleConeList3D)
177  {
178  const CartesianVector coneBaseCentre3D(simpleCone3D.GetConeApex() + simpleCone3D.GetConeDirection() * coneLength3D);
179  const CartesianVector coneApex2D(LArGeometryHelper::ProjectPosition(this->GetPandora(), simpleCone3D.GetConeApex(), hitType));
180  const CartesianVector coneBaseCentre2D(LArGeometryHelper::ProjectPosition(this->GetPandora(), coneBaseCentre3D, hitType));
181 
182  const CartesianVector apexToBase2D(coneBaseCentre2D - coneApex2D);
183  const SimpleCone simpleCone2D(coneApex2D, apexToBase2D.GetUnitVector(), apexToBase2D.GetMagnitude(), m_coneTanHalfAngle);
184  const ClusterMerge clusterMerge(pShowerCluster, simpleCone2D.GetBoundedHitFraction(pNearbyCluster2D), simpleCone2D.GetMeanRT(pNearbyCluster2D));
185 
186  if (clusterMerge < bestClusterMerge)
187  bestClusterMerge = clusterMerge;
188  }
189 
190  if (bestClusterMerge.GetParentCluster() && (bestClusterMerge.GetBoundedFraction() > m_coneBoundedFraction))
191  clusterMergeMap[pNearbyCluster2D].push_back(bestClusterMerge);
192  }
193  }
194 
195  for (ClusterMergeMap::value_type &mapEntry : clusterMergeMap)
196  std::sort(mapEntry.second.begin(), mapEntry.second.end());
197 }
198 
199 //------------------------------------------------------------------------------------------------------------------------------------------
200 
201 void SlidingConeClusterMopUpAlgorithm::MakeClusterMerges(const ClusterToPfoMap &clusterToPfoMap, const ClusterMergeMap &clusterMergeMap) const
202 {
203  ClusterVector daughterClusters;
204  for (const ClusterMergeMap::value_type &mapEntry : clusterMergeMap) daughterClusters.push_back(mapEntry.first);
205  std::sort(daughterClusters.begin(), daughterClusters.end(), LArClusterHelper::SortByNHits);
206 
207  for (ClusterVector::const_reverse_iterator rIter = daughterClusters.rbegin(), rIterEnd = daughterClusters.rend(); rIter != rIterEnd; ++rIter)
208  {
209  const Cluster *const pDaughterCluster(*rIter);
210  const HitType daughterHitType(LArClusterHelper::GetClusterHitType(pDaughterCluster));
211  const Cluster *const pParentCluster3D(clusterMergeMap.at(pDaughterCluster).at(0).GetParentCluster());
212 
213  const Pfo *const pParentPfo(clusterToPfoMap.at(pParentCluster3D));
214  const Cluster *const pParentCluster(this->GetParentCluster(pParentPfo->GetClusterList(), daughterHitType));
215 
216  if (pParentCluster)
217  {
218  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::MergeAndDeleteClusters(*this, pParentCluster, pDaughterCluster,
219  this->GetListName(pParentCluster), this->GetListName(pDaughterCluster)));
220  }
221  else
222  {
223  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToPfo(*this, pParentPfo, pDaughterCluster));
224  }
225  }
226 }
227 
228 //------------------------------------------------------------------------------------------------------------------------------------------
229 //------------------------------------------------------------------------------------------------------------------------------------------
230 
232 {
233  if (!this->GetParentCluster() && !rhs.GetParentCluster())
234  return false;
235 
236  if (this->GetParentCluster() && !rhs.GetParentCluster())
237  return true;
238 
239  if (!this->GetParentCluster() && rhs.GetParentCluster())
240  return false;
241 
242  if (std::fabs(this->GetBoundedFraction() - rhs.GetBoundedFraction()) > std::numeric_limits<float>::epsilon())
243  return (this->GetBoundedFraction() > rhs.GetBoundedFraction());
244 
245  if (std::fabs(this->GetMeanRT() - rhs.GetMeanRT()) > std::numeric_limits<float>::epsilon())
246  return (this->GetMeanRT() < rhs.GetMeanRT());
247 
249 }
250 
251 //------------------------------------------------------------------------------------------------------------------------------------------
252 //------------------------------------------------------------------------------------------------------------------------------------------
253 
254 StatusCode SlidingConeClusterMopUpAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
255 {
256  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadVectorOfValues(xmlHandle,
257  "InputPfoListNames", m_inputPfoListNames));
258 
259  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
260  "UseVertex", m_useVertex));
261 
262  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
263  "MaxIterations", m_maxIterations));
264 
265  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
266  "MaxHitsToConsider3DTrack", m_maxHitsToConsider3DTrack));
267 
268  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
269  "MinHitsToConsider3DShower", m_minHitsToConsider3DShower));
270 
271  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
272  "MaxHitsToConsider2DCluster", m_maxHitsToConsider2DCluster));
273 
274  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
275  "SlidingFitHalfWindow", m_halfWindowLayers));
276 
277  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
278  "NConeFitLayers", m_nConeFitLayers));
279 
280  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
281  "NConeFits", m_nConeFits));
282 
283  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
284  "ConeLengthMultiplier", m_coneLengthMultiplier));
285 
286  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
287  "MaxConeLength", m_maxConeLength));
288 
289  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
290  "ConeTanHalfAngle", m_coneTanHalfAngle));
291 
292  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
293  "ConeBoundedFraction", m_coneBoundedFraction));
294 
295  return PfoMopUpBaseAlgorithm::ReadSettings(xmlHandle);
296 }
297 
298 } // 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.
bool m_useVertex
Whether to use the interaction vertex to select useful cone directions.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
Header file for the pfo helper class.
Header file for the sliding cone cluster mop up algorithm class.
std::unordered_map< const pandora::Cluster *, ClusterMergeList > ClusterMergeMap
void GetThreeDClusters(pandora::ClusterVector &clusters3D, ClusterToPfoMap &clusterToPfoMap) const
Get all 3d clusters contained in the input pfo lists and a mapping from clusters to pfos...
unsigned int m_nConeFits
The number of cone fits to perform, spread roughly uniformly along the shower length.
const pandora::Cluster * GetParentCluster() const
Get the address of the candidate parent (shower) cluster.
std::vector< SimpleCone > SimpleConeList
unsigned int m_maxIterations
The maximum allowed number of algorithm iterations.
float GetMeanRT() const
Get the mean transverse distance of all hits (whether contained or not)
pandora::StringVector m_daughterListNames
The list of potential daughter object list names.
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
float GetBoundedFraction() const
Get the bounded fraction for algorithm-specified cone angle.
unsigned int m_nConeFitLayers
The number of layers over which to sum fitted direction to obtain cone fit.
const ThreeDSlidingFitResult & GetSlidingFitResult() const
Get the sliding fit result for the full cluster.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
std::unordered_map< const pandora::Cluster *, const pandora::ParticleFlowObject * > ClusterToPfoMap
static bool IsTrack(const pandora::ParticleFlowObject *const pPfo)
Return track flag based on Pfo Particle ID.
unsigned int m_maxHitsToConsider2DCluster
The maximum number of hits in a 2d cluster to allow pick-up via sliding cone fits.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
Header file for the lar three dimensional sliding cone fit result class.
void GetSimpleConeList(const unsigned int nLayersForConeFit, const unsigned int nCones, const ConeSelection coneSelection, SimpleConeList &simpleConeList) const
Get the list of simple cones fitted to the three dimensional cluster.
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
const std::string GetListName(const T *const pT) const
Find the name of the list hosting a specific object.
static bool IsShower(const pandora::ParticleFlowObject *const pPfo)
Return shower flag based on Pfo Particle ID.
void GetClusterMergeMap(const pandora::Vertex *const pVertex, const pandora::ClusterVector &clusters3D, const pandora::ClusterVector &availableClusters2D, ClusterMergeMap &clusterMergeMap) const
Get the cluster merge map describing all potential 3d cluster merges.
Header file for the cluster helper class.
float m_coneLengthMultiplier
The cone length multiplier to use when calculating bounded cluster fractions.
unsigned int m_halfWindowLayers
The number of layers to use for half-window of sliding fit.
unsigned int m_minHitsToConsider3DShower
The minimum number of hits in a 3d shower cluster to attempt cone fits.
unsigned int m_maxHitsToConsider3DTrack
The maximum number of hits in a 3d track cluster to warrant inclusion in algorithm.
void MakeClusterMerges(const ClusterToPfoMap &clusterToPfoMap, const ClusterMergeMap &clusterMergeMap) const
Make cluster merges based on the provided cluster merge map.
static const pandora::Cluster * GetParentCluster(const pandora::ClusterList &clusterList, const pandora::HitType hitType)
Select the parent cluster (same hit type and most hits) using a provided cluster list and hit type...
static void GetThreeDClusterList(const pandora::ParticleFlowObject *const pPfo, pandora::ClusterList &clusterList)
Get the list of 3D clusters from an input pfo.
const pandora::CartesianVector & GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
Int_t min
Definition: plot.C:26
float m_coneBoundedFraction
The minimum cluster bounded fraction for association.
pandora::StringVector m_inputPfoListNames
The input pfo list names.
float m_coneTanHalfAngle
The cone tan half angle to use when calculating bounded cluster fractions.
ConeSelection
ConeSelection enum.
void GetAvailableTwoDClusters(pandora::ClusterVector &availableClusters2D) const
Get all available 2d clusters contained in the input cluster lists.
void GetInteractionVertex(const pandora::Vertex *&pVertex) const
Get the neutrino interaction vertex if it is available and if the algorithm is configured to do so...
float m_maxConeLength
The maximum allowed cone length to use when calculating bounded cluster fractions.
std::list< Vertex > VertexList
Definition: DCEL.h:178
const pandora::CartesianVector & GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.