LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
trkf::SeedFinderAlgorithm Class Reference

#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::SeedGetSeedsFromUnSortedHits (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &, std::vector< art::PtrVector< recob::Hit >> &, unsigned int StopAfter=0) const
 
SpacePointAlgGetSpacePointAlg () const
 

Private Member Functions

std::vector< recob::SeedFindSeeds (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

SpacePointAlgfSptalg
 
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 30 of file SeedFinderAlgorithm.h.

Constructor & Destructor Documentation

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

Definition at line 28 of file SeedFinderAlgorithm.cxx.

References reconfigure().

29  {
30  reconfigure(pset);
31  }
void reconfigure(fhicl::ParameterSet const &pset)

Member Function Documentation

void trkf::SeedFinderAlgorithm::CalculateGeometricalElements ( )
private

Definition at line 931 of file SeedFinderAlgorithm.cxx.

References fNChannels, fPitchDir, fPitches, fWireDir, fWireZeroOffset, fXDir, fYDir, fZDir, geo::kU, geo::kV, geo::kW, n, geo::GeometryCore::Nchannels(), geo::GeometryCore::WireEndPoints(), and geo::GeometryCore::WirePitch().

Referenced by reconfigure().

932  {
934 
935  // Total number of channels in the detector
936  fNChannels = geom->Nchannels();
937 
938  // Find pitch of each wireplane
939  fPitches.resize(3);
940  fPitches.at(0) = fabs(geom->WirePitch(geo::kU));
941  fPitches.at(1) = fabs(geom->WirePitch(geo::kV));
942  fPitches.at(2) = fabs(geom->WirePitch(geo::kW));
943 
944  // Setup basis vectors
945  fXDir = TVector3(1, 0, 0);
946  fYDir = TVector3(0, 1, 0);
947  fZDir = TVector3(0, 0, 1);
948 
949  fWireDir.resize(3);
950  fPitchDir.resize(3);
951  fWireZeroOffset.resize(3);
952 
953  double xyzStart1[3], xyzStart2[3];
954  double xyzEnd1[3], xyzEnd2[3];
955 
956  // Calculate wire coordinate systems
957  for (unsigned int n = 0; n != 3u; ++n) {
958  geom->WireEndPoints(geo::WireID{0, 0, n, 0}, xyzStart1, xyzEnd1);
959  geom->WireEndPoints(geo::WireID{0, 0, n, 1}, xyzStart2, xyzEnd2);
960  fWireDir[n] =
961  TVector3(xyzEnd1[0] - xyzStart1[0], xyzEnd1[1] - xyzStart1[1], xyzEnd1[2] - xyzStart1[2])
962  .Unit();
963  fPitchDir[n] = fWireDir[n].Cross(fXDir).Unit();
964  if (fPitchDir[n].Dot(TVector3(
965  xyzEnd2[0] - xyzEnd1[0], xyzEnd2[1] - xyzEnd1[1], xyzEnd2[2] - xyzEnd1[2])) < 0)
966  fPitchDir[n] = -fPitchDir[n];
967 
968  fWireZeroOffset[n] =
969  xyzEnd1[0] * fPitchDir[n][0] + xyzEnd1[1] * fPitchDir[n][1] + xyzEnd1[2] * fPitchDir[n][2];
970  }
971  }
std::vector< double > fPitches
Planes which measure V.
Definition: geo_types.h:136
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
Planes which measure U.
Definition: geo_types.h:135
std::vector< TVector3 > fWireDir
Char_t n[5]
std::vector< TVector3 > fPitchDir
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:137
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
std::vector< double > fWireZeroOffset
void WireEndPoints(WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
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 375 of file SeedFinderAlgorithm.cxx.

References art::PtrVector< T >::at(), dir, larg4::dist(), fAllow2DSeeds, fHitResolution, fNChannels, fOccupancyCut, recob::Seed::GetDirection(), GetHitDistAndProj(), recob::Seed::GetPoint(), n, pt, recob::Seed::SetDirection(), recob::Seed::SetPoint(), recob::Seed::SetValidity(), and util::size().

Referenced by FindSeeds().

381  {
382 
383  bool ThrowOutSeed = false;
384 
385  // This will keep track of what hits are in this seed
386  std::map<geo::View_t, std::map<uint32_t, std::vector<int>>> HitsInThisSeed;
387 
388  int NHitsThisSeed = 0;
389 
390  double MinS = 1000, MaxS = -1000;
391  for (size_t i = 0; i != HitStatus.size(); ++i) {
392  if (HitStatus.at(i) == 2) {
393  double disp, s;
394  GetHitDistAndProj(detProp, TheSeed, HitsFlat.at(i), disp, s);
395  if (fabs(s) > 1.2) {
396  // This hit is not rightfully part of this seed, toss it.
397  HitStatus[i] = 0;
398  }
399  else {
400  NHitsThisSeed++;
401 
402  if (s < MinS) MinS = s;
403  if (s > MaxS) MaxS = s;
404  HitsInThisSeed[HitsFlat.at(i)->View()][HitsFlat.at(i)->Channel()].push_back(i);
405  }
406  }
407  }
408 
409  double LengthRescale = (MaxS - MinS) / 2.;
410  double PositionShift = (MaxS + MinS) / 2.;
411 
412  double pt[3], dir[3], err[3];
413  TheSeed.GetPoint(pt, err);
414  TheSeed.GetDirection(dir, err);
415 
416  for (size_t n = 0; n != 3; ++n) {
417  pt[n] += dir[n] * PositionShift;
418  dir[n] *= LengthRescale;
419  }
420 
421  TheSeed.SetPoint(pt, err);
422  TheSeed.SetDirection(dir, err);
423 
424  // Run through checking if we missed any hits
425  for (auto itP = HitsInThisSeed.begin(); itP != HitsInThisSeed.end(); ++itP) {
426  double dist, s;
427  geo::View_t View = itP->first;
428  uint32_t LowestChan = itP->second.begin()->first;
429  uint32_t HighestChan = itP->second.rbegin()->first;
430  for (uint32_t c = LowestChan; c != HighestChan; ++c) {
431  for (size_t h = 0; h != OrgHits[View][c].size(); ++h) {
432  if (HitStatus[OrgHits[View][c].at(h)] == 0) {
433  GetHitDistAndProj(detProp, TheSeed, HitsFlat[OrgHits[View][c].at(h)], dist, s);
434  if (dist < fHitResolution) {
435  NHitsThisSeed++;
436 
437  HitStatus[OrgHits[View][c].at(h)] = 2;
438 
439  HitsInThisSeed[View][c].push_back(OrgHits[View][c].at(h));
440  }
441  else
442  HitStatus[OrgHits[View][c].at(h)] = 0;
443  }
444  }
445  }
446  }
447 
448  if (NHitsThisSeed == 0) ThrowOutSeed = true;
449 
450  // Check seed occupancy
451 
452  uint32_t LowestChanInSeed[3], HighestChanInSeed[3];
453  double Occupancy[] = {0., 0., 0.};
454  int nHitsPerView[] = {0, 0, 0};
455 
456  for (auto itP = HitsInThisSeed.begin(); itP != HitsInThisSeed.end(); ++itP) {
457 
458  geo::View_t View = itP->first;
459 
460  LowestChanInSeed[View] = itP->second.begin()->first;
461  HighestChanInSeed[View] = itP->second.rbegin()->first;
462 
463  nHitsPerView[View]++;
464 
465  int FilledChanCount = 0;
466 
467  for (size_t c = LowestChanInSeed[View]; c != HighestChanInSeed[View]; ++c) {
468  if (itP->second[c].size() > 0) ++FilledChanCount;
469  }
470 
471  Occupancy[View] =
472  float(FilledChanCount) / float(HighestChanInSeed[View] - LowestChanInSeed[View]);
473  }
474 
475  int nBelowCut(0);
476  int nViewsWithHits(0);
477  for (size_t n = 0; n != 3; ++n) {
478  if (Occupancy[n] < fOccupancyCut) nBelowCut++;
479  if (nHitsPerView[n] > 0) nViewsWithHits++;
480  }
481 
482  int belowCut(0);
483 
484  if (fAllow2DSeeds && nViewsWithHits < 3) belowCut = 1;
485 
486  if (nBelowCut > belowCut) ThrowOutSeed = true;
487 
488  if ((Extend) && (!ThrowOutSeed)) {
489  std::vector<std::vector<double>> ToAddNegativeS(3, std::vector<double>());
490  std::vector<std::vector<double>> ToAddPositiveS(3, std::vector<double>());
491  std::vector<std::vector<int>> ToAddNegativeH(3, std::vector<int>());
492  std::vector<std::vector<int>> ToAddPositiveH(3, std::vector<int>());
493 
494  for (auto itP = HitsInThisSeed.begin(); itP != HitsInThisSeed.end(); ++itP) {
495  double dist, s;
496 
497  geo::View_t View = itP->first;
498 
499  if (LowestChanInSeed[View] > 0) {
500  for (uint32_t c = LowestChanInSeed[View] - 1; c != 0; --c) {
501  bool GotOneThisChannel = false;
502  for (size_t h = 0; h != OrgHits[View][c].size(); ++h) {
503  if (HitStatus[OrgHits[View][c][h]] == 0) {
504  GetHitDistAndProj(detProp, TheSeed, HitsFlat[OrgHits[View][c].at(h)], dist, s);
505  if (dist < fHitResolution) {
506  GotOneThisChannel = true;
507  if (s < 0) {
508  ToAddNegativeS[View].push_back(s);
509  ToAddNegativeH[View].push_back(OrgHits[View][c].at(h));
510  }
511  else {
512  ToAddPositiveS[View].push_back(s);
513  ToAddPositiveH[View].push_back(OrgHits[View][c].at(h));
514  }
515  }
516  }
517  }
518  if (GotOneThisChannel == false) break;
519  }
520  }
521  if (HighestChanInSeed[View] < fNChannels)
522 
523  for (uint32_t c = HighestChanInSeed[View] + 1; c != fNChannels; ++c) {
524  bool GotOneThisChannel = false;
525  for (size_t h = 0; h != OrgHits[View][c].size(); ++h) {
526  if (HitStatus[OrgHits[View][c][h]] == 0) {
527  GetHitDistAndProj(detProp, TheSeed, HitsFlat[OrgHits[View][c].at(h)], dist, s);
528  if (dist < fHitResolution) {
529  GotOneThisChannel = true;
530  if (s < 0) {
531 
532  ToAddNegativeS[View].push_back(s);
533  ToAddNegativeH[View].push_back(OrgHits[View][c].at(h));
534  }
535  else {
536  ToAddPositiveS[View].push_back(s);
537  ToAddPositiveH[View].push_back(OrgHits[View][c].at(h));
538  }
539  }
540  }
541  }
542  if (GotOneThisChannel == false) break;
543  }
544  }
545 
546  double ExtendPositiveS = 0, ExtendNegativeS = 0;
547 
548  if ((ToAddPositiveS[0].size() > 0) && (ToAddPositiveS[1].size() > 0) &&
549  (ToAddPositiveS[2].size() > 0)) {
550  for (size_t n = 0; n != 3; ++n) {
551  int n1 = (n + 1) % 3;
552  int n2 = (n + 2) % 3;
553 
554  if ((ToAddPositiveS[n].back() <= ToAddPositiveS[n1].back()) &&
555  (ToAddPositiveS[n].back() <= ToAddPositiveS[n2].back())) {
556  ExtendPositiveS = ToAddPositiveS[n].back();
557  }
558  }
559  }
560 
561  if ((ToAddNegativeS[0].size() > 0) && (ToAddNegativeS[1].size() > 0) &&
562  (ToAddNegativeS[2].size() > 0)) {
563  for (size_t n = 0; n != 3; ++n) {
564  int n1 = (n + 1) % 3;
565  int n2 = (n + 2) % 3;
566  if ((ToAddNegativeS[n].back() >= ToAddNegativeS[n1].back()) &&
567  (ToAddNegativeS[n].back() >= ToAddNegativeS[n2].back())) {
568  ExtendNegativeS = ToAddNegativeS[n].back();
569  }
570  }
571  }
572 
573  if (fabs(ExtendNegativeS) < 1.) ExtendNegativeS = -1.;
574  if (fabs(ExtendPositiveS) < 1.) ExtendPositiveS = 1.;
575 
576  LengthRescale = (ExtendPositiveS - ExtendNegativeS) / 2.;
577  PositionShift = (ExtendPositiveS + ExtendNegativeS) / 2.;
578 
579  for (size_t n = 0; n != 3; ++n) {
580  pt[n] += dir[n] * PositionShift;
581  dir[n] *= LengthRescale;
582 
583  for (size_t i = 0; i != ToAddPositiveS[n].size(); ++i) {
584  if (ToAddPositiveS[n].at(i) < ExtendPositiveS)
585  HitStatus[ToAddPositiveH[n].at(i)] = 2;
586  else
587  HitStatus[ToAddPositiveH[n].at(i)] = 0;
588  }
589 
590  for (size_t i = 0; i != ToAddNegativeS[n].size(); ++i) {
591  if (ToAddNegativeS[n].at(i) > ExtendNegativeS)
592  HitStatus[ToAddNegativeH[n].at(i)] = 2;
593  else
594  HitStatus[ToAddNegativeH[n].at(i)] = 0;
595  }
596  }
597 
598  TheSeed.SetPoint(pt, err);
599  TheSeed.SetDirection(dir, err);
600  }
601 
602  if (ThrowOutSeed) TheSeed.SetValidity(false);
603  }
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
void GetPoint(double *Pt, double *Err) const
Definition: Seed.cxx:89
void GetHitDistAndProj(detinfo::DetectorPropertiesData const &detProp, recob::Seed const &ASeed, art::Ptr< recob::Hit > const &AHit, double &disp, double &s) const
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
void SetPoint(double *Pt, double *Err)
Definition: Seed.cxx:122
TMarker * pt
Definition: egs.C:25
reference at(size_type n)
Definition: PtrVector.h:359
TDirectory * dir
Definition: macro.C:5
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Char_t n[5]
void SetValidity(bool Validity)
Definition: Seed.cxx:74
void GetDirection(double *Dir, double *Err) const
Definition: Seed.cxx:80
void SetDirection(double *Dir, double *Err)
Definition: Seed.cxx:112
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 648 of file SeedFinderAlgorithm.cxx.

References art::PtrVector< T >::begin(), util::counter(), art::PtrVector< T >::end(), fAllow2DSeeds, fInitSeedLength, fMaxViewRMS, fMinPointsInSeed, fSptalg, trkf::SpacePointAlg::getAssociatedHits(), GetCenterAndDirection(), and n.

Referenced by FindSeeds().

655  {
656  // This pointer will be returned later
657  recob::Seed ReturnSeed;
658 
659  // Keep track of spacepoints we used, not just their IDs
660  std::vector<recob::SpacePoint> PointsUsed;
661 
662  // Clear output vector
663  PointsInRange.clear();
664 
665  // Loop through hits looking for highest Z seedable point
666  TVector3 HighestZPoint;
667  bool NoPointFound = true;
668  int counter = Points.size() - 1;
669  while ((NoPointFound == true) && (counter >= 0)) {
670  if (PointStatus[counter] == 0) {
671  HighestZPoint = TVector3(
672  Points.at(counter).XYZ()[0], Points.at(counter).XYZ()[1], Points.at(counter).XYZ()[2]);
673  NoPointFound = false;
674  }
675  else
676  counter--;
677  }
678  if (NoPointFound) {
679  // We didn't find a high point at all
680  // - let the algorithm know to give up.
681  PointStatus[0] = 3;
682  }
683 
684  // Now we have the high Z point, loop through collecting
685  // near enough hits. We look 2 seed lengths away, since
686  // the seed is bidirectional from the central point
687 
688  double TwiceLength = 2.0 * fInitSeedLength;
689 
690  for (int index = Points.size() - 1; index != -1; --index) {
691  if (PointStatus[index] == 0) {
692  // first check z, then check total distance
693  // (much faster, since most will be out of range in z anyway)
694  if ((HighestZPoint[2] - Points.at(index).XYZ()[2]) < TwiceLength) {
695  double DistanceToHighZ = pow(pow(HighestZPoint[1] - Points.at(index).XYZ()[1], 2) +
696  pow(HighestZPoint[2] - Points.at(index).XYZ()[2], 2),
697  0.5);
698  if (DistanceToHighZ < TwiceLength) {
699  PointsInRange.push_back(index);
700  PointsUsed.push_back(Points.at(index));
701  }
702  }
703  else
704  break;
705  }
706  }
707 
708  TVector3 SeedCenter(0, 0, 0);
709  TVector3 SeedDirection(0, 0, 0);
710 
711  // Check we have enough points in here to form a seed,
712  // otherwise return a dud
713  int NPoints = PointsInRange.size();
714 
715  if (NPoints < fMinPointsInSeed) return recob::Seed();
716 
717  std::map<int, bool> HitMap;
718  std::vector<int> HitList;
719 
720  for (unsigned int i = 0; i != PointsInRange.size(); i++) {
721  std::vector<art::PtrVector<recob::Hit>> HitsInThisCollection(3);
722 
723  art::PtrVector<recob::Hit> HitsThisSP =
724  fSptalg->getAssociatedHits((Points.at(PointsInRange.at(i))));
725  for (art::PtrVector<recob::Hit>::const_iterator itHit = HitsThisSP.begin();
726  itHit != HitsThisSP.end();
727  ++itHit) {
728  uint32_t Channel = (*itHit)->Channel();
729  geo::View_t View = (*itHit)->View();
730 
731  double eta = 0.01;
732  for (size_t iH = 0; iH != OrgHits[View][Channel].size(); ++iH) {
733  if (fabs(HitsFlat[OrgHits[View][Channel][iH]]->PeakTime() - (*itHit)->PeakTime()) < eta) {
734  HitMap[OrgHits[View][Channel][iH]] = true;
735  }
736  }
737  }
738  }
739 
740  for (auto itH = HitMap.begin(); itH != HitMap.end(); ++itH) {
741  HitList.push_back(itH->first);
742  }
743 
744  std::vector<double> ViewRMS;
745  std::vector<int> HitsPerView;
746 
748  detProp, HitsFlat, HitList, SeedCenter, SeedDirection, ViewRMS, HitsPerView);
749 
750  HitMap.clear();
751  HitList.clear();
752 
753  // See if seed points have some linearity
754 
755  bool ThrowOutSeed = false;
756 
757  double PtArray[3], DirArray[3];
758  double AngleFactor =
759  pow(pow(SeedDirection.Y(), 2) + pow(SeedDirection.Z(), 2), 0.5) / SeedDirection.Mag();
760 
761  int nViewsWithHits(0);
762 
763  for (size_t n = 0; n != 3; ++n) {
764  DirArray[n] = SeedDirection[n] * fInitSeedLength / AngleFactor;
765  PtArray[n] = SeedCenter[n];
766  if (HitsPerView[n] > 0) nViewsWithHits++;
767  }
768 
769  if (nViewsWithHits < 2 || (nViewsWithHits < 3 && !fAllow2DSeeds)) ThrowOutSeed = true;
770 
771  ReturnSeed = recob::Seed(PtArray, DirArray);
772 
773  if (fMaxViewRMS.at(0) > 0) {
774 
775  for (size_t j = 0; j != fMaxViewRMS.size(); j++) {
776  if (fMaxViewRMS.at(j) < ViewRMS.at(j)) { ThrowOutSeed = true; }
777  // mf::LogVerbatim("SeedFinderAlgorithm") << RMS.at(j);
778  }
779  }
780 
781  // If the seed is marked as bad, return a dud, otherwise
782  // return the ReturnSeed pointer
783  if (!ThrowOutSeed)
784  return ReturnSeed;
785  else
786  return recob::Seed();
787  }
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
iterator begin()
Definition: PtrVector.h:217
const art::PtrVector< recob::Hit > & getAssociatedHits(const recob::SpacePoint &spt) const
std::map< int, art::Ptr< recob::Hit > > HitMap
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:295
iterator end()
Definition: PtrVector.h:231
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
Char_t n[5]
std::set< art::Ptr< recob::Hit > > HitList
std::vector< double > fMaxViewRMS
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 62 of file SeedFinderAlgorithm.cxx.

References art::PtrVector< T >::at(), art::PtrVector< T >::clear(), trkf::SpacePointAlg::clearHitMap(), ConsolidateSeed(), dir, fAllow2DSeeds, fExtendSeeds, FindSeedAtEnd(), fLengthCut, fMaxViewRMS, fMinPointsInSeed, fNChannels, fPitches, fRefits, fSptalg, trkf::SpacePointAlg::getAssociatedHits(), GetCenterAndDirection(), recob::Seed::GetDirection(), recob::Seed::GetLength(), recob::Seed::GetPoint(), recob::Seed::IsValid(), tca::Length(), trkf::SpacePointAlg::makeSpacePoints(), n, pt, art::PtrVector< T >::push_back(), r, recob::Seed::SetDirection(), recob::Seed::SetPoint(), recob::Seed::SetValidity(), art::PtrVector< T >::size(), and util::size().

Referenced by GetSeedsFromSortedHits(), and GetSeedsFromUnSortedHits().

68  {
69  // Vector of seeds found to return
70  std::vector<recob::Seed> ReturnVector;
71 
72  // First check if there is any overlap
73  std::vector<recob::SpacePoint> spts;
74 
76  fSptalg->makeSpacePoints(clockData, detProp, HitsFlat, spts);
77 
78  if (int(spts.size()) < fMinPointsInSeed) { return std::vector<recob::Seed>(); }
79 
80  // If we got here, we have enough spacepoints to potentially make seeds from these hits.
81 
82  // This is table will let us quickly look up which hits are in a given view / channel.
83  // structure is OrgHits[View][Channel] = {index1,index2...}
84  // where the indices are of HitsFlat[index].
85 
86  std::vector<std::vector<std::vector<int>>> OrgHits(3);
87  for (size_t n = 0; n != 3; ++n)
88  OrgHits[n].resize(fNChannels);
89 
90  // These two tables contain the hit to spacepoint mappings
91 
92  std::vector<std::vector<int>> SpacePointsPerHit(HitsFlat.size(), std::vector<int>());
93  std::vector<std::vector<int>> HitsPerSpacePoint(spts.size(), std::vector<int>());
94 
95  // This vector records the status of each hit.
96  // The possible values are:
97  // 0. hit unused
98  // 1. hit used in a seed
99  // Initially none are used.
100 
101  std::vector<char> HitStatus(HitsFlat.size(), 0);
102 
103  // Fill the organizing map
104 
105  for (size_t i = 0; i != HitsFlat.size(); ++i) {
106  OrgHits[HitsFlat.at(i)->View()][HitsFlat.at(i)->Channel()].push_back(i);
107  }
108 
109  // Fill the spacepoint / hit lookups
110 
111  for (size_t iSP = 0; iSP != spts.size(); ++iSP) {
112  art::PtrVector<recob::Hit> HitsThisSP = fSptalg->getAssociatedHits(spts.at(iSP));
113 
114  for (size_t iH = 0; iH != HitsThisSP.size(); ++iH) {
115  geo::View_t ThisView = HitsThisSP.at(iH)->View();
116  uint32_t ThisChannel = HitsThisSP.at(iH)->Channel();
117  float ThisTime = HitsThisSP.at(iH)->PeakTime();
118  float eta = 0.001;
119 
120  for (size_t iOrg = 0; iOrg != OrgHits[ThisView][ThisChannel].size(); ++iOrg) {
121  if (fabs(ThisTime - HitsFlat.at(OrgHits[ThisView][ThisChannel][iOrg])->PeakTime()) <
122  eta) {
123  SpacePointsPerHit.at(OrgHits[ThisView][ThisChannel][iOrg]).push_back(iSP);
124  HitsPerSpacePoint.at(iSP).push_back(OrgHits[ThisView][ThisChannel][iOrg]);
125  }
126  }
127  }
128  }
129 
130  // The general idea will be:
131  // A. Make spacepoints from remaining hits
132  // B. Look for 1 seed in that vector
133  // C. Discard used hits
134  // D. Repeat
135 
136  // This vector keeps track of the status of each space point.
137  // The key is the position in the AllSpacePoints vector.
138  // The value is
139  // 0: point unused
140  // 1: point used in seed
141  // 2: point thrown but unused
142  // 3: flag to terminate seed finding
143 
144  std::vector<char> PointStatus(spts.size(), 0);
145 
146  std::vector<std::map<geo::View_t, std::vector<int>>> WhichHitsPerSeed;
147 
148  bool KeepChopping = true;
149 
150  while (KeepChopping) {
151  // This vector keeps a list of the points used in this seed
152  std::vector<int> PointsUsed;
153 
154  // Find exactly one seed, starting at high Z
155  recob::Seed TheSeed =
156  FindSeedAtEnd(detProp, spts, PointStatus, PointsUsed, HitsFlat, OrgHits);
157 
158  // If it was a good seed, collect up the relevant spacepoints
159  // and add the seed to the return vector
160 
161  if (TheSeed.IsValid()) {
162 
163  for (size_t iP = 0; iP != PointsUsed.size(); ++iP) {
164  for (size_t iH = 0; iH != HitsPerSpacePoint.at(PointsUsed.at(iP)).size(); ++iH) {
165  int UsedHitID = HitsPerSpacePoint.at(PointsUsed.at(iP)).at(iH);
166  HitStatus[UsedHitID] = 2;
167  }
168  }
169  PointStatus[PointsUsed.at(0)] = 1;
170  ConsolidateSeed(detProp, TheSeed, HitsFlat, HitStatus, OrgHits, false);
171  }
172 
173  if (TheSeed.IsValid()) {
174  if (fRefits > 0) {
175  std::vector<char> HitStatusGood;
176  recob::Seed SeedGood;
177  for (size_t r = 0; r != (unsigned int)fRefits; ++r) {
178  double PrevLength = TheSeed.GetLength();
179 
180  SeedGood = TheSeed;
181  HitStatusGood = HitStatus;
182 
183  std::vector<int> PresentHitList;
184  for (size_t iH = 0; iH != HitStatus.size(); ++iH) {
185  if (HitStatus[iH] == 2) { PresentHitList.push_back(iH); }
186  }
187  double pt[3], dir[3], err[3];
188 
189  TheSeed.GetPoint(pt, err);
190  TheSeed.GetDirection(dir, err);
191 
192  TVector3 Center, Direction;
193  std::vector<double> ViewRMS;
194  std::vector<int> HitsPerView;
196  detProp, HitsFlat, PresentHitList, Center, Direction, ViewRMS, HitsPerView);
197 
198  Direction = Direction.Unit() * TheSeed.GetLength();
199 
200  int nViewsWithHits(0);
201  for (size_t n = 0; n != 3; ++n) {
202  pt[n] = Center[n];
203  dir[n] = Direction[n];
204 
205  TheSeed.SetPoint(pt, err);
206  TheSeed.SetDirection(dir, err);
207 
208  if (HitsPerView[n] > 0) nViewsWithHits++;
209  }
210 
211  if (nViewsWithHits < 2) TheSeed.SetValidity(false);
212 
213  if (TheSeed.IsValid())
214  ConsolidateSeed(detProp, TheSeed, HitsFlat, HitStatus, OrgHits, fExtendSeeds);
215 
216  // if we accidentally invalidated the seed, go back to the old one and escape
217  else {
218  // If we invalidated the seed, go back one interaction
219  // and kill the loop
220  HitStatus = HitStatusGood;
221  TheSeed = SeedGood;
222  break;
223  }
224  if ((r > 0) && (fabs(PrevLength - TheSeed.GetLength()) < fPitches.at(0))) {
225  // If we didn't change much, kill the loop
226  break;
227  }
228  }
229  }
230  }
231 
232  if (TheSeed.IsValid()) {
233  WhichHitsPerSeed.push_back(std::map<geo::View_t, std::vector<int>>());
234 
235  art::PtrVector<recob::Hit> HitsWithThisSeed;
236  for (size_t iH = 0; iH != HitStatus.size(); ++iH) {
237  if (HitStatus.at(iH) == 2) {
238  WhichHitsPerSeed.at(WhichHitsPerSeed.size() - 1)[HitsFlat[iH]->View()].push_back(iH);
239  HitsWithThisSeed.push_back(HitsFlat.at(iH));
240  HitStatus.at(iH) = 1;
241 
242  for (size_t iSP = 0; iSP != SpacePointsPerHit.at(iH).size(); ++iSP) {
243  PointStatus[SpacePointsPerHit.at(iH).at(iSP)] = 1;
244  }
245  }
246  }
247 
248  // Record that we used this set of hits with this seed in the return catalogue
249  ReturnVector.push_back(TheSeed);
250  CataloguedHits.push_back(HitsWithThisSeed);
251 
252  // Tidy up
253  HitsWithThisSeed.clear();
254  }
255  else {
256  // If it was not a good seed, throw out the top SP and try again
257  PointStatus.at(PointsUsed.at(0)) = 2;
258  }
259 
260  int TotalSPsUsed = 0;
261  for (size_t i = 0; i != PointStatus.size(); ++i) {
262  if (PointStatus[i] != 0) TotalSPsUsed++;
263  }
264 
265  if ((int(spts.size()) - TotalSPsUsed) < fMinPointsInSeed) KeepChopping = false;
266 
267  if ((PointStatus[0] == 3) || (PointStatus.size() == 0)) KeepChopping = false;
268 
269  PointsUsed.clear();
270 
271  if ((ReturnVector.size() >= StopAfter) && (StopAfter > 0)) break;
272 
273  } // end spt-level loop
274 
275  // If we didn't find any seeds, we can make one last ditch attempt. See
276  // if the whole collection is colinear enough to make one megaseed
277  // (good for patching up high angle tracks)
278 
279  if (ReturnVector.size() == 0) {
280  std::vector<int> ListAllHits;
281  for (size_t i = 0; i != HitsFlat.size(); ++i) {
282  ListAllHits.push_back(i);
283  HitStatus[i] = 2;
284  }
285  TVector3 SeedCenter(0, 0, 0);
286  TVector3 SeedDirection(0, 0, 0);
287 
288  std::vector<double> ViewRMS;
289  std::vector<int> HitsPerView;
290 
291  std::vector<art::PtrVector<recob::Hit>> HitsInThisCollection(3);
292 
294  detProp, HitsFlat, ListAllHits, SeedCenter, SeedDirection, ViewRMS, HitsPerView);
295 
296  bool ThrowOutSeed = false;
297 
298  double PtArray[3], DirArray[3];
299  int nViewsWithHits(0);
300  for (size_t n = 0; n != 3; ++n) {
301  PtArray[n] = SeedCenter[n];
302  DirArray[n] = SeedDirection[n];
303  if (HitsPerView[n] > 0) nViewsWithHits++;
304  }
305  recob::Seed TheSeed(PtArray, DirArray);
306 
307  if (nViewsWithHits < 2 || (nViewsWithHits < 3 && !fAllow2DSeeds)) ThrowOutSeed = true;
308 
309  if (!ThrowOutSeed) {
310  ConsolidateSeed(detProp, TheSeed, HitsFlat, HitStatus, OrgHits, false);
311 
312  // Now we have consolidated, grab the right
313  // hits to find the RMS and refitted direction
314  ListAllHits.clear();
315  for (size_t i = 0; i != HitStatus.size(); ++i) {
316  if (HitStatus.at(i) == 2) ListAllHits.push_back(i);
317  }
318  std::vector<int> HitsPerView;
320  detProp, HitsFlat, ListAllHits, SeedCenter, SeedDirection, ViewRMS, HitsPerView);
321 
322  int nViewsWithHits(0);
323  for (size_t n = 0; n != 3; ++n) {
324  PtArray[n] = SeedCenter[n];
325  DirArray[n] = SeedDirection[n];
326  // ThrowOutSeed = true;
327  if (HitsPerView[n] > 0) nViewsWithHits++;
328  }
329 
330  if (nViewsWithHits < 2 || (nViewsWithHits < 3 && !fAllow2DSeeds)) ThrowOutSeed = true;
331 
332  TheSeed = recob::Seed(PtArray, DirArray);
333 
334  if (fMaxViewRMS.at(0) > 0) {
335  for (size_t j = 0; j != fMaxViewRMS.size(); j++) {
336  if (fMaxViewRMS.at(j) < ViewRMS.at(j)) { ThrowOutSeed = true; }
337  }
338  }
339  }
340 
341  if ((!ThrowOutSeed) && (TheSeed.IsValid())) {
342  ReturnVector.push_back(TheSeed);
343  art::PtrVector<recob::Hit> HitsThisSeed;
344  for (size_t i = 0; i != ListAllHits.size(); ++i) {
345  HitsThisSeed.push_back(HitsFlat.at(ListAllHits.at(i)));
346  }
347  CataloguedHits.push_back(HitsThisSeed);
348  }
349  }
350 
351  // Tidy up
352  SpacePointsPerHit.clear();
353  HitsPerSpacePoint.clear();
354  PointStatus.clear();
355  OrgHits.clear();
356  HitStatus.clear();
357 
358  for (size_t i = 0; i != ReturnVector.size(); ++i) {
359  double CrazyValue = 1000000;
360  double Length = ReturnVector.at(i).GetLength();
361  if (!((Length > fLengthCut) && (Length < CrazyValue))) {
362  ReturnVector.erase(ReturnVector.begin() + i);
363  CataloguedHits.erase(CataloguedHits.begin() + i);
364  --i;
365  }
366  }
367 
368  return ReturnVector;
369  }
TRandom r
Definition: spectrum.C:23
std::vector< double > fPitches
float Length(const PFPStruct &pfp)
Definition: PFPUtils.cxx:3269
bool IsValid() const
Definition: Seed.cxx:58
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
void GetPoint(double *Pt, double *Err) const
Definition: Seed.cxx:89
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
const art::PtrVector< recob::Hit > & getAssociatedHits(const recob::SpacePoint &spt) 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
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
void SetPoint(double *Pt, double *Err)
Definition: Seed.cxx:122
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
void makeSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
TMarker * pt
Definition: egs.C:25
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
reference at(size_type n)
Definition: PtrVector.h:359
size_type size() const
Definition: PtrVector.h:302
void clearHitMap() const
TDirectory * dir
Definition: macro.C:5
Char_t n[5]
Direction
Definition: types.h:12
void clear()
Definition: PtrVector.h:533
void SetValidity(bool Validity)
Definition: Seed.cxx:74
void GetDirection(double *Dir, double *Err) const
Definition: Seed.cxx:80
double GetLength() const
Definition: Seed.cxx:132
void SetDirection(double *Dir, double *Err)
Definition: Seed.cxx:112
std::vector< double > fMaxViewRMS
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 791 of file SeedFinderAlgorithm.cxx.

References art::PtrVector< T >::begin(), detinfo::DetectorPropertiesData::ConvertTicksToX(), fPitchDir, fPitches, fWireDir, fWireZeroOffset, fXDir, geo::kU, geo::kV, geo::kW, art::errors::LogicError, n, geo::PlaneGeo::ViewName(), x, xx, and y.

Referenced by FindSeedAtEnd(), and FindSeeds().

798  {
799  N.resize(3);
800 
801  std::map<uint32_t, bool> HitsClaimed;
802 
803  // We'll store hit coordinates in each view into this vector
804 
805  std::vector<std::vector<double>> HitTimes(3);
806  std::vector<std::vector<double>> HitWires(3);
807  std::vector<std::vector<double>> HitWidths(3);
808  std::vector<double> MeanWireCoord(3, 0);
809  std::vector<double> MeanTimeCoord(3, 0);
810 
811  // Run through the collection getting hit info for these spacepoints
812 
813  std::vector<double> x(3, 0), y(3, 0), xx(3, 0), xy(3, 0), yy(3, 0), sig(3, 0);
814 
815  for (size_t i = 0; i != HitsToUse.size(); ++i) {
816  auto itHit = HitsFlat.begin() + HitsToUse[i];
817 
818  size_t ViewIndex;
819 
820  auto const hitView = (*itHit)->View();
821  if (hitView == geo::kU)
822  ViewIndex = 0;
823  else if (hitView == geo::kV)
824  ViewIndex = 1;
825  else if (hitView == geo::kW)
826  ViewIndex = 2;
827  else {
829  << "SpacePointAlg does not support view " << geo::PlaneGeo::ViewName(hitView) << " (#"
830  << hitView << ")\n";
831  }
832  double WireCoord = (*itHit)->WireID().Wire * fPitches.at(ViewIndex);
833  double TimeCoord = detProp.ConvertTicksToX((*itHit)->PeakTime(), ViewIndex, 0, 0);
834  double TimeUpper = detProp.ConvertTicksToX((*itHit)->PeakTimePlusRMS(), ViewIndex, 0, 0);
835  double TimeLower = detProp.ConvertTicksToX((*itHit)->PeakTimeMinusRMS(), ViewIndex, 0, 0);
836  double Width = fabs(0.5 * (TimeUpper - TimeLower));
837  double Width2 = pow(Width, 2);
838 
839  HitWires.at(ViewIndex).push_back(WireCoord);
840  HitTimes.at(ViewIndex).push_back(TimeCoord);
841  HitWidths.at(ViewIndex).push_back(fabs(0.5 * (TimeUpper - TimeLower)));
842 
843  MeanWireCoord.at(ViewIndex) += WireCoord;
844  MeanTimeCoord.at(ViewIndex) += TimeCoord;
845 
846  // Elements for LS
847  x.at(ViewIndex) += WireCoord / Width2;
848  y.at(ViewIndex) += TimeCoord / Width2;
849  xy.at(ViewIndex) += (TimeCoord * WireCoord) / Width2;
850  xx.at(ViewIndex) += (WireCoord * WireCoord) / Width2;
851  yy.at(ViewIndex) += (TimeCoord * TimeCoord) / Width2;
852  sig.at(ViewIndex) += 1. / Width2;
853  N.at(ViewIndex)++;
854  }
855 
856  ViewRMS.clear();
857  ViewRMS.resize(3);
858  std::vector<double> ViewGrad(3);
859  std::vector<double> ViewOffset(3);
860 
861  for (size_t n = 0; n != 3; ++n) {
862  MeanWireCoord[n] /= N[n];
863  MeanTimeCoord[n] /= N[n];
864 
865  double BigN = 1000000;
866  double SmallN = 1. / BigN;
867 
868  if (N[n] > 2) {
869  double Numerator = (y[n] / sig[n] - xy[n] / x[n]);
870  double Denominator = (x[n] / sig[n] - xx[n] / x[n]);
871  if (fabs(Denominator) > SmallN)
872  ViewGrad.at(n) = Numerator / Denominator;
873  else
874  ViewGrad[n] = BigN;
875  }
876  else if (N[n] == 2)
877  ViewGrad[n] = xy[n] / xx[n];
878  else
879  ViewGrad[n] = BigN;
880 
881  ViewOffset.at(n) = (y[n] - ViewGrad[n] * x[n]) / sig[n];
882  ViewRMS.at(n) = pow((yy[n] + pow(ViewGrad[n], 2) * xx[n] + pow(ViewOffset[n], 2) * sig[n] -
883  2 * ViewGrad[n] * xy[n] - 2 * ViewOffset[n] * y[n] +
884  2 * ViewGrad[n] * ViewOffset[n] * x[n]) /
885  N[n],
886  0.5);
887  // Make RMS rotation perp to track
888  if (ViewGrad.at(n) != 0) ViewRMS[n] *= sin(atan(1. / ViewGrad.at(n)));
889  }
890 
891  for (size_t n = 0; n != 3; ++n) {
892  size_t n1 = (n + 1) % 3;
893  size_t n2 = (n + 2) % 3;
894  if ((N[n] <= N[n1]) && (N[n] <= N[n2])) {
895 
896  if (N[n1] < N[n2]) { std::swap(n1, n2); }
897  if ((N[n1] == 0) || (N[n2] == 0)) continue;
898 
899  Direction =
900  (fXDir + fPitchDir[n1] * (1. / ViewGrad[n1]) +
901  fWireDir[n1] *
902  (((1. / ViewGrad[n2]) - fPitchDir[n1].Dot(fPitchDir[n2]) * (1. / ViewGrad[n1])) /
903  fWireDir[n1].Dot(fPitchDir[n2])))
904  .Unit();
905 
906  /*
907  Center2D[n] =
908  fXDir * 0.5 * (MeanTimeCoord[n1]+MeanTimeCoord[n2])
909  + fPitchDir[n1] * (MeanWireCoord[n1] + fWireZeroOffset[n1])
910  + fWireDir[n1] * ( ((MeanWireCoord[n2] + fWireZeroOffset[n2]) - ( MeanWireCoord[n1] + fWireZeroOffset[n1] )*fPitchDir[n1].Dot(fPitchDir[n2]))/(fPitchDir[n2].Dot(fWireDir[n1])) );
911  */
912 
913  double TimeCoord = 0.5 * (MeanTimeCoord[n1] + MeanTimeCoord[n2]);
914  double WireCoordIn1 = (TimeCoord - ViewOffset[n1]) / ViewGrad[n1] + fWireZeroOffset[n1];
915  double WireCoordIn2 = (TimeCoord - ViewOffset[n2]) / ViewGrad[n2] + fWireZeroOffset[n2];
916 
917  Center = fXDir * TimeCoord + fPitchDir[n1] * WireCoordIn1 +
918  fWireDir[n1] * ((WireCoordIn2 - WireCoordIn1 * fPitchDir[n1].Dot(fPitchDir[n2])) /
919  (fPitchDir[n2].Dot(fWireDir[n1])));
920 
921  ViewRMS[n] = -fabs(ViewRMS[n]);
922  ViewRMS[n1] = fabs(ViewRMS[n1]);
923  ViewRMS[n2] = fabs(ViewRMS[n2]);
924 
925  break;
926  }
927  }
928  }
Float_t x
Definition: compare.C:6
std::vector< double > fPitches
static std::string ViewName(geo::View_t view)
Returns the name of the specified view.
Definition: PlaneGeo.cxx:682
Double_t xx
Definition: macro.C:12
Planes which measure V.
Definition: geo_types.h:136
iterator begin()
Definition: PtrVector.h:217
Float_t y
Definition: compare.C:6
Planes which measure U.
Definition: geo_types.h:135
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::vector< TVector3 > fWireDir
void swap(lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &a, lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &b)
Char_t n[5]
Direction
Definition: types.h:12
std::vector< TVector3 > fPitchDir
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:137
std::vector< double > fWireZeroOffset
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 607 of file SeedFinderAlgorithm.cxx.

References detinfo::DetectorPropertiesData::ConvertTicksToX(), dir, recob::Seed::GetDirection(), recob::Seed::GetPoint(), recob::Hit::PeakTime(), recob::Hit::PeakTimeMinusRMS(), recob::Hit::PeakTimePlusRMS(), geo::PlaneID::Plane, pt, geo::GeometryCore::WireEndPoints(), and recob::Hit::WireID().

Referenced by ConsolidateSeed().

612  {
614 
615  double xyzStart[3], xyzEnd[3];
616 
617  constexpr geo::TPCID tpcid{0, 0};
618  geo::PlaneID const planeID{tpcid, AHit->WireID().Plane};
619 
620  geom->WireEndPoints(geo::WireID{planeID, 0}, xyzStart, xyzEnd);
621 
622  double HitX = detProp.ConvertTicksToX(AHit->PeakTime(), planeID);
623  double HitXHigh = detProp.ConvertTicksToX(AHit->PeakTimePlusRMS(), planeID);
624  double HitXLow = detProp.ConvertTicksToX(AHit->PeakTimeMinusRMS(), planeID);
625 
626  double HitWidth = HitXHigh - HitXLow;
627 
628  double pt[3], dir[3], err[3];
629 
630  ASeed.GetDirection(dir, err);
631  ASeed.GetPoint(pt, err);
632 
633  TVector3 sPt(pt[0], pt[1], pt[2]);
634  TVector3 sDir(dir[0], dir[1], dir[2]);
635  TVector3 hPt(HitX, xyzStart[1], xyzStart[2]);
636  TVector3 hDir(0, xyzStart[1] - xyzEnd[1], xyzStart[2] - xyzEnd[2]);
637 
638  s = (sPt - hPt).Dot(hDir * (hDir.Dot(sDir)) - sDir * (hDir.Dot(hDir))) /
639  (hDir.Dot(hDir) * sDir.Dot(sDir) - pow(hDir.Dot(sDir), 2));
640 
641  disp = fabs((sPt - hPt).Dot(sDir.Cross(hDir)) / (sDir.Cross(hDir)).Mag()) / HitWidth;
642  }
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
TMarker * pt
Definition: egs.C:25
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
float PeakTimeMinusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:290
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
TDirectory * dir
Definition: macro.C:5
float PeakTimePlusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:285
void WireEndPoints(WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
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 987 of file SeedFinderAlgorithm.cxx.

References trkf::SpacePointAlg::enableU(), trkf::SpacePointAlg::enableV(), trkf::SpacePointAlg::enableW(), FindSeeds(), fSptalg, geo::kU, geo::kV, geo::kW, art::PtrVector< T >::push_back(), and lar::dump::vector().

Referenced by trkf::SeedFinderModule::produce().

993  {
994 
995  std::vector<std::vector<recob::Seed>> ReturnVec;
996 
997  // This piece of code looks detector specific, but its not -
998  // it also works for 2 planes, but one vector is empty.
999 
1000  if (!(fSptalg->enableU() && fSptalg->enableV() && fSptalg->enableW()))
1001  mf::LogWarning("BezierTrackerModule")
1002  << "Warning: SpacePointAlg is does not have three views enabled. This may cause unexpected "
1003  "behaviour in the bezier tracker.";
1004 
1005  try {
1007  SortedHits.at(geo::kU).begin();
1008  itU != SortedHits.at(geo::kU).end();
1009  ++itU)
1011  SortedHits.at(geo::kV).begin();
1012  itV != SortedHits.at(geo::kV).end();
1013  ++itV)
1015  SortedHits.at(geo::kW).begin();
1016  itW != SortedHits.at(geo::kW).end();
1017  ++itW) {
1018  art::PtrVector<recob::Hit> HitsThisComboFlat;
1019 
1020  if (fSptalg->enableU())
1021  for (size_t i = 0; i != itU->size(); ++i)
1022  HitsThisComboFlat.push_back(itU->at(i));
1023 
1024  if (fSptalg->enableV())
1025  for (size_t i = 0; i != itV->size(); ++i)
1026  HitsThisComboFlat.push_back(itV->at(i));
1027 
1028  if (fSptalg->enableW())
1029  for (size_t i = 0; i != itW->size(); ++i)
1030  HitsThisComboFlat.push_back(itW->at(i));
1031 
1032  std::vector<art::PtrVector<recob::Hit>> CataloguedHits;
1033 
1034  std::vector<recob::Seed> Seeds =
1035  FindSeeds(clockData, detProp, HitsThisComboFlat, CataloguedHits, StopAfter);
1036 
1037  // Add this harvest to return vectors
1038  HitsPerSeed.push_back(CataloguedHits);
1039  ReturnVec.push_back(Seeds);
1040 
1041  // Tidy up
1042  CataloguedHits.clear();
1043  Seeds.clear();
1044  }
1045  }
1046  catch (...) {
1047  mf::LogWarning("SeedFinderTrackerModule")
1048  << " bailed during hit map lookup - have you enabled all 3 planes?";
1049  ReturnVec.push_back(std::vector<recob::Seed>());
1050  }
1051 
1052  return ReturnVec;
1053  }
bool enableW() const noexcept
Definition: SpacePointAlg.h:93
bool enableV() const noexcept
Definition: SpacePointAlg.h:92
Planes which measure V.
Definition: geo_types.h:136
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
Planes which measure U.
Definition: geo_types.h:135
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
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
bool enableU() const noexcept
Definition: SpacePointAlg.h:91
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:137
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.

References FindSeeds().

Referenced by trkf::Track3DKalmanHitAlg::makeTracks(), and trkf::SeedFinderModule::produce().

981  {
982  return FindSeeds(clockData, detProp, Hits, HitCatalogue, StopAfter);
983  }
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
SpacePointAlg* trkf::SeedFinderAlgorithm::GetSpacePointAlg ( ) const
inline

Definition at line 66 of file SeedFinderAlgorithm.h.

References lar::dump::vector().

66 { return fSptalg; }
void trkf::SeedFinderAlgorithm::reconfigure ( fhicl::ParameterSet const &  pset)

Definition at line 34 of file SeedFinderAlgorithm.cxx.

References CalculateGeometricalElements(), fAllow2DSeeds, fExtendSeeds, fHitResolution, fInitSeedLength, fLengthCut, fMaxViewRMS, fMinPointsInSeed, fOccupancyCut, fRefits, fSptalg, and fhicl::ParameterSet::get().

Referenced by SeedFinderAlgorithm().

35  {
36 
37  fSptalg = new trkf::SpacePointAlg(pset.get<fhicl::ParameterSet>("SpacePointAlg"));
38 
39  fInitSeedLength = pset.get<double>("InitSeedLength");
40  fMinPointsInSeed = pset.get<int>("MinPointsInSeed");
41 
42  fRefits = pset.get<double>("Refits");
43 
44  fHitResolution = pset.get<double>("HitResolution");
45 
46  fOccupancyCut = pset.get<double>("OccupancyCut");
47  fLengthCut = pset.get<double>("LengthCut");
48  fExtendSeeds = pset.get<bool>("ExtendSeeds");
49 
50  fAllow2DSeeds = pset.get<bool>("Allow2DSeeds", false);
51 
52  fMaxViewRMS.resize(3);
53  fMaxViewRMS = pset.get<std::vector<double>>("MaxViewRMS");
54 
56  }
std::vector< double > fMaxViewRMS

Member Data Documentation

bool trkf::SeedFinderAlgorithm::fAllow2DSeeds
private

Definition at line 139 of file SeedFinderAlgorithm.h.

Referenced by ConsolidateSeed(), FindSeedAtEnd(), FindSeeds(), and reconfigure().

bool trkf::SeedFinderAlgorithm::fExtendSeeds
private

Definition at line 137 of file SeedFinderAlgorithm.h.

Referenced by FindSeeds(), and reconfigure().

float trkf::SeedFinderAlgorithm::fHitResolution
private

Definition at line 131 of file SeedFinderAlgorithm.h.

Referenced by ConsolidateSeed(), and reconfigure().

double trkf::SeedFinderAlgorithm::fInitSeedLength
private

Definition at line 123 of file SeedFinderAlgorithm.h.

Referenced by FindSeedAtEnd(), and reconfigure().

double trkf::SeedFinderAlgorithm::fLengthCut
private

Definition at line 135 of file SeedFinderAlgorithm.h.

Referenced by FindSeeds(), and reconfigure().

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

Definition at line 129 of file SeedFinderAlgorithm.h.

Referenced by FindSeedAtEnd(), FindSeeds(), and reconfigure().

int trkf::SeedFinderAlgorithm::fMinPointsInSeed
private

Definition at line 125 of file SeedFinderAlgorithm.h.

Referenced by FindSeedAtEnd(), FindSeeds(), and reconfigure().

size_t trkf::SeedFinderAlgorithm::fNChannels
private

Definition at line 146 of file SeedFinderAlgorithm.h.

Referenced by CalculateGeometricalElements(), ConsolidateSeed(), and FindSeeds().

float trkf::SeedFinderAlgorithm::fOccupancyCut
private

Definition at line 133 of file SeedFinderAlgorithm.h.

Referenced by ConsolidateSeed(), and reconfigure().

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

Definition at line 142 of file SeedFinderAlgorithm.h.

Referenced by CalculateGeometricalElements(), and GetCenterAndDirection().

std::vector<double> trkf::SeedFinderAlgorithm::fPitches
private
int trkf::SeedFinderAlgorithm::fRefits
private

Definition at line 127 of file SeedFinderAlgorithm.h.

Referenced by FindSeeds(), and reconfigure().

SpacePointAlg* trkf::SeedFinderAlgorithm::fSptalg
private
std::vector<TVector3> trkf::SeedFinderAlgorithm::fWireDir
private

Definition at line 143 of file SeedFinderAlgorithm.h.

Referenced by CalculateGeometricalElements(), and GetCenterAndDirection().

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

Definition at line 144 of file SeedFinderAlgorithm.h.

Referenced by CalculateGeometricalElements(), and GetCenterAndDirection().

TVector3 trkf::SeedFinderAlgorithm::fXDir
private

Definition at line 145 of file SeedFinderAlgorithm.h.

Referenced by CalculateGeometricalElements(), and GetCenterAndDirection().

TVector3 trkf::SeedFinderAlgorithm::fYDir
private

Definition at line 145 of file SeedFinderAlgorithm.h.

Referenced by CalculateGeometricalElements().

TVector3 trkf::SeedFinderAlgorithm::fZDir
private

Definition at line 145 of file SeedFinderAlgorithm.h.

Referenced by CalculateGeometricalElements().


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