LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DBScanAlg_tool.cc
Go to the documentation of this file.
1 
8 // Framework Includes
12 #include "cetlib/cpu_timer.h"
13 #include "fhiclcpp/ParameterSet.h"
15 
16 // LArSoft includes
22 
23 // std includes
24 #include <memory>
25 
26 //------------------------------------------------------------------------------------------------------------------------------------------
27 // implementation follows
28 
29 namespace lar_cluster3d {
30 
34  class DBScanAlg : virtual public IClusterAlg {
35  public:
41  explicit DBScanAlg(fhicl::ParameterSet const& pset);
42 
46  ~DBScanAlg();
47 
48  void configure(const fhicl::ParameterSet&) override;
49 
56  void Cluster3DHits(reco::HitPairList& hitPairList,
57  reco::ClusterParametersList& clusterParametersList) const override;
58 
65  void Cluster3DHits(reco::HitPairListPtr& hitPairList,
66  reco::ClusterParametersList& clusterParametersList) const override;
67 
71  float getTimeToExecute(IClusterAlg::TimeValues index) const override
72  {
73  return m_timeVector[index];
74  }
75 
76  private:
83  size_t) const;
84 
89  size_t m_minPairPts;
90  mutable std::vector<float> m_timeVector;
91 
92  std::unique_ptr<lar_cluster3d::IClusterParametersBuilder>
94  kdTree m_kdTree; // For the kdTree
95  };
96 
98  {
99  this->configure(pset);
100  }
101 
102  //------------------------------------------------------------------------------------------------------------------------------------------
103 
105 
106  //------------------------------------------------------------------------------------------------------------------------------------------
107 
109  {
110  m_enableMonitoring = pset.get<bool>("EnableMonitoring", true);
111  m_minPairPts = pset.get<size_t>("MinPairPts", 2);
112 
113  m_clusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
114  pset.get<fhicl::ParameterSet>("ClusterParamsBuilder"));
115 
116  // Recover the parameter set for the kdTree
117  fhicl::ParameterSet kdTreeParams(pset.get<fhicl::ParameterSet>("kdTree"));
118 
119  // Now work out the maximum wire pitch
121 
122  // Returns the wire pitch per plane assuming they will be the same for all TPCs
123  constexpr geo::TPCID tpcid{0, 0};
124  std::vector<double> const wirePitchVec{geometry->WirePitch(geo::PlaneID{tpcid, 0}),
125  geometry->WirePitch(geo::PlaneID{tpcid, 1}),
126  geometry->WirePitch(geo::PlaneID{tpcid, 2})};
127 
128  float maxBestDist = 1.99 * *std::max_element(wirePitchVec.begin(), wirePitchVec.end());
129 
130  kdTreeParams.put_or_replace<float>("RefLeafBestDist", maxBestDist);
131 
132  m_kdTree = kdTree(kdTreeParams);
133  }
134 
136  reco::ClusterParametersList& clusterParametersList) const
137  {
142  cet::cpu_timer theClockDBScan;
143 
144  m_timeVector.resize(NUMTIMEVALUES, 0.);
145 
146  // DBScan is driven of its "epsilon neighborhood". Computing adjacency within DBScan can be time
147  // consuming so the idea is the prebuild the adjaceny map and then run DBScan.
148  // We'll employ a kdTree to implement this scheme
149  kdTree::KdTreeNodeList kdTreeNodeContainer;
150  kdTree::KdTreeNode topNode = m_kdTree.BuildKdTree(hitPairList, kdTreeNodeContainer);
151 
153 
154  if (m_enableMonitoring) theClockDBScan.start();
155 
156  // Ok, here we go!
157  // The idea is to loop through all of the input 3D hits and do the clustering
158  for (const auto& hit : hitPairList) {
159  // Check if the hit has already been visited
160  if (hit.getStatusBits() & reco::ClusterHit3D::CLUSTERVISITED) continue;
161 
162  // Mark as visited
164 
165  // Find the neighborhood for this hit
166  kdTree::CandPairList candPairList;
167  float bestDistance(std::numeric_limits<float>::max());
168 
169  m_kdTree.FindNearestNeighbors(&hit, topNode, candPairList, bestDistance);
170 
171  if (candPairList.size() < m_minPairPts) {
173  }
174  else {
175  // "Create" a new cluster and get a reference to it
176  clusterParametersList.push_back(reco::ClusterParameters());
177 
178  reco::ClusterParameters& curCluster = clusterParametersList.back();
179 
181  curCluster.addHit3D(&hit);
182 
183  // expand the cluster
184  expandCluster(topNode, candPairList, curCluster, m_minPairPts);
185  }
186  }
187 
188  if (m_enableMonitoring) {
189  theClockDBScan.stop();
190 
191  m_timeVector[RUNDBSCAN] = theClockDBScan.accumulated_real_time();
192  }
193 
194  // Initial clustering is done, now trim the list and get output parameters
195  cet::cpu_timer theClockBuildClusters;
196 
197  // Start clocks if requested
198  if (m_enableMonitoring) theClockBuildClusters.start();
199 
200  m_clusterBuilder->BuildClusterInfo(clusterParametersList);
201 
202  if (m_enableMonitoring) {
203  theClockBuildClusters.stop();
204 
205  m_timeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
206  }
207 
208  mf::LogDebug("Cluster3D") << ">>>>> DBScan done, found " << clusterParametersList.size()
209  << " clusters" << std::endl;
210 
211  return;
212  }
213 
215  reco::ClusterParametersList& clusterParametersList) const
216  {
221  cet::cpu_timer theClockDBScan;
222 
223  m_timeVector.resize(NUMTIMEVALUES, 0.);
224 
225  // DBScan is driven of its "epsilon neighborhood". Computing adjacency within DBScan can be time
226  // consuming so the idea is the prebuild the adjaceny map and then run DBScan.
227  // We'll employ a kdTree to implement this scheme
228  kdTree::KdTreeNodeList kdTreeNodeContainer;
229  kdTree::KdTreeNode topNode = m_kdTree.BuildKdTree(hitPairList, kdTreeNodeContainer);
230 
232 
233  if (m_enableMonitoring) theClockDBScan.start();
234 
235  // Ok, here we go!
236  // The idea is to loop through all of the input 3D hits and do the clustering
237  for (const auto& hit : hitPairList) {
238  // Check if the hit has already been visited
239  if (hit->getStatusBits() & reco::ClusterHit3D::CLUSTERVISITED) continue;
240 
241  // Mark as visited
243 
244  // Find the neighborhood for this hit
245  kdTree::CandPairList candPairList;
246  float bestDistance(std::numeric_limits<float>::max());
247 
248  m_kdTree.FindNearestNeighbors(hit, topNode, candPairList, bestDistance);
249 
250  if (candPairList.size() < m_minPairPts) {
251  hit->setStatusBit(reco::ClusterHit3D::CLUSTERNOISE);
252  }
253  else {
254  // "Create" a new cluster and get a reference to it
255  clusterParametersList.push_back(reco::ClusterParameters());
256 
257  reco::ClusterParameters& curCluster = clusterParametersList.back();
258 
260  curCluster.addHit3D(hit);
261 
262  // expand the cluster
263  expandCluster(topNode, candPairList, curCluster, m_minPairPts);
264  }
265  }
266 
267  if (m_enableMonitoring) {
268  theClockDBScan.stop();
269 
270  m_timeVector[RUNDBSCAN] = theClockDBScan.accumulated_real_time();
271  }
272 
273  // Initial clustering is done, now trim the list and get output parameters
274  cet::cpu_timer theClockBuildClusters;
275 
276  // Start clocks if requested
277  if (m_enableMonitoring) theClockBuildClusters.start();
278 
279  m_clusterBuilder->BuildClusterInfo(clusterParametersList);
280 
281  if (m_enableMonitoring) {
282  theClockBuildClusters.stop();
283 
284  m_timeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
285  }
286 
287  mf::LogDebug("Cluster3D") << ">>>>> DBScan done, found " << clusterParametersList.size()
288  << " clusters" << std::endl;
289 
290  return;
291  }
292 
294  kdTree::CandPairList& candPairList,
296  size_t minPts) const
297  {
298  // This is the main inside loop for the DBScan based clustering algorithm
299 
300  // Loop over added hits until list has been exhausted
301  while (!candPairList.empty()) {
302  // Dereference the point so we can see in the debugger...
303  const reco::ClusterHit3D* neighborHit = candPairList.front().second;
304 
305  // Process if we've not been here before
306  if (!(neighborHit->getStatusBits() & reco::ClusterHit3D::CLUSTERVISITED)) {
307  // set as visited
309 
310  // get the neighborhood around this point
311  kdTree::CandPairList neighborCandPairList;
312  float bestDistance(std::numeric_limits<float>::max());
313 
314  m_kdTree.FindNearestNeighbors(neighborHit, topNode, neighborCandPairList, bestDistance);
315 
316  // If the epsilon neighborhood of this point is large enough then add its points to our list
317  if (neighborCandPairList.size() >= minPts) {
318  std::copy(neighborCandPairList.begin(),
319  neighborCandPairList.end(),
320  std::back_inserter(candPairList));
321  }
322  }
323 
324  // If the point is not yet in a cluster then we now add
325  if (!(neighborHit->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)) {
327  cluster.addHit3D(neighborHit);
328  }
329 
330  candPairList.pop_front();
331  }
332 
333  return;
334  }
335 
336  //------------------------------------------------------------------------------------------------------------------------------------------
337 
339 } // namespace lar_cluster3d
void expandCluster(const kdTree::KdTreeNode &, kdTree::CandPairList &, reco::ClusterParameters &, size_t) const
the main routine for DBScan
std::list< reco::ClusterHit3D > HitPairList
Definition: Cluster3D.h:330
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
kdTree class definiton
Definition: kdTree.h:31
void configure(const fhicl::ParameterSet &) override
Interface for configuring the particular algorithm tool.
bool m_enableMonitoring
Data members to follow.
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
std::list< KdTreeNode > KdTreeNodeList
Definition: kdTree.h:69
Implements a kdTree for use in clustering.
size_t FindNearestNeighbors(const reco::ClusterHit3D *, const KdTreeNode &, CandPairList &, float &) const
Definition: kdTree.cxx:183
Labelled "noise" by a clustering algorithm.
Definition: Cluster3D.h:105
Cluster finding and building.
IClusterAlg interface class definiton.
Definition: IClusterAlg.h:26
unsigned int getStatusBits() const
Definition: Cluster3D.h:154
define a kd tree node
Definition: kdTree.h:127
T get(std::string const &key) const
Definition: ParameterSet.h:314
TimeValues
enumerate the possible values for time checking if monitoring timing
Definition: IClusterAlg.h:61
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
This provides an art tool interface definition for 3D Cluster algorithms.
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
std::unique_ptr< lar_cluster3d::IClusterParametersBuilder > m_clusterBuilder
Common cluster builder tool.
Detector simulation of raw signals on wires.
std::list< CandPair > CandPairList
Definition: kdTree.h:84
void Cluster3DHits(reco::HitPairList &hitPairList, reco::ClusterParametersList &clusterParametersList) const override
Given a set of recob hits, run DBscan to form 3D clusters.
DBScanAlg(fhicl::ParameterSet const &pset)
Constructor.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
DBScanAlg class definiton.
float getTimeToExecute(IClusterAlg::TimeValues index) const override
If monitoring, recover the time to execute a particular function.
float getTimeToExecute() const
Definition: kdTree.h:104
"visited" by a clustering algorithm
Definition: Cluster3D.h:104
std::vector< float > m_timeVector
void put_or_replace(std::string const &key)
KdTreeNode & BuildKdTree(Hit3DVec::iterator, Hit3DVec::iterator, KdTreeNodeList &, int depth=0) const
Given an input set of ClusterHit3D objects, build a kd tree structure.
Definition: kdTree.cxx:109
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
art framework interface to geometry description
std::list< ClusterParameters > ClusterParametersList
Definition: Cluster3D.h:393
attached to a cluster
Definition: Cluster3D.h:106
void addHit3D(const reco::ClusterHit3D *hit3D)
Definition: Cluster3D.h:445
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:177