Classes | Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
lar_cluster3d::HoughSeedFinderAlg Class Reference

HoughSeedFinderAlg class. More...

#include <HoughSeedFinderAlg.h>

Inheritance diagram for lar_cluster3d::HoughSeedFinderAlg:
lar_cluster3d::SeedFinderAlgBase lar_cluster3d::SeedFinderAlgBase

Classes

class  AccumulatorBin
 
struct  SortBinIndexList
 
class  SortHoughClusterList
 

Public Member Functions

 HoughSeedFinderAlg (fhicl::ParameterSet const &pset)
 Constructor. More...
 
virtual ~HoughSeedFinderAlg ()
 Destructor. More...
 
virtual void reconfigure (fhicl::ParameterSet const &pset)
 a handler for the case where the algorithm control parameters are to be reset More...
 
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. More...
 
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. More...
 
 HoughSeedFinderAlg (fhicl::ParameterSet const &pset)
 Constructor. More...
 
bool findTrackSeeds (reco::HitPairListPtr &hitPairListPtr, reco::PrincipalComponents &inputPCA, SeedHitPairListPairVec &seedHitPairVec) const override
 Given the list of hits this will search for candidate Seed objects and return them. More...
 
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. More...
 

Private Types

typedef std::pair< int, int > BinIndex
 
typedef std::map< BinIndex, AccumulatorBinRhoThetaAccumulatorBinMap
 
typedef std::list< BinIndexHoughCluster
 
typedef std::list< HoughClusterHoughClusterList
 
typedef std::pair< int, int > BinIndex
 
typedef std::map< BinIndex, AccumulatorBinRhoThetaAccumulatorBinMap
 
typedef std::list< BinIndexHoughCluster
 
typedef std::list< HoughClusterHoughClusterList
 

Private Member Functions

void findHitGaps (reco::HitPairListPtr &inputHitList, reco::HitPairListPtr &outputList) const
 Using Principal Components Axis, look for gaps in a list of 3D hits. More...
 
void HoughRegionQuery (BinIndex &curBin, RhoThetaAccumulatorBinMap &rhoThetaAccumulatorBinMap, HoughCluster &neighborPts, size_t threshold) const
 
void expandHoughCluster (BinIndex &curBin, HoughCluster &neighborPts, HoughCluster &houghCluster, RhoThetaAccumulatorBinMap &rhoThetaAccumulatorBinMap, size_t threshold) const
 
void findHoughClusters (const reco::HitPairListPtr &inputHits, reco::PrincipalComponents &pca, int &nLoops, RhoThetaAccumulatorBinMap &rhoThetaMap, HoughClusterList &clusterList) const
 
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. More...
 
void LineFit2DHits (std::set< const reco::ClusterHit2D * > &hitList, double XOrigin, TVector3 &Pos, TVector3 &Dir, double &ChiDOF) const
 
void findHitGaps (reco::HitPairListPtr &inputHitList, reco::HitPairListPtr &outputList) const
 Using Principal Components Axis, look for gaps in a list of 3D hits. More...
 
void HoughRegionQuery (BinIndex &curBin, RhoThetaAccumulatorBinMap &rhoThetaAccumulatorBinMap, HoughCluster &neighborPts, size_t threshold) const
 
void expandHoughCluster (BinIndex &curBin, HoughCluster &neighborPts, HoughCluster &houghCluster, RhoThetaAccumulatorBinMap &rhoThetaAccumulatorBinMap, size_t threshold) const
 
void findHoughClusters (const reco::HitPairListPtr &inputHits, reco::PrincipalComponents &pca, int &nLoops, RhoThetaAccumulatorBinMap &rhoThetaMap, HoughClusterList &clusterList) const
 
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. More...
 
void LineFit2DHits (std::set< const reco::ClusterHit2D * > &hitList, double XOrigin, TVector3 &Pos, TVector3 &Dir, double &ChiDOF) const
 

Private Attributes

size_t m_minimum3DHits
 
int m_thetaBins
 
int m_rhoBins
 
size_t m_hiThresholdMin
 
double m_hiThresholdFrac
 
double m_loThresholdFrac
 
size_t m_numSeed2DHits
 
double m_numAveDocas
 
int m_numSkippedHits
 
int m_maxLoopsPerCluster
 
double m_maximumGap
 
geo::Geometrym_geometry
 
PrincipalComponentsAlg m_pcaAlg
 
bool m_displayHist
 
std::vector< std::unique_ptr< TCanvas > > m_Canvases
 Graphical trace canvases. More...
 
std::vector< TVirtualPad * > m_Pads
 View pads in current canvas. More...
 
geo::Geometry const * m_geometry
 

Detailed Description

HoughSeedFinderAlg class.

Definition at line 31 of file HoughSeedFinderAlg.h.

Member Typedef Documentation

typedef std::pair<int, int> lar_cluster3d::HoughSeedFinderAlg::BinIndex
private

Definition at line 73 of file HoughSeedFinderAlg.h.

typedef std::pair<int, int> lar_cluster3d::HoughSeedFinderAlg::BinIndex
private

Definition at line 82 of file HoughSeedFinderAlg.h.

Definition at line 79 of file HoughSeedFinderAlg.h.

Definition at line 88 of file HoughSeedFinderAlg.h.

Definition at line 80 of file HoughSeedFinderAlg.h.

Definition at line 89 of file HoughSeedFinderAlg.h.

Definition at line 78 of file HoughSeedFinderAlg.h.

Definition at line 87 of file HoughSeedFinderAlg.h.

Constructor & Destructor Documentation

lar_cluster3d::HoughSeedFinderAlg::HoughSeedFinderAlg ( fhicl::ParameterSet const &  pset)

Constructor.

Parameters
pset

Definition at line 49 of file HoughSeedFinderAlg.cxx.

49  :
50  m_minimum3DHits(5),
51  m_thetaBins(360),
52  m_rhoBins(75),
54  m_hiThresholdFrac(.05),
55  m_loThresholdFrac(0.85),
56  m_numSeed2DHits(80),
57  m_numAveDocas(6.),
58  m_numSkippedHits(10),
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 }
virtual void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
lar_cluster3d::HoughSeedFinderAlg::~HoughSeedFinderAlg ( )
virtual

Destructor.

Definition at line 74 of file HoughSeedFinderAlg.cxx.

75 {
76 }
lar_cluster3d::HoughSeedFinderAlg::HoughSeedFinderAlg ( fhicl::ParameterSet const &  pset)

Constructor.

Parameters
pset

Member Function Documentation

bool lar_cluster3d::HoughSeedFinderAlg::buildSeed ( reco::HitPairListPtr seed3DHits,
SeedHitPairListPair seedHitPair 
) const
private

Given a list of candidate "seed" 3D hits, build the seed and get associated unique 2D hits.

bool lar_cluster3d::HoughSeedFinderAlg::buildSeed ( reco::HitPairListPtr seed3DHits,
SeedHitPairListPair seedHitPair 
) const
private

Given a list of candidate "seed" 3D hits, build the seed and get associated unique 2D hits.

Definition at line 552 of file HoughSeedFinderAlg.cxx.

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());
567  seed3DHits.sort(SeedFinderAlgBase::Sort3DHitsByAbsArcLen3D());
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
661  seed3DHits.sort(SeedFinderAlgBase::Sort3DHitsByAbsArcLen3D());
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 }
void LineFit2DHits(std::set< const reco::ClusterHit2D * > &hitList, double XOrigin, TVector3 &Pos, TVector3 &Dir, double &ChiDOF) const
intermediate_table::iterator iterator
void PCAAnalysis_calc3DDocas(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca) const
bool getSvdOK() const
Definition: Cluster3D.h:244
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
std::pair< recob::Seed, reco::HitPairListPtr > SeedHitPairListPair
float getDocaToAxis() const
Definition: Cluster3D.h:169
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:248
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
Detector simulation of raw signals on wires.
const float getAveHitDoca() const
Definition: Cluster3D.h:249
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:171
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:247
void lar_cluster3d::HoughSeedFinderAlg::expandHoughCluster ( BinIndex curBin,
HoughCluster neighborPts,
HoughCluster houghCluster,
RhoThetaAccumulatorBinMap rhoThetaAccumulatorBinMap,
size_t  threshold 
) const
private
void lar_cluster3d::HoughSeedFinderAlg::expandHoughCluster ( BinIndex curBin,
HoughCluster neighborPts,
HoughCluster houghCluster,
RhoThetaAccumulatorBinMap rhoThetaAccumulatorBinMap,
size_t  threshold 
) const
private

The workhorse routine for a DBScan like clustering routine to identify peak bins in Hough Space

The workhorse routine for a DBScan like clustering routine to identify peak bins in Hough Space

Definition at line 228 of file HoughSeedFinderAlg.cxx.

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 }
void HoughRegionQuery(BinIndex &curBin, RhoThetaAccumulatorBinMap &rhoThetaAccumulatorBinMap, HoughCluster &neighborPts, size_t threshold) const
void lar_cluster3d::HoughSeedFinderAlg::findHitGaps ( reco::HitPairListPtr inputHitList,
reco::HitPairListPtr outputList 
) const
private

Using Principal Components Axis, look for gaps in a list of 3D hits.

Parameters
inputHitList- input list of 3D hits to check
pca- Principal Components Axis to use
hitListList- output list of hit lists which are gap free
void lar_cluster3d::HoughSeedFinderAlg::findHitGaps ( reco::HitPairListPtr inputHitList,
reco::HitPairListPtr outputList 
) const
private

Using Principal Components Axis, look for gaps in a list of 3D hits.

Parameters
inputHitList- input list of 3D hits to check
pca- Principal Components Axis to use
hitListList- output list of hit lists which are gap free

Definition at line 306 of file HoughSeedFinderAlg.cxx.

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 }
void PCAAnalysis_calc3DDocas(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca) const
bool getSvdOK() const
Definition: Cluster3D.h:244
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
T copy(T const &v)
void lar_cluster3d::HoughSeedFinderAlg::findHoughClusters ( const reco::HitPairListPtr inputHits,
reco::PrincipalComponents pca,
int &  nLoops,
RhoThetaAccumulatorBinMap rhoThetaMap,
HoughClusterList clusterList 
) const
private
void lar_cluster3d::HoughSeedFinderAlg::findHoughClusters ( const reco::HitPairListPtr inputHits,
reco::PrincipalComponents pca,
int &  nLoops,
RhoThetaAccumulatorBinMap rhoThetaMap,
HoughClusterList clusterList 
) const
private

Definition at line 374 of file HoughSeedFinderAlg.cxx.

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 }
intermediate_table::iterator iterator
intermediate_table::const_iterator const_iterator
std::vector< TVirtualPad * > m_Pads
View pads in current canvas.
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
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
void HoughRegionQuery(BinIndex &curBin, RhoThetaAccumulatorBinMap &rhoThetaAccumulatorBinMap, HoughCluster &neighborPts, size_t threshold) const
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:248
p
Definition: test.py:223
static int max(int a, int b)
#define M_PI
Definition: includeROOT.h:54
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:247
bool lar_cluster3d::HoughSeedFinderAlg::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.

bool lar_cluster3d::HoughSeedFinderAlg::findTrackHits ( reco::HitPairListPtr hitPairListPtr,
reco::PrincipalComponents inputPCA,
reco::HitPairListPtrList hitPairListPtrList 
) const
virtual

Given the list of hits this will return the sets of hits which belong on the same line.

Definition at line 872 of file HoughSeedFinderAlg.cxx.

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 }
intermediate_table::const_iterator const_iterator
void findHoughClusters(const reco::HitPairListPtr &inputHits, reco::PrincipalComponents &pca, int &nLoops, RhoThetaAccumulatorBinMap &rhoThetaMap, HoughClusterList &clusterList) const
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:246
std::list< HoughCluster > HoughClusterList
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
std::map< BinIndex, AccumulatorBin > RhoThetaAccumulatorBinMap
T copy(T const &v)
bool lar_cluster3d::HoughSeedFinderAlg::findTrackSeeds ( reco::HitPairListPtr hitPairListPtr,
reco::PrincipalComponents inputPCA,
SeedHitPairListPairVec seedHitPairVec 
) const
overridevirtual

Given the list of hits this will search for candidate Seed objects and return them.

Implements lar_cluster3d::SeedFinderAlgBase.

bool lar_cluster3d::HoughSeedFinderAlg::findTrackSeeds ( reco::HitPairListPtr hitPairListPtr,
reco::PrincipalComponents inputPCA,
SeedHitPairListPairVec seedHitPairVec 
) const
virtual

Given the list of hits this will search for candidate Seed objects and return them.

Implements lar_cluster3d::SeedFinderAlgBase.

Definition at line 676 of file HoughSeedFinderAlg.cxx.

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 }
intermediate_table::iterator iterator
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
intermediate_table::const_iterator const_iterator
void findHoughClusters(const reco::HitPairListPtr &inputHits, reco::PrincipalComponents &pca, int &nLoops, RhoThetaAccumulatorBinMap &rhoThetaMap, HoughClusterList &clusterList) const
std::pair< recob::Seed, reco::HitPairListPtr > SeedHitPairListPair
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:246
std::list< HoughCluster > HoughClusterList
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
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...
std::map< BinIndex, AccumulatorBin > RhoThetaAccumulatorBinMap
void lar_cluster3d::HoughSeedFinderAlg::HoughRegionQuery ( BinIndex curBin,
RhoThetaAccumulatorBinMap rhoThetaAccumulatorBinMap,
HoughCluster neighborPts,
size_t  threshold 
) const
private
void lar_cluster3d::HoughSeedFinderAlg::HoughRegionQuery ( BinIndex curBin,
RhoThetaAccumulatorBinMap rhoThetaAccumulatorBinMap,
HoughCluster neighborPts,
size_t  threshold 
) const
private

Does a query of nearest neighbors to look for matching bins

Does a query of nearest neighbors to look for matching bins

Definition at line 192 of file HoughSeedFinderAlg.cxx.

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 }
intermediate_table::iterator iterator
void lar_cluster3d::HoughSeedFinderAlg::LineFit2DHits ( std::set< const reco::ClusterHit2D * > &  hitList,
double  XOrigin,
TVector3 &  Pos,
TVector3 &  Dir,
double &  ChiDOF 
) const
private
void lar_cluster3d::HoughSeedFinderAlg::LineFit2DHits ( std::set< const reco::ClusterHit2D * > &  hitList,
double  XOrigin,
TVector3 &  Pos,
TVector3 &  Dir,
double &  ChiDOF 
) const
private

Definition at line 989 of file HoughSeedFinderAlg.cxx.

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()
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
auto norm(Vector const &v)
Return norm of the specified vector.
Detector simulation of raw signals on wires.
#define A
Definition: memgrp.cpp:38
list x
Definition: train.py:276
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
void lar_cluster3d::HoughSeedFinderAlg::reconfigure ( fhicl::ParameterSet const &  pset)
virtual

a handler for the case where the algorithm control parameters are to be reset

Implements lar_cluster3d::SeedFinderAlgBase.

Definition at line 78 of file HoughSeedFinderAlg.cxx.

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 }
void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset

Member Data Documentation

std::vector< std::unique_ptr< TCanvas > > lar_cluster3d::HoughSeedFinderAlg::m_Canvases
mutableprivate

Graphical trace canvases.

Definition at line 130 of file HoughSeedFinderAlg.h.

bool lar_cluster3d::HoughSeedFinderAlg::m_displayHist
private

Definition at line 129 of file HoughSeedFinderAlg.h.

geo::Geometry const* lar_cluster3d::HoughSeedFinderAlg::m_geometry
private

Definition at line 122 of file HoughSeedFinderAlg.h.

geo::Geometry* lar_cluster3d::HoughSeedFinderAlg::m_geometry
private

Definition at line 124 of file HoughSeedFinderAlg.h.

double lar_cluster3d::HoughSeedFinderAlg::m_hiThresholdFrac
private

Definition at line 116 of file HoughSeedFinderAlg.h.

size_t lar_cluster3d::HoughSeedFinderAlg::m_hiThresholdMin
private

Definition at line 115 of file HoughSeedFinderAlg.h.

double lar_cluster3d::HoughSeedFinderAlg::m_loThresholdFrac
private

Definition at line 117 of file HoughSeedFinderAlg.h.

double lar_cluster3d::HoughSeedFinderAlg::m_maximumGap
private

Definition at line 122 of file HoughSeedFinderAlg.h.

int lar_cluster3d::HoughSeedFinderAlg::m_maxLoopsPerCluster
private

Definition at line 121 of file HoughSeedFinderAlg.h.

size_t lar_cluster3d::HoughSeedFinderAlg::m_minimum3DHits
private

Definition at line 112 of file HoughSeedFinderAlg.h.

double lar_cluster3d::HoughSeedFinderAlg::m_numAveDocas
private

Definition at line 119 of file HoughSeedFinderAlg.h.

size_t lar_cluster3d::HoughSeedFinderAlg::m_numSeed2DHits
private

Definition at line 118 of file HoughSeedFinderAlg.h.

int lar_cluster3d::HoughSeedFinderAlg::m_numSkippedHits
private

Definition at line 120 of file HoughSeedFinderAlg.h.

std::vector< TVirtualPad * > lar_cluster3d::HoughSeedFinderAlg::m_Pads
mutableprivate

View pads in current canvas.

Definition at line 131 of file HoughSeedFinderAlg.h.

PrincipalComponentsAlg lar_cluster3d::HoughSeedFinderAlg::m_pcaAlg
private

Definition at line 127 of file HoughSeedFinderAlg.h.

int lar_cluster3d::HoughSeedFinderAlg::m_rhoBins
private

Definition at line 114 of file HoughSeedFinderAlg.h.

int lar_cluster3d::HoughSeedFinderAlg::m_thetaBins
private

Definition at line 113 of file HoughSeedFinderAlg.h.


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