HoughSeedFinderAlg.cxx
Go to the documentation of this file.
1 /**
2  * @file HoughSeedFinderAlg.cxx
3  *
4  * @brief Implementation of the Seed Finder Algorithm based on a Hough Transform
5  *
6  * The intent of this algorithm is to take an input list of 3D space points and from those
7  * to find candidate track start points and directions. It does so by first performing a
8  * Principal Components Analysis of the input 3D hits and then projects them to the plane
9  * of largest spread. A standard Hough Transform method is then applied to attempt to identify
10  * straight line segments which can be used as seeds to the kalman filter tracker.
11  */
12 
13 // The main include
15 // Framework Includes
16 
17 // LArSoft includes
22 #include "lardata/RecoObjects/Cluster3D.h"
25 
26 // ROOT includes
27 #include "TTree.h"
28 #include "TCanvas.h"
29 #include "TFrame.h"
30 #include "TH2D.h"
31 #include "TVectorD.h"
32 #include "TMatrixD.h"
33 #include "TDecompSVD.h"
34 
35 
36 // std includes
37 #include <string>
38 #include <functional>
39 #include <iostream>
40 #include <memory>
41 
42 //------------------------------------------------------------------------------------------------------------------------------------------
43 
44 //------------------------------------------------------------------------------------------------------------------------------------------
45 // implementation follows
46 
47 namespace lar_cluster3d {
48 
50  m_minimum3DHits(5),
51  m_thetaBins(360),
52  m_rhoBins(75),
53  m_hiThresholdMin(5),
54  m_hiThresholdFrac(.05),
55  m_loThresholdFrac(0.85),
56  m_numSeed2DHits(80),
57  m_numAveDocas(6.),
58  m_numSkippedHits(10),
59  m_maxLoopsPerCluster(3),
60  m_maximumGap(5.),
61  m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg")),
62  m_displayHist(false)
63 {
64  this->reconfigure(pset);
65 
67 
68  m_geometry = &*geometry;
69  // m_detector = detectorProperties->provider();
70 }
71 
72 //------------------------------------------------------------------------------------------------------------------------------------------
73 
75 {
76 }
77 
79 {
80  m_minimum3DHits = pset.get<size_t>("Minimum3DHits", 5);
81  m_thetaBins = pset.get<int> ("ThetaBins", 360);
82  m_rhoBins = pset.get<int> ("HalfRhoBins", 75);
83  m_hiThresholdMin = pset.get<size_t>("HiThresholdMin", 5);
84  m_hiThresholdFrac = pset.get<double>("HiThresholdFrac", 0.05);
85  m_loThresholdFrac = pset.get<double>("LoThresholdFrac", 0.85);
86  m_numSeed2DHits = pset.get<size_t>("NumSeed2DHits", 80);
87  m_numAveDocas = pset.get<double>("NumAveDocas", 6.);
88  m_numSkippedHits = pset.get<int> ("NumSkippedHits", 10);
89  m_maxLoopsPerCluster = pset.get<int>("MaxLoopsPerCluster", 3);
90  m_maximumGap = pset.get<double>("MaximumGap", 5.);
91  m_displayHist = pset.get<bool> ("DisplayHoughHist", false);
92 
93  m_pcaAlg.reconfigure(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"));
94 }
95 
97 {
98  /**
99  * @brief A utility class to contain the values of a given "bin" in Hough Space
100  *
101  * Specifically, this is keeping track of the projected x,y coordinates of a given
102  * 3D hit projected to the plane of largest spread in PCA and an interator to
103  * that hit in the input container
104  */
105 public:
106  AccumulatorValues() : m_position(0.,0.,0.) {}
108  m_position(position), m_hit3DIterator(itr) {}
109 
110  const TVector3& getPosition() const {return m_position;}
111  reco::HitPairListPtr::const_iterator getHitIterator() const {return m_hit3DIterator;}
112 
113 private:
114  TVector3 m_position; ///< We really only need the x,y coordinates here but keep all three for now
115  reco::HitPairListPtr::const_iterator m_hit3DIterator; ///< This will be used to take us back to our 3D hit
116 };
117 
118 typedef std::vector<AccumulatorValues> AccumulatorValuesVec;
119 
121 {
122  /**
123  * @brief A utility class used to accumulate the above values
124  *
125  * One of these objects will exist for each "bin" in rho-theta space and this will
126  * be used to accumulate the 3D hits which contribute to this bin
127  */
128 public:
129  AccumulatorBin() : m_visited(false), m_noise(false), m_inCluster(false) {}
130  AccumulatorBin(AccumulatorValues& values) : m_visited(false), m_noise(false), m_inCluster(false)
131  {
132  m_accumulatorValuesVec.push_back(values);
133  }
134 
135  void setVisited() {m_visited = true;}
136  void setNoise() {m_noise = true;}
137  void setInCluster() {m_inCluster = true;}
138 
139  void addAccumulatorValue(AccumulatorValues& value) {m_accumulatorValuesVec.push_back(value);}
140 
141  bool isVisited() const {return m_visited;}
142  bool isNoise() const {return m_noise;}
143  bool isInCluster() const {return m_inCluster;}
144 
145  const AccumulatorValuesVec& getAccumulatorValues() const {return m_accumulatorValuesVec;}
146 private:
147  bool m_visited;
148  bool m_noise;
150  AccumulatorValuesVec m_accumulatorValuesVec;
151 };
152 
154 {
155  /**
156  * @brief This is used to sort "Hough Clusters" by the maximum entries in a bin
157  */
158 public:
160 
162  {
163  size_t peakCountLeft(0);
164  size_t peakCountRight(0);
165 
166  for(const auto& binIndex : left)
167  peakCountLeft = std::max(peakCountLeft, m_accMap[binIndex].getAccumulatorValues().size());
168  for(const auto& binIndex : right)
169  peakCountRight = std::max(peakCountRight, m_accMap[binIndex].getAccumulatorValues().size());
170 
171  return peakCountLeft > peakCountRight;
172  }
173 private:
175 };
176 
178 {
179  /**
180  * @brief This is used to sort bins in Hough Space
181  */
183  {
184  size_t leftSize = left->second.getAccumulatorValues().size();
185  size_t rightSize = right->second.getAccumulatorValues().size();
186 
187  return leftSize > rightSize;
188  }
189 };
190 
191 
193  RhoThetaAccumulatorBinMap& rhoThetaAccumulatorBinMap,
194  HoughCluster& neighborPts,
195  size_t threshold) const
196 {
197  /**
198  * @brief Does a query of nearest neighbors to look for matching bins
199  */
200 
201  // We simply loop over the nearest indices and see if we have any friends over threshold
202  for(int rhoIdx = curBin.first - 1; rhoIdx <= curBin.first + 1; rhoIdx++)
203  {
204  for(int jdx = curBin.second - 1; jdx <= curBin.second + 1; jdx++)
205  {
206  // Skip the self reference
207  if (rhoIdx == curBin.first && jdx == curBin.second) continue;
208 
209  // Theta bin needs to handle the wrap.
210  int thetaIdx(jdx);
211 
212  if (thetaIdx < 0) thetaIdx = m_thetaBins - 1;
213  else if (thetaIdx > m_thetaBins -1) thetaIdx = 0;
214 
215  BinIndex binIndex(rhoIdx,thetaIdx);
216  RhoThetaAccumulatorBinMap::iterator mapItr = rhoThetaAccumulatorBinMap.find(binIndex);
217 
218  if (mapItr != rhoThetaAccumulatorBinMap.end())
219  {
220  if (mapItr->second.getAccumulatorValues().size() >= threshold) neighborPts.push_back(binIndex);
221  }
222  }
223  }
224 
225  return;
226 }
227 
229  HoughCluster& neighborPts,
230  HoughCluster& houghCluster,
231  RhoThetaAccumulatorBinMap& rhoThetaAccumulatorBinMap,
232  size_t threshold) const
233 {
234  /**
235  * @brief The workhorse routine for a DBScan like clustering routine to identify peak bins in Hough Space
236  */
237 
238  // Start by adding the input point to our Hough Cluster
239  houghCluster.push_back(curBin);
240 
241  for(auto& binIndex : neighborPts)
242  {
243  AccumulatorBin& accBin = rhoThetaAccumulatorBinMap[binIndex];
244 
245  if (!accBin.isVisited())
246  {
247  accBin.setVisited();
248 
249  HoughCluster nextNeighborPts;
250 
251  HoughRegionQuery(binIndex, rhoThetaAccumulatorBinMap, nextNeighborPts, threshold);
252 
253  if(!nextNeighborPts.empty())
254  {
255  for(auto& neighborIdx : nextNeighborPts)
256  {
257  neighborPts.push_back(neighborIdx);
258  }
259  }
260  }
261 
262  if (!accBin.isInCluster())
263  {
264  houghCluster.push_back(binIndex);
265  accBin.setInCluster();
266  }
267  }
268 
269  return;
270 }
271 
273 {
274  return *left < *right;
275 }
276 
278 {
279  bool operator() (const reco::ClusterHit3D* left, const reco::ClusterHit3D* right) const {return Hit3DCompare(left, right);}
280 };
281 
283 {
284 public:
285  OrderHitsAlongWire(int view = 0) : m_view(view) {}
286 
288  {
289  int viewToCheck = (m_view + 1) % 3;
290 
291  return left->getHits()[viewToCheck]->getHit().WireID().Wire < right->getHits()[viewToCheck]->getHit().WireID().Wire;
292  }
293 private:
294  int m_view;
295 };
296 
298 {
299  bool operator()(const std::pair<size_t,size_t>& left, const std::pair<size_t,size_t>& right)
300  {
301  return left.second < right.second;
302  }
303 
304 };
305 
307  reco::HitPairListPtr& outputList) const
308 {
309  // Intention is to try to find the largest "contiguous" block of hits in the input list
310  // In a nutshell, the idea is to order input hits according to the pca, then
311  // loop through the hits and store them in a new hit list. If a gap is detected then
312  // we terminate the previous list, create a new one and continue. After the loop over
313  // hits is complete then simply pick the biggest list.
314  // Step I is to run the pca and order the hits
316 
317  m_pcaAlg.PCAAnalysis_3D(inputHitList, pca, true);
318 
319  // It would seem that the analysis can fail!
320  if (!pca.getSvdOK())
321  {
322  outputList = inputHitList;
323  return;
324  }
325 
326  m_pcaAlg.PCAAnalysis_calc3DDocas(inputHitList, pca);
327 
328  inputHitList.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
329  outputList.clear();
330 
331  // Create containers to hold our hit lists
332  reco::HitPairListPtr continuousHitList;
333 
334  std::map<size_t, reco::HitPairListPtr > gapHitListMap;
335 
336  // Remember the distance arc length of the last hit
337  double arcLenLastHit = inputHitList.front()->getArclenToPoca();
338 
339  // Loop through the input hits
340  for(const auto& hit3D : inputHitList)
341  {
342  // Hits in order, delta arclength should always be positive
343  double arcLen = hit3D->getArclenToPoca();
344  double deltaArcLen = arcLen - arcLenLastHit;
345 
346  if (deltaArcLen > m_maximumGap)
347  {
348  gapHitListMap[continuousHitList.size()] = continuousHitList;
349  continuousHitList.clear();
350  }
351 
352  continuousHitList.emplace_back(hit3D);
353 
354  arcLenLastHit = arcLen;
355  }
356 
357  if (!continuousHitList.empty()) gapHitListMap[continuousHitList.size()] = continuousHitList;
358 
359  // Get the largest list of hits
360  std::map<size_t, reco::HitPairListPtr >::reverse_iterator longestListItr = gapHitListMap.rbegin();
361 
362  if (longestListItr != gapHitListMap.rend())
363  {
364  size_t nContinuousHits = longestListItr->first;
365  reco::HitPairListPtr longestList = longestListItr->second;
366 
367  outputList.resize(nContinuousHits);
368  std::copy(longestList.begin(), longestList.end(), outputList.begin());
369  }
370 
371  return;
372 }
373 
376  int& nLoops,
377  RhoThetaAccumulatorBinMap& rhoThetaAccumulatorBinMap,
378  HoughClusterList& houghClusters) const
379 {
380  // The goal of this function is to do a basic Hough Transform on the input list of 3D hits.
381  // In order to transform this to a 2D problem, the 3D hits are projected to the plane of the two
382  // largest eigen values from the input principal components analysis axis.
383  // There are two basic steps to the job here:
384  // 1) Build and accumlate a rho-theta map
385  // 2) Go through that rho-theta map and find candidate Hough "clusters"
386  // Unfortunately, the following may not be suitable viewing for those who may be feint of heart
387  //
388  // Define some constants
389  static int histCount(0);
390  const double maxTheta(M_PI); // Obviously, 180 degrees
391  const double thetaBinSize(maxTheta/double(m_thetaBins)); // around 4 degrees (45 bins)
392  const double rhoBinSizeMin(m_geometry->WirePitch()); // Wire spacing gives a natural bin size?
393 
394  // Recover the parameters from the Principal Components Analysis that we need to project and accumulate
395  TVector3 pcaCenter(pca.getAvePosition()[0],pca.getAvePosition()[1],pca.getAvePosition()[2]);
396  TVector3 planeVec0(pca.getEigenVectors()[0][0],pca.getEigenVectors()[0][1],pca.getEigenVectors()[0][2]);
397  TVector3 planeVec1(pca.getEigenVectors()[1][0],pca.getEigenVectors()[1][1],pca.getEigenVectors()[1][2]);
398  TVector3 pcaPlaneNrml(pca.getEigenVectors()[2][0],pca.getEigenVectors()[2][1],pca.getEigenVectors()[2][2]);
399  double eigenVal0 = 3. * sqrt(pca.getEigenValues()[0]);
400  double eigenVal1 = 3. * sqrt(pca.getEigenValues()[1]);
401  double maxRho = std::sqrt(eigenVal0*eigenVal0 + eigenVal1*eigenVal1) * 2. / 3.;
402  double rhoBinSize = maxRho / double(m_rhoBins);
403 
404  if (rhoBinSize < rhoBinSizeMin) rhoBinSize = rhoBinSizeMin;
405 
406  // **********************************************************************
407  // Part I: Accumulate values in the rho-theta map
408  // **********************************************************************
409 
410  size_t maxBinCount(0);
411  int nAccepted3DHits(0);
412 
413  // Commence looping over the input 3D hits and fill our accumulator bins
414  for(reco::HitPairListPtr::const_iterator hit3DItr = hitPairListPtr.begin();
415  hit3DItr != hitPairListPtr.end();
416  hit3DItr++)
417  {
418  // Recover the hit
419  const reco::ClusterHit3D* hit3D(*hit3DItr);
420 
421  // Skip hits which are not skeleton points
422  if (!(hit3D->getStatusBits() & 0x10000000)) continue;
423 
424  nAccepted3DHits++;
425 
426  TVector3 hit3DPosition(hit3D->getPosition()[0], hit3D->getPosition()[1], hit3D->getPosition()[2]);
427  TVector3 pcaToHitVec = hit3DPosition - pcaCenter;
428  TVector3 pcaToHitPlaneVec(pcaToHitVec.Dot(planeVec0), pcaToHitVec.Dot(planeVec1), 0.);
429 
430  double xPcaToHit = pcaToHitPlaneVec[0];
431  double yPcaToHit = pcaToHitPlaneVec[1];
432 
433  // Create an accumulator value
434  AccumulatorValues accValue(pcaToHitPlaneVec, hit3DItr);
435 
436  // Commence loop over theta to fill accumulator bins
437  // Note that with theta in the range 0-pi then we can have negative values for rho
438  for(int thetaIdx = 0; thetaIdx < m_thetaBins; thetaIdx++)
439  {
440  // We need to convert our theta index to an angle
441  double theta = thetaBinSize * double(thetaIdx);
442 
443  // calculate rho for this angle
444  double rho = xPcaToHit * std::cos(theta) + yPcaToHit * std::sin(theta);
445 
446  // Get the rho index
447  int rhoIdx = std::round(rho / rhoBinSize);
448 
449  // Accumulate
450  BinIndex binIndex(rhoIdx, thetaIdx);
451 
452  rhoThetaAccumulatorBinMap[binIndex].addAccumulatorValue(accValue);
453 
454  if (rhoThetaAccumulatorBinMap[binIndex].getAccumulatorValues().size() > maxBinCount)
455  maxBinCount = rhoThetaAccumulatorBinMap[binIndex].getAccumulatorValues().size();
456  }
457  }
458 
459  // Accumulation done, if asked now display the hist
460  if (m_displayHist)
461  {
462  std::ostringstream ostr;
463  ostr << "Hough Histogram " << histCount++;
464  m_Canvases.emplace_back(new TCanvas(ostr.str().c_str(), ostr.str().c_str(), 1000, 1000));
465 
466  std::ostringstream ostr2;
467  ostr2 << "Plane";
468 
469  m_Canvases.back()->GetFrame()->SetFillColor(46);
470  m_Canvases.back()->SetFillColor(19);
471  m_Canvases.back()->SetBorderMode(19);
472  m_Canvases.back()->cd(1);
473 
474  double zmin = 0.06;
475  double zmax = 0.94;
476  double xmin = 0.04;
477  double xmax = 0.95;
478  TPad* p = new TPad(ostr2.str().c_str(), ostr2.str().c_str(), zmin, xmin, zmax, xmax);
479  p->SetBit(kCanDelete); // Give away ownership.
480  p->Range(zmin, xmin, zmax, xmax);
481  p->SetFillStyle(4000); // Transparent.
482  p->Draw();
483  m_Pads.push_back(p);
484 
485  TH2D* houghHist = new TH2D("HoughHist", "Hough Space", 2*m_rhoBins, -m_rhoBins+0.5, m_rhoBins+0.5, m_thetaBins, 0., m_thetaBins);
486 
487  for(const auto& rhoThetaMap : rhoThetaAccumulatorBinMap)
488  {
489  houghHist->Fill(rhoThetaMap.first.first, rhoThetaMap.first.second+0.5, rhoThetaMap.second.getAccumulatorValues().size());
490  }
491 
492  houghHist->SetBit(kCanDelete);
493  houghHist->Draw();
494  m_Canvases.back()->Update();
495  }
496 
497  // **********************************************************************
498  // Part II: Use DBScan (or a slight variation) to find clusters of bins
499  // **********************************************************************
500 
501  size_t thresholdLo = std::max(size_t(m_hiThresholdFrac*nAccepted3DHits), m_hiThresholdMin);
502  size_t thresholdHi = m_loThresholdFrac * maxBinCount;
503 
504  std::list<RhoThetaAccumulatorBinMap::iterator> binIndexList;
505 
506  for(RhoThetaAccumulatorBinMap::iterator mapItr = rhoThetaAccumulatorBinMap.begin();
507  mapItr != rhoThetaAccumulatorBinMap.end();
508  mapItr++)
509  binIndexList.push_back(mapItr);
510 
511  binIndexList.sort(SortBinIndexList());
512 
513  for(auto& mapItr : binIndexList)
514  {
515  // If we have been here before we skip
516  //if (mapItr.second.isVisited()) continue;
517  if (mapItr->second.isInCluster()) continue;
518 
519  // Mark this bin as visited
520  // Actually, don't mark it since we are double thresholding and don't want it missed
521  //mapItr.second.setVisited();
522 
523  // Make sure over threshold
524  if (mapItr->second.getAccumulatorValues().size() < thresholdLo)
525  {
526  mapItr->second.setNoise();
527  continue;
528  }
529 
530  // Set the low threshold to make sure we merge bins that might be either side of a boundary trajectory
531  thresholdHi = std::max(size_t(m_loThresholdFrac * mapItr->second.getAccumulatorValues().size()), m_hiThresholdMin);
532 
533  // Recover our neighborhood
534  HoughCluster neighborhood;
535  BinIndex curBin(mapItr->first);
536 
537  HoughRegionQuery(curBin, rhoThetaAccumulatorBinMap, neighborhood, thresholdHi);
538 
539  houghClusters.push_back(HoughCluster());
540 
541  HoughCluster& houghCluster = houghClusters.back();
542 
543  expandHoughCluster(curBin, neighborhood, houghCluster, rhoThetaAccumulatorBinMap, thresholdHi);
544  }
545 
546  // Sort the clusters using the SortHoughClusterList metric
547  if (!houghClusters.empty()) houghClusters.sort(SortHoughClusterList(rhoThetaAccumulatorBinMap));
548 
549  return;
550 }
551 
553 {
554  if (seed3DHits.size() < m_minimum3DHits) return false;
555 
556  reco::PrincipalComponents seedFullPca;
557 
558  m_pcaAlg.PCAAnalysis_3D(seed3DHits, seedFullPca, true);
559 
560  if (!seedFullPca.getSvdOK()) return false;
561 
562  // Use the following to set the 3D doca and arclength for each hit
563  m_pcaAlg.PCAAnalysis_calc3DDocas(seed3DHits, seedFullPca);
564 
565  // Use this info to sort the hits along the principle axis
566  //seed3DHits.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
568 
569  // The idea here is to search for the first hit that lies "close" to the principle axis
570  // At that point we count out n hits to use as the seed
571  reco::HitPairListPtr seedHit3DList;
572  std::set<const reco::ClusterHit2D*> seedHitSet;
573  double aveDocaToAxis = seedFullPca.getAveHitDoca();
574  int gapCount(0);
575 
576  // Now loop through hits to search for a "continuous" block of at least m_numSeed2DHits
577  // We'll arrive at that number by collecting 2D hits in an stl set which will keep track of unique occurances
578  for(reco::HitPairListPtr::iterator peakBinItr = seed3DHits.begin(); peakBinItr != seed3DHits.end(); peakBinItr++)
579  {
580  const reco::ClusterHit3D* hit3D = *peakBinItr;
581 
582  if (hit3D->getDocaToAxis() < m_numAveDocas*aveDocaToAxis)
583  {
584  // Check if we need to reset because of gap count
585  if (gapCount > m_numSkippedHits)
586  {
587  seedHit3DList.clear();
588  seedHitSet.clear();
589  }
590 
591  seedHit3DList.push_back(hit3D);
592 
593  for(const auto& hit : hit3D->getHits()) seedHitSet.insert(hit);
594 
595  gapCount = 0;
596  }
597  else gapCount++;
598 
599  if (seedHitSet.size() > m_numSeed2DHits) break;
600  }
601 
602  // If not enough hits then we are done
603  if (seedHit3DList.size() < m_minimum3DHits) return false;
604 
606 
607  // Use only the "seed" 3D hits to get new PCA axes
608  m_pcaAlg.PCAAnalysis_3D(seedHit3DList, seedPca, true);
609 
610  if (!seedPca.getSvdOK()) return false;
611 
612  m_pcaAlg.PCAAnalysis_calc3DDocas(seedHit3DList, seedPca);
613  //seedHit3DList.sort(SeedFinderAlgBase::Sort3DHitsByAbsArcLen3D());
614  seedHit3DList.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
615 
616  // Now translate the seedCenter by the arc len to the first hit
617  double seedCenter[3] = {seedPca.getAvePosition()[0], seedPca.getAvePosition()[1], seedPca.getAvePosition()[2]};
618  double seedDir[3] = {seedPca.getEigenVectors()[0][0], seedPca.getEigenVectors()[0][1], seedPca.getEigenVectors()[0][2]};
619 
620  double arcLen = seedHit3DList.front()->getArclenToPoca();
621  double seedStart[3] = {seedCenter[0]+arcLen*seedDir[0], seedCenter[1]+arcLen*seedDir[1], seedCenter[2]+arcLen*seedDir[2]};
622 
623  //seedStart[0] = seedHit3DList.front()->getX();
624  //seedStart[1] = seedHit3DList.front()->getY();
625  //seedStart[2] = seedHit3DList.front()->getZ();
626 
627  if (seedHitSet.size() >= 10)
628  {
629  TVector3 newSeedPos;
630  TVector3 newSeedDir;
631  double chiDOF;
632 
633  LineFit2DHits(seedHitSet, seedStart[0], newSeedPos, newSeedDir, chiDOF);
634 
635  if (chiDOF > 0.)
636  {
637  // check angles between new/old directions
638  double cosAng = seedDir[0]*newSeedDir[0]+seedDir[1]*newSeedDir[1]+seedDir[2]*newSeedDir[2];
639 
640  if (cosAng < 0.) newSeedDir *= -1.;
641 
642  seedStart[0] = newSeedPos[0];
643  seedStart[1] = newSeedPos[1];
644  seedStart[2] = newSeedPos[2];
645  seedDir[0] = newSeedDir[0];
646  seedDir[1] = newSeedDir[1];
647  seedDir[2] = newSeedDir[2];
648  }
649  }
650 
651  // Keep track of this seed and the 3D hits that make it up
652  seedHitPair = SeedHitPairListPair(recob::Seed(seedStart, seedDir), seedHit3DList);
653 
654  // We are going to drop a few hits off the ends in the hope this facilitates finding more track like objects, provided there are enough hits
655  if (seed3DHits.size() > 100)
656  {
657  // Need to reset the doca/arclen first
658  m_pcaAlg.PCAAnalysis_calc3DDocas(seed3DHits, seedFullPca);
659 
660  // Now resort the hits
662 
663  size_t numToKeep = seed3DHits.size() - 10;
664 
665  reco::HitPairListPtr::iterator endItr = seed3DHits.begin();
666 
667  std::advance(endItr, numToKeep);
668 
669  seed3DHits.erase(endItr, seed3DHits.end());
670  }
671 
672  return true;
673 }
674 
675 
677  reco::PrincipalComponents& inputPCA,
678  SeedHitPairListPairVec& seedHitPairVec) const
679 {
680  // This will be a busy routine... the basic tasks are:
681  // 1) loop through hits and project to the plane defined by the two largest eigen values, accumulate in Hough space
682  // 2) "Cluster" the Hough space to associate hits which are common to a line
683  // 3) Process these clusters (still to be defined exactly)
684 
685  // Create an interim data structure which will allow us to sort our seeds by "best"
686  // before we return them in the seedHitPairVec
687  typedef std::map<size_t, SeedHitPairListPairVec > SizeToSeedToHitMap;
688 
689  SizeToSeedToHitMap seedHitPairMap;
690  SeedHitPairListPair seedHitPair;
691 
692  // Make sure we are using the right pca
693  reco::HitPairListPtr hitPairListPtr = inputHitPairListPtr;
694 
695  int nLoops(0);
696 
697  // Make a local copy of the input PCA
698  reco::PrincipalComponents pca = inputPCA;
699 
700  // We loop over hits in our list until there are no more
701  while(!hitPairListPtr.empty())
702  {
703  // We also require that there be some spread in the data, otherwise not worth running?
704  double eigenVal0 = 3. * sqrt(pca.getEigenValues()[0]);
705  double eigenVal1 = 3. * sqrt(pca.getEigenValues()[1]);
706 
707  if (eigenVal0 > 5. && eigenVal1 > 0.001)
708  {
709  // **********************************************************************
710  // Part I: Build Hough space and find Hough clusters
711  // **********************************************************************
712  RhoThetaAccumulatorBinMap rhoThetaAccumulatorBinMap;
713  HoughClusterList houghClusters;
714 
715  findHoughClusters(hitPairListPtr, pca, nLoops, rhoThetaAccumulatorBinMap, houghClusters);
716 
717  // If no clusters then done
718  if (houghClusters.empty()) break;
719 
720  // **********************************************************************
721  // Part II: Go through the clusters to find the peak bins
722  // **********************************************************************
723 
724  // We need to use a set so we can be sure to have unique hits
725  reco::HitPairListPtr clusterHitsList;
726  std::set<const reco::ClusterHit3D*> masterHitPtrList;
727  std::set<const reco::ClusterHit3D*> peakBinPtrList;
728 
729  size_t firstPeakCount(0);
730 
731  // Loop through the list of all clusters found above
732  for(auto& houghCluster : houghClusters)
733  {
734  BinIndex peakBin = houghCluster.front();
735  size_t peakCount = 0;
736  size_t totalHits = 0;
737 
738  // Make a local (to this cluster) set of of hits
739  std::set<const reco::ClusterHit3D*> localHitPtrList;
740 
741  // Now loop through the bins that were attached to this cluster
742  for(auto& binIndex : houghCluster)
743  {
744  // An even more local list so we can keep track of peak values
745  std::set<const reco::ClusterHit3D*> tempHitPtrList;
746 
747  // Recover the hits associated to this cluster
748  for(auto& hitItr : rhoThetaAccumulatorBinMap[binIndex].getAccumulatorValues())
749  {
750  reco::HitPairListPtr::const_iterator hit3DItr = hitItr.getHitIterator();
751 
752  tempHitPtrList.insert(*hit3DItr);
753  }
754 
755  // count hits before we remove any
756  totalHits += tempHitPtrList.size();
757 
758  // Trim out any hits already used by a bigger/better cluster
759  std::set<const reco::ClusterHit3D*> tempHit3DSet;
760 
761  std::set_difference(tempHitPtrList.begin(), tempHitPtrList.end(),
762  masterHitPtrList.begin(), masterHitPtrList.end(),
763  std::inserter(tempHit3DSet, tempHit3DSet.end()) );
764 
765  tempHitPtrList = tempHit3DSet;
766 
767  size_t binCount = tempHitPtrList.size();
768 
769  if (peakCount < binCount)
770  {
771  peakCount = binCount;
772  peakBin = binIndex;
773  peakBinPtrList = tempHitPtrList;
774  }
775 
776  // Add this to our local list
777  localHitPtrList.insert(tempHitPtrList.begin(),tempHitPtrList.end());
778  }
779 
780  if (localHitPtrList.size() < m_minimum3DHits) continue;
781 
782  if (!firstPeakCount) firstPeakCount = peakCount;
783 
784  // If the peak counts are significantly less than the first cluster's peak then skip
785  if (peakCount < firstPeakCount / 10) continue;
786 
787  // **********************************************************************
788  // Part III: Make a Seed from the peak bin hits
789  // **********************************************************************
790 
791  reco::HitPairListPtr allPeakBinHits;
792 
793  for(const auto& hit3D : localHitPtrList) allPeakBinHits.push_back(hit3D);
794 
795  reco::HitPairListPtr peakBinHits;
796 
797  // Find longest "continuous" set of hits and use these for the seed
798  findHitGaps(allPeakBinHits, peakBinHits);
799 
800  if (peakBinHits.size() < m_minimum3DHits) continue;
801 
802  // We now build the actual seed.
803  if (buildSeed(peakBinHits, seedHitPair))
804  {
805  // Keep track of this in our map (which will do ordering for us)
806  seedHitPairMap[peakBinHits.size()].push_back(seedHitPair);
807 
808  // For visual testing in event display, mark all the hits in the first seed so we can see them
809  if (seedHitPairMap.size() == 1)
810  {
811  for(const auto& hit3D : peakBinHits) hit3D->setStatusBit(0x40000000);
812  }
813 
814 // for(const auto& hit3D : seedHitPair.second) hit3D->setStatusBit(0x40000000);
815  }
816 
817  // Our peakBinHits collection will most likely be a subset of the localHitPtrList collection
818  // We want to remove only the "pure" hits which are those in the peakBinHits collection
819  // So, sort them and then add to our master list
820  peakBinHits.sort();
821 
822  masterHitPtrList.insert(peakBinHits.begin(),peakBinHits.end());
823 
824  if (hitPairListPtr.size() - masterHitPtrList.size() < m_minimum3DHits) break;
825  } // end loop over hough clusters
826 
827  // If the masterHitPtrList is empty then nothing happened and we're done
828  if (masterHitPtrList.empty()) break;
829 
830  // **********************************************************************
831  // Part IV: Remove remaining peak bin hits from HitPairPtrList
832  // **********************************************************************
833 
834  hitPairListPtr.sort();
835 
836  reco::HitPairListPtr::iterator newListEnd =
837  std::set_difference(hitPairListPtr.begin(), hitPairListPtr.end(),
838  masterHitPtrList.begin(), masterHitPtrList.end(),
839  hitPairListPtr.begin() );
840 
841  hitPairListPtr.erase(newListEnd, hitPairListPtr.end());
842 
843  if (hitPairListPtr.size() < m_minimum3DHits) break;
844 
845  if (nLoops++ > m_maxLoopsPerCluster) break;
846 
847  // ********************************************************
848  }
849  else break; // eigen values not in range
850 
851  // At this point run the PCA on the remaining hits
852  m_pcaAlg.PCAAnalysis_3D(hitPairListPtr, pca, true);
853 
854  if (!pca.getSvdOK()) break;
855  }
856 
857  // The final task before returning is to transfer the stored seeds into the output seed vector
858  // What we want to do is make sure the first seeds are the "best" seeds which is defined as the
859  // seeds which were associated to the most hits by the Hough Transform. Our seed map will have
860  // the reverse of this ordering so we simply iterate through it "backwards"
861  for(SizeToSeedToHitMap::reverse_iterator seedMapItr = seedHitPairMap.rbegin(); seedMapItr != seedHitPairMap.rend(); seedMapItr++)
862  {
863  for(const auto& seedHitPair : seedMapItr->second)
864  {
865  seedHitPairVec.emplace_back(seedHitPair);
866  }
867  }
868 
869  return true;
870 }
871 
873  reco::PrincipalComponents& inputPCA,
874  reco::HitPairListPtrList& hitPairListPtrList) const
875 {
876  // The goal of this routine is run the Hough Transform on the input set of hits
877  // and then to return a list of lists of hits which are associated to a given line
878 
879  // Make sure we are using the right pca
880  reco::HitPairListPtr hitPairListPtr = inputHitPairListPtr;
881 
882  int nLoops(0);
883 
884  // Make a local copy of the input PCA
885  reco::PrincipalComponents pca = inputPCA;
886 
887  // We also require that there be some spread in the data, otherwise not worth running?
888  double eigenVal0 = 3. * sqrt(pca.getEigenValues()[0]);
889  double eigenVal1 = 3. * sqrt(pca.getEigenValues()[1]);
890 
891  if (eigenVal0 > 5. && eigenVal1 > 0.001)
892  {
893  // **********************************************************************
894  // Part I: Build Hough space and find Hough clusters
895  // **********************************************************************
896  RhoThetaAccumulatorBinMap rhoThetaAccumulatorBinMap;
897  HoughClusterList houghClusters;
898 
899  findHoughClusters(hitPairListPtr, pca, nLoops, rhoThetaAccumulatorBinMap, houghClusters);
900 
901  // **********************************************************************
902  // Part II: Go through the clusters to find the peak bins
903  // **********************************************************************
904 
905  // We need to use a set so we can be sure to have unique hits
906  reco::HitPairListPtr clusterHitsList;
907  std::set<const reco::ClusterHit3D*> masterHitPtrList;
908  std::set<const reco::ClusterHit3D*> peakBinPtrList;
909 
910  size_t firstPeakCount(0);
911 
912  // Loop through the list of all clusters found above
913  for(auto& houghCluster : houghClusters)
914  {
915  BinIndex peakBin = houghCluster.front();
916  size_t peakCount = 0;
917  size_t totalHits = 0;
918 
919  // Make a local (to this cluster) set of of hits
920  std::set<const reco::ClusterHit3D*> localHitPtrList;
921 
922  // Now loop through the bins that were attached to this cluster
923  for(auto& binIndex : houghCluster)
924  {
925  // An even more local list so we can keep track of peak values
926  std::set<const reco::ClusterHit3D*> tempHitPtrList;
927 
928  // Recover the hits associated to this cluster
929  for(auto& hitItr : rhoThetaAccumulatorBinMap[binIndex].getAccumulatorValues())
930  {
931  reco::HitPairListPtr::const_iterator hit3DItr = hitItr.getHitIterator();
932 
933  tempHitPtrList.insert(*hit3DItr);
934  }
935 
936  // count hits before we remove any
937  totalHits += tempHitPtrList.size();
938 
939  // Trim out any hits already used by a bigger/better cluster
940  std::set<const reco::ClusterHit3D*> tempHit3DSet;
941 
942  std::set_difference(tempHitPtrList.begin(), tempHitPtrList.end(),
943  masterHitPtrList.begin(), masterHitPtrList.end(),
944  std::inserter(tempHit3DSet, tempHit3DSet.end()) );
945 
946  tempHitPtrList = tempHit3DSet;
947 
948  size_t binCount = tempHitPtrList.size();
949 
950  if (peakCount < binCount)
951  {
952  peakCount = binCount;
953  peakBin = binIndex;
954  peakBinPtrList = tempHitPtrList;
955  }
956 
957  // Add this to our local list
958  localHitPtrList.insert(tempHitPtrList.begin(),tempHitPtrList.end());
959  }
960 
961  if (localHitPtrList.size() < m_minimum3DHits) continue;
962 
963  if (!firstPeakCount) firstPeakCount = peakCount;
964 
965  // If the peak counts are significantly less than the first cluster's peak then skip
966  if (peakCount < firstPeakCount / 10) continue;
967 
968  // **********************************************************************
969  // Part III: Make a list of hits from the total number associated
970  // **********************************************************************
971 
972  hitPairListPtrList.push_back(reco::HitPairListPtr());
973 
974  hitPairListPtrList.back().resize(localHitPtrList.size());
975  std::copy(localHitPtrList.begin(), localHitPtrList.end(), hitPairListPtrList.back().begin());
976 
977  // We want to remove the hits which have been used from further contention
978  masterHitPtrList.insert(localHitPtrList.begin(),localHitPtrList.end());
979 
980  if (hitPairListPtr.size() - masterHitPtrList.size() < m_minimum3DHits) break;
981  } // end loop over hough clusters
982  }
983 
984  return true;
985 }
986 
987 
988 //------------------------------------------------------------------------------
989 void HoughSeedFinderAlg::LineFit2DHits(std::set<const reco::ClusterHit2D*>& hit2DSet,
990  double XOrigin,
991  TVector3& Pos,
992  TVector3& Dir,
993  double& ChiDOF) const
994 {
995  // The following is lifted from Bruce Baller to try to get better
996  // initial parameters for a candidate Seed. It is slightly reworked
997  // which is why it is included here instead of used as is.
998  //
999  // Linear fit using X as the independent variable. Hits to be fitted
1000  // are passed in the hits vector in a pair form (X, WireID). The
1001  // fitted track position at XOrigin is returned in the Pos vector.
1002  // The direction cosines are returned in the Dir vector.
1003  //
1004  // SVD fit adapted from $ROOTSYS/tutorials/matrix/solveLinear.C
1005  // Fit equation is w = A(X)v, where w is a vector of hit wires, A is
1006  // a matrix to calculate a track projected to a point at X, and v is
1007  // a vector (Yo, Zo, dY/dX, dZ/dX).
1008  //
1009  // Note: The covariance matrix should also be returned
1010  // B. Baller August 2014
1011 
1012  // assume failure
1013  ChiDOF = -1;
1014 
1015  if(hit2DSet.size() < 4) return;
1016 
1017  const unsigned int nvars = 4;
1018  unsigned int npts = hit2DSet.size();
1019 
1020  TMatrixD A(npts, nvars);
1021  // vector holding the Wire number
1022  TVectorD w(npts);
1023  unsigned short ninpl[3] = {0};
1024  unsigned short nok = 0;
1025  unsigned short iht(0), cstat, tpc, ipl;
1026  double x, cw, sw, off;
1027 
1028  // Loop over unique 2D hits from the input list of 3D hits
1029  for (const auto& hit : hit2DSet)
1030  {
1031  geo::WireID wireID = hit->getHit().WireID();
1032 
1033  cstat = wireID.Cryostat;
1034  tpc = wireID.TPC;
1035  ipl = wireID.Plane;
1036 
1037  // get the wire plane offset
1038  off = m_geometry->WireCoordinate(0, 0, ipl, tpc, cstat);
1039 
1040  // get the "cosine-like" component
1041  cw = m_geometry->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
1042 
1043  // the "sine-like" component
1044  sw = m_geometry->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
1045 
1046  x = hit->getXPosition() - XOrigin;
1047 
1048  A[iht][0] = cw;
1049  A[iht][1] = sw;
1050  A[iht][2] = cw * x;
1051  A[iht][3] = sw * x;
1052  w[iht] = wireID.Wire - off;
1053 
1054  ++ninpl[ipl];
1055 
1056  // need at least two points in a plane
1057  if(ninpl[ipl] == 2) ++nok;
1058 
1059  iht++;
1060  }
1061 
1062  // need at least 2 planes with at least two points
1063  if(nok < 2) return;
1064 
1065  TDecompSVD svd(A);
1066  bool ok;
1067  TVectorD tVec = svd.Solve(w, ok);
1068 
1069  ChiDOF = 0;
1070 
1071  // not enough points to calculate Chisq
1072  if(npts <= 4) return;
1073 
1074  double ypr, zpr, diff;
1075 
1076  for (const auto& hit : hit2DSet)
1077  {
1078  geo::WireID wireID = hit->getHit().WireID();
1079 
1080  cstat = wireID.Cryostat;
1081  tpc = wireID.TPC;
1082  ipl = wireID.Plane;
1083  off = m_geometry->WireCoordinate(0, 0, ipl, tpc, cstat);
1084  cw = m_geometry->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
1085  sw = m_geometry->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
1086  x = hit->getXPosition() - XOrigin;
1087  ypr = tVec[0] + tVec[2] * x;
1088  zpr = tVec[1] + tVec[3] * x;
1089  diff = ypr * cw + zpr * sw - (wireID.Wire - off);
1090  ChiDOF += diff * diff;
1091  }
1092 
1093 
1094  float werr2 = m_geometry->WirePitch() * m_geometry->WirePitch();
1095  ChiDOF /= werr2;
1096  ChiDOF /= (float)(npts - 4);
1097 
1098  double norm = sqrt(1 + tVec[2] * tVec[2] + tVec[3] * tVec[3]);
1099  Dir[0] = 1 / norm;
1100  Dir[1] = tVec[2] / norm;
1101  Dir[2] = tVec[3] / norm;
1102 
1103  Pos[0] = XOrigin;
1104  Pos[1] = tVec[0];
1105  Pos[2] = tVec[1];
1106 
1107 } // TrkLineFit()
1108 
1109 
1110 } // namespace lar_cluster3d
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
void LineFit2DHits(std::set< const reco::ClusterHit2D * > &hitList, double XOrigin, TVector3 &Pos, TVector3 &Dir, double &ChiDOF) const
intermediate_table::iterator iterator
Define a comparator which will sort hits by arc length along a PCA axis.
void PCAAnalysis_calc3DDocas(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca) const
void findHitGaps(reco::HitPairListPtr &inputHitList, reco::HitPairListPtr &outputList) const
Using Principal Components Axis, look for gaps in a list of 3D hits.
bool getSvdOK() const
Definition: Cluster3D.h:244
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:158
void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
virtual void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
std::vector< AccumulatorValues > AccumulatorValuesVec
HoughSeedFinderAlg(fhicl::ParameterSet const &pset)
Constructor.
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
intermediate_table::const_iterator const_iterator
std::list< HitPairListPtr > HitPairListPtrList
Definition: Cluster3D.h:337
std::vector< TVirtualPad * > m_Pads
View pads in current canvas.
void findHoughClusters(const reco::HitPairListPtr &inputHits, reco::PrincipalComponents &pca, int &nLoops, RhoThetaAccumulatorBinMap &rhoThetaMap, HoughClusterList &clusterList) const
HoughSeedFinderAlg::RhoThetaAccumulatorBinMap & m_accMap
unsigned int getStatusBits() const
Definition: Cluster3D.h:157
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
std::pair< recob::Seed, reco::HitPairListPtr > SeedHitPairListPair
float getDocaToAxis() const
Definition: Cluster3D.h:169
art framework interface to geometry description
Define a comparator which will sort hits by the absolute value of arc length so hits are ordered clos...
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
std::vector< std::unique_ptr< TCanvas > > m_Canvases
Graphical trace canvases.
void expandHoughCluster(BinIndex &curBin, HoughCluster &neighborPts, HoughCluster &houghCluster, RhoThetaAccumulatorBinMap &rhoThetaAccumulatorBinMap, size_t threshold) const
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:246
bool operator()(const std::pair< size_t, size_t > &left, const std::pair< size_t, size_t > &right)
void HoughRegionQuery(BinIndex &curBin, RhoThetaAccumulatorBinMap &rhoThetaAccumulatorBinMap, HoughCluster &neighborPts, size_t threshold) const
T get(std::string const &key) const
Definition: ParameterSet.h:271
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:248
std::list< HoughCluster > HoughClusterList
p
Definition: test.py:223
const AccumulatorValuesVec & getAccumulatorValues() const
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
bool operator()(const HoughSeedFinderAlg::HoughCluster &left, const HoughSeedFinderAlg::HoughCluster &right)
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right)
Q_UINT16 values[128]
static int max(int a, int b)
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
bool buildSeed(reco::HitPairListPtr &seed3DHits, SeedHitPairListPair &seedHitPair) const
Given a list of candidate "seed" 3D hits, build the seed and get associated unique 2D hits...
bool operator()(const HoughSeedFinderAlg::RhoThetaAccumulatorBinMap::iterator &left, const HoughSeedFinderAlg::RhoThetaAccumulatorBinMap::iterator &right)
This is used to sort bins in Hough Space.
std::map< BinIndex, AccumulatorBin > RhoThetaAccumulatorBinMap
auto norm(Vector const &v)
Return norm of the specified vector.
Detector simulation of raw signals on wires.
Encapsulate the geometry of a wire.
reco::HitPairListPtr::const_iterator m_hit3DIterator
This will be used to take us back to our 3D hit.
#define M_PI
Definition: includeROOT.h:54
AccumulatorValues(const TVector3 &position, const reco::HitPairListPtr::const_iterator &itr)
virtual bool findTrackHits(reco::HitPairListPtr &hitPairListPtr, reco::PrincipalComponents &inputPCA, reco::HitPairListPtrList &hitPairListPtrList) const
Given the list of hits this will return the sets of hits which belong on the same line...
Declaration of signal hit object.
Encapsulate the construction of a single detector plane.
const float getAveHitDoca() const
Definition: Cluster3D.h:249
virtual ~HoughSeedFinderAlg()
Destructor.
TVector3 m_position
We really only need the x,y coordinates here but keep all three for now.
T copy(T const &v)
#define A
Definition: memgrp.cpp:38
const TVector3 & getPosition() const
list x
Definition: train.py:276
virtual bool findTrackSeeds(reco::HitPairListPtr &hitPairListPtr, reco::PrincipalComponents &inputPCA, SeedHitPairListPairVec &seedHitPairVec) const
Given the list of hits this will search for candidate Seed objects and return them.
SortHoughClusterList(HoughSeedFinderAlg::RhoThetaAccumulatorBinMap &accMap)
This is used to sort "Hough Clusters" by the maximum entries in a bin.
std::vector< SeedHitPairListPair > SeedHitPairListPairVec
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
AccumulatorValues()
A utility class to contain the values of a given "bin" in Hough Space.
AccumulatorBin()
A utility class used to accumulate the above values.
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:171
bool Hit3DCompare(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right)
reco::HitPairListPtr::const_iterator getHitIterator() const
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:247