#include <SeedFinderAlgorithm.h>

Public Member Functions
	SeedFinderAlgorithm (const fhicl::ParameterSet &pset)

void	reconfigure (fhicl::ParameterSet const &pset)

std::vector< std::vector< recob::Seed > >	GetSeedsFromSortedHits (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< std::vector< art::PtrVector< recob::Hit >>> const &SortedHits, std::vector< std::vector< art::PtrVector< recob::Hit >>> &HitsPerSeed, unsigned int StopAfter=0) const

std::vector< recob::Seed >	GetSeedsFromUnSortedHits (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &, std::vector< art::PtrVector< recob::Hit >> &, unsigned int StopAfter=0) const

SpacePointAlg *	GetSpacePointAlg () const

Private Member Functions
std::vector< recob::Seed >	FindSeeds (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &HitsFlat, std::vector< art::PtrVector< recob::Hit >> &CataloguedHits, unsigned int StopAfter) const

recob::Seed	FindSeedAtEnd (detinfo::DetectorPropertiesData const &detProp, std::vector< recob::SpacePoint > const &, std::vector< char > &, std::vector< int > &, art::PtrVector< recob::Hit > const &HitsFlat, std::vector< std::vector< std::vector< int >>> &OrgHits) const

void	GetCenterAndDirection (detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &HitsFlat, std::vector< int > &HitsToUse, TVector3 &Center, TVector3 &Direction, std::vector< double > &ViewRMS, std::vector< int > &HitsPerView) const

void	ConsolidateSeed (detinfo::DetectorPropertiesData const &detProp, recob::Seed &TheSeed, art::PtrVector< recob::Hit > const &, std::vector< char > &HitStatus, std::vector< std::vector< std::vector< int >>> &OrgHits, bool Extend) const

void	GetHitDistAndProj (detinfo::DetectorPropertiesData const &detProp, recob::Seed const &ASeed, art::Ptr< recob::Hit > const &AHit, double &disp, double &s) const

void	CalculateGeometricalElements ()

Private Attributes
SpacePointAlg *	fSptalg

double	fInitSeedLength

int	fMinPointsInSeed

int	fRefits

std::vector< double >	fMaxViewRMS

float	fHitResolution

float	fOccupancyCut

double	fLengthCut

bool	fExtendSeeds

bool	fAllow2DSeeds

std::vector< double >	fPitches

std::vector< TVector3 >	fPitchDir

std::vector< TVector3 >	fWireDir

std::vector< double >	fWireZeroOffset

TVector3	fXDir

TVector3	fYDir

TVector3	fZDir

size_t	fNChannels

Detailed Description

Definition at line 29 of file SeedFinderAlgorithm.h.

Constructor & Destructor Documentation

trkf::SeedFinderAlgorithm::SeedFinderAlgorithm ( const fhicl::ParameterSet & pset )

Definition at line 27 of file SeedFinderAlgorithm.cxx.

27 { reconfigure(pset); }

trkf::SeedFinderAlgorithm::reconfigure

void reconfigure(fhicl::ParameterSet const &pset)

Definition: SeedFinderAlgorithm.cxx:31

Member Function Documentation

void trkf::SeedFinderAlgorithm::CalculateGeometricalElements ( )

private

Definition at line 930 of file SeedFinderAlgorithm.cxx.

   {
     art::ServiceHandle<geo::Geometry const> geom;
 
     // Total number of channels in the detector
     fNChannels = geom->Nchannels();
 
     // Find pitch of each wireplane
     fPitches.resize(3);
     fPitches.at(0) = fabs(geom->WirePitch(geo::kU));
     fPitches.at(1) = fabs(geom->WirePitch(geo::kV));
     fPitches.at(2) = fabs(geom->WirePitch(geo::kW));
 
     // Setup basis vectors
     fXDir = TVector3(1, 0, 0);
     fYDir = TVector3(0, 1, 0);
     fZDir = TVector3(0, 0, 1);
 
     fWireDir.resize(3);
     fPitchDir.resize(3);
     fWireZeroOffset.resize(3);
 
     double xyzStart1[3], xyzStart2[3];
     double xyzEnd1[3], xyzEnd2[3];
 
     // Calculate wire coordinate systems
     for (size_t n = 0; n != 3; ++n) {
       geom->WireEndPoints(0, 0, n, 0, xyzStart1, xyzEnd1);
       geom->WireEndPoints(0, 0, n, 1, xyzStart2, xyzEnd2);
       fWireDir[n] =
         TVector3(xyzEnd1[0] - xyzStart1[0], xyzEnd1[1] - xyzStart1[1], xyzEnd1[2] - xyzStart1[2])
           .Unit();
       fPitchDir[n] = fWireDir[n].Cross(fXDir).Unit();
       if (fPitchDir[n].Dot(TVector3(
             xyzEnd2[0] - xyzEnd1[0], xyzEnd2[1] - xyzEnd1[1], xyzEnd2[2] - xyzEnd1[2])) < 0)
         fPitchDir[n] = -fPitchDir[n];
 
       fWireZeroOffset[n] =
         xyzEnd1[0] * fPitchDir[n][0] + xyzEnd1[1] * fPitchDir[n][1] + xyzEnd1[2] * fPitchDir[n][2];
     }
   }

void trkf::SeedFinderAlgorithm::ConsolidateSeed	(	detinfo::DetectorPropertiesData const &	detProp,
		recob::Seed &	TheSeed,
		art::PtrVector< recob::Hit > const &	HitsFlat,
		std::vector< char > &	HitStatus,
		std::vector< std::vector< std::vector< int >>> &	OrgHits,
		bool	Extend
	)		const

private

Definition at line 373 of file SeedFinderAlgorithm.cxx.

   {
 
     bool ThrowOutSeed = false;
 
     // This will keep track of what hits are in this seed
     std::map<geo::View_t, std::map<uint32_t, std::vector<int>>> HitsInThisSeed;
 
     int NHitsThisSeed = 0;
 
     double MinS = 1000, MaxS = -1000;
     for (size_t i = 0; i != HitStatus.size(); ++i) {
       if (HitStatus.at(i) == 2) {
         double disp, s;
         GetHitDistAndProj(detProp, TheSeed, HitsFlat.at(i), disp, s);
         if (fabs(s) > 1.2) {
           // This hit is not rightfully part of this seed, toss it.
           HitStatus[i] = 0;
         }
         else {
           NHitsThisSeed++;
 
           if (s < MinS) MinS = s;
           if (s > MaxS) MaxS = s;
           HitsInThisSeed[HitsFlat.at(i)->View()][HitsFlat.at(i)->Channel()].push_back(i);
         }
       }
     }
 
     double LengthRescale = (MaxS - MinS) / 2.;
     double PositionShift = (MaxS + MinS) / 2.;
 
     double pt[3], dir[3], err[3];
     TheSeed.GetPoint(pt, err);
     TheSeed.GetDirection(dir, err);
 
     for (size_t n = 0; n != 3; ++n) {
       pt[n] += dir[n] * PositionShift;
       dir[n] *= LengthRescale;
     }
 
     TheSeed.SetPoint(pt, err);
     TheSeed.SetDirection(dir, err);
 
     // Run through checking if we missed any hits
     for (auto itP = HitsInThisSeed.begin(); itP != HitsInThisSeed.end(); ++itP) {
       double dist, s;
       geo::View_t View = itP->first;
       uint32_t LowestChan = itP->second.begin()->first;
       uint32_t HighestChan = itP->second.rbegin()->first;
       for (uint32_t c = LowestChan; c != HighestChan; ++c) {
         for (size_t h = 0; h != OrgHits[View][c].size(); ++h) {
           if (HitStatus[OrgHits[View][c].at(h)] == 0) {
             GetHitDistAndProj(detProp, TheSeed, HitsFlat[OrgHits[View][c].at(h)], dist, s);
             if (dist < fHitResolution) {
               NHitsThisSeed++;
 
               HitStatus[OrgHits[View][c].at(h)] = 2;
 
               HitsInThisSeed[View][c].push_back(OrgHits[View][c].at(h));
             }
             else
               HitStatus[OrgHits[View][c].at(h)] = 0;
           }
         }
       }
     }
 
     if (NHitsThisSeed == 0) ThrowOutSeed = true;
 
     // Check seed occupancy
 
     uint32_t LowestChanInSeed[3], HighestChanInSeed[3];
     double Occupancy[] = {0., 0., 0.};
     int nHitsPerView[] = {0, 0, 0};
 
     for (auto itP = HitsInThisSeed.begin(); itP != HitsInThisSeed.end(); ++itP) {
 
       geo::View_t View = itP->first;
 
       LowestChanInSeed[View] = itP->second.begin()->first;
       HighestChanInSeed[View] = itP->second.rbegin()->first;
 
       nHitsPerView[View]++;
 
       int FilledChanCount = 0;
 
       for (size_t c = LowestChanInSeed[View]; c != HighestChanInSeed[View]; ++c) {
         if (itP->second[c].size() > 0) ++FilledChanCount;
       }
 
       Occupancy[View] =
         float(FilledChanCount) / float(HighestChanInSeed[View] - LowestChanInSeed[View]);
     }
 
     int nBelowCut(0);
     int nViewsWithHits(0);
     for (size_t n = 0; n != 3; ++n) {
       if (Occupancy[n] < fOccupancyCut) nBelowCut++;
       if (nHitsPerView[n] > 0) nViewsWithHits++;
     }
 
     int belowCut(0);
 
     if (fAllow2DSeeds && nViewsWithHits < 3) belowCut = 1;
 
     if (nBelowCut > belowCut) ThrowOutSeed = true;
 
     if ((Extend) && (!ThrowOutSeed)) {
       std::vector<std::vector<double>> ToAddNegativeS(3, std::vector<double>());
       std::vector<std::vector<double>> ToAddPositiveS(3, std::vector<double>());
       std::vector<std::vector<int>> ToAddNegativeH(3, std::vector<int>());
       std::vector<std::vector<int>> ToAddPositiveH(3, std::vector<int>());
 
       for (auto itP = HitsInThisSeed.begin(); itP != HitsInThisSeed.end(); ++itP) {
         double dist, s;
 
         geo::View_t View = itP->first;
 
         if (LowestChanInSeed[View] > 0) {
           for (uint32_t c = LowestChanInSeed[View] - 1; c != 0; --c) {
             bool GotOneThisChannel = false;
             for (size_t h = 0; h != OrgHits[View][c].size(); ++h) {
               if (HitStatus[OrgHits[View][c][h]] == 0) {
                 GetHitDistAndProj(detProp, TheSeed, HitsFlat[OrgHits[View][c].at(h)], dist, s);
                 if (dist < fHitResolution) {
                   GotOneThisChannel = true;
                   if (s < 0) {
                     ToAddNegativeS[View].push_back(s);
                     ToAddNegativeH[View].push_back(OrgHits[View][c].at(h));
                   }
                   else {
                     ToAddPositiveS[View].push_back(s);
                     ToAddPositiveH[View].push_back(OrgHits[View][c].at(h));
                   }
                 }
               }
             }
             if (GotOneThisChannel == false) break;
           }
         }
         if (HighestChanInSeed[View] < fNChannels)
 
           for (uint32_t c = HighestChanInSeed[View] + 1; c != fNChannels; ++c) {
             bool GotOneThisChannel = false;
             for (size_t h = 0; h != OrgHits[View][c].size(); ++h) {
               if (HitStatus[OrgHits[View][c][h]] == 0) {
                 GetHitDistAndProj(detProp, TheSeed, HitsFlat[OrgHits[View][c].at(h)], dist, s);
                 if (dist < fHitResolution) {
                   GotOneThisChannel = true;
                   if (s < 0) {
 
                     ToAddNegativeS[View].push_back(s);
                     ToAddNegativeH[View].push_back(OrgHits[View][c].at(h));
                   }
                   else {
                     ToAddPositiveS[View].push_back(s);
                     ToAddPositiveH[View].push_back(OrgHits[View][c].at(h));
                   }
                 }
               }
             }
             if (GotOneThisChannel == false) break;
           }
       }
 
       double ExtendPositiveS = 0, ExtendNegativeS = 0;
 
       if ((ToAddPositiveS[0].size() > 0) && (ToAddPositiveS[1].size() > 0) &&
           (ToAddPositiveS[2].size() > 0)) {
         for (size_t n = 0; n != 3; ++n) {
           int n1 = (n + 1) % 3;
           int n2 = (n + 2) % 3;
 
           if ((ToAddPositiveS[n].back() <= ToAddPositiveS[n1].back()) &&
               (ToAddPositiveS[n].back() <= ToAddPositiveS[n2].back())) {
             ExtendPositiveS = ToAddPositiveS[n].back();
           }
         }
       }
 
       if ((ToAddNegativeS[0].size() > 0) && (ToAddNegativeS[1].size() > 0) &&
           (ToAddNegativeS[2].size() > 0)) {
         for (size_t n = 0; n != 3; ++n) {
           int n1 = (n + 1) % 3;
           int n2 = (n + 2) % 3;
           if ((ToAddNegativeS[n].back() >= ToAddNegativeS[n1].back()) &&
               (ToAddNegativeS[n].back() >= ToAddNegativeS[n2].back())) {
             ExtendNegativeS = ToAddNegativeS[n].back();
           }
         }
       }
 
       if (fabs(ExtendNegativeS) < 1.) ExtendNegativeS = -1.;
       if (fabs(ExtendPositiveS) < 1.) ExtendPositiveS = 1.;
 
       LengthRescale = (ExtendPositiveS - ExtendNegativeS) / 2.;
       PositionShift = (ExtendPositiveS + ExtendNegativeS) / 2.;
 
       for (size_t n = 0; n != 3; ++n) {
         pt[n] += dir[n] * PositionShift;
         dir[n] *= LengthRescale;
 
         for (size_t i = 0; i != ToAddPositiveS[n].size(); ++i) {
           if (ToAddPositiveS[n].at(i) < ExtendPositiveS)
             HitStatus[ToAddPositiveH[n].at(i)] = 2;
           else
             HitStatus[ToAddPositiveH[n].at(i)] = 0;
         }
 
         for (size_t i = 0; i != ToAddNegativeS[n].size(); ++i) {
           if (ToAddNegativeS[n].at(i) > ExtendNegativeS)
             HitStatus[ToAddNegativeH[n].at(i)] = 2;
           else
             HitStatus[ToAddNegativeH[n].at(i)] = 0;
         }
       }
 
       TheSeed.SetPoint(pt, err);
       TheSeed.SetDirection(dir, err);
     }
 
     if (ThrowOutSeed) TheSeed.SetValidity(false);
   }

recob::Seed trkf::SeedFinderAlgorithm::FindSeedAtEnd	(	detinfo::DetectorPropertiesData const &	detProp,
		std::vector< recob::SpacePoint > const &	Points,
		std::vector< char > &	PointStatus,
		std::vector< int > &	PointsInRange,
		art::PtrVector< recob::Hit > const &	HitsFlat,
		std::vector< std::vector< std::vector< int >>> &	OrgHits
	)		const

private

Definition at line 646 of file SeedFinderAlgorithm.cxx.

   {
     // This pointer will be returned later
     recob::Seed ReturnSeed;
 
     // Keep track of spacepoints we used, not just their IDs
     std::vector<recob::SpacePoint> PointsUsed;
 
     // Clear output vector
     PointsInRange.clear();
 
     // Loop through hits looking for highest Z seedable point
     TVector3 HighestZPoint;
     bool NoPointFound = true;
     int counter = Points.size() - 1;
     while ((NoPointFound == true) && (counter >= 0)) {
       if (PointStatus[counter] == 0) {
         HighestZPoint = TVector3(
           Points.at(counter).XYZ()[0], Points.at(counter).XYZ()[1], Points.at(counter).XYZ()[2]);
         NoPointFound = false;
       }
       else
         counter--;
     }
     if (NoPointFound) {
       // We didn't find a high point at all
       //  - let the algorithm know to give up.
       PointStatus[0] = 3;
     }
 
     // Now we have the high Z point, loop through collecting
     // near enough hits.  We look 2 seed lengths away, since
     // the seed is bidirectional from the central point
 
     double TwiceLength = 2.0 * fInitSeedLength;
 
     for (int index = Points.size() - 1; index != -1; --index) {
       if (PointStatus[index] == 0) {
         // first check z, then check total distance
         //  (much faster, since most will be out of range in z anyway)
         if ((HighestZPoint[2] - Points.at(index).XYZ()[2]) < TwiceLength) {
           double DistanceToHighZ = pow(pow(HighestZPoint[1] - Points.at(index).XYZ()[1], 2) +
                                          pow(HighestZPoint[2] - Points.at(index).XYZ()[2], 2),
                                        0.5);
           if (DistanceToHighZ < TwiceLength) {
             PointsInRange.push_back(index);
             PointsUsed.push_back(Points.at(index));
           }
         }
         else
           break;
       }
     }
 
     TVector3 SeedCenter(0, 0, 0);
     TVector3 SeedDirection(0, 0, 0);
 
     // Check we have enough points in here to form a seed,
     // otherwise return a dud
     int NPoints = PointsInRange.size();
 
     if (NPoints < fMinPointsInSeed) return recob::Seed();
 
     std::map<int, bool> HitMap;
     std::vector<int> HitList;
 
     for (unsigned int i = 0; i != PointsInRange.size(); i++) {
       std::vector<art::PtrVector<recob::Hit>> HitsInThisCollection(3);
 
       art::PtrVector<recob::Hit> HitsThisSP =
         fSptalg->getAssociatedHits((Points.at(PointsInRange.at(i))));
       for (art::PtrVector<recob::Hit>::const_iterator itHit = HitsThisSP.begin();
            itHit != HitsThisSP.end();
            ++itHit) {
         uint32_t Channel = (*itHit)->Channel();
         geo::View_t View = (*itHit)->View();
 
         double eta = 0.01;
         for (size_t iH = 0; iH != OrgHits[View][Channel].size(); ++iH) {
           if (fabs(HitsFlat[OrgHits[View][Channel][iH]]->PeakTime() - (*itHit)->PeakTime()) < eta) {
             HitMap[OrgHits[View][Channel][iH]] = true;
           }
         }
       }
     }
 
     for (auto itH = HitMap.begin(); itH != HitMap.end(); ++itH) {
       HitList.push_back(itH->first);
     }
 
     std::vector<double> ViewRMS;
     std::vector<int> HitsPerView;
 
     GetCenterAndDirection(
       detProp, HitsFlat, HitList, SeedCenter, SeedDirection, ViewRMS, HitsPerView);
 
     HitMap.clear();
     HitList.clear();
 
     // See if seed points have some linearity
 
     bool ThrowOutSeed = false;
 
     double PtArray[3], DirArray[3];
     double AngleFactor =
       pow(pow(SeedDirection.Y(), 2) + pow(SeedDirection.Z(), 2), 0.5) / SeedDirection.Mag();
 
     int nViewsWithHits(0);
 
     for (size_t n = 0; n != 3; ++n) {
       DirArray[n] = SeedDirection[n] * fInitSeedLength / AngleFactor;
       PtArray[n] = SeedCenter[n];
       if (HitsPerView[n] > 0) nViewsWithHits++;
     }
 
     if (nViewsWithHits < 2 || (nViewsWithHits < 3 && !fAllow2DSeeds)) ThrowOutSeed = true;
 
     ReturnSeed = recob::Seed(PtArray, DirArray);
 
     if (fMaxViewRMS.at(0) > 0) {
 
       for (size_t j = 0; j != fMaxViewRMS.size(); j++) {
         if (fMaxViewRMS.at(j) < ViewRMS.at(j)) { ThrowOutSeed = true; }
         //   mf::LogVerbatim("SeedFinderAlgorithm") << RMS.at(j);
       }
     }
 
     // If the seed is marked as bad, return a dud, otherwise
     //  return the ReturnSeed pointer
     if (!ThrowOutSeed)
       return ReturnSeed;
     else
       return recob::Seed();
   }

std::vector< recob::Seed > trkf::SeedFinderAlgorithm::FindSeeds	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		art::PtrVector< recob::Hit > const &	HitsFlat,
		std::vector< art::PtrVector< recob::Hit >> &	CataloguedHits,
		unsigned int	StopAfter
	)		const

private

Definition at line 60 of file SeedFinderAlgorithm.cxx.

   {
     // Vector of seeds found to return
     std::vector<recob::Seed> ReturnVector;
 
     // First check if there is any overlap
     std::vector<recob::SpacePoint> spts;
 
     fSptalg->clearHitMap();
     fSptalg->makeSpacePoints(clockData, detProp, HitsFlat, spts);
 
     if (int(spts.size()) < fMinPointsInSeed) { return std::vector<recob::Seed>(); }
 
     // If we got here, we have enough spacepoints to potentially make seeds from these hits.
 
     // This is table will let us quickly look up which hits are in a given view / channel.
     //  structure is OrgHits[View][Channel] = {index1,index2...}
     // where the indices are of HitsFlat[index].
 
     std::vector<std::vector<std::vector<int>>> OrgHits(3);
     for (size_t n = 0; n != 3; ++n)
       OrgHits[n].resize(fNChannels);
 
     // These two tables contain the hit to spacepoint mappings
 
     std::vector<std::vector<int>> SpacePointsPerHit(HitsFlat.size(), std::vector<int>());
     std::vector<std::vector<int>> HitsPerSpacePoint(spts.size(), std::vector<int>());
 
     // This vector records the status of each hit.
     // The possible values are:
     //  0. hit unused
     //  1. hit used in a seed
     // Initially none are used.
 
     std::vector<char> HitStatus(HitsFlat.size(), 0);
 
     // Fill the organizing map
 
     for (size_t i = 0; i != HitsFlat.size(); ++i) {
       OrgHits[HitsFlat.at(i)->View()][HitsFlat.at(i)->Channel()].push_back(i);
     }
 
     // Fill the spacepoint / hit lookups
 
     for (size_t iSP = 0; iSP != spts.size(); ++iSP) {
       art::PtrVector<recob::Hit> HitsThisSP = fSptalg->getAssociatedHits(spts.at(iSP));
 
       for (size_t iH = 0; iH != HitsThisSP.size(); ++iH) {
         geo::View_t ThisView = HitsThisSP.at(iH)->View();
         uint32_t ThisChannel = HitsThisSP.at(iH)->Channel();
         float ThisTime = HitsThisSP.at(iH)->PeakTime();
         float eta = 0.001;
 
         for (size_t iOrg = 0; iOrg != OrgHits[ThisView][ThisChannel].size(); ++iOrg) {
           if (fabs(ThisTime - HitsFlat.at(OrgHits[ThisView][ThisChannel][iOrg])->PeakTime()) <
               eta) {
             SpacePointsPerHit.at(OrgHits[ThisView][ThisChannel][iOrg]).push_back(iSP);
             HitsPerSpacePoint.at(iSP).push_back(OrgHits[ThisView][ThisChannel][iOrg]);
           }
         }
       }
     }
 
     // The general idea will be:
     //  A. Make spacepoints from remaining hits
     //  B. Look for 1 seed in that vector
     //  C. Discard used hits
     //  D. Repeat
 
     // This vector keeps track of the status of each space point.
     // The key is the position in the AllSpacePoints vector.
     // The value is
     //    0: point unused
     //    1: point used in seed
     //    2: point thrown but unused
     //    3: flag to terminate seed finding
 
     std::vector<char> PointStatus(spts.size(), 0);
 
     std::vector<std::map<geo::View_t, std::vector<int>>> WhichHitsPerSeed;
 
     bool KeepChopping = true;
 
     while (KeepChopping) {
       // This vector keeps a list of the points used in this seed
       std::vector<int> PointsUsed;
 
       // Find exactly one seed, starting at high Z
       recob::Seed TheSeed =
         FindSeedAtEnd(detProp, spts, PointStatus, PointsUsed, HitsFlat, OrgHits);
 
       // If it was a good seed, collect up the relevant spacepoints
       // and add the seed to the return vector
 
       if (TheSeed.IsValid()) {
 
         for (size_t iP = 0; iP != PointsUsed.size(); ++iP) {
           for (size_t iH = 0; iH != HitsPerSpacePoint.at(PointsUsed.at(iP)).size(); ++iH) {
             int UsedHitID = HitsPerSpacePoint.at(PointsUsed.at(iP)).at(iH);
             HitStatus[UsedHitID] = 2;
           }
         }
         PointStatus[PointsUsed.at(0)] = 1;
         ConsolidateSeed(detProp, TheSeed, HitsFlat, HitStatus, OrgHits, false);
       }
 
       if (TheSeed.IsValid()) {
         if (fRefits > 0) {
           std::vector<char> HitStatusGood;
           recob::Seed SeedGood;
           for (size_t r = 0; r != (unsigned int)fRefits; ++r) {
             double PrevLength = TheSeed.GetLength();
 
             SeedGood = TheSeed;
             HitStatusGood = HitStatus;
 
             std::vector<int> PresentHitList;
             for (size_t iH = 0; iH != HitStatus.size(); ++iH) {
               if (HitStatus[iH] == 2) { PresentHitList.push_back(iH); }
             }
             double pt[3], dir[3], err[3];
 
             TheSeed.GetPoint(pt, err);
             TheSeed.GetDirection(dir, err);
 
             TVector3 Center, Direction;
             std::vector<double> ViewRMS;
             std::vector<int> HitsPerView;
             GetCenterAndDirection(
               detProp, HitsFlat, PresentHitList, Center, Direction, ViewRMS, HitsPerView);
 
             Direction = Direction.Unit() * TheSeed.GetLength();
 
             int nViewsWithHits(0);
             for (size_t n = 0; n != 3; ++n) {
               pt[n] = Center[n];
               dir[n] = Direction[n];
 
               TheSeed.SetPoint(pt, err);
               TheSeed.SetDirection(dir, err);
 
               if (HitsPerView[n] > 0) nViewsWithHits++;
             }
 
             if (nViewsWithHits < 2) TheSeed.SetValidity(false);
 
             if (TheSeed.IsValid())
               ConsolidateSeed(detProp, TheSeed, HitsFlat, HitStatus, OrgHits, fExtendSeeds);
 
             // if we accidentally invalidated the seed, go back to the old one and escape
             else {
               // If we invalidated the seed, go back one interaction
               //  and kill the loop
               HitStatus = HitStatusGood;
               TheSeed = SeedGood;
               break;
             }
             if ((r > 0) && (fabs(PrevLength - TheSeed.GetLength()) < fPitches.at(0))) {
               // If we didn't change much, kill the loop
               break;
             }
           }
         }
       }
 
       if (TheSeed.IsValid()) {
         WhichHitsPerSeed.push_back(std::map<geo::View_t, std::vector<int>>());
 
         art::PtrVector<recob::Hit> HitsWithThisSeed;
         for (size_t iH = 0; iH != HitStatus.size(); ++iH) {
           if (HitStatus.at(iH) == 2) {
             WhichHitsPerSeed.at(WhichHitsPerSeed.size() - 1)[HitsFlat[iH]->View()].push_back(iH);
             HitsWithThisSeed.push_back(HitsFlat.at(iH));
             HitStatus.at(iH) = 1;
 
             for (size_t iSP = 0; iSP != SpacePointsPerHit.at(iH).size(); ++iSP) {
               PointStatus[SpacePointsPerHit.at(iH).at(iSP)] = 1;
             }
           }
         }
 
         // Record that we used this set of hits with this seed in the return catalogue
         ReturnVector.push_back(TheSeed);
         CataloguedHits.push_back(HitsWithThisSeed);
 
         // Tidy up
         HitsWithThisSeed.clear();
       }
       else {
         // If it was not a good seed, throw out the top SP and try again
         PointStatus.at(PointsUsed.at(0)) = 2;
       }
 
       int TotalSPsUsed = 0;
       for (size_t i = 0; i != PointStatus.size(); ++i) {
         if (PointStatus[i] != 0) TotalSPsUsed++;
       }
 
       if ((int(spts.size()) - TotalSPsUsed) < fMinPointsInSeed) KeepChopping = false;
 
       if ((PointStatus[0] == 3) || (PointStatus.size() == 0)) KeepChopping = false;
 
       PointsUsed.clear();
 
       if ((ReturnVector.size() >= StopAfter) && (StopAfter > 0)) break;
 
     } // end spt-level loop
 
     // If we didn't find any seeds, we can make one last ditch attempt. See
     //  if the whole collection is colinear enough to make one megaseed
     //  (good for patching up high angle tracks)
 
     if (ReturnVector.size() == 0) {
       std::vector<int> ListAllHits;
       for (size_t i = 0; i != HitsFlat.size(); ++i) {
         ListAllHits.push_back(i);
         HitStatus[i] = 2;
       }
       TVector3 SeedCenter(0, 0, 0);
       TVector3 SeedDirection(0, 0, 0);
 
       std::vector<double> ViewRMS;
       std::vector<int> HitsPerView;
 
       std::vector<art::PtrVector<recob::Hit>> HitsInThisCollection(3);
 
       GetCenterAndDirection(
         detProp, HitsFlat, ListAllHits, SeedCenter, SeedDirection, ViewRMS, HitsPerView);
 
       bool ThrowOutSeed = false;
 
       double PtArray[3], DirArray[3];
       int nViewsWithHits(0);
       for (size_t n = 0; n != 3; ++n) {
         PtArray[n] = SeedCenter[n];
         DirArray[n] = SeedDirection[n];
         if (HitsPerView[n] > 0) nViewsWithHits++;
       }
       recob::Seed TheSeed(PtArray, DirArray);
 
       if (nViewsWithHits < 2 || (nViewsWithHits < 3 && !fAllow2DSeeds)) ThrowOutSeed = true;
 
       if (!ThrowOutSeed) {
         ConsolidateSeed(detProp, TheSeed, HitsFlat, HitStatus, OrgHits, false);
 
         // Now we have consolidated, grab the right
         //  hits to find the RMS and refitted direction
         ListAllHits.clear();
         for (size_t i = 0; i != HitStatus.size(); ++i) {
           if (HitStatus.at(i) == 2) ListAllHits.push_back(i);
         }
         std::vector<int> HitsPerView;
         GetCenterAndDirection(
           detProp, HitsFlat, ListAllHits, SeedCenter, SeedDirection, ViewRMS, HitsPerView);
 
         int nViewsWithHits(0);
         for (size_t n = 0; n != 3; ++n) {
           PtArray[n] = SeedCenter[n];
           DirArray[n] = SeedDirection[n];
           //                ThrowOutSeed = true;
           if (HitsPerView[n] > 0) nViewsWithHits++;
         }
 
         if (nViewsWithHits < 2 || (nViewsWithHits < 3 && !fAllow2DSeeds)) ThrowOutSeed = true;
 
         TheSeed = recob::Seed(PtArray, DirArray);
 
         if (fMaxViewRMS.at(0) > 0) {
           for (size_t j = 0; j != fMaxViewRMS.size(); j++) {
             if (fMaxViewRMS.at(j) < ViewRMS.at(j)) { ThrowOutSeed = true; }
           }
         }
       }
 
       if ((!ThrowOutSeed) && (TheSeed.IsValid())) {
         ReturnVector.push_back(TheSeed);
         art::PtrVector<recob::Hit> HitsThisSeed;
         for (size_t i = 0; i != ListAllHits.size(); ++i) {
           HitsThisSeed.push_back(HitsFlat.at(ListAllHits.at(i)));
         }
         CataloguedHits.push_back(HitsThisSeed);
       }
     }
 
     // Tidy up
     SpacePointsPerHit.clear();
     HitsPerSpacePoint.clear();
     PointStatus.clear();
     OrgHits.clear();
     HitStatus.clear();
 
     for (size_t i = 0; i != ReturnVector.size(); ++i) {
       double CrazyValue = 1000000;
       double Length = ReturnVector.at(i).GetLength();
       if (!((Length > fLengthCut) && (Length < CrazyValue))) {
         ReturnVector.erase(ReturnVector.begin() + i);
         CataloguedHits.erase(CataloguedHits.begin() + i);
         --i;
       }
     }
 
     return ReturnVector;
   }

void trkf::SeedFinderAlgorithm::GetCenterAndDirection	(	detinfo::DetectorPropertiesData const &	detProp,
		art::PtrVector< recob::Hit > const &	HitsFlat,
		std::vector< int > &	HitsToUse,
		TVector3 &	Center,
		TVector3 &	Direction,
		std::vector< double > &	ViewRMS,
		std::vector< int > &	HitsPerView
	)		const

private

Definition at line 789 of file SeedFinderAlgorithm.cxx.

   {
     N.resize(3);
 
     std::map<uint32_t, bool> HitsClaimed;
 
     // We'll store hit coordinates in each view into this vector
 
     std::vector<std::vector<double>> HitTimes(3);
     std::vector<std::vector<double>> HitWires(3);
     std::vector<std::vector<double>> HitWidths(3);
     std::vector<double> MeanWireCoord(3, 0);
     std::vector<double> MeanTimeCoord(3, 0);
 
     // Run through the collection getting hit info for these spacepoints
 
     std::vector<double> x(3, 0), y(3, 0), xx(3, 0), xy(3, 0), yy(3, 0), sig(3, 0);
 
     for (size_t i = 0; i != HitsToUse.size(); ++i) {
       auto itHit = HitsFlat.begin() + HitsToUse[i];
 
       size_t ViewIndex;
 
       auto const hitView = (*itHit)->View();
       if (hitView == geo::kU)
         ViewIndex = 0;
       else if (hitView == geo::kV)
         ViewIndex = 1;
       else if (hitView == geo::kW)
         ViewIndex = 2;
       else {
         throw art::Exception(art::errors::LogicError)
           << "SpacePointAlg does not support view " << geo::PlaneGeo::ViewName(hitView) << " (#"
           << hitView << ")\n";
       }
       double WireCoord = (*itHit)->WireID().Wire * fPitches.at(ViewIndex);
       double TimeCoord = detProp.ConvertTicksToX((*itHit)->PeakTime(), ViewIndex, 0, 0);
       double TimeUpper = detProp.ConvertTicksToX((*itHit)->PeakTimePlusRMS(), ViewIndex, 0, 0);
       double TimeLower = detProp.ConvertTicksToX((*itHit)->PeakTimeMinusRMS(), ViewIndex, 0, 0);
       double Width = fabs(0.5 * (TimeUpper - TimeLower));
       double Width2 = pow(Width, 2);
 
       HitWires.at(ViewIndex).push_back(WireCoord);
       HitTimes.at(ViewIndex).push_back(TimeCoord);
       HitWidths.at(ViewIndex).push_back(fabs(0.5 * (TimeUpper - TimeLower)));
 
       MeanWireCoord.at(ViewIndex) += WireCoord;
       MeanTimeCoord.at(ViewIndex) += TimeCoord;
 
       // Elements for LS
       x.at(ViewIndex) += WireCoord / Width2;
       y.at(ViewIndex) += TimeCoord / Width2;
       xy.at(ViewIndex) += (TimeCoord * WireCoord) / Width2;
       xx.at(ViewIndex) += (WireCoord * WireCoord) / Width2;
       yy.at(ViewIndex) += (TimeCoord * TimeCoord) / Width2;
       sig.at(ViewIndex) += 1. / Width2;
       N.at(ViewIndex)++;
     }
 
     ViewRMS.clear();
     ViewRMS.resize(3);
     std::vector<double> ViewGrad(3);
     std::vector<double> ViewOffset(3);
 
     for (size_t n = 0; n != 3; ++n) {
       MeanWireCoord[n] /= N[n];
       MeanTimeCoord[n] /= N[n];
 
       double BigN = 1000000;
       double SmallN = 1. / BigN;
 
       if (N[n] > 2) {
         double Numerator = (y[n] / sig[n] - xy[n] / x[n]);
         double Denominator = (x[n] / sig[n] - xx[n] / x[n]);
         if (fabs(Denominator) > SmallN)
           ViewGrad.at(n) = Numerator / Denominator;
         else
           ViewGrad[n] = BigN;
       }
       else if (N[n] == 2)
         ViewGrad[n] = xy[n] / xx[n];
       else
         ViewGrad[n] = BigN;
 
       ViewOffset.at(n) = (y[n] - ViewGrad[n] * x[n]) / sig[n];
       ViewRMS.at(n) = pow((yy[n] + pow(ViewGrad[n], 2) * xx[n] + pow(ViewOffset[n], 2) * sig[n] -
                            2 * ViewGrad[n] * xy[n] - 2 * ViewOffset[n] * y[n] +
                            2 * ViewGrad[n] * ViewOffset[n] * x[n]) /
                             N[n],
                           0.5);
       // Make RMS rotation perp to track
       if (ViewGrad.at(n) != 0) ViewRMS[n] *= sin(atan(1. / ViewGrad.at(n)));
     }
 
     for (size_t n = 0; n != 3; ++n) {
       size_t n1 = (n + 1) % 3;
       size_t n2 = (n + 2) % 3;
       if ((N[n] <= N[n1]) && (N[n] <= N[n2])) {
 
         if (N[n1] < N[n2]) { std::swap(n1, n2); }
         if ((N[n1] == 0) || (N[n2] == 0)) continue;
 
         Direction =
           (fXDir + fPitchDir[n1] * (1. / ViewGrad[n1]) +
            fWireDir[n1] *
              (((1. / ViewGrad[n2]) - fPitchDir[n1].Dot(fPitchDir[n2]) * (1. / ViewGrad[n1])) /
               fWireDir[n1].Dot(fPitchDir[n2])))
             .Unit();
 
         /*
                  Center2D[n] =
                  fXDir * 0.5 * (MeanTimeCoord[n1]+MeanTimeCoord[n2])
                  + fPitchDir[n1] * (MeanWireCoord[n1] + fWireZeroOffset[n1])
                  + fWireDir[n1] *  ( ((MeanWireCoord[n2] + fWireZeroOffset[n2]) - ( MeanWireCoord[n1] + fWireZeroOffset[n1] )*fPitchDir[n1].Dot(fPitchDir[n2]))/(fPitchDir[n2].Dot(fWireDir[n1])) );
                  */
 
         double TimeCoord = 0.5 * (MeanTimeCoord[n1] + MeanTimeCoord[n2]);
         double WireCoordIn1 = (TimeCoord - ViewOffset[n1]) / ViewGrad[n1] + fWireZeroOffset[n1];
         double WireCoordIn2 = (TimeCoord - ViewOffset[n2]) / ViewGrad[n2] + fWireZeroOffset[n2];
 
         Center = fXDir * TimeCoord + fPitchDir[n1] * WireCoordIn1 +
                  fWireDir[n1] * ((WireCoordIn2 - WireCoordIn1 * fPitchDir[n1].Dot(fPitchDir[n2])) /
                                  (fPitchDir[n2].Dot(fWireDir[n1])));
 
         ViewRMS[n] = -fabs(ViewRMS[n]);
         ViewRMS[n1] = fabs(ViewRMS[n1]);
         ViewRMS[n2] = fabs(ViewRMS[n2]);
 
         break;
       }
     }
   }

void trkf::SeedFinderAlgorithm::GetHitDistAndProj	(	detinfo::DetectorPropertiesData const &	detProp,
		recob::Seed const &	ASeed,
		art::Ptr< recob::Hit > const &	AHit,
		double &	disp,
		double &	s
	)		const

private

Definition at line 606 of file SeedFinderAlgorithm.cxx.

   {
     art::ServiceHandle<geo::Geometry const> geom;
 
     double xyzStart[3], xyzEnd[3];
 
     geom->WireEndPoints(0, 0, AHit->WireID().Plane, AHit->WireID().Wire, xyzStart, xyzEnd);
 
     double HitX = detProp.ConvertTicksToX(AHit->PeakTime(), AHit->WireID().Plane, 0, 0);
 
     double HitXHigh = detProp.ConvertTicksToX(AHit->PeakTimePlusRMS(), AHit->WireID().Plane, 0, 0);
     double HitXLow = detProp.ConvertTicksToX(AHit->PeakTimeMinusRMS(), AHit->WireID().Plane, 0, 0);
 
     double HitWidth = HitXHigh - HitXLow;
 
     double pt[3], dir[3], err[3];
 
     ASeed.GetDirection(dir, err);
     ASeed.GetPoint(pt, err);
 
     TVector3 sPt(pt[0], pt[1], pt[2]);
     TVector3 sDir(dir[0], dir[1], dir[2]);
     TVector3 hPt(HitX, xyzStart[1], xyzStart[2]);
     TVector3 hDir(0, xyzStart[1] - xyzEnd[1], xyzStart[2] - xyzEnd[2]);
 
     s = (sPt - hPt).Dot(hDir * (hDir.Dot(sDir)) - sDir * (hDir.Dot(hDir))) /
         (hDir.Dot(hDir) * sDir.Dot(sDir) - pow(hDir.Dot(sDir), 2));
 
     disp = fabs((sPt - hPt).Dot(sDir.Cross(hDir)) / (sDir.Cross(hDir)).Mag()) / HitWidth;
   }

std::vector< std::vector< recob::Seed > > trkf::SeedFinderAlgorithm::GetSeedsFromSortedHits	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		std::vector< std::vector< art::PtrVector< recob::Hit >>> const &	SortedHits,
		std::vector< std::vector< art::PtrVector< recob::Hit >>> &	HitsPerSeed,
		unsigned int	StopAfter = `0`
	)		const

Definition at line 989 of file SeedFinderAlgorithm.cxx.

   {
 
     std::vector<std::vector<recob::Seed>> ReturnVec;
 
     // This piece of code looks detector specific, but its not -
     //   it also works for 2 planes, but one vector is empty.
 
     if (!(fSptalg->enableU() && fSptalg->enableV() && fSptalg->enableW()))
       mf::LogWarning("BezierTrackerModule")
         << "Warning: SpacePointAlg is does not have three views enabled. This may cause unexpected "
            "behaviour in the bezier tracker.";
 
     try {
       for (std::vector<art::PtrVector<recob::Hit>>::const_iterator itU =
              SortedHits.at(geo::kU).begin();
            itU != SortedHits.at(geo::kU).end();
            ++itU)
         for (std::vector<art::PtrVector<recob::Hit>>::const_iterator itV =
                SortedHits.at(geo::kV).begin();
              itV != SortedHits.at(geo::kV).end();
              ++itV)
           for (std::vector<art::PtrVector<recob::Hit>>::const_iterator itW =
                  SortedHits.at(geo::kW).begin();
                itW != SortedHits.at(geo::kW).end();
                ++itW) {
             art::PtrVector<recob::Hit> HitsThisComboFlat;
 
             if (fSptalg->enableU())
               for (size_t i = 0; i != itU->size(); ++i)
                 HitsThisComboFlat.push_back(itU->at(i));
 
             if (fSptalg->enableV())
               for (size_t i = 0; i != itV->size(); ++i)
                 HitsThisComboFlat.push_back(itV->at(i));
 
             if (fSptalg->enableW())
               for (size_t i = 0; i != itW->size(); ++i)
                 HitsThisComboFlat.push_back(itW->at(i));
 
             std::vector<art::PtrVector<recob::Hit>> CataloguedHits;
 
             std::vector<recob::Seed> Seeds =
               FindSeeds(clockData, detProp, HitsThisComboFlat, CataloguedHits, StopAfter);
 
             // Add this harvest to return vectors
             HitsPerSeed.push_back(CataloguedHits);
             ReturnVec.push_back(Seeds);
 
             // Tidy up
             CataloguedHits.clear();
             Seeds.clear();
           }
     }
     catch (...) {
       mf::LogWarning("SeedFinderTrackerModule")
         << " bailed during hit map lookup - have you enabled all 3 planes?";
       ReturnVec.push_back(std::vector<recob::Seed>());
     }
 
     return ReturnVec;
   }

std::vector< recob::Seed > trkf::SeedFinderAlgorithm::GetSeedsFromUnSortedHits	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		art::PtrVector< recob::Hit > const &	Hits,
		std::vector< art::PtrVector< recob::Hit >> &	HitCatalogue,
		unsigned int	StopAfter = `0`
	)		const

Definition at line 975 of file SeedFinderAlgorithm.cxx.

   {
     std::vector<recob::Seed> ReturnVec;
     return FindSeeds(clockData, detProp, Hits, HitCatalogue, StopAfter);
   }

SpacePointAlg* trkf::SeedFinderAlgorithm::GetSpacePointAlg ( ) const

inline

Definition at line 66 of file SeedFinderAlgorithm.h.

     {
       return fSptalg;
     }

void trkf::SeedFinderAlgorithm::reconfigure ( fhicl::ParameterSet const & pset )

Definition at line 31 of file SeedFinderAlgorithm.cxx.

   {
 
     fSptalg = new trkf::SpacePointAlg(pset.get<fhicl::ParameterSet>("SpacePointAlg"));
 
     fInitSeedLength = pset.get<double>("InitSeedLength");
     fMinPointsInSeed = pset.get<int>("MinPointsInSeed");
 
     fRefits = pset.get<double>("Refits");
 
     fHitResolution = pset.get<double>("HitResolution");
 
     fOccupancyCut = pset.get<double>("OccupancyCut");
     fLengthCut = pset.get<double>("LengthCut");
     fExtendSeeds = pset.get<bool>("ExtendSeeds");
 
     fAllow2DSeeds = pset.get<bool>("Allow2DSeeds", false);
 
     fMaxViewRMS.resize(3);
     fMaxViewRMS = pset.get<std::vector<double>>("MaxViewRMS");
 
     CalculateGeometricalElements();
   }

Member Data Documentation

bool trkf::SeedFinderAlgorithm::fAllow2DSeeds

private

Definition at line 142 of file SeedFinderAlgorithm.h.

bool trkf::SeedFinderAlgorithm::fExtendSeeds

private

Definition at line 140 of file SeedFinderAlgorithm.h.

float trkf::SeedFinderAlgorithm::fHitResolution

private

Definition at line 134 of file SeedFinderAlgorithm.h.

double trkf::SeedFinderAlgorithm::fInitSeedLength

private

Definition at line 126 of file SeedFinderAlgorithm.h.

double trkf::SeedFinderAlgorithm::fLengthCut

private

Definition at line 138 of file SeedFinderAlgorithm.h.

std::vector<double> trkf::SeedFinderAlgorithm::fMaxViewRMS

private

Definition at line 132 of file SeedFinderAlgorithm.h.

int trkf::SeedFinderAlgorithm::fMinPointsInSeed

private

Definition at line 128 of file SeedFinderAlgorithm.h.

size_t trkf::SeedFinderAlgorithm::fNChannels

private

Definition at line 149 of file SeedFinderAlgorithm.h.

float trkf::SeedFinderAlgorithm::fOccupancyCut

private

Definition at line 136 of file SeedFinderAlgorithm.h.

std::vector<TVector3> trkf::SeedFinderAlgorithm::fPitchDir

private

Definition at line 145 of file SeedFinderAlgorithm.h.

std::vector<double> trkf::SeedFinderAlgorithm::fPitches

private

Definition at line 144 of file SeedFinderAlgorithm.h.

int trkf::SeedFinderAlgorithm::fRefits

private

Definition at line 130 of file SeedFinderAlgorithm.h.

SpacePointAlg* trkf::SeedFinderAlgorithm::fSptalg

private

Definition at line 124 of file SeedFinderAlgorithm.h.

std::vector<TVector3> trkf::SeedFinderAlgorithm::fWireDir

private

Definition at line 146 of file SeedFinderAlgorithm.h.

std::vector<double> trkf::SeedFinderAlgorithm::fWireZeroOffset

private

Definition at line 147 of file SeedFinderAlgorithm.h.

TVector3 trkf::SeedFinderAlgorithm::fXDir

private

Definition at line 148 of file SeedFinderAlgorithm.h.

TVector3 trkf::SeedFinderAlgorithm::fYDir

private

Definition at line 148 of file SeedFinderAlgorithm.h.

TVector3 trkf::SeedFinderAlgorithm::fZDir

private

Definition at line 148 of file SeedFinderAlgorithm.h.

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

larreco/larreco/RecoAlg/SeedFinderAlgorithm.h
larreco/larreco/RecoAlg/SeedFinderAlgorithm.cxx

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation