lar_cluster3d::ClusterParamsBuilder Class Reference

ClusterParamsBuilder class definiton. More...

#include "ClusterParamsBuilder.h"

Public Member Functions
	ClusterParamsBuilder (fhicl::ParameterSet const &pset)
	Constructor. More...
virtual	~ClusterParamsBuilder ()
	Destructor. More...
void	configure (const fhicl::ParameterSet &)
void	BuildClusterInfo (reco::ClusterParametersList &clusterParametersList) const
	Given the results of running DBScan, format the clusters so that they can be easily transferred back to the larsoft world. More...
void	FillClusterParams (reco::ClusterParameters &, reco::Hit2DToClusterMap &, double minUniqueFrac=0., double maxLostFrac=1.) const
	Fill the cluster parameters (expose to outside world for case of splitting/merging clusters) More...

Private Member Functions
void	removeUsedHitsFromMap (reco::ClusterParameters &, reco::HitPairListPtr &, reco::Hit2DToClusterMap &) const

Private Attributes
size_t	m_clusterMinHits
	Data members to follow. More...
double	m_clusterMinUniqueFraction
double	m_clusterMaxLostFraction
PrincipalComponentsAlg	m_pcaAlg

Detailed Description

ClusterParamsBuilder class definiton.

Definition at line 34 of file ClusterParamsBuilder.h.

Constructor & Destructor Documentation

lar_cluster3d::ClusterParamsBuilder::ClusterParamsBuilder

(

fhicl::ParameterSet const &

pset

)

Constructor.

Parameters

pset

Definition at line 30 of file ClusterParamsBuilder.cxx.

References configure().

                                                                        :
    m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
{
    this->configure(pset);
}

lar_cluster3d::ClusterParamsBuilder::~ClusterParamsBuilder ( )

virtual

Destructor.

Definition at line 38 of file ClusterParamsBuilder.cxx.

{
}

Member Function Documentation

void lar_cluster3d::ClusterParamsBuilder::BuildClusterInfo

(

reco::ClusterParametersList &

clusterParametersList

)

const

Given the results of running DBScan, format the clusters so that they can be easily transferred back to the larsoft world.

Parameters

hitPairClusterMap	map between view and a list of 3D hits
clusterParametersList	a container for our candidate 3D clusters
rejectionFraction	Used for determine "hit purity" when rejecting clusters The last two parameters are passed through to the FillClusterParams method

Given a list of a list of candidate cluster hits, build these out into the intermediate 3D cluster objects to pass to the final stage

Note that this routine will also reject unworthy clusters, in particular those that share too many hits with other clusters. The criteria is that a larger cluster (more hits) will be superior to a smaller one, if the smaller one shares too many hits with the larger it is zapped. *** THIS IS AN AREA FOR CONTINUED STUDY ***

Definition at line 53 of file ClusterParamsBuilder.cxx.

References FillClusterParams(), reco::ClusterParameters::getClusterParams(), reco::ClusterParameters::getFullPCA(), reco::PrincipalComponents::getSvdOK(), m_clusterMaxLostFraction, m_clusterMinHits, and m_clusterMinUniqueFraction.

Referenced by lar_cluster3d::DBScanAlg::Cluster3DHits(), and lar_cluster3d::MinSpanTreeAlg::Cluster3DHits().

{
    // This is a remote possibility but why not check?
    if (!clusterParametersList.empty())
    {
        // We want to order our clusters on by largest (most number hits) to smallest. So, we'll loop through the clusters,
        // weeding out the unwanted ones and keep track of things in a set of "good" clusters which we'll order
        // by cluster size.
        clusterParametersList.sort();
        
        // The smallest clusters are now at the end, drop those off the back that are less than the mininum necessary
        while(!clusterParametersList.empty() && clusterParametersList.back().getHitPairListPtr().size() < m_clusterMinHits) clusterParametersList.pop_back();
        
        // The next step is to build out a mapping of all 2D hits to clusters
        // Keep track of where the hits get distributed...
        reco::Hit2DToClusterMap hit2DToClusterMap;
        
        reco::ClusterParametersList::iterator clusterItr = clusterParametersList.begin();
        
        for(auto& clusterParams : clusterParametersList)
        {
            for(const auto& hit3D : clusterParams.getHitPairListPtr())
            {
                for(const auto& hit2D : hit3D->getHits())
                {
                    if (!hit2D) continue;
                    
                    hit2DToClusterMap[hit2D][&clusterParams].insert(hit3D);
                }
            }
        }
        
        // Ok, spin through again to remove ambiguous hits
        //        for(auto& clusterParams : clusterParametersList) PruneAmbiguousHits(clusterParams,hit2DToClusterMap);
        
        // What remains is an order set of clusters, largest first
        // Now go through and obtain cluster parameters
        clusterItr = clusterParametersList.begin();
        
        while(clusterItr != clusterParametersList.end())
        {
            // Dereference for ease...
            reco::ClusterParameters& clusterParams = *clusterItr;
            
            // Do the actual work of filling the parameters
            FillClusterParams(clusterParams, hit2DToClusterMap, m_clusterMinUniqueFraction, m_clusterMaxLostFraction);
            
            // If this cluster is rejected then the parameters will be empty
            if (clusterParams.getClusterParams().empty() || !clusterParams.getFullPCA().getSvdOK())
            {
                clusterItr = clusterParametersList.erase(clusterItr);
            }
            else clusterItr++;
        }
    }
    
    return;
}

void lar_cluster3d::ClusterParamsBuilder::configure

(

const fhicl::ParameterSet &

)

Definition at line 44 of file ClusterParamsBuilder.cxx.

References fhicl::ParameterSet::get(), m_clusterMaxLostFraction, m_clusterMinHits, and m_clusterMinUniqueFraction.

Referenced by ClusterParamsBuilder().

{
    m_clusterMinHits           = pset.get<size_t>("ClusterMinHits",            3     );
    m_clusterMinUniqueFraction = pset.get<double>("ClusterMinUniqueFraction",  0.5   );
    m_clusterMaxLostFraction   = pset.get<double>("ClusterMaxLostFraction",    0.5   );
    
    return;
}

void lar_cluster3d::ClusterParamsBuilder::FillClusterParams	(	reco::ClusterParameters &	clusterParams,
		reco::Hit2DToClusterMap &	hit2DToClusterMap,
		double	minUniqueFrac = 0.,
		double	maxLostFrac = 1.
	)		const

Fill the cluster parameters (expose to outside world for case of splitting/merging clusters)

Parameters

ClusterParameters	The cluster parameters container to be modified
Hit2DToClusterMap	Map to keep track of 2D hit to cluster association
double	minimum fraction of unique hits
double	maximum fraction of "lost" hits

Given a list of hits fill out the remaining parameters for this cluster and evaluate the candidate's worthiness to achieve stardom in the event display

Definition at line 121 of file ClusterParamsBuilder.cxx.

References reco::ClusterParameters::getClusterParams(), reco::ClusterParameters::getFullPCA(), reco::ClusterParameters::getHit2DToHit3DListMap(), reco::ClusterParameters::getHit3DToEdgeMap(), reco::ClusterParameters::getHitPairListPtr(), reco::ClusterParameters::getSkeletonPCA(), reco::PrincipalComponents::getSvdOK(), art::left(), m_pcaAlg, max, min, lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_3D(), reco::ClusterHit3D::REJECTEDHIT, art::right(), reco::ClusterParameters::UpdateParameters(), and reco::ClusterHit2D::USED.

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

{
    // Recover the HitPairListPtr from the input clusterParams (which will be the
    // only thing that has been provided)
    reco::HitPairListPtr& hitPairVector = clusterParams.getHitPairListPtr();
    
    // To be sure, we should clear the other data members
    clusterParams.getClusterParams().clear();
    clusterParams.getFullPCA() = reco::PrincipalComponents();
    
    // A test of the emergency broadcast system...
    //    FindBestPathInCluster(clusterParams);
    //    CheckHitSorting(clusterParams);
    
    // See if we can avoid duplicates by temporarily transferring to a set
    std::set<const reco::ClusterHit2D*> hitSet;
    
    // Ultimately we want to keep track of the number of unique 2D hits in this cluster
    // So use a vector (by plane) of sets of hits
    std::vector<size_t> planeHit2DVec;
    std::vector<size_t> planeUniqueHit2DVec;

    planeHit2DVec.resize(3);
    planeUniqueHit2DVec.resize(3);
    
    // Map from 2D hits to associated 3D hits
    reco::Hit2DToHit3DListMap& hit2DToHit3DListMap = clusterParams.getHit2DToHit3DListMap();

    // The map from 2D to 3D hits will contain unique entries for 2D hits so we can do some quick accounting here
    for(const auto& hitMapPair : hit2DToHit3DListMap)
    {
        size_t plane = hitMapPair.first->getHit().WireID().Plane;
        
        planeHit2DVec[plane] += hitMapPair.second.size();
        if (!(hitMapPair.first->getStatusBits() & reco::ClusterHit2D::USED)) planeUniqueHit2DVec[plane] += hitMapPair.second.size();
    }

    // Get totals
    int numTotalHits(0);
    int numUniqueHits(0);
    
    // Also consider the number of hits shared on a given view...
    std::vector<float> uniqueHitFracVec(3,0.);
    int                nPlanesWithHits(0);
    int                nPlanesWithUniqueHits(0);
    size_t             minPlane(0);
    size_t             minPlaneCnt = planeUniqueHit2DVec.at(0);
    
    // Loop through the planes
    for(int idx = 0; idx < 3; idx++)
    {
        // numerology
        numTotalHits  += planeHit2DVec.at(idx);
        numUniqueHits += planeUniqueHit2DVec.at(idx);
        
        if (planeHit2DVec.at(idx)       > 0) nPlanesWithHits++;
        if (planeUniqueHit2DVec.at(idx) > 0) nPlanesWithUniqueHits++;
        
        // Compute the fraction of unique hits in this plane
        uniqueHitFracVec[idx] = float(planeUniqueHit2DVec.at(idx)) / std::max(float(planeHit2DVec.at(idx)),float(1.));
        
        // Finding the plane with the fewest hits
        if (planeHit2DVec.at(idx) < minPlaneCnt)
        {
            minPlaneCnt = planeHit2DVec.at(idx);
            minPlane    = idx;
        }
    }

    // If we have something left then at this point we make one more check
    // This check is intended to weed out clusters made from isolated groups of ambiguous hits which
    // really belong to a larger cluster
    if (numUniqueHits > 0.25 * numTotalHits && nPlanesWithHits > 1 && nPlanesWithUniqueHits > 1)
    {
        // Sorts lowest to highest
        std::sort(uniqueHitFracVec.begin(),uniqueHitFracVec.end());
        
        float acceptRatio = 0.;

        if(uniqueHitFracVec[0] * uniqueHitFracVec[1] > 0.25) acceptRatio = 1.;
        
        float uniqueFraction = uniqueHitFracVec[0] * uniqueHitFracVec[1] * uniqueHitFracVec[2];
        
        // Arbitrary rejection criteria... need to understand
        if (uniqueFraction > 0.6 && acceptRatio > 0.)
        {
            // Create a list to hold 3D hits which are already in use (criteria below)
            reco::HitPairListPtr usedHitPairList;
            
//            std::cout << "--------> Starting 3D hit removal, # 2D hits/plane: " << planeHit2DVec[0] << "/" << planeHit2DVec[1] << "/" << planeHit2DVec[2] << std::endl;
            
            // Have survived laugh test, do final processing...
            // In this first loop go through all the 2D hits and identify the 3D hits that are candidates for deletion
            for(auto& pair : hit2DToHit3DListMap)
            {
                // Check to see if this can happen
                if (pair.second.empty())
                {
                    std::cout << "<<<<<< no matching 3D hits for reco hit in final hit processing >>>>>>" << std::endl;
                    continue;
                }
                
                // Which plane for this hit?
                size_t hitPlane = pair.first->getHit().WireID().Plane;
                
                // Only reject hits on the planes not the fewest 2D hits and really only do this if more than a couple
                if (hitPlane != minPlane && pair.second.size() > 2)
                {
                    // If this hit is associated to a number of 3D hits then do some arbitration
                    // Start by sorting the 3D hits by "significance"
                    // --> Really should do this by the significance of adding the hit we are looking at?
                    //pair.second.sort([hitPlane](const auto& left, const auto& right){return left->getHitDelTSigVec()[hitPlane] < right->getHitDelTSigVec()[hitPlane];});
                    pair.second.sort([](const auto& left, const auto& right){return left->getHitChiSquare() < right->getHitChiSquare();});

                    //std::cout << "~~~~> Checking hit removal, # matches: " << pair.second.size() << ", first params: " << pair.second.front()->getHitChiSquare() << ", last params: "<< pair.second.back()->getHitChiSquare();

                    // From sorted list, determine a rejection value to eliminate bad hits
                    //float cutDeltaTSig = std::min(2.0,std::max(0.5, double((pair.second.front()->getHitDelTSigVec()[hitPlane]))));
                    float cutDeltaTSig = std::min(2.0,std::max(0.5, double(pair.second.front()->getHitChiSquare())));

                    //std::cout << ", cutDeltaTSig: " << cutDeltaTSig;
                    
                    cutDeltaTSig = 10.;
   
                    // And here go through the process of eliminating it
                    //reco::HitPairListPtr::iterator firstBadHitItr = std::find_if(pair.second.begin(),pair.second.end(),[hitPlane,cutDeltaTSig](const auto& hitPtr){return hitPtr->getHitDelTSigVec()[hitPlane] > cutDeltaTSig;});
                    reco::HitPairListPtr::iterator firstBadHitItr = std::find_if(pair.second.begin(),pair.second.end(),[cutDeltaTSig](const auto& hitPtr){return hitPtr->getHitChiSquare() > cutDeltaTSig;});

                    // We need to worry about cutting too many hits... use this loop to try to expand the range in a reasonable fashion
//                    while(std::distance(pair.second.begin(),firstBadHitItr) < int(pair.second.size()/3) && cutDeltaTSig < 0.5)
//                    {
//                        float candDeltaTSig = (*firstBadHitItr)->getHitDelTSigVec()[hitPlane];
//
//                        if (candDeltaTSig > 2. * cutDeltaTSig) break;
//
//                        firstBadHitItr++;
//                        cutDeltaTSig = candDeltaTSig;
//                    }
                
                    reco::HitPairListPtr rejectCandList;
                
                    std::copy(firstBadHitItr,pair.second.end(),std::back_inserter(rejectCandList));
                
                    //std::cout << ", bad hits: " << rejectCandList.size() << std::endl;
                
                    // Remove the 3D hits from all the lists
                    for(const auto& hit3D : rejectCandList)
                    {
                        bool rejectThisHit(true);
                        std::vector<std::pair<reco::HitPairListPtr&,reco::HitPairListPtr::iterator>> deleteVec;
                    
                        for(const auto& hit2D : hit3D->getHits())
                        {
                            // Watch for null hit (dead channels)
                            if (!hit2D) continue;
                        
                            reco::HitPairListPtr& removeHitList = hit2DToHit3DListMap.at(hit2D);

                            // Don't allow all the 3D hits associated to this 2D hit to be rejected?
                            if (removeHitList.size() < 2)
                            {
                                //std::cout << "   ---> remove list too small, size: " << removeHitList.size() << " for hit: " << hit2D << ", pair.first: " << pair.first << std::endl;
                                rejectThisHit = false;
                                break;
                            }
                        
                            reco::HitPairListPtr::iterator removeItr = std::find(removeHitList.begin(),removeHitList.end(),hit3D);
                        
                            if (removeItr != removeHitList.end()) deleteVec.emplace_back(removeHitList,removeItr);
                            //else std::cout << "======>> Did not find 3D hit to remove from list for 2D hit! <<+++++++++" << std::endl;
                        }
                    
                        if (rejectThisHit)
                        {
                            for(auto& rejectPair : deleteVec) rejectPair.first.erase(rejectPair.second);
                        
                            usedHitPairList.push_back(hit3D);
                        }
                    }
                }

                hitSet.insert(pair.first);
            }
            
            // Now we go through the list of candidates and delete those which are unworthy of being processed...
            if (!usedHitPairList.empty())
            {
                // Loop through the hits watching out for double counting
                const reco::ClusterHit3D* lastHit3D = 0;
                
                for(const auto& hit3D : usedHitPairList)
                {
                    if (hit3D == lastHit3D) continue;
                    
                    reco::HitPairListPtr::iterator hit3DItr = std::find(hitPairVector.begin(),hitPairVector.end(),hit3D);
                    
                    if (hit3DItr != hitPairVector.end())
                    {
                        // Mark the hit
                        hit3D->setStatusBit(reco::ClusterHit3D::REJECTEDHIT);
                        
                        // Remove from the cluster's hit container
                        hitPairVector.erase(hit3DItr);
                        
                        // If the clustering algorithm includes edges then need to get rid of those as well
                        if (!clusterParams.getHit3DToEdgeMap().empty())
                        {
                            reco::Hit3DToEdgeMap& edgeMap = clusterParams.getHit3DToEdgeMap();
                            
                            edgeMap.erase(edgeMap.find(hit3D));
                        }
                    }
                    
                    lastHit3D = hit3D;
                }
//                removeUsedHitsFromMap(clusterParams, usedHitPairList, hit2DToClusterMap);
            }

            // First stage of feature extraction runs here
            m_pcaAlg.PCAAnalysis_3D(hitPairVector, clusterParams.getFullPCA());
            
            // Must have a valid pca
            if (clusterParams.getFullPCA().getSvdOK())
            {
                // Set the skeleton PCA to make sure it has some value
                clusterParams.getSkeletonPCA() = clusterParams.getFullPCA();
                
                // Add the "good" hits to our cluster parameters
                for(const auto& hit2D : hitSet)
                {
                    hit2D->setStatusBit(reco::ClusterHit2D::USED);
                    clusterParams.UpdateParameters(hit2D);
                }
            }
        }
    }
    
    return;
}

void lar_cluster3d::ClusterParamsBuilder::removeUsedHitsFromMap ( reco::ClusterParameters &  clusterParams, reco::HitPairListPtr &  usedHitPairList, reco::Hit2DToClusterMap &  hit2DToClusterMap  ) const

private

Definition at line 369 of file ClusterParamsBuilder.cxx.

{
    // Clean up our hit to cluster map
    for(const auto& hit3D : usedHitPairList)
    {
        for(const auto& hit2D : hit3D->getHits())
        {
            if (!hit2D) continue;
            
            reco::Hit2DToClusterMap::iterator hitToClusMapItr = hit2DToClusterMap.find(hit2D);
            
            // I am pretty sure this can't happen but let's check anyway...
            if (hitToClusMapItr == hit2DToClusterMap.end())
            {
                std::cout << "*********** COULD NOT FIND ENTRY FOR 2D HIT! **************" << std::endl;
                break;
            }
            
            // Ok, the same hit can be shared in the same cluster so must be careful here
            // First loop over clusters looking for match
            reco::ClusterToHitPairSetMap::iterator clusToHit3DMapItr = hitToClusMapItr->second.find(&clusterParams);
            
            // This also can't happen
            if (clusToHit3DMapItr == hitToClusMapItr->second.end())
            {
                std::cout << "************ DUCK! THE SKY HAS FALLEN!! *********" << std::endl;
                break;
            }
            
            // If this hit is shared by more than one 3D hit then pick the right one
            if (clusToHit3DMapItr->second.size() > 1)
            {
                reco::HitPairSetPtr::iterator hit3DItr = clusToHit3DMapItr->second.find(hit3D);
                
                clusToHit3DMapItr->second.erase(hit3DItr);
            }
            else
                hitToClusMapItr->second.erase(clusToHit3DMapItr);
            
        }
    }

    return;
}

Member Data Documentation

double lar_cluster3d::ClusterParamsBuilder::m_clusterMaxLostFraction

private

Definition at line 86 of file ClusterParamsBuilder.h.

Referenced by BuildClusterInfo(), and configure().

size_t lar_cluster3d::ClusterParamsBuilder::m_clusterMinHits

private

Data members to follow.

Definition at line 84 of file ClusterParamsBuilder.h.

Referenced by BuildClusterInfo(), and configure().

double lar_cluster3d::ClusterParamsBuilder::m_clusterMinUniqueFraction

private

Definition at line 85 of file ClusterParamsBuilder.h.

Referenced by BuildClusterInfo(), and configure().

PrincipalComponentsAlg lar_cluster3d::ClusterParamsBuilder::m_pcaAlg

private

Definition at line 88 of file ClusterParamsBuilder.h.

Referenced by FillClusterParams().

---

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

• /cvmfs/larsoft.opensciencegrid.org/products/larreco/v07_10_02/source/larreco/RecoAlg/Cluster3DAlgs/ClusterParamsBuilder.h
• /cvmfs/larsoft.opensciencegrid.org/products/larreco/v07_10_02/source/larreco/RecoAlg/Cluster3DAlgs/ClusterParamsBuilder.cxx