Cluster3D class. More...

#include "SkeletonAlg.h"

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

	~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 28 of file SkeletonAlg.h.

Constructor & Destructor Documentation

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

Constructor.

Parameters

pset

Definition at line 24 of file SkeletonAlg.cxx.

References reconfigure().

   {
     reconfigure(pset);
   }

lar_cluster3d::SkeletonAlg::~SkeletonAlg ( )

Destructor.

Definition at line 31 of file SkeletonAlg.cxx.

31 {}

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 309 of file SkeletonAlg.cxx.

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

   {
     // 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->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->WireID().Plane, hit2DToHit3DMap[hit2D->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]];
 
         Eigen::Vector3f avePosition(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->getOverlapFraction(),
                                                         hit3D->getChargeAsymmetry(),
                                                         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 41 of file SkeletonAlg.cxx.

References m_maximumDeltaTicks, and m_minimumDeltaTicks.

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]->WireID().Wire;
     lastWire = hitVec.back()->getHits()[planeToCheck]->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]->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 140 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->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]->WireID().Wire),
                                int(hitPair->getHits()[1]->WireID().Wire),
                                int(hitPair->getHits()[2]->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 300 of file SkeletonAlg.cxx.

References reco::ClusterHit3D::SKELETONHIT.

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

   {
     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 35 of file SkeletonAlg.cxx.

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

Referenced by 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 85 of file SkeletonAlg.h.

Referenced by FindFirstAndLastWires(), and reconfigure().

double lar_cluster3d::SkeletonAlg::m_minimumDeltaTicks

private

Definition at line 84 of file SkeletonAlg.h.

Referenced by FindFirstAndLastWires(), and reconfigure().

fhicl::ParameterSet lar_cluster3d::SkeletonAlg::m_pset

private

Definition at line 87 of file SkeletonAlg.h.

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

larreco/v09_25_00/source/larreco/RecoAlg/Cluster3DAlgs/SkeletonAlg.h
larreco/v09_25_00/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