#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.

61 {

62 }

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"));
 }    