#include "HoughSeedFinderAlg.h"

Inheritance diagram for lar_cluster3d::HoughSeedFinderAlg:

Classes
class	AccumulatorBin

struct	SortBinIndexList

class	SortHoughClusterList

Public Member Functions
	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, AccumulatorBin >	RhoThetaAccumulatorBinMap

typedef std::list< BinIndex >	HoughCluster

typedef std::list< HoughCluster >	HoughClusterList

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, 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::Geometry const *	m_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...

Detailed Description

HoughSeedFinderAlg class.

Definition at line 36 of file HoughSeedFinderAlg.h.

Member Typedef Documentation

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

private

Definition at line 75 of file HoughSeedFinderAlg.h.

typedef std::list<BinIndex> lar_cluster3d::HoughSeedFinderAlg::HoughCluster

private

Definition at line 81 of file HoughSeedFinderAlg.h.

typedef std::list<HoughCluster> lar_cluster3d::HoughSeedFinderAlg::HoughClusterList

private

Definition at line 82 of file HoughSeedFinderAlg.h.

typedef std::map<BinIndex, AccumulatorBin> lar_cluster3d::HoughSeedFinderAlg::RhoThetaAccumulatorBinMap

private

Definition at line 80 of file HoughSeedFinderAlg.h.

Constructor & Destructor Documentation

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

Constructor.

Parameters

pset

Definition at line 41 of file HoughSeedFinderAlg.cxx.

References art::ServiceHandle< T, SCOPE >::get(), fhicl::ParameterSet::get(), m_displayHist, m_geometry, m_hiThresholdFrac, m_hiThresholdMin, m_loThresholdFrac, m_maximumGap, m_maxLoopsPerCluster, m_minimum3DHits, m_numAveDocas, m_numSeed2DHits, m_numSkippedHits, m_rhoBins, and m_thetaBins.

     : m_minimum3DHits(5)
     , m_thetaBins(360)
     , m_rhoBins(75)
     , m_hiThresholdMin(5)
     , m_hiThresholdFrac(.05)
     , m_loThresholdFrac(0.85)
     , m_numSeed2DHits(80)
     , m_numAveDocas(6.)
     , m_numSkippedHits(10)
     , m_maxLoopsPerCluster(3)
     , m_maximumGap(5.)
     , m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
     , m_displayHist(false)
   {
     m_minimum3DHits = pset.get<size_t>("Minimum3DHits", 5);
     m_thetaBins = pset.get<int>("ThetaBins", 360);
     m_rhoBins = pset.get<int>("HalfRhoBins", 75);
     m_hiThresholdMin = pset.get<size_t>("HiThresholdMin", 5);
     m_hiThresholdFrac = pset.get<double>("HiThresholdFrac", 0.05);
     m_loThresholdFrac = pset.get<double>("LoThresholdFrac", 0.85);
     m_numSeed2DHits = pset.get<size_t>("NumSeed2DHits", 80);
     m_numAveDocas = pset.get<double>("NumAveDocas", 6.);
     m_numSkippedHits = pset.get<int>("NumSkippedHits", 10);
     m_maxLoopsPerCluster = pset.get<int>("MaxLoopsPerCluster", 3);
     m_maximumGap = pset.get<double>("MaximumGap", 5.);
     m_displayHist = pset.get<bool>("DisplayHoughHist", false);
     m_geometry = art::ServiceHandle<geo::Geometry const>{}.get();
   }

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.

Definition at line 531 of file HoughSeedFinderAlg.cxx.

References reco::PrincipalComponents::getAveHitDoca(), reco::PrincipalComponents::getAvePosition(), reco::ClusterHit3D::getDocaToAxis(), reco::PrincipalComponents::getEigenVectors(), reco::ClusterHit3D::getHits(), reco::PrincipalComponents::getSvdOK(), LineFit2DHits(), m_minimum3DHits, m_numAveDocas, m_numSeed2DHits, m_numSkippedHits, m_pcaAlg, lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_3D(), and lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_calc3DDocas().

Referenced by findTrackSeeds().

   {
     if (seed3DHits.size() < m_minimum3DHits) return false;
 
     reco::PrincipalComponents seedFullPca;
 
     m_pcaAlg.PCAAnalysis_3D(seed3DHits, seedFullPca, true);
 
     if (!seedFullPca.getSvdOK()) return false;
 
     // Use the following to set the 3D doca and arclength for each hit
     m_pcaAlg.PCAAnalysis_calc3DDocas(seed3DHits, seedFullPca);
 
     // Use this info to sort the hits along the principle axis
     //seed3DHits.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
     seed3DHits.sort(SeedFinderAlgBase::Sort3DHitsByAbsArcLen3D());
 
     // The idea here is to search for the first hit that lies "close" to the principle axis
     // At that point we count out n hits to use as the seed
     reco::HitPairListPtr seedHit3DList;
     std::set<const reco::ClusterHit2D*> seedHitSet;
     double aveDocaToAxis = seedFullPca.getAveHitDoca();
     int gapCount(0);
 
     // Now loop through hits to search for a "continuous" block of at least m_numSeed2DHits
     // We'll arrive at that number by collecting 2D hits in an stl set which will keep track of unique occurances
     for (reco::HitPairListPtr::iterator peakBinItr = seed3DHits.begin();
          peakBinItr != seed3DHits.end();
          peakBinItr++) {
       const reco::ClusterHit3D* hit3D = *peakBinItr;
 
       if (hit3D->getDocaToAxis() < m_numAveDocas * aveDocaToAxis) {
         // Check if we need to reset because of gap count
         if (gapCount > m_numSkippedHits) {
           seedHit3DList.clear();
           seedHitSet.clear();
         }
 
         seedHit3DList.push_back(hit3D);
 
         for (const auto& hit : hit3D->getHits())
           seedHitSet.insert(hit);
 
         gapCount = 0;
       }
       else
         gapCount++;
 
       if (seedHitSet.size() > m_numSeed2DHits) break;
     }
 
     // If not enough hits then we are done
     if (seedHit3DList.size() < m_minimum3DHits) return false;
 
     reco::PrincipalComponents seedPca;
 
     // Use only the "seed" 3D hits to get new PCA axes
     m_pcaAlg.PCAAnalysis_3D(seedHit3DList, seedPca, true);
 
     if (!seedPca.getSvdOK()) return false;
 
     m_pcaAlg.PCAAnalysis_calc3DDocas(seedHit3DList, seedPca);
     //seedHit3DList.sort(SeedFinderAlgBase::Sort3DHitsByAbsArcLen3D());
     seedHit3DList.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
 
     // Now translate the seedCenter by the arc len to the first hit
     double seedCenter[3] = {
       seedPca.getAvePosition()[0], seedPca.getAvePosition()[1], seedPca.getAvePosition()[2]};
     double seedDir[3] = {seedPca.getEigenVectors().row(2)[0],
                          seedPca.getEigenVectors().row(2)[1],
                          seedPca.getEigenVectors().row(2)[2]};
 
     double arcLen = seedHit3DList.front()->getArclenToPoca();
     double seedStart[3] = {seedCenter[0] + arcLen * seedDir[0],
                            seedCenter[1] + arcLen * seedDir[1],
                            seedCenter[2] + arcLen * seedDir[2]};
 
     //seedStart[0] = seedHit3DList.front()->getX();
     //seedStart[1] = seedHit3DList.front()->getY();
     //seedStart[2] = seedHit3DList.front()->getZ();
 
     if (seedHitSet.size() >= 10) {
       TVector3 newSeedPos;
       TVector3 newSeedDir;
       double chiDOF;
 
       LineFit2DHits(seedHitSet, seedStart[0], newSeedPos, newSeedDir, chiDOF);
 
       if (chiDOF > 0.) {
         // check angles between new/old directions
         double cosAng =
           seedDir[0] * newSeedDir[0] + seedDir[1] * newSeedDir[1] + seedDir[2] * newSeedDir[2];
 
         if (cosAng < 0.) newSeedDir *= -1.;
 
         seedStart[0] = newSeedPos[0];
         seedStart[1] = newSeedPos[1];
         seedStart[2] = newSeedPos[2];
         seedDir[0] = newSeedDir[0];
         seedDir[1] = newSeedDir[1];
         seedDir[2] = newSeedDir[2];
       }
     }
 
     // Keep track of this seed and the 3D hits that make it up
     seedHitPair = SeedHitPairListPair(recob::Seed(seedStart, seedDir), seedHit3DList);
 
     // 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
     if (seed3DHits.size() > 100) {
       // Need to reset the doca/arclen first
       m_pcaAlg.PCAAnalysis_calc3DDocas(seed3DHits, seedFullPca);
 
       // Now resort the hits
       seed3DHits.sort(SeedFinderAlgBase::Sort3DHitsByAbsArcLen3D());
 
       size_t numToKeep = seed3DHits.size() - 10;
 
       reco::HitPairListPtr::iterator endItr = seed3DHits.begin();
 
       std::advance(endItr, numToKeep);
 
       seed3DHits.erase(endItr, seed3DHits.end());
     }
 
     return true;
   }

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

Definition at line 209 of file HoughSeedFinderAlg.cxx.

References HoughRegionQuery(), lar_cluster3d::HoughSeedFinderAlg::AccumulatorBin::isInCluster(), lar_cluster3d::HoughSeedFinderAlg::AccumulatorBin::isVisited(), lar_cluster3d::HoughSeedFinderAlg::AccumulatorBin::setInCluster(), and lar_cluster3d::HoughSeedFinderAlg::AccumulatorBin::setVisited().

Referenced by findHoughClusters().

   {
     // Start by adding the input point to our Hough Cluster
     houghCluster.push_back(curBin);
 
     for (auto& binIndex : neighborPts) {
       AccumulatorBin& accBin = rhoThetaAccumulatorBinMap[binIndex];
 
       if (!accBin.isVisited()) {
         accBin.setVisited();
 
         HoughCluster nextNeighborPts;
 
         HoughRegionQuery(binIndex, rhoThetaAccumulatorBinMap, nextNeighborPts, threshold);
 
         if (!nextNeighborPts.empty()) {
           for (auto& neighborIdx : nextNeighborPts) {
             neighborPts.push_back(neighborIdx);
           }
         }
       }
 
       if (!accBin.isInCluster()) {
         houghCluster.push_back(binIndex);
         accBin.setInCluster();
       }
     }
 
     return;
   }

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 283 of file HoughSeedFinderAlg.cxx.

References reco::PrincipalComponents::getSvdOK(), m_maximumGap, m_pcaAlg, lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_3D(), and lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_calc3DDocas().

Referenced by findTrackSeeds().

   {
     // Intention is to try to find the largest "contiguous" block of hits in the input list
     // In a nutshell, the idea is to order input hits according to the pca, then
     // loop through the hits and store them in a new hit list. If a gap is detected then
     // we terminate the previous list, create a new one and continue. After the loop over
     // hits is complete then simply pick the biggest list.
     // Step I is to run the pca and order the hits
     reco::PrincipalComponents pca;
 
     m_pcaAlg.PCAAnalysis_3D(inputHitList, pca, true);
 
     // It would seem that the analysis can fail!
     if (!pca.getSvdOK()) {
       outputList = inputHitList;
       return;
     }
 
     m_pcaAlg.PCAAnalysis_calc3DDocas(inputHitList, pca);
 
     inputHitList.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
     outputList.clear();
 
     // Create containers to hold our hit lists
     reco::HitPairListPtr continuousHitList;
 
     std::map<size_t, reco::HitPairListPtr> gapHitListMap;
 
     // Remember the distance arc length of the last hit
     double arcLenLastHit = inputHitList.front()->getArclenToPoca();
 
     // Loop through the input hits
     for (const auto& hit3D : inputHitList) {
       // Hits in order, delta arclength should always be positive
       double arcLen = hit3D->getArclenToPoca();
       double deltaArcLen = arcLen - arcLenLastHit;
 
       if (deltaArcLen > m_maximumGap) {
         gapHitListMap[continuousHitList.size()] = continuousHitList;
         continuousHitList.clear();
       }
 
       continuousHitList.emplace_back(hit3D);
 
       arcLenLastHit = arcLen;
     }
 
     if (!continuousHitList.empty()) gapHitListMap[continuousHitList.size()] = continuousHitList;
 
     // Get the largest list of hits
     std::map<size_t, reco::HitPairListPtr>::reverse_iterator longestListItr =
       gapHitListMap.rbegin();
 
     if (longestListItr != gapHitListMap.rend()) {
       size_t nContinuousHits = longestListItr->first;
       reco::HitPairListPtr longestList = longestListItr->second;
 
       outputList.resize(nContinuousHits);
       std::copy(longestList.begin(), longestList.end(), outputList.begin());
     }
 
     return;
   }

void lar_cluster3d::HoughSeedFinderAlg::findHoughClusters	(	const reco::HitPairListPtr &	inputHits,
		reco::PrincipalComponents &	pca,
		RhoThetaAccumulatorBinMap &	rhoThetaMap,
		HoughClusterList &	clusterList
	)		const

private

Definition at line 348 of file HoughSeedFinderAlg.cxx.

References expandHoughCluster(), reco::PrincipalComponents::getAvePosition(), reco::PrincipalComponents::getEigenValues(), reco::PrincipalComponents::getEigenVectors(), reco::ClusterHit3D::getPosition(), reco::ClusterHit3D::getStatusBits(), HoughRegionQuery(), m_Canvases, m_displayHist, m_geometry, m_hiThresholdFrac, m_hiThresholdMin, m_loThresholdFrac, m_Pads, m_rhoBins, m_thetaBins, util::size(), and geo::GeometryCore::WirePitch().

Referenced by findTrackHits(), and findTrackSeeds().

   {
     // The goal of this function is to do a basic Hough Transform on the input list of 3D hits.
     // In order to transform this to a 2D problem, the 3D hits are projected to the plane of the two
     // largest eigen values from the input principal components analysis axis.
     // There are two basic steps to the job here:
     // 1) Build and accumlate a rho-theta map
     // 2) Go through that rho-theta map and find candidate Hough "clusters"
     // Unfortunately, the following may not be suitable viewing for those who may be feint of heart
     //
     // Define some constants
     static int histCount(0);
     const double maxTheta(M_PI);                               // Obviously, 180 degrees
     const double thetaBinSize(maxTheta / double(m_thetaBins)); // around 4 degrees (45 bins)
     const double rhoBinSizeMin(m_geometry->WirePitch()); // Wire spacing gives a natural bin size?
 
     // Recover the parameters from the Principal Components Analysis that we need to project and accumulate
     Eigen::Vector3f pcaCenter(
       pca.getAvePosition()[0], pca.getAvePosition()[1], pca.getAvePosition()[2]);
     Eigen::Vector3f planeVec0(pca.getEigenVectors().row(2));
     Eigen::Vector3f planeVec1(pca.getEigenVectors().row(1));
     Eigen::Vector3f pcaPlaneNrml(pca.getEigenVectors().row(0));
     double eigenVal0 = 3. * sqrt(pca.getEigenValues()[2]);
     double eigenVal1 = 3. * sqrt(pca.getEigenValues()[1]);
     double maxRho = std::sqrt(eigenVal0 * eigenVal0 + eigenVal1 * eigenVal1) * 2. / 3.;
     double rhoBinSize = maxRho / double(m_rhoBins);
 
     if (rhoBinSize < rhoBinSizeMin) rhoBinSize = rhoBinSizeMin;
 
     // **********************************************************************
     // Part I: Accumulate values in the rho-theta map
     // **********************************************************************
 
     size_t maxBinCount(0);
     int nAccepted3DHits(0);
 
     // Commence looping over the input 3D hits and fill our accumulator bins
     for (reco::HitPairListPtr::const_iterator hit3DItr = hitPairListPtr.begin();
          hit3DItr != hitPairListPtr.end();
          hit3DItr++) {
       // Recover the hit
       const reco::ClusterHit3D* hit3D(*hit3DItr);
 
       // Skip hits which are not skeleton points
       if (!(hit3D->getStatusBits() & 0x10000000)) continue;
 
       nAccepted3DHits++;
 
       Eigen::Vector3f hit3DPosition(
         hit3D->getPosition()[0], hit3D->getPosition()[1], hit3D->getPosition()[2]);
       Eigen::Vector3f pcaToHitVec = hit3DPosition - pcaCenter;
       Eigen::Vector3f pcaToHitPlaneVec(pcaToHitVec.dot(planeVec0), pcaToHitVec.dot(planeVec1), 0.);
       double xPcaToHit = pcaToHitPlaneVec[0];
       double yPcaToHit = pcaToHitPlaneVec[1];
 
       // Create an accumulator value
       AccumulatorValues accValue(pcaToHitPlaneVec, hit3DItr);
 
       // Commence loop over theta to fill accumulator bins
       // Note that with theta in the range 0-pi then we can have negative values for rho
       for (int thetaIdx = 0; thetaIdx < m_thetaBins; thetaIdx++) {
         // We need to convert our theta index to an angle
         double theta = thetaBinSize * double(thetaIdx);
 
         // calculate rho for this angle
         double rho = xPcaToHit * std::cos(theta) + yPcaToHit * std::sin(theta);
 
         // Get the rho index
         int rhoIdx = std::round(rho / rhoBinSize);
 
         // Accumulate
         BinIndex binIndex(rhoIdx, thetaIdx);
 
         rhoThetaAccumulatorBinMap[binIndex].addAccumulatorValue(accValue);
 
         if (rhoThetaAccumulatorBinMap[binIndex].getAccumulatorValues().size() > maxBinCount)
           maxBinCount = rhoThetaAccumulatorBinMap[binIndex].getAccumulatorValues().size();
       }
     }
 
     // Accumulation done, if asked now display the hist
     if (m_displayHist) {
       std::ostringstream ostr;
       ostr << "Hough Histogram " << histCount++;
       m_Canvases.emplace_back(new TCanvas(ostr.str().c_str(), ostr.str().c_str(), 1000, 1000));
 
       std::ostringstream ostr2;
       ostr2 << "Plane";
 
       m_Canvases.back()->GetFrame()->SetFillColor(46);
       m_Canvases.back()->SetFillColor(19);
       m_Canvases.back()->SetBorderMode(19);
       m_Canvases.back()->cd(1);
 
       double zmin = 0.06;
       double zmax = 0.94;
       double xmin = 0.04;
       double xmax = 0.95;
       TPad* p = new TPad(ostr2.str().c_str(), ostr2.str().c_str(), zmin, xmin, zmax, xmax);
       p->SetBit(kCanDelete); // Give away ownership.
       p->Range(zmin, xmin, zmax, xmax);
       p->SetFillStyle(4000); // Transparent.
       p->Draw();
       m_Pads.push_back(p);
 
       TH2D* houghHist = new TH2D("HoughHist",
                                  "Hough Space",
                                  2 * m_rhoBins,
                                  -m_rhoBins + 0.5,
                                  m_rhoBins + 0.5,
                                  m_thetaBins,
                                  0.,
                                  m_thetaBins);
 
       for (const auto& rhoThetaMap : rhoThetaAccumulatorBinMap) {
         houghHist->Fill(rhoThetaMap.first.first,
                         rhoThetaMap.first.second + 0.5,
                         rhoThetaMap.second.getAccumulatorValues().size());
       }
 
       houghHist->SetBit(kCanDelete);
       houghHist->Draw();
       m_Canvases.back()->Update();
     }
 
     // **********************************************************************
     // Part II: Use DBScan (or a slight variation) to find clusters of bins
     // **********************************************************************
 
     size_t thresholdLo = std::max(size_t(m_hiThresholdFrac * nAccepted3DHits), m_hiThresholdMin);
     size_t thresholdHi = m_loThresholdFrac * maxBinCount;
 
     std::list<RhoThetaAccumulatorBinMap::iterator> binIndexList;
 
     for (RhoThetaAccumulatorBinMap::iterator mapItr = rhoThetaAccumulatorBinMap.begin();
          mapItr != rhoThetaAccumulatorBinMap.end();
          mapItr++)
       binIndexList.push_back(mapItr);
 
     binIndexList.sort(SortBinIndexList());
 
     for (auto& mapItr : binIndexList) {
       // If we have been here before we skip
       //if (mapItr.second.isVisited()) continue;
       if (mapItr->second.isInCluster()) continue;
 
       // Mark this bin as visited
       // Actually, don't mark it since we are double thresholding and don't want it missed
       //mapItr.second.setVisited();
 
       // Make sure over threshold
       if (mapItr->second.getAccumulatorValues().size() < thresholdLo) {
         mapItr->second.setNoise();
         continue;
       }
 
       // Set the low threshold to make sure we merge bins that might be either side of a boundary trajectory
       thresholdHi = std::max(
         size_t(m_loThresholdFrac * mapItr->second.getAccumulatorValues().size()), m_hiThresholdMin);
 
       // Recover our neighborhood
       HoughCluster neighborhood;
       BinIndex curBin(mapItr->first);
 
       HoughRegionQuery(curBin, rhoThetaAccumulatorBinMap, neighborhood, thresholdHi);
 
       houghClusters.push_back(HoughCluster());
 
       HoughCluster& houghCluster = houghClusters.back();
 
       expandHoughCluster(
         curBin, neighborhood, houghCluster, rhoThetaAccumulatorBinMap, thresholdHi);
     }
 
     // Sort the clusters using the SortHoughClusterList metric
     if (!houghClusters.empty()) houghClusters.sort(SortHoughClusterList(rhoThetaAccumulatorBinMap));
 
     return;
   }

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.

Definition at line 849 of file HoughSeedFinderAlg.cxx.

References findHoughClusters(), reco::PrincipalComponents::getEigenValues(), and m_minimum3DHits.

Referenced by lar_cluster3d::Cluster3D::splitClustersWithHough().

   {
     // The goal of this routine is run the Hough Transform on the input set of hits
     // and then to return a list of lists of hits which are associated to a given line
 
     // Make sure we are using the right pca
     reco::HitPairListPtr hitPairListPtr = inputHitPairListPtr;
 
     // Make a local copy of the input PCA
     reco::PrincipalComponents pca = inputPCA;
 
     // We also require that there be some spread in the data, otherwise not worth running?
     double eigenVal0 = 3. * sqrt(pca.getEigenValues()[2]);
     double eigenVal1 = 3. * sqrt(pca.getEigenValues()[1]);
 
     if (eigenVal0 > 5. && eigenVal1 > 0.001) {
       // **********************************************************************
       // Part I: Build Hough space and find Hough clusters
       // **********************************************************************
       RhoThetaAccumulatorBinMap rhoThetaAccumulatorBinMap;
       HoughClusterList houghClusters;
 
       findHoughClusters(hitPairListPtr, pca, rhoThetaAccumulatorBinMap, houghClusters);
 
       // **********************************************************************
       // Part II: Go through the clusters to find the peak bins
       // **********************************************************************
 
       // We need to use a set so we can be sure to have unique hits
       reco::HitPairListPtr clusterHitsList;
       std::set<const reco::ClusterHit3D*> masterHitPtrList;
       std::set<const reco::ClusterHit3D*> peakBinPtrList;
 
       size_t firstPeakCount(0);
 
       // Loop through the list of all clusters found above
       for (auto& houghCluster : houghClusters) {
         BinIndex peakBin = houghCluster.front();
         size_t peakCount = 0;
 
         // Make a local (to this cluster) set of of hits
         std::set<const reco::ClusterHit3D*> localHitPtrList;
 
         // Now loop through the bins that were attached to this cluster
         for (auto& binIndex : houghCluster) {
           // An even more local list so we can keep track of peak values
           std::set<const reco::ClusterHit3D*> tempHitPtrList;
 
           // Recover the hits associated to this cluster
           for (auto& hitItr : rhoThetaAccumulatorBinMap[binIndex].getAccumulatorValues()) {
             reco::HitPairListPtr::const_iterator hit3DItr = hitItr.getHitIterator();
 
             tempHitPtrList.insert(*hit3DItr);
           }
 
           // Trim out any hits already used by a bigger/better cluster
           std::set<const reco::ClusterHit3D*> tempHit3DSet;
 
           std::set_difference(tempHitPtrList.begin(),
                               tempHitPtrList.end(),
                               masterHitPtrList.begin(),
                               masterHitPtrList.end(),
                               std::inserter(tempHit3DSet, tempHit3DSet.end()));
 
           tempHitPtrList = tempHit3DSet;
 
           size_t binCount = tempHitPtrList.size();
 
           if (peakCount < binCount) {
             peakCount = binCount;
             peakBin = binIndex;
             peakBinPtrList = tempHitPtrList;
           }
 
           // Add this to our local list
           localHitPtrList.insert(tempHitPtrList.begin(), tempHitPtrList.end());
         }
 
         if (localHitPtrList.size() < m_minimum3DHits) continue;
 
         if (!firstPeakCount) firstPeakCount = peakCount;
 
         // If the peak counts are significantly less than the first cluster's peak then skip
         if (peakCount < firstPeakCount / 10) continue;
 
         // **********************************************************************
         // Part III: Make a list of hits from the total number associated
         // **********************************************************************
 
         hitPairListPtrList.push_back(reco::HitPairListPtr());
 
         hitPairListPtrList.back().resize(localHitPtrList.size());
         std::copy(
           localHitPtrList.begin(), localHitPtrList.end(), hitPairListPtrList.back().begin());
 
         // We want to remove the hits which have been used from further contention
         masterHitPtrList.insert(localHitPtrList.begin(), localHitPtrList.end());
 
         if (hitPairListPtr.size() - masterHitPtrList.size() < m_minimum3DHits) break;
       } // end loop over hough clusters
     }
 
     return true;
   }

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.

Definition at line 659 of file HoughSeedFinderAlg.cxx.

References buildSeed(), findHitGaps(), findHoughClusters(), reco::PrincipalComponents::getEigenValues(), reco::PrincipalComponents::getSvdOK(), m_maxLoopsPerCluster, m_minimum3DHits, m_pcaAlg, and lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_3D().

Referenced by lar_cluster3d::Cluster3D::findTrackSeeds().

   {
     // This will be a busy routine... the basic tasks are:
     // 1) loop through hits and project to the plane defined by the two largest eigen values, accumulate in Hough space
     // 2) "Cluster" the Hough space to associate hits which are common to a line
     // 3) Process these clusters (still to be defined exactly)
 
     // Create an interim data structure which will allow us to sort our seeds by "best"
     // before we return them in the seedHitPairVec
     typedef std::map<size_t, SeedHitPairListPairVec> SizeToSeedToHitMap;
 
     SizeToSeedToHitMap seedHitPairMap;
     SeedHitPairListPair seedHitPair;
 
     // Make sure we are using the right pca
     reco::HitPairListPtr hitPairListPtr = inputHitPairListPtr;
 
     int nLoops(0);
 
     // Make a local copy of the input PCA
     reco::PrincipalComponents pca = inputPCA;
 
     // We loop over hits in our list until there are no more
     while (!hitPairListPtr.empty()) {
       // We also require that there be some spread in the data, otherwise not worth running?
       double eigenVal0 = 3. * sqrt(pca.getEigenValues()[2]);
       double eigenVal1 = 3. * sqrt(pca.getEigenValues()[1]);
 
       if (eigenVal0 > 5. && eigenVal1 > 0.001) {
         // **********************************************************************
         // Part I: Build Hough space and find Hough clusters
         // **********************************************************************
         RhoThetaAccumulatorBinMap rhoThetaAccumulatorBinMap;
         HoughClusterList houghClusters;
 
         findHoughClusters(hitPairListPtr, pca, rhoThetaAccumulatorBinMap, houghClusters);
 
         // If no clusters then done
         if (houghClusters.empty()) break;
 
         // **********************************************************************
         // Part II: Go through the clusters to find the peak bins
         // **********************************************************************
 
         // We need to use a set so we can be sure to have unique hits
         reco::HitPairListPtr clusterHitsList;
         std::set<const reco::ClusterHit3D*> masterHitPtrList;
         std::set<const reco::ClusterHit3D*> peakBinPtrList;
 
         size_t firstPeakCount(0);
 
         // Loop through the list of all clusters found above
         for (auto& houghCluster : houghClusters) {
           BinIndex peakBin = houghCluster.front();
           size_t peakCount = 0;
 
           // Make a local (to this cluster) set of of hits
           std::set<const reco::ClusterHit3D*> localHitPtrList;
 
           // Now loop through the bins that were attached to this cluster
           for (auto& binIndex : houghCluster) {
             // An even more local list so we can keep track of peak values
             std::set<const reco::ClusterHit3D*> tempHitPtrList;
 
             // Recover the hits associated to this cluster
             for (auto& hitItr : rhoThetaAccumulatorBinMap[binIndex].getAccumulatorValues()) {
               reco::HitPairListPtr::const_iterator hit3DItr = hitItr.getHitIterator();
 
               tempHitPtrList.insert(*hit3DItr);
             }
 
             // Trim out any hits already used by a bigger/better cluster
             std::set<const reco::ClusterHit3D*> tempHit3DSet;
 
             std::set_difference(tempHitPtrList.begin(),
                                 tempHitPtrList.end(),
                                 masterHitPtrList.begin(),
                                 masterHitPtrList.end(),
                                 std::inserter(tempHit3DSet, tempHit3DSet.end()));
 
             tempHitPtrList = tempHit3DSet;
 
             size_t binCount = tempHitPtrList.size();
 
             if (peakCount < binCount) {
               peakCount = binCount;
               peakBin = binIndex;
               peakBinPtrList = tempHitPtrList;
             }
 
             // Add this to our local list
             localHitPtrList.insert(tempHitPtrList.begin(), tempHitPtrList.end());
           }
 
           if (localHitPtrList.size() < m_minimum3DHits) continue;
 
           if (!firstPeakCount) firstPeakCount = peakCount;
 
           // If the peak counts are significantly less than the first cluster's peak then skip
           if (peakCount < firstPeakCount / 10) continue;
 
           // **********************************************************************
           // Part III: Make a Seed from the peak bin hits
           // **********************************************************************
 
           reco::HitPairListPtr allPeakBinHits;
 
           for (const auto& hit3D : localHitPtrList)
             allPeakBinHits.push_back(hit3D);
 
           reco::HitPairListPtr peakBinHits;
 
           // Find longest "continuous" set of hits and use these for the seed
           findHitGaps(allPeakBinHits, peakBinHits);
 
           if (peakBinHits.size() < m_minimum3DHits) continue;
 
           // We now build the actual seed.
           if (buildSeed(peakBinHits, seedHitPair)) {
             // Keep track of this in our map (which will do ordering for us)
             seedHitPairMap[peakBinHits.size()].push_back(seedHitPair);
 
             // For visual testing in event display, mark all the hits in the first seed so we can see them
             if (seedHitPairMap.size() == 1) {
               for (const auto& hit3D : peakBinHits)
                 hit3D->setStatusBit(0x40000000);
             }
 
             //                    for(const auto& hit3D : seedHitPair.second) hit3D->setStatusBit(0x40000000);
           }
 
           // Our peakBinHits collection will most likely be a subset of the localHitPtrList collection
           // We want to remove only the "pure" hits which are those in the peakBinHits collection
           // So, sort them and then add to our master list
           peakBinHits.sort();
 
           masterHitPtrList.insert(peakBinHits.begin(), peakBinHits.end());
 
           if (hitPairListPtr.size() - masterHitPtrList.size() < m_minimum3DHits) break;
         } // end loop over hough clusters
 
         // If the masterHitPtrList is empty then nothing happened and we're done
         if (masterHitPtrList.empty()) break;
 
         // **********************************************************************
         // Part IV: Remove remaining peak bin hits from HitPairPtrList
         // **********************************************************************
 
         hitPairListPtr.sort();
 
         reco::HitPairListPtr::iterator newListEnd = std::set_difference(hitPairListPtr.begin(),
                                                                         hitPairListPtr.end(),
                                                                         masterHitPtrList.begin(),
                                                                         masterHitPtrList.end(),
                                                                         hitPairListPtr.begin());
 
         hitPairListPtr.erase(newListEnd, hitPairListPtr.end());
 
         if (hitPairListPtr.size() < m_minimum3DHits) break;
 
         if (nLoops++ > m_maxLoopsPerCluster) break;
 
         // ********************************************************
       }
       else
         break; // eigen values not in range
 
       // At this point run the PCA on the remaining hits
       m_pcaAlg.PCAAnalysis_3D(hitPairListPtr, pca, true);
 
       if (!pca.getSvdOK()) break;
     }
 
     // The final task before returning is to transfer the stored seeds into the output seed vector
     // What we want to do is make sure the first seeds are the "best" seeds which is defined as the
     // seeds which were associated to the most hits by the Hough Transform. Our seed map will have
     // the reverse of this ordering so we simply iterate through it "backwards"
     for (SizeToSeedToHitMap::reverse_iterator seedMapItr = seedHitPairMap.rbegin();
          seedMapItr != seedHitPairMap.rend();
          seedMapItr++) {
       for (const auto& seedHitPair : seedMapItr->second) {
         seedHitPairVec.emplace_back(seedHitPair);
       }
     }
 
     return true;
   }

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

Definition at line 173 of file HoughSeedFinderAlg.cxx.

References m_thetaBins.

Referenced by expandHoughCluster(), and findHoughClusters().

   {
     // We simply loop over the nearest indices and see if we have any friends over threshold
     for (int rhoIdx = curBin.first - 1; rhoIdx <= curBin.first + 1; rhoIdx++) {
       for (int jdx = curBin.second - 1; jdx <= curBin.second + 1; jdx++) {
         // Skip the self reference
         if (rhoIdx == curBin.first && jdx == curBin.second) continue;
 
         // Theta bin needs to handle the wrap.
         int thetaIdx(jdx);
 
         if (thetaIdx < 0)
           thetaIdx = m_thetaBins - 1;
         else if (thetaIdx > m_thetaBins - 1)
           thetaIdx = 0;
 
         BinIndex binIndex(rhoIdx, thetaIdx);
         RhoThetaAccumulatorBinMap::iterator mapItr = rhoThetaAccumulatorBinMap.find(binIndex);
 
         if (mapItr != rhoThetaAccumulatorBinMap.end()) {
           if (mapItr->second.getAccumulatorValues().size() >= threshold)
             neighborPts.push_back(binIndex);
         }
       }
     }
 
     return;
   }

void lar_cluster3d::HoughSeedFinderAlg::LineFit2DHits	(	std::set< const reco::ClusterHit2D * > &	hitList,
		double	XOrigin,
		TVector3 &	Pos,
		TVector3 &	Dir,
		double &	ChiDOF
	)		const

private

Definition at line 957 of file HoughSeedFinderAlg.cxx.

References m_geometry, norm, geo::PlaneID::Plane, sw, w, geo::WireID::Wire, geo::GeometryCore::WireCoordinate(), geo::GeometryCore::WirePitch(), and x.

Referenced by buildSeed().

   {
     // The following is lifted from Bruce Baller to try to get better
     // initial parameters for a candidate Seed. It is slightly reworked
     // which is why it is included here instead of used as is.
     //
     // Linear fit using X as the independent variable. Hits to be fitted
     // are passed in the hits vector in a pair form (X, WireID). The
     // fitted track position at XOrigin is returned in the Pos vector.
     // The direction cosines are returned in the Dir vector.
     //
     // SVD fit adapted from $ROOTSYS/tutorials/matrix/solveLinear.C
     // Fit equation is w = A(X)v, where w is a vector of hit wires, A is
     // a matrix to calculate a track projected to a point at X, and v is
     // a vector (Yo, Zo, dY/dX, dZ/dX).
     //
     // Note: The covariance matrix should also be returned
     // B. Baller August 2014
 
     // assume failure
     ChiDOF = -1;
 
     if (hit2DSet.size() < 4) return;
 
     const unsigned int nvars = 4;
     unsigned int npts = hit2DSet.size();
 
     TMatrixD A(npts, nvars);
     // vector holding the Wire number
     TVectorD w(npts);
     unsigned short ninpl[3] = {0};
     unsigned short nok = 0;
     unsigned short iht(0);
 
     // Loop over unique 2D hits from the input list of 3D hits
     for (const auto& hit : hit2DSet) {
       geo::WireID const& wireID = hit->WireID();
       unsigned int const ipl = wireID.Plane;
 
       // get the wire plane offset
       double const off = m_geometry->WireCoordinate(geo::Point_t{}, wireID);
 
       // get the "cosine-like" component
       double const cw = m_geometry->WireCoordinate(geo::Point_t{0., 1., 0.}, wireID) - off;
 
       // the "sine-like" component
       double const sw = m_geometry->WireCoordinate(geo::Point_t{0., 0., 1.}, wireID) - off;
 
       double const x = hit->getXPosition() - XOrigin;
 
       A[iht][0] = cw;
       A[iht][1] = sw;
       A[iht][2] = cw * x;
       A[iht][3] = sw * x;
       w[iht] = wireID.Wire - off;
 
       ++ninpl[ipl];
 
       // need at least two points in a plane
       if (ninpl[ipl] == 2) ++nok;
 
       ++iht;
     }
 
     // need at least 2 planes with at least two points
     if (nok < 2) return;
 
     TDecompSVD svd(A);
     bool ok;
     TVectorD tVec = svd.Solve(w, ok);
 
     ChiDOF = 0;
 
     // not enough points to calculate Chisq
     if (npts <= 4) return;
 
     for (const auto& hit : hit2DSet) {
       geo::WireID const& wireID = hit->WireID();
 
       double const off = m_geometry->WireCoordinate(geo::Point_t{0., 0., 0.}, wireID);
       double const cw = m_geometry->WireCoordinate(geo::Point_t{0., 1., 0.}, wireID) - off;
       double const sw = m_geometry->WireCoordinate(geo::Point_t{0., 0., 1.}, wireID) - off;
       double const x = hit->getXPosition() - XOrigin;
       double const ypr = tVec[0] + tVec[2] * x;
       double const zpr = tVec[1] + tVec[3] * x;
       double const diff = ypr * cw + zpr * sw - (wireID.Wire - off);
       ChiDOF += diff * diff;
     }
 
     float werr2 = m_geometry->WirePitch() * m_geometry->WirePitch();
     ChiDOF /= werr2;
     ChiDOF /= (float)(npts - 4);
 
     double norm = sqrt(1 + tVec[2] * tVec[2] + tVec[3] * tVec[3]);
     Dir[0] = 1 / norm;
     Dir[1] = tVec[2] / norm;
     Dir[2] = tVec[3] / norm;
 
     Pos[0] = XOrigin;
     Pos[1] = tVec[0];
     Pos[2] = tVec[1];
 
   } // TrkLineFit()