Cluster3D class. More...

#include "SkeletonAlg.h"

Public Member Functions
	SkeletonAlg (fhicl::ParameterSet const &pset)
	Constructor. More...

virtual	~SkeletonAlg ()
	Destructor. More...

void	reconfigure (fhicl::ParameterSet const &pset)
	a handler for the case where the algorithm control parameters are to be reset More...

int	FindMedialSkeleton (reco::HitPairListPtr &hitPairList) const
	This is intended to find the medial skeleton given a list of input hit pairs. More...

void	GetSkeletonHits (const reco::HitPairListPtr &inputHitList, reco::HitPairListPtr &skeletonHitList) const
	Return the skeleton hits from the input list. More...

void	AverageSkeletonPositions (reco::HitPairListPtr &skeletonHitList) const
	Modifies the position of input skeleton hits by averaging along the "best" wire direction. More...

Private Member Functions
double	FindFirstAndLastWires (std::vector< const reco::ClusterHit3D * > &hitVec, int planeToCheck, int referenceWire, double referenceTicks, int &firstWire, int &lastWire) const
	A function to find the bounding wires in a given view. More...

Private Attributes
double	m_minimumDeltaTicks

double	m_maximumDeltaTicks

fhicl::ParameterSet	m_pset

Detailed Description

Cluster3D class.

Definition at line 33 of file SkeletonAlg.h.

Constructor & Destructor Documentation

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

Constructor.

Parameters

pset

Definition at line 30 of file SkeletonAlg.cxx.

References reconfigure().

 {
     reconfigure(pset);
 }

lar_cluster3d::SkeletonAlg::~SkeletonAlg ( )

virtual

Destructor.

Definition at line 37 of file SkeletonAlg.cxx.

38 {

39 }

Member Function Documentation

void lar_cluster3d::SkeletonAlg::AverageSkeletonPositions ( reco::HitPairListPtr & skeletonHitList ) const

Modifies the position of input skeleton hits by averaging along the "best" wire direction.

Parameters

skeletonHitList - input list of skeleton hits

Definition at line 331 of file SkeletonAlg.cxx.

Referenced by lar_cluster3d::Cluster3D::findTrackSeeds(), lar_cluster3d::Cluster3D::splitClustersWithHough(), and lar_cluster3d::Cluster3D::splitClustersWithMST().

 {
     // NOTE: This method assumes the list being given to it is comprised of skeleton hits
     //       YMMV if you send in a complete hit collection!
 
     // Define a mapping between 2D hits and the 3D hits they make
     typedef std::map<const reco::ClusterHit2D*, std::vector<const reco::ClusterHit3D*> > Hit2DtoHit3DMap;
     
     // We want to keep track of this mapping by view
     Hit2DtoHit3DMap hit2DToHit3DMap[3];
     
     // Keep count of the number of skeleton hits selected
     unsigned int nSkeletonHits(0);
     
     // Execute the first loop through the hits to build the map
     for(const auto& hitPair : skeletonHitList)
     {
         // Don't consider points "rejected" earlier
         if (hitPair->bitsAreSet(reco::ClusterHit3D::REJECTEDHIT)) continue;
         
         // Count only those skeleton hits which have not been averaged
         if (!hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONPOSAVE)) nSkeletonHits++;
         
         for(const auto& hit2D : hitPair->getHits())
         {
             size_t plane = hit2D->getHit().WireID().Plane;
             
             hit2DToHit3DMap[plane][hit2D].push_back(hitPair);
         }
     }
     
     // Exit early if no skeleton hits to process
     if (!nSkeletonHits) return;
     
     // The list of 3D hits associated to any given 2D hit is most useful to us if it is ordered
     // So this will loop through entries in the map and do the ordering
     for(size_t idx = 0; idx < 3; idx++)
     {
         for(auto& mapPair : hit2DToHit3DMap[idx])
         {
             size_t numHitPairs = mapPair.second.size();
             int    planeToCheck = (idx + 1) % 3;    // have forgotten reasoning here...
             
             if (numHitPairs > 1) std::sort(mapPair.second.begin(), mapPair.second.end(), OrderHitsAlongWire(planeToCheck));
         }
     }
     
     // Ok, so the basic strategy is not entirely different from that used to build the skeleton hits in the first
     // place. The idea is to loop through all skeleton hits and then determine the average position for the hit
     // based on the wire direction with the fewest hits to the edge.
     // Currently this is a pretty simple minded approach and it would appear that some effort could be spent
     // here to improve this.
     // NOTE: the averaging being performed is ONLY in the Y-Z plane since this is where the hit ambiguity
     //       issue arises.
     reco::HitPairListPtr::iterator hitPairItr = skeletonHitList.begin();
     
     for(int bestPlaneVecIdx = 0; bestPlaneVecIdx < 2; bestPlaneVecIdx++)
     {
         std::list<reco::ClusterHit3D> tempHitPairList;
         reco::HitPairListPtr          tempHitPairListPtr;
         
         std::map<const reco::ClusterHit3D*, const reco::ClusterHit3D*> hit3DToHit3DMap;
         
         while(hitPairItr != skeletonHitList.end())
         {
             const reco::ClusterHit3D* hit3D = *hitPairItr++;
             
             std::vector<std::pair<size_t,size_t> > bestPlaneVec;
             
             for(const auto& hit2D : hit3D->getHits())
             {
                 bestPlaneVec.push_back(std::pair<size_t,size_t>(hit2D->getHit().WireID().Plane,hit2DToHit3DMap[hit2D->getHit().WireID().Plane][hit2D].size()));
             }
             
             std::sort(bestPlaneVec.begin(), bestPlaneVec.end(), OrderBestPlanes());
             
             size_t bestPlaneIdx = bestPlaneVec[bestPlaneVecIdx].first;
             size_t bestPlaneCnt = bestPlaneVec[bestPlaneVecIdx].second;
             
             if (bestPlaneCnt > 5) continue;
             
             std::vector<const reco::ClusterHit3D*>& hit3DVec = hit2DToHit3DMap[bestPlaneIdx][hit3D->getHits()[bestPlaneIdx]];
             
             float avePosition[3] = {hit3D->getPosition()[0],0.,0.};
             
             for(const auto& tempHit3D : hit3DVec)
             {
                 avePosition[1] += tempHit3D->getPosition()[1];
                 avePosition[2] += tempHit3D->getPosition()[2];
             }
             
             avePosition[1] *= 1./double(hit3DVec.size());
             avePosition[2] *= 1./double(hit3DVec.size());
             
             tempHitPairList.emplace_back(reco::ClusterHit3D(hit3D->getID(),
                                                             hit3D->getStatusBits(),
                                                             avePosition,
                                                             hit3D->getTotalCharge(),
                                                             hit3D->getAvePeakTime(),
                                                             hit3D->getDeltaPeakTime(),
                                                             hit3D->getSigmaPeakTime(),
                                                             hit3D->getHitChiSquare(),
                                                             hit3D->getDocaToAxis(),
                                                             hit3D->getArclenToPoca(),
                                                             hit3D->getHits(),
                                                             hit3D->getHitDelTSigVec(),
                                                             hit3D->getWireIDs()));
             
             tempHitPairListPtr.push_back(&tempHitPairList.back());
             
             hit3DToHit3DMap[tempHitPairListPtr.back()] = hit3D;
         }
         
         for(const auto& pair : hit3DToHit3DMap)
         {
             pair.second->setPosition(pair.first->getPosition());
             pair.second->setStatusBit(reco::ClusterHit3D::SKELETONPOSAVE);
         }
     }
     
     return;
 }

double lar_cluster3d::SkeletonAlg::FindFirstAndLastWires	(	std::vector< const reco::ClusterHit3D * > &	hitVec,
		int	planeToCheck,
		int	referenceWire,
		double	referenceTicks,
		int &	firstWire,
		int &	lastWire
	)		const

private

A function to find the bounding wires in a given view.

Definition at line 49 of file SkeletonAlg.cxx.

References m_maximumDeltaTicks, m_minimumDeltaTicks, max, and min.

Referenced by FindMedialSkeleton().

 {
     // In the simple case the first and last wires are simply the front and back of the input vector
     firstWire = hitVec.front()->getHits()[planeToCheck]->getHit().WireID().Wire;
     lastWire  = hitVec.back()->getHits()[planeToCheck]->getHit().WireID().Wire;
     
     double maxDeltaTicks = referenceTicks - hitVec.front()->getHits()[planeToCheck]->getTimeTicks();
     double minDeltaTicks = referenceTicks - hitVec.back()->getHits()[planeToCheck]->getTimeTicks();
     
     if (minDeltaTicks > maxDeltaTicks) std::swap(maxDeltaTicks, minDeltaTicks);
 
     double bestDeltaTicks = 1000000.;
     
     // Can't have a gap if only one element
     if (hitVec.size() > 1)
     {
         // The issue is that there may be a gap in wires which we need to find.
         // Reset the last wire
         lastWire = firstWire;
         
         // Keep track of all the deltas
         double nextBestDeltaTicks = bestDeltaTicks;
         
         for(const auto& hitPair : hitVec)
         {
             int    curWire    = hitPair->getHits()[planeToCheck]->getHit().WireID().Wire;
             double deltaTicks = referenceTicks - hitPair->getHits()[planeToCheck]->getTimeTicks();
             
             maxDeltaTicks = std::max(maxDeltaTicks, deltaTicks);
             minDeltaTicks = std::min(minDeltaTicks, deltaTicks);
             
             if (bestDeltaTicks > fabs(deltaTicks))
             {
                 // By definition, the "next best" is now the current best
                 nextBestDeltaTicks = bestDeltaTicks;
                 bestDeltaTicks     = fabs(deltaTicks);
             }
             else if (nextBestDeltaTicks > fabs(deltaTicks))
             {
                 nextBestDeltaTicks = fabs(deltaTicks);
             }
         
             // if gap detected, take action depending on where the input wire is
             // But remember we need to be willing to accept a 1 wire gap... for efficiency
             if (fabs(curWire - lastWire) > 2000)
             {
                 if (referenceWire <= lastWire) break;
             
                 // Stepped over gap, reset everything
                 firstWire          = curWire;
                 maxDeltaTicks      = deltaTicks;
                 minDeltaTicks      = deltaTicks;
                 bestDeltaTicks     = fabs(deltaTicks);
                 nextBestDeltaTicks = bestDeltaTicks;
             }
             
             lastWire = curWire;
         }
         
         bestDeltaTicks = nextBestDeltaTicks;
     }
     
     bestDeltaTicks = std::max(std::min(bestDeltaTicks,m_maximumDeltaTicks), m_minimumDeltaTicks);
     
     if (minDeltaTicks * maxDeltaTicks > 0. && bestDeltaTicks > 30.) bestDeltaTicks = 0.;
     
     return bestDeltaTicks;
 }

int lar_cluster3d::SkeletonAlg::FindMedialSkeleton ( reco::HitPairListPtr & hitPairList ) const

This is intended to find the medial skeleton given a list of input hit pairs.

Parameters

hitPairList - input list of pointers to internal Cluster3D 3D hits

Definition at line 159 of file SkeletonAlg.cxx.

References reco::ClusterHit3D::EDGEHIT, FindFirstAndLastWires(), reco::ClusterHit2D::getTimeTicks(), reco::ClusterHit3D::REJECTEDHIT, and reco::ClusterHit3D::SKELETONHIT.

 {
     // Our mission is to try to find the medial skeletion of the input list of hits
     // We define that as the set of hit pairs where the pairs share the same hit in a given direction
     // and the selected medial hit is equal distance from the edges.
     // The first step in trying to do this is to build a map which relates a give 2D hit to an ordered list of hitpairs using it
     typedef std::map<const reco::ClusterHit2D*, std::vector<const reco::ClusterHit3D*> > Hit2DtoHit3DMapU;
     
     Hit2DtoHit3DMapU hit2DToHit3DMap[3];
     
     for(const auto& hitPair : hitPairList)
     {
         // Don't consider points "rejected" earlier
         if (hitPair->bitsAreSet(reco::ClusterHit3D::REJECTEDHIT)) continue;
         
         // Don't consider hitPair's with less than 3 hits
         unsigned status = hitPair->getStatusBits();
         
         if ((status & 0x7) != 0x7) continue;
         
         for(const auto& hit2D : hitPair->getHits())
         {
             size_t plane = hit2D->getHit().WireID().Plane;
             
             hit2DToHit3DMap[plane][hit2D].push_back(hitPair);
         }
     }
     
     // Do an explicit loop through to sort?
     for(size_t idx = 0; idx < 3; idx++)
     {
         for(auto& mapPair : hit2DToHit3DMap[idx])
         {
             size_t numHitPairs = mapPair.second.size();
             int    planeToCheck = (idx + 1) % 3;    // have forgotten reasoning here...
             
             if (numHitPairs > 1) std::sort(mapPair.second.begin(), mapPair.second.end(), OrderHitsAlongWire(planeToCheck));
         }
     }
     
     // Keep a count of the number of skeleton points to be returned
     int nSkeletonPoints(0);
     
     // The idea is go through all the hits again and determine if they could be "skeleton" elements
     for(const auto& hitPair : hitPairList)
     {
         // Don't consider points "rejected" earlier
         if (hitPair->bitsAreSet(reco::ClusterHit3D::REJECTEDHIT)) continue;
         
         // If a hit pair we skip for now
         if ((hitPair->getStatusBits() & 0x7) != 7) continue;
         
         // Hopefully I am not confusing myself here.
         // The goal is to know, for a given 3D hit, how many other 3D hits share the 2D hits that it is made of
         // We want this to be a bit more precise in that we don't count other 3D hits if there is a gap between
         // the current hit and the one we are comparing.
         // How do we do this?
         // We consider each 2D hit comprising the 3D hit in turn... for each 2D hit we have a list of the 3D
         // hits which share this hit. Note that for the 3D hits to be unique, there must be an equal number of 2D
         // hits from the other two views. So, to count "wires" to the boundary, we can simply pick one of the other views
         
         // so the idea is to through each of the views to determine
         // a) the difference in distance to the last hit on this wire (in each direction)
         // b) the number of 3D hits using the 2D hit in the given 3D hit
         size_t numHitPairs[3]    = {0,  0,  0};    // This keeps track of the number of 3D hits a given 2D hit is associated with
         int    deltaWires[3]     = {0,  0,  0};
         double planeDeltaT[3]    = {0., 0., 0.};
         double bestDeltaTicks[3] = {0., 0., 0.};
         int    wireNumByPlane[3] = {int(hitPair->getHits()[0]->getHit().WireID().Wire),
                                     int(hitPair->getHits()[1]->getHit().WireID().Wire),
                                     int(hitPair->getHits()[2]->getHit().WireID().Wire)};
         
         size_t bestPlaneIdx(0);
         
         // Initiate the loop over views
         for(size_t planeIdx = 0; planeIdx < 3; planeIdx++)
         {
             const reco::ClusterHit2D* hit2D          = hitPair->getHits()[planeIdx];
             double                    hit2DTimeTicks = hit2D->getTimeTicks();
             
             std::vector<const reco::ClusterHit3D*>& hitVec(hit2DToHit3DMap[planeIdx][hit2D]);
             
             numHitPairs[planeIdx] = hitVec.size();
             
             if (numHitPairs[planeIdx] > 1)
             {
                 int planeToCheck = (planeIdx + 1) % 3;
                 int firstWire;
                 int lastWire;
                 
                 bestDeltaTicks[planeIdx] = FindFirstAndLastWires(hitVec,
                                                                  planeToCheck,
                                                                  wireNumByPlane[planeToCheck],
                                                                  hit2DTimeTicks,
                                                                  firstWire,
                                                                  lastWire);
                 
                 int deltaFirst = wireNumByPlane[planeToCheck] - firstWire;
                 int deltaLast  = wireNumByPlane[planeToCheck] - lastWire;
                 
                 // If either distance to the first or last wire is zero then we are an edge point
                 if (deltaFirst == 0 || deltaLast == 0) hitPair->setStatusBit(reco::ClusterHit3D::EDGEHIT);
                 
                 deltaWires[planeIdx]  = deltaFirst + deltaLast;
                 numHitPairs[planeIdx] = lastWire - firstWire + 1;
                 planeDeltaT[planeIdx]  = fabs(hit2DTimeTicks - hitPair->getHits()[planeToCheck]->getTimeTicks());
             }
             // Otherwise, by definition, it is both a skeleton point and an edge point
             else hitPair->setStatusBit(reco::ClusterHit3D::SKELETONHIT | reco::ClusterHit3D::EDGEHIT);
             
             if (numHitPairs[planeIdx] < numHitPairs[bestPlaneIdx])
             {
                 bestPlaneIdx = planeIdx;
             }
         }
         
         // Make sure a real index was found (which should always happen)
         if (bestPlaneIdx < 3 && numHitPairs[bestPlaneIdx] > 0)
         {
             // Make a sliding value for the width of the core region
             int maxBestWires = 0.05 * numHitPairs[bestPlaneIdx] + 1;
             
             maxBestWires = 5000;
             
             // Check condition for a "skeleton" point
             if (fabs(deltaWires[bestPlaneIdx]) < maxBestWires + 1)
             {
                 
                 // If exactly in the middle then we are done
                 // Set a bit in the ClusterHit3D word to signify
 //                if (fabs(deltaWires[bestViewIdx]) < maxBestWires || numHitPairs[bestViewIdx] < 3)
 //                    hitPair->setStatusBit(reco::ClusterHit3D::SKELETONHIT);
 
                 // Otherwise we can try to look at the other view
 //                else
                 {
                     // Find the next best view
                     int nextBestIdx = (bestPlaneIdx + 1) % 3;
                     
                     for(size_t idx = 0; idx < 3; idx++)
                     {
                         if (idx == bestPlaneIdx) continue;
                         
                         if (numHitPairs[idx] < numHitPairs[nextBestIdx]) nextBestIdx = idx;
                     }
                     
                     if (nextBestIdx > 2)
                     {
                         std::cout << "***** invalid next best plane: " << nextBestIdx << " *******" << std::endl;
                         continue;
                     }
                     
                     if (planeDeltaT[bestPlaneIdx] < 1.01*bestDeltaTicks[bestPlaneIdx] && planeDeltaT[nextBestIdx] < 6.01*bestDeltaTicks[nextBestIdx])
                         hitPair->setStatusBit(reco::ClusterHit3D::SKELETONHIT);
                 }
             }
         }
         
         // We want to keep count of "pure" skeleton points only
         if (hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONHIT) && !hitPair->bitsAreSet(reco::ClusterHit3D::EDGEHIT)) nSkeletonPoints++;
     }
     
     return nSkeletonPoints;
 }

void lar_cluster3d::SkeletonAlg::GetSkeletonHits	(	const reco::HitPairListPtr &	inputHitList,
		reco::HitPairListPtr &	skeletonHitList
	)		const

Return the skeleton hits from the input list.

note that this presumes the skeleton hits have been found already

Parameters

inputHitList	- input list of pointers to internal Cluster3D 3D hits
skeletonHitList	- output list of skeleton hits

Definition at line 324 of file SkeletonAlg.cxx.

References reco::ClusterHit3D::SKELETONHIT.

Referenced by lar_cluster3d::Cluster3D::findTrackSeeds(), lar_cluster3d::Cluster3D::splitClustersWithHough(), and lar_cluster3d::Cluster3D::splitClustersWithMST().

 {
     for(const auto& hit3D : inputHitList) if (hit3D->bitsAreSet(reco::ClusterHit3D::SKELETONHIT)) skeletonHitList.emplace_back(hit3D);
     
     return;
 }

void lar_cluster3d::SkeletonAlg::reconfigure ( fhicl::ParameterSet const & pset )

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

Definition at line 43 of file SkeletonAlg.cxx.

References fhicl::ParameterSet::get(), m_maximumDeltaTicks, and m_minimumDeltaTicks.

Referenced by lar_cluster3d::Cluster3D::reconfigure(), and SkeletonAlg().

 {
     m_minimumDeltaTicks = pset.get<double>("MinimumDeltaTicks",  0.05);
     m_maximumDeltaTicks = pset.get<double>("MaximumDeltaTicks", 10.0 );
 }

Member Data Documentation

double lar_cluster3d::SkeletonAlg::m_maximumDeltaTicks

private

Definition at line 91 of file SkeletonAlg.h.

Referenced by FindFirstAndLastWires(), and reconfigure().

double lar_cluster3d::SkeletonAlg::m_minimumDeltaTicks

private

Definition at line 90 of file SkeletonAlg.h.

Referenced by FindFirstAndLastWires(), and reconfigure().

fhicl::ParameterSet lar_cluster3d::SkeletonAlg::m_pset

private

Definition at line 93 of file SkeletonAlg.h.

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

larreco/v06_64_02/source/larreco/RecoAlg/Cluster3DAlgs/SkeletonAlg.h
larreco/v06_64_02/source/larreco/RecoAlg/Cluster3DAlgs/SkeletonAlg.cxx

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation