lar_cluster3d::PCASeedFinderAlg Class Reference

PCASeedFinderAlg class. More...

#include "PCASeedFinderAlg.h"

Inheritance diagram for lar_cluster3d::PCASeedFinderAlg:

Public Member Functions
	PCASeedFinderAlg (fhicl::ParameterSet const &pset)
	Constructor. More...
virtual	~PCASeedFinderAlg ()
	Destructor. More...
virtual void	reconfigure (fhicl::ParameterSet const &pset)
	a handler for the case where the algorithm control parameters are to be reset More...
virtual bool	findTrackSeeds (reco::HitPairListPtr &hitPairListPtr, reco::PrincipalComponents &inputPCA, SeedHitPairListPairVec &seedHitMap) const
	Given the list of hits this will search for candidate Seed objects and return them. More...

Private Member Functions
bool	getHitsAtEnd (reco::HitPairListPtr &hit3DList, reco::PrincipalComponents &seedPca) const
	Separate function to find hits at the ends of the input hits. More...
void	LineFit2DHits (const reco::HitPairListPtr &hitList, double XOrigin, TVector3 &Pos, TVector3 &Dir, double &ChiDOF) const

Private Attributes
geo::Geometry *	m_geometry
double	m_gapDistance
size_t	m_numSeed2DHits
double	m_minAllowedCosAng
	The minimum cos(ang) between input and seed axes. More...
PrincipalComponentsAlg	m_pcaAlg

Detailed Description

PCASeedFinderAlg class.

Definition at line 36 of file PCASeedFinderAlg.h.

Constructor & Destructor Documentation

lar_cluster3d::PCASeedFinderAlg::PCASeedFinderAlg

(

fhicl::ParameterSet const &

pset

)

Constructor.

Parameters

pset

Definition at line 43 of file PCASeedFinderAlg.cxx.

References m_geometry, and reconfigure().

                                                                :
     m_gapDistance(5.),
     m_numSeed2DHits(80),
     m_minAllowedCosAng(0.7),
     m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
{
    this->reconfigure(pset);
    
    art::ServiceHandle<geo::Geometry>            geometry;
    //    auto const* detectorProperties = lar::providerFrom<detinfo::DetectorPropertiesService>();
    
    m_geometry = &*geometry;
    //    m_detector = detectorProperties->provider();
}

lar_cluster3d::PCASeedFinderAlg::~PCASeedFinderAlg ( )

virtual

Destructor.

Definition at line 60 of file PCASeedFinderAlg.cxx.

{
}

Member Function Documentation

bool lar_cluster3d::PCASeedFinderAlg::findTrackSeeds ( reco::HitPairListPtr &  hitPairListPtr, reco::PrincipalComponents &  inputPCA, SeedHitPairListPairVec &  seedHitMap  ) const

virtual

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

Implements lar_cluster3d::SeedFinderAlgBase.

Definition at line 72 of file PCASeedFinderAlg.cxx.

References reco::PrincipalComponents::flipAxis(), reco::PrincipalComponents::getEigenValues(), reco::PrincipalComponents::getEigenVectors(), getHitsAtEnd(), LineFit2DHits(), m_pcaAlg, and lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_calc3DDocas().

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

{
    bool foundGoodSeed(false);
    
    // 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()[0]);
    double eigenVal1 = 3. * sqrt(pca.getEigenValues()[1]);
    
    if (eigenVal0 > 5. && eigenVal1 > 0.001)
    {
        // Presume CR muons will be "downward going"...
        if (pca.getEigenVectors()[0][1] > 0.) pca.flipAxis(0);
        
        // Use the following to set the 3D doca and arclength for each hit
        m_pcaAlg.PCAAnalysis_calc3DDocas(hitPairListPtr, pca);
        
        // Use this info to sort the hits along the principle axis
        hitPairListPtr.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
        //hitPairListPtr.sort(PcaSort3DHitsByAbsArcLen3D());
        
        // Make a local copy of the 3D hits
        reco::HitPairListPtr hit3DList;
        
        hit3DList.resize(hitPairListPtr.size());
        
        std::copy(hitPairListPtr.begin(), hitPairListPtr.end(), hit3DList.begin());
        
        reco::PrincipalComponents seedPca = pca;
        
        if (getHitsAtEnd(hit3DList, seedPca))
        {
            // We can use the same routine to check hits at the opposite end to make sure
            // we have consistency between both ends of the track.
            // So... follow the same general program as above
            reco::HitPairListPtr checkList;
            
            checkList.resize(hitPairListPtr.size());
            
            std::copy(hitPairListPtr.begin(), hitPairListPtr.end(), checkList.begin());
            
            std::reverse(checkList.begin(), checkList.end());
            
            reco::PrincipalComponents checkPca = pca;
            
            if (getHitsAtEnd(checkList, checkPca))
            {
                // We want our seed to be in the middle of our seed points
                // First step to getting there is to reorder the hits along the
                // new pca axis
                m_pcaAlg.PCAAnalysis_calc3DDocas(hit3DList, seedPca);
            
                // Use this info to sort the hits along the principle axis - note by absolute value of arc length
                hit3DList.sort(SeedFinderAlgBase::Sort3DHitsByAbsArcLen3D());
            
                // Now translate the seedCenter by the arc len to the first hit
                double seedDir[3]   = {seedPca.getEigenVectors()[0][0], seedPca.getEigenVectors()[0][1], seedPca.getEigenVectors()[0][2]};
                double seedStart[3] = {hit3DList.front()->getX(), hit3DList.front()->getY(), hit3DList.front()->getZ()};
                
                if (hit3DList.size() > 10)
                {
                    TVector3 newSeedPos;
                    TVector3 newSeedDir;
                    double   chiDOF;
                
                    LineFit2DHits(hit3DList, 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];
                    }
                }
            
                for(const auto& hit3D : hit3DList) hit3D->setStatusBit(0x40000000);
            
                seedHitPairVec.emplace_back(std::pair<recob::Seed, reco::HitPairListPtr>(recob::Seed(seedStart, seedDir), hit3DList));
                
                foundGoodSeed = true;
            }
        }
    }
    
    return foundGoodSeed;
}

bool lar_cluster3d::PCASeedFinderAlg::getHitsAtEnd ( reco::HitPairListPtr &  hit3DList, reco::PrincipalComponents &  seedPca  ) const

private

Separate function to find hits at the ends of the input hits.

Definition at line 174 of file PCASeedFinderAlg.cxx.

References reco::PrincipalComponents::flipAxis(), reco::ClusterHit3D::getArclenToPoca(), reco::PrincipalComponents::getEigenVectors(), reco::ClusterHit3D::getHits(), reco::PrincipalComponents::getSvdOK(), m_gapDistance, m_minAllowedCosAng, m_numSeed2DHits, m_pcaAlg, and lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_3D().

Referenced by findTrackSeeds().

{
    bool foundGoodSeedHits(false);
    
    // Use a set to count 2D hits by keeping only a single copy
    std::set<const reco::ClusterHit2D*> hit2DSet;
    
    // Look for numSeed2DHits which are continuous
    double lastArcLen = hit3DList.front()->getArclenToPoca();
    
    reco::HitPairListPtr::iterator startItr = hit3DList.begin();
    reco::HitPairListPtr::iterator lastItr  = hit3DList.begin();
    
    while(++lastItr != hit3DList.end())
    {
        const reco::ClusterHit3D* hit3D  = *lastItr;
        double                    arcLen = hit3D->getArclenToPoca();
        
        if (fabs(arcLen-lastArcLen) > m_gapDistance)
        {
            startItr = lastItr;
            hit2DSet.clear();
        }
        
        for(const auto& hit2D : hit3D->getHits())
        {
            hit2DSet.insert(hit2D);
        }
        
        if (hit2DSet.size() > m_numSeed2DHits) break;
        
        lastArcLen = arcLen;
    }
    
    if (hit2DSet.size() > m_numSeed2DHits)
    {
        if (startItr != hit3DList.begin()) hit3DList.erase(hit3DList.begin(), startItr);
        if (lastItr  != hit3DList.end()  ) hit3DList.erase(lastItr,           hit3DList.end());
        
        // On input, the seedPca will contain the original values so we can recover the original axis now
        TVector3 planeVec0(seedPca.getEigenVectors()[0][0],seedPca.getEigenVectors()[0][1],seedPca.getEigenVectors()[0][2]);
        
        m_pcaAlg.PCAAnalysis_3D(hit3DList, seedPca, true);
        
        if (seedPca.getSvdOK())
        {
            // Still looking to point "down"
            if (seedPca.getEigenVectors()[0][1] > 0.) seedPca.flipAxis(0);
            
            // Check that the seed PCA we have found is consistent with the input PCA
            TVector3 primarySeedAxis(seedPca.getEigenVectors()[0][0],seedPca.getEigenVectors()[0][1],seedPca.getEigenVectors()[0][2]);
            
            double cosAng = primarySeedAxis.Dot(planeVec0);
            
            // If the proposed seed axis is not relatively aligned with the input PCA then
            // we should not be using this method to return seeds. Check that here
            if (cosAng > m_minAllowedCosAng) foundGoodSeedHits = true;
        }
    }
    
    return foundGoodSeedHits;
}

void lar_cluster3d::PCASeedFinderAlg::LineFit2DHits ( const reco::HitPairListPtr &  hitList, double  XOrigin, TVector3 &  Pos, TVector3 &  Dir, double &  ChiDOF  ) const

private

Definition at line 238 of file PCASeedFinderAlg.cxx.

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

Referenced by findTrackSeeds().

{
    // 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;
    
    // Make a set of 2D hits so we consider only unique hits
    std::set<const reco::ClusterHit2D*> hit2DSet;
    
    for(const auto& hit3D : hit3DList)
    {
        for(const auto& hit2D : hit3D->getHits()) hit2DSet.insert(hit2D);
    }
    
    if(hit2DSet.size() < 4) return;
    
    const unsigned int nvars = 4;
    unsigned int       npts  = 3 * hit3DList.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), cstat, tpc, ipl;
    double x, cw, sw, off;
    
    // Loop over the 2D hits in the above
    for (const auto& hit : hit2DSet)
    {
        geo::WireID wireID = hit->getHit().WireID();
        
        cstat = wireID.Cryostat;
        tpc   = wireID.TPC;
        ipl   = wireID.Plane;
        
        // get the wire plane offset
        off = m_geometry->WireCoordinate(0, 0, ipl, tpc, cstat);
        
        // get the "cosine-like" component
        cw  = m_geometry->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
        
        // the "sine-like" component
        sw  = m_geometry->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
        
        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;
    
    double ypr, zpr, diff;
    
    for (const auto& hit : hit2DSet)
    {
        geo::WireID wireID = hit->getHit().WireID();
        
        cstat = wireID.Cryostat;
        tpc   = wireID.TPC;
        ipl   = wireID.Plane;
        off   = m_geometry->WireCoordinate(0, 0, ipl, tpc, cstat);
        cw    = m_geometry->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
        sw    = m_geometry->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
        x     = hit->getXPosition() - XOrigin;
        ypr   = tVec[0] + tVec[2] * x;
        zpr   = tVec[1] + tVec[3] * x;
        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()

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

virtual

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

Implements lar_cluster3d::SeedFinderAlgBase.

Definition at line 64 of file PCASeedFinderAlg.cxx.

References fhicl::ParameterSet::get(), m_gapDistance, m_minAllowedCosAng, m_numSeed2DHits, m_pcaAlg, and lar_cluster3d::PrincipalComponentsAlg::reconfigure().

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

{
    m_gapDistance      = pset.get<double>("GapDistance",       5.);
    m_numSeed2DHits    = pset.get<size_t>("NumSeed2DHits",     80);
    m_minAllowedCosAng = pset.get<double>("MinAllowedCosAng", 0.7);
    m_pcaAlg.reconfigure(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"));
}    

Member Data Documentation

double lar_cluster3d::PCASeedFinderAlg::m_gapDistance

private

Definition at line 75 of file PCASeedFinderAlg.h.

Referenced by getHitsAtEnd(), and reconfigure().

geo::Geometry* lar_cluster3d::PCASeedFinderAlg::m_geometry

private

Definition at line 72 of file PCASeedFinderAlg.h.

Referenced by LineFit2DHits(), and PCASeedFinderAlg().

double lar_cluster3d::PCASeedFinderAlg::m_minAllowedCosAng

private

The minimum cos(ang) between input and seed axes.

Definition at line 77 of file PCASeedFinderAlg.h.

Referenced by getHitsAtEnd(), and reconfigure().

size_t lar_cluster3d::PCASeedFinderAlg::m_numSeed2DHits

private

Definition at line 76 of file PCASeedFinderAlg.h.

Referenced by getHitsAtEnd(), and reconfigure().

PrincipalComponentsAlg lar_cluster3d::PCASeedFinderAlg::m_pcaAlg

private

Definition at line 79 of file PCASeedFinderAlg.h.

Referenced by findTrackSeeds(), getHitsAtEnd(), and reconfigure().

---

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

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