LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
SeedFinderAlgorithm.cxx
Go to the documentation of this file.
1 //
2 // Name: SeedFinderAlgorithm.cxx
3 //
4 // Purpose: Implementation file for module SeedFinderAlgorithm.
5 //
6 // Ben Jones, MIT
7 //
8 
9 #include <stdint.h>
10 
15 
25 
26 namespace trkf {
27 
28  //----------------------------------------------------------------------------
30  {
31  fSptalg = new trkf::SpacePointAlg(pset.get<fhicl::ParameterSet>("SpacePointAlg"));
32 
33  fInitSeedLength = pset.get<double>("InitSeedLength");
34  fMinPointsInSeed = pset.get<int>("MinPointsInSeed");
35 
36  fRefits = pset.get<double>("Refits");
37 
38  fHitResolution = pset.get<double>("HitResolution");
39 
40  fOccupancyCut = pset.get<double>("OccupancyCut");
41  fLengthCut = pset.get<double>("LengthCut");
42  fExtendSeeds = pset.get<bool>("ExtendSeeds");
43 
44  fAllow2DSeeds = pset.get<bool>("Allow2DSeeds", false);
45 
46  fMaxViewRMS.resize(3);
47  fMaxViewRMS = pset.get<std::vector<double>>("MaxViewRMS");
48 
50  }
51 
52  //------------------------------------------------------------
53  // Given a set of spacepoints, find seeds, and catalogue
54  // spacepoints by the seeds they formed
55  //
56  std::vector<recob::Seed> SeedFinderAlgorithm::FindSeeds(
57  detinfo::DetectorClocksData const& clockData,
58  detinfo::DetectorPropertiesData const& detProp,
59  art::PtrVector<recob::Hit> const& HitsFlat,
60  std::vector<art::PtrVector<recob::Hit>>& CataloguedHits,
61  unsigned int StopAfter) const
62  {
63  // Vector of seeds found to return
64  std::vector<recob::Seed> ReturnVector;
65 
66  // First check if there is any overlap
67  std::vector<recob::SpacePoint> spts;
68 
70  fSptalg->makeSpacePoints(clockData, detProp, HitsFlat, spts);
71 
72  if (int(spts.size()) < fMinPointsInSeed) { return std::vector<recob::Seed>(); }
73 
74  // If we got here, we have enough spacepoints to potentially make seeds from these hits.
75 
76  // This is table will let us quickly look up which hits are in a given view / channel.
77  // structure is OrgHits[View][Channel] = {index1,index2...}
78  // where the indices are of HitsFlat[index].
79 
80  std::vector<std::vector<std::vector<int>>> OrgHits(3);
81  for (size_t n = 0; n != 3; ++n)
82  OrgHits[n].resize(fNChannels);
83 
84  // These two tables contain the hit to spacepoint mappings
85 
86  std::vector<std::vector<int>> SpacePointsPerHit(HitsFlat.size(), std::vector<int>());
87  std::vector<std::vector<int>> HitsPerSpacePoint(spts.size(), std::vector<int>());
88 
89  // This vector records the status of each hit.
90  // The possible values are:
91  // 0. hit unused
92  // 1. hit used in a seed
93  // Initially none are used.
94 
95  std::vector<char> HitStatus(HitsFlat.size(), 0);
96 
97  // Fill the organizing map
98 
99  for (size_t i = 0; i != HitsFlat.size(); ++i) {
100  OrgHits[HitsFlat.at(i)->View()][HitsFlat.at(i)->Channel()].push_back(i);
101  }
102 
103  // Fill the spacepoint / hit lookups
104 
105  for (size_t iSP = 0; iSP != spts.size(); ++iSP) {
106  art::PtrVector<recob::Hit> HitsThisSP = fSptalg->getAssociatedHits(spts.at(iSP));
107 
108  for (size_t iH = 0; iH != HitsThisSP.size(); ++iH) {
109  geo::View_t ThisView = HitsThisSP.at(iH)->View();
110  uint32_t ThisChannel = HitsThisSP.at(iH)->Channel();
111  float ThisTime = HitsThisSP.at(iH)->PeakTime();
112  float eta = 0.001;
113 
114  for (size_t iOrg = 0; iOrg != OrgHits[ThisView][ThisChannel].size(); ++iOrg) {
115  if (fabs(ThisTime - HitsFlat.at(OrgHits[ThisView][ThisChannel][iOrg])->PeakTime()) <
116  eta) {
117  SpacePointsPerHit.at(OrgHits[ThisView][ThisChannel][iOrg]).push_back(iSP);
118  HitsPerSpacePoint.at(iSP).push_back(OrgHits[ThisView][ThisChannel][iOrg]);
119  }
120  }
121  }
122  }
123 
124  // The general idea will be:
125  // A. Make spacepoints from remaining hits
126  // B. Look for 1 seed in that vector
127  // C. Discard used hits
128  // D. Repeat
129 
130  // This vector keeps track of the status of each space point.
131  // The key is the position in the AllSpacePoints vector.
132  // The value is
133  // 0: point unused
134  // 1: point used in seed
135  // 2: point thrown but unused
136  // 3: flag to terminate seed finding
137 
138  std::vector<char> PointStatus(spts.size(), 0);
139 
140  std::vector<std::map<geo::View_t, std::vector<int>>> WhichHitsPerSeed;
141 
142  bool KeepChopping = true;
143 
144  while (KeepChopping) {
145  // This vector keeps a list of the points used in this seed
146  std::vector<int> PointsUsed;
147 
148  // Find exactly one seed, starting at high Z
149  recob::Seed TheSeed =
150  FindSeedAtEnd(detProp, spts, PointStatus, PointsUsed, HitsFlat, OrgHits);
151 
152  // If it was a good seed, collect up the relevant spacepoints
153  // and add the seed to the return vector
154 
155  if (TheSeed.IsValid()) {
156 
157  for (size_t iP = 0; iP != PointsUsed.size(); ++iP) {
158  for (size_t iH = 0; iH != HitsPerSpacePoint.at(PointsUsed.at(iP)).size(); ++iH) {
159  int UsedHitID = HitsPerSpacePoint.at(PointsUsed.at(iP)).at(iH);
160  HitStatus[UsedHitID] = 2;
161  }
162  }
163  PointStatus[PointsUsed.at(0)] = 1;
164  ConsolidateSeed(detProp, TheSeed, HitsFlat, HitStatus, OrgHits, false);
165  }
166 
167  if (TheSeed.IsValid()) {
168  if (fRefits > 0) {
169  std::vector<char> HitStatusGood;
170  recob::Seed SeedGood;
171  for (size_t r = 0; r != (unsigned int)fRefits; ++r) {
172  double PrevLength = TheSeed.GetLength();
173 
174  SeedGood = TheSeed;
175  HitStatusGood = HitStatus;
176 
177  std::vector<int> PresentHitList;
178  for (size_t iH = 0; iH != HitStatus.size(); ++iH) {
179  if (HitStatus[iH] == 2) { PresentHitList.push_back(iH); }
180  }
181  double pt[3], dir[3], err[3];
182 
183  TheSeed.GetPoint(pt, err);
184  TheSeed.GetDirection(dir, err);
185 
186  TVector3 Center, Direction;
187  std::vector<double> ViewRMS;
188  std::vector<int> HitsPerView;
190  detProp, HitsFlat, PresentHitList, Center, Direction, ViewRMS, HitsPerView);
191 
192  Direction = Direction.Unit() * TheSeed.GetLength();
193 
194  int nViewsWithHits(0);
195  for (size_t n = 0; n != 3; ++n) {
196  pt[n] = Center[n];
197  dir[n] = Direction[n];
198 
199  TheSeed.SetPoint(pt, err);
200  TheSeed.SetDirection(dir, err);
201 
202  if (HitsPerView[n] > 0) nViewsWithHits++;
203  }
204 
205  if (nViewsWithHits < 2) TheSeed.SetValidity(false);
206 
207  if (TheSeed.IsValid())
208  ConsolidateSeed(detProp, TheSeed, HitsFlat, HitStatus, OrgHits, fExtendSeeds);
209 
210  // if we accidentally invalidated the seed, go back to the old one and escape
211  else {
212  // If we invalidated the seed, go back one interaction
213  // and kill the loop
214  HitStatus = HitStatusGood;
215  TheSeed = SeedGood;
216  break;
217  }
218  if ((r > 0) && (fabs(PrevLength - TheSeed.GetLength()) < fPitches.at(0))) {
219  // If we didn't change much, kill the loop
220  break;
221  }
222  }
223  }
224  }
225 
226  if (TheSeed.IsValid()) {
227  WhichHitsPerSeed.push_back(std::map<geo::View_t, std::vector<int>>());
228 
229  art::PtrVector<recob::Hit> HitsWithThisSeed;
230  for (size_t iH = 0; iH != HitStatus.size(); ++iH) {
231  if (HitStatus.at(iH) == 2) {
232  WhichHitsPerSeed.at(WhichHitsPerSeed.size() - 1)[HitsFlat[iH]->View()].push_back(iH);
233  HitsWithThisSeed.push_back(HitsFlat.at(iH));
234  HitStatus.at(iH) = 1;
235 
236  for (size_t iSP = 0; iSP != SpacePointsPerHit.at(iH).size(); ++iSP) {
237  PointStatus[SpacePointsPerHit.at(iH).at(iSP)] = 1;
238  }
239  }
240  }
241 
242  // Record that we used this set of hits with this seed in the return catalogue
243  ReturnVector.push_back(TheSeed);
244  CataloguedHits.push_back(HitsWithThisSeed);
245 
246  // Tidy up
247  HitsWithThisSeed.clear();
248  }
249  else {
250  // If it was not a good seed, throw out the top SP and try again
251  PointStatus.at(PointsUsed.at(0)) = 2;
252  }
253 
254  int TotalSPsUsed = 0;
255  for (size_t i = 0; i != PointStatus.size(); ++i) {
256  if (PointStatus[i] != 0) TotalSPsUsed++;
257  }
258 
259  if ((int(spts.size()) - TotalSPsUsed) < fMinPointsInSeed) KeepChopping = false;
260 
261  if ((PointStatus[0] == 3) || (PointStatus.size() == 0)) KeepChopping = false;
262 
263  PointsUsed.clear();
264 
265  if ((ReturnVector.size() >= StopAfter) && (StopAfter > 0)) break;
266 
267  } // end spt-level loop
268 
269  // If we didn't find any seeds, we can make one last ditch attempt. See
270  // if the whole collection is colinear enough to make one megaseed
271  // (good for patching up high angle tracks)
272 
273  if (ReturnVector.size() == 0) {
274  std::vector<int> ListAllHits;
275  for (size_t i = 0; i != HitsFlat.size(); ++i) {
276  ListAllHits.push_back(i);
277  HitStatus[i] = 2;
278  }
279  TVector3 SeedCenter(0, 0, 0);
280  TVector3 SeedDirection(0, 0, 0);
281 
282  std::vector<double> ViewRMS;
283  std::vector<int> HitsPerView;
284 
285  std::vector<art::PtrVector<recob::Hit>> HitsInThisCollection(3);
286 
288  detProp, HitsFlat, ListAllHits, SeedCenter, SeedDirection, ViewRMS, HitsPerView);
289 
290  bool ThrowOutSeed = false;
291 
292  double PtArray[3], DirArray[3];
293  int nViewsWithHits(0);
294  for (size_t n = 0; n != 3; ++n) {
295  PtArray[n] = SeedCenter[n];
296  DirArray[n] = SeedDirection[n];
297  if (HitsPerView[n] > 0) nViewsWithHits++;
298  }
299  recob::Seed TheSeed(PtArray, DirArray);
300 
301  if (nViewsWithHits < 2 || (nViewsWithHits < 3 && !fAllow2DSeeds)) ThrowOutSeed = true;
302 
303  if (!ThrowOutSeed) {
304  ConsolidateSeed(detProp, TheSeed, HitsFlat, HitStatus, OrgHits, false);
305 
306  // Now we have consolidated, grab the right
307  // hits to find the RMS and refitted direction
308  ListAllHits.clear();
309  for (size_t i = 0; i != HitStatus.size(); ++i) {
310  if (HitStatus.at(i) == 2) ListAllHits.push_back(i);
311  }
312  std::vector<int> HitsPerView;
314  detProp, HitsFlat, ListAllHits, SeedCenter, SeedDirection, ViewRMS, HitsPerView);
315 
316  int nViewsWithHits(0);
317  for (size_t n = 0; n != 3; ++n) {
318  PtArray[n] = SeedCenter[n];
319  DirArray[n] = SeedDirection[n];
320  // ThrowOutSeed = true;
321  if (HitsPerView[n] > 0) nViewsWithHits++;
322  }
323 
324  if (nViewsWithHits < 2 || (nViewsWithHits < 3 && !fAllow2DSeeds)) ThrowOutSeed = true;
325 
326  TheSeed = recob::Seed(PtArray, DirArray);
327 
328  if (fMaxViewRMS.at(0) > 0) {
329  for (size_t j = 0; j != fMaxViewRMS.size(); j++) {
330  if (fMaxViewRMS.at(j) < ViewRMS.at(j)) { ThrowOutSeed = true; }
331  }
332  }
333  }
334 
335  if ((!ThrowOutSeed) && (TheSeed.IsValid())) {
336  ReturnVector.push_back(TheSeed);
337  art::PtrVector<recob::Hit> HitsThisSeed;
338  for (size_t i = 0; i != ListAllHits.size(); ++i) {
339  HitsThisSeed.push_back(HitsFlat.at(ListAllHits.at(i)));
340  }
341  CataloguedHits.push_back(HitsThisSeed);
342  }
343  }
344 
345  // Tidy up
346  SpacePointsPerHit.clear();
347  HitsPerSpacePoint.clear();
348  PointStatus.clear();
349  OrgHits.clear();
350  HitStatus.clear();
351 
352  for (size_t i = 0; i != ReturnVector.size(); ++i) {
353  double CrazyValue = 1000000;
354  double Length = ReturnVector.at(i).GetLength();
355  if (!((Length > fLengthCut) && (Length < CrazyValue))) {
356  ReturnVector.erase(ReturnVector.begin() + i);
357  CataloguedHits.erase(CataloguedHits.begin() + i);
358  --i;
359  }
360  }
361 
362  return ReturnVector;
363  }
364 
365  //------------------------------------------------------------
366  // Latest extendseed method
367  //
368 
370  recob::Seed& TheSeed,
371  art::PtrVector<recob::Hit> const& HitsFlat,
372  std::vector<char>& HitStatus,
373  std::vector<std::vector<std::vector<int>>>& OrgHits,
374  bool Extend) const
375  {
376 
377  bool ThrowOutSeed = false;
378 
379  // This will keep track of what hits are in this seed
380  std::map<geo::View_t, std::map<uint32_t, std::vector<int>>> HitsInThisSeed;
381 
382  int NHitsThisSeed = 0;
383 
384  double MinS = 1000, MaxS = -1000;
385  for (size_t i = 0; i != HitStatus.size(); ++i) {
386  if (HitStatus.at(i) == 2) {
387  double disp, s;
388  GetHitDistAndProj(detProp, TheSeed, HitsFlat.at(i), disp, s);
389  if (fabs(s) > 1.2) {
390  // This hit is not rightfully part of this seed, toss it.
391  HitStatus[i] = 0;
392  }
393  else {
394  NHitsThisSeed++;
395 
396  if (s < MinS) MinS = s;
397  if (s > MaxS) MaxS = s;
398  HitsInThisSeed[HitsFlat.at(i)->View()][HitsFlat.at(i)->Channel()].push_back(i);
399  }
400  }
401  }
402 
403  double LengthRescale = (MaxS - MinS) / 2.;
404  double PositionShift = (MaxS + MinS) / 2.;
405 
406  double pt[3], dir[3], err[3];
407  TheSeed.GetPoint(pt, err);
408  TheSeed.GetDirection(dir, err);
409 
410  for (size_t n = 0; n != 3; ++n) {
411  pt[n] += dir[n] * PositionShift;
412  dir[n] *= LengthRescale;
413  }
414 
415  TheSeed.SetPoint(pt, err);
416  TheSeed.SetDirection(dir, err);
417 
418  // Run through checking if we missed any hits
419  for (auto itP = HitsInThisSeed.begin(); itP != HitsInThisSeed.end(); ++itP) {
420  double dist, s;
421  geo::View_t View = itP->first;
422  uint32_t LowestChan = itP->second.begin()->first;
423  uint32_t HighestChan = itP->second.rbegin()->first;
424  for (uint32_t c = LowestChan; c != HighestChan; ++c) {
425  for (size_t h = 0; h != OrgHits[View][c].size(); ++h) {
426  if (HitStatus[OrgHits[View][c].at(h)] == 0) {
427  GetHitDistAndProj(detProp, TheSeed, HitsFlat[OrgHits[View][c].at(h)], dist, s);
428  if (dist < fHitResolution) {
429  NHitsThisSeed++;
430 
431  HitStatus[OrgHits[View][c].at(h)] = 2;
432 
433  HitsInThisSeed[View][c].push_back(OrgHits[View][c].at(h));
434  }
435  else
436  HitStatus[OrgHits[View][c].at(h)] = 0;
437  }
438  }
439  }
440  }
441 
442  if (NHitsThisSeed == 0) ThrowOutSeed = true;
443 
444  // Check seed occupancy
445 
446  uint32_t LowestChanInSeed[3], HighestChanInSeed[3];
447  double Occupancy[] = {0., 0., 0.};
448  int nHitsPerView[] = {0, 0, 0};
449 
450  for (auto itP = HitsInThisSeed.begin(); itP != HitsInThisSeed.end(); ++itP) {
451 
452  geo::View_t View = itP->first;
453 
454  LowestChanInSeed[View] = itP->second.begin()->first;
455  HighestChanInSeed[View] = itP->second.rbegin()->first;
456 
457  nHitsPerView[View]++;
458 
459  int FilledChanCount = 0;
460 
461  for (size_t c = LowestChanInSeed[View]; c != HighestChanInSeed[View]; ++c) {
462  if (itP->second[c].size() > 0) ++FilledChanCount;
463  }
464 
465  Occupancy[View] =
466  float(FilledChanCount) / float(HighestChanInSeed[View] - LowestChanInSeed[View]);
467  }
468 
469  int nBelowCut(0);
470  int nViewsWithHits(0);
471  for (size_t n = 0; n != 3; ++n) {
472  if (Occupancy[n] < fOccupancyCut) nBelowCut++;
473  if (nHitsPerView[n] > 0) nViewsWithHits++;
474  }
475 
476  int belowCut(0);
477 
478  if (fAllow2DSeeds && nViewsWithHits < 3) belowCut = 1;
479 
480  if (nBelowCut > belowCut) ThrowOutSeed = true;
481 
482  if ((Extend) && (!ThrowOutSeed)) {
483  std::vector<std::vector<double>> ToAddNegativeS(3, std::vector<double>());
484  std::vector<std::vector<double>> ToAddPositiveS(3, std::vector<double>());
485  std::vector<std::vector<int>> ToAddNegativeH(3, std::vector<int>());
486  std::vector<std::vector<int>> ToAddPositiveH(3, std::vector<int>());
487 
488  for (auto itP = HitsInThisSeed.begin(); itP != HitsInThisSeed.end(); ++itP) {
489  double dist, s;
490 
491  geo::View_t View = itP->first;
492 
493  if (LowestChanInSeed[View] > 0) {
494  for (uint32_t c = LowestChanInSeed[View] - 1; c != 0; --c) {
495  bool GotOneThisChannel = false;
496  for (size_t h = 0; h != OrgHits[View][c].size(); ++h) {
497  if (HitStatus[OrgHits[View][c][h]] == 0) {
498  GetHitDistAndProj(detProp, TheSeed, HitsFlat[OrgHits[View][c].at(h)], dist, s);
499  if (dist < fHitResolution) {
500  GotOneThisChannel = true;
501  if (s < 0) {
502  ToAddNegativeS[View].push_back(s);
503  ToAddNegativeH[View].push_back(OrgHits[View][c].at(h));
504  }
505  else {
506  ToAddPositiveS[View].push_back(s);
507  ToAddPositiveH[View].push_back(OrgHits[View][c].at(h));
508  }
509  }
510  }
511  }
512  if (GotOneThisChannel == false) break;
513  }
514  }
515  if (HighestChanInSeed[View] < fNChannels)
516 
517  for (uint32_t c = HighestChanInSeed[View] + 1; c != fNChannels; ++c) {
518  bool GotOneThisChannel = false;
519  for (size_t h = 0; h != OrgHits[View][c].size(); ++h) {
520  if (HitStatus[OrgHits[View][c][h]] == 0) {
521  GetHitDistAndProj(detProp, TheSeed, HitsFlat[OrgHits[View][c].at(h)], dist, s);
522  if (dist < fHitResolution) {
523  GotOneThisChannel = true;
524  if (s < 0) {
525 
526  ToAddNegativeS[View].push_back(s);
527  ToAddNegativeH[View].push_back(OrgHits[View][c].at(h));
528  }
529  else {
530  ToAddPositiveS[View].push_back(s);
531  ToAddPositiveH[View].push_back(OrgHits[View][c].at(h));
532  }
533  }
534  }
535  }
536  if (GotOneThisChannel == false) break;
537  }
538  }
539 
540  double ExtendPositiveS = 0, ExtendNegativeS = 0;
541 
542  if ((ToAddPositiveS[0].size() > 0) && (ToAddPositiveS[1].size() > 0) &&
543  (ToAddPositiveS[2].size() > 0)) {
544  for (size_t n = 0; n != 3; ++n) {
545  int n1 = (n + 1) % 3;
546  int n2 = (n + 2) % 3;
547 
548  if ((ToAddPositiveS[n].back() <= ToAddPositiveS[n1].back()) &&
549  (ToAddPositiveS[n].back() <= ToAddPositiveS[n2].back())) {
550  ExtendPositiveS = ToAddPositiveS[n].back();
551  }
552  }
553  }
554 
555  if ((ToAddNegativeS[0].size() > 0) && (ToAddNegativeS[1].size() > 0) &&
556  (ToAddNegativeS[2].size() > 0)) {
557  for (size_t n = 0; n != 3; ++n) {
558  int n1 = (n + 1) % 3;
559  int n2 = (n + 2) % 3;
560  if ((ToAddNegativeS[n].back() >= ToAddNegativeS[n1].back()) &&
561  (ToAddNegativeS[n].back() >= ToAddNegativeS[n2].back())) {
562  ExtendNegativeS = ToAddNegativeS[n].back();
563  }
564  }
565  }
566 
567  if (fabs(ExtendNegativeS) < 1.) ExtendNegativeS = -1.;
568  if (fabs(ExtendPositiveS) < 1.) ExtendPositiveS = 1.;
569 
570  LengthRescale = (ExtendPositiveS - ExtendNegativeS) / 2.;
571  PositionShift = (ExtendPositiveS + ExtendNegativeS) / 2.;
572 
573  for (size_t n = 0; n != 3; ++n) {
574  pt[n] += dir[n] * PositionShift;
575  dir[n] *= LengthRescale;
576 
577  for (size_t i = 0; i != ToAddPositiveS[n].size(); ++i) {
578  if (ToAddPositiveS[n].at(i) < ExtendPositiveS)
579  HitStatus[ToAddPositiveH[n].at(i)] = 2;
580  else
581  HitStatus[ToAddPositiveH[n].at(i)] = 0;
582  }
583 
584  for (size_t i = 0; i != ToAddNegativeS[n].size(); ++i) {
585  if (ToAddNegativeS[n].at(i) > ExtendNegativeS)
586  HitStatus[ToAddNegativeH[n].at(i)] = 2;
587  else
588  HitStatus[ToAddNegativeH[n].at(i)] = 0;
589  }
590  }
591 
592  TheSeed.SetPoint(pt, err);
593  TheSeed.SetDirection(dir, err);
594  }
595 
596  if (ThrowOutSeed) TheSeed.SetValidity(false);
597  }
598 
599  //------------------------------------------------------------
600 
602  recob::Seed const& ASeed,
603  art::Ptr<recob::Hit> const& AHit,
604  double& disp,
605  double& s) const
606  {
607  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
608 
609  constexpr geo::TPCID tpcid{0, 0};
610  geo::PlaneID const planeID{tpcid, AHit->WireID().Plane};
611 
612  auto const [start, end] = wireReadoutGeom.WireEndPoints(geo::WireID{planeID, 0});
613 
614  double HitX = detProp.ConvertTicksToX(AHit->PeakTime(), planeID);
615  double HitXHigh = detProp.ConvertTicksToX(AHit->PeakTimePlusRMS(), planeID);
616  double HitXLow = detProp.ConvertTicksToX(AHit->PeakTimeMinusRMS(), planeID);
617 
618  double HitWidth = HitXHigh - HitXLow;
619 
620  double pt[3], dir[3], err[3];
621 
622  ASeed.GetDirection(dir, err);
623  ASeed.GetPoint(pt, err);
624 
625  TVector3 sPt(pt[0], pt[1], pt[2]);
626  TVector3 sDir(dir[0], dir[1], dir[2]);
627  TVector3 hPt(HitX, start.Y(), start.Z());
628  TVector3 hDir(0, start.Y() - end.Y(), start.Z() - end.Z());
629 
630  s = (sPt - hPt).Dot(hDir * (hDir.Dot(sDir)) - sDir * (hDir.Dot(hDir))) /
631  (hDir.Dot(hDir) * sDir.Dot(sDir) - pow(hDir.Dot(sDir), 2));
632 
633  disp = fabs((sPt - hPt).Dot(sDir.Cross(hDir)) / (sDir.Cross(hDir)).Mag()) / HitWidth;
634  }
635 
636  //------------------------------------------------------------
637  // Try to find one seed at the high Z end of a set of spacepoints
638  //
639 
641  detinfo::DetectorPropertiesData const& detProp,
642  std::vector<recob::SpacePoint> const& Points,
643  std::vector<char>& PointStatus,
644  std::vector<int>& PointsInRange,
645  art::PtrVector<recob::Hit> const& HitsFlat,
646  std::vector<std::vector<std::vector<int>>>& OrgHits) const
647  {
648  // This pointer will be returned later
649  recob::Seed ReturnSeed;
650 
651  // Keep track of spacepoints we used, not just their IDs
652  std::vector<recob::SpacePoint> PointsUsed;
653 
654  // Clear output vector
655  PointsInRange.clear();
656 
657  // Loop through hits looking for highest Z seedable point
658  TVector3 HighestZPoint;
659  bool NoPointFound = true;
660  int counter = Points.size() - 1;
661  while ((NoPointFound == true) && (counter >= 0)) {
662  if (PointStatus[counter] == 0) {
663  HighestZPoint = TVector3(
664  Points.at(counter).XYZ()[0], Points.at(counter).XYZ()[1], Points.at(counter).XYZ()[2]);
665  NoPointFound = false;
666  }
667  else
668  counter--;
669  }
670  if (NoPointFound) {
671  // We didn't find a high point at all
672  // - let the algorithm know to give up.
673  PointStatus[0] = 3;
674  }
675 
676  // Now we have the high Z point, loop through collecting
677  // near enough hits. We look 2 seed lengths away, since
678  // the seed is bidirectional from the central point
679 
680  double TwiceLength = 2.0 * fInitSeedLength;
681 
682  for (int index = Points.size() - 1; index != -1; --index) {
683  if (PointStatus[index] == 0) {
684  // first check z, then check total distance
685  // (much faster, since most will be out of range in z anyway)
686  if ((HighestZPoint[2] - Points.at(index).XYZ()[2]) < TwiceLength) {
687  double DistanceToHighZ = pow(pow(HighestZPoint[1] - Points.at(index).XYZ()[1], 2) +
688  pow(HighestZPoint[2] - Points.at(index).XYZ()[2], 2),
689  0.5);
690  if (DistanceToHighZ < TwiceLength) {
691  PointsInRange.push_back(index);
692  PointsUsed.push_back(Points.at(index));
693  }
694  }
695  else
696  break;
697  }
698  }
699 
700  TVector3 SeedCenter(0, 0, 0);
701  TVector3 SeedDirection(0, 0, 0);
702 
703  // Check we have enough points in here to form a seed,
704  // otherwise return a dud
705  int NPoints = PointsInRange.size();
706 
707  if (NPoints < fMinPointsInSeed) return recob::Seed();
708 
709  std::map<int, bool> HitMap;
710  std::vector<int> HitList;
711 
712  for (unsigned int i = 0; i != PointsInRange.size(); i++) {
713  std::vector<art::PtrVector<recob::Hit>> HitsInThisCollection(3);
714 
715  art::PtrVector<recob::Hit> HitsThisSP =
716  fSptalg->getAssociatedHits((Points.at(PointsInRange.at(i))));
717  for (art::PtrVector<recob::Hit>::const_iterator itHit = HitsThisSP.begin();
718  itHit != HitsThisSP.end();
719  ++itHit) {
720  uint32_t Channel = (*itHit)->Channel();
721  geo::View_t View = (*itHit)->View();
722 
723  double eta = 0.01;
724  for (size_t iH = 0; iH != OrgHits[View][Channel].size(); ++iH) {
725  if (fabs(HitsFlat[OrgHits[View][Channel][iH]]->PeakTime() - (*itHit)->PeakTime()) < eta) {
726  HitMap[OrgHits[View][Channel][iH]] = true;
727  }
728  }
729  }
730  }
731 
732  for (auto itH = HitMap.begin(); itH != HitMap.end(); ++itH) {
733  HitList.push_back(itH->first);
734  }
735 
736  std::vector<double> ViewRMS;
737  std::vector<int> HitsPerView;
738 
740  detProp, HitsFlat, HitList, SeedCenter, SeedDirection, ViewRMS, HitsPerView);
741 
742  HitMap.clear();
743  HitList.clear();
744 
745  // See if seed points have some linearity
746 
747  bool ThrowOutSeed = false;
748 
749  double PtArray[3], DirArray[3];
750  double AngleFactor =
751  pow(pow(SeedDirection.Y(), 2) + pow(SeedDirection.Z(), 2), 0.5) / SeedDirection.Mag();
752 
753  int nViewsWithHits(0);
754 
755  for (size_t n = 0; n != 3; ++n) {
756  DirArray[n] = SeedDirection[n] * fInitSeedLength / AngleFactor;
757  PtArray[n] = SeedCenter[n];
758  if (HitsPerView[n] > 0) nViewsWithHits++;
759  }
760 
761  if (nViewsWithHits < 2 || (nViewsWithHits < 3 && !fAllow2DSeeds)) ThrowOutSeed = true;
762 
763  ReturnSeed = recob::Seed(PtArray, DirArray);
764 
765  if (fMaxViewRMS.at(0) > 0) {
766 
767  for (size_t j = 0; j != fMaxViewRMS.size(); j++) {
768  if (fMaxViewRMS.at(j) < ViewRMS.at(j)) { ThrowOutSeed = true; }
769  // mf::LogVerbatim("SeedFinderAlgorithm") << RMS.at(j);
770  }
771  }
772 
773  // If the seed is marked as bad, return a dud, otherwise
774  // return the ReturnSeed pointer
775  if (!ThrowOutSeed)
776  return ReturnSeed;
777  else
778  return recob::Seed();
779  }
780 
781  //-----------------------------------------------------------
782 
784  art::PtrVector<recob::Hit> const& HitsFlat,
785  std::vector<int>& HitsToUse,
786  TVector3& Center,
787  TVector3& Direction,
788  std::vector<double>& ViewRMS,
789  std::vector<int>& N) const
790  {
791  N.resize(3);
792 
793  std::map<uint32_t, bool> HitsClaimed;
794 
795  // We'll store hit coordinates in each view into this vector
796 
797  std::vector<std::vector<double>> HitTimes(3);
798  std::vector<std::vector<double>> HitWires(3);
799  std::vector<std::vector<double>> HitWidths(3);
800  std::vector<double> MeanWireCoord(3, 0);
801  std::vector<double> MeanTimeCoord(3, 0);
802 
803  // Run through the collection getting hit info for these spacepoints
804 
805  std::vector<double> x(3, 0), y(3, 0), xx(3, 0), xy(3, 0), yy(3, 0), sig(3, 0);
806 
807  for (size_t i = 0; i != HitsToUse.size(); ++i) {
808  auto itHit = HitsFlat.begin() + HitsToUse[i];
809 
810  size_t ViewIndex;
811 
812  auto const hitView = (*itHit)->View();
813  if (hitView == geo::kU)
814  ViewIndex = 0;
815  else if (hitView == geo::kV)
816  ViewIndex = 1;
817  else if (hitView == geo::kW)
818  ViewIndex = 2;
819  else {
821  << "SpacePointAlg does not support view " << geo::PlaneGeo::ViewName(hitView) << " (#"
822  << hitView << ")\n";
823  }
824  double WireCoord = (*itHit)->WireID().Wire * fPitches.at(ViewIndex);
825  double TimeCoord = detProp.ConvertTicksToX((*itHit)->PeakTime(), ViewIndex, 0, 0);
826  double TimeUpper = detProp.ConvertTicksToX((*itHit)->PeakTimePlusRMS(), ViewIndex, 0, 0);
827  double TimeLower = detProp.ConvertTicksToX((*itHit)->PeakTimeMinusRMS(), ViewIndex, 0, 0);
828  double Width = fabs(0.5 * (TimeUpper - TimeLower));
829  double Width2 = pow(Width, 2);
830 
831  HitWires.at(ViewIndex).push_back(WireCoord);
832  HitTimes.at(ViewIndex).push_back(TimeCoord);
833  HitWidths.at(ViewIndex).push_back(fabs(0.5 * (TimeUpper - TimeLower)));
834 
835  MeanWireCoord.at(ViewIndex) += WireCoord;
836  MeanTimeCoord.at(ViewIndex) += TimeCoord;
837 
838  // Elements for LS
839  x.at(ViewIndex) += WireCoord / Width2;
840  y.at(ViewIndex) += TimeCoord / Width2;
841  xy.at(ViewIndex) += (TimeCoord * WireCoord) / Width2;
842  xx.at(ViewIndex) += (WireCoord * WireCoord) / Width2;
843  yy.at(ViewIndex) += (TimeCoord * TimeCoord) / Width2;
844  sig.at(ViewIndex) += 1. / Width2;
845  N.at(ViewIndex)++;
846  }
847 
848  ViewRMS.clear();
849  ViewRMS.resize(3);
850  std::vector<double> ViewGrad(3);
851  std::vector<double> ViewOffset(3);
852 
853  for (size_t n = 0; n != 3; ++n) {
854  MeanWireCoord[n] /= N[n];
855  MeanTimeCoord[n] /= N[n];
856 
857  double BigN = 1000000;
858  double SmallN = 1. / BigN;
859 
860  if (N[n] > 2) {
861  double Numerator = (y[n] / sig[n] - xy[n] / x[n]);
862  double Denominator = (x[n] / sig[n] - xx[n] / x[n]);
863  if (fabs(Denominator) > SmallN)
864  ViewGrad.at(n) = Numerator / Denominator;
865  else
866  ViewGrad[n] = BigN;
867  }
868  else if (N[n] == 2)
869  ViewGrad[n] = xy[n] / xx[n];
870  else
871  ViewGrad[n] = BigN;
872 
873  ViewOffset.at(n) = (y[n] - ViewGrad[n] * x[n]) / sig[n];
874  ViewRMS.at(n) = pow((yy[n] + pow(ViewGrad[n], 2) * xx[n] + pow(ViewOffset[n], 2) * sig[n] -
875  2 * ViewGrad[n] * xy[n] - 2 * ViewOffset[n] * y[n] +
876  2 * ViewGrad[n] * ViewOffset[n] * x[n]) /
877  N[n],
878  0.5);
879  // Make RMS rotation perp to track
880  if (ViewGrad.at(n) != 0) ViewRMS[n] *= sin(atan(1. / ViewGrad.at(n)));
881  }
882 
883  for (size_t n = 0; n != 3; ++n) {
884  size_t n1 = (n + 1) % 3;
885  size_t n2 = (n + 2) % 3;
886  if ((N[n] <= N[n1]) && (N[n] <= N[n2])) {
887 
888  if (N[n1] < N[n2]) { std::swap(n1, n2); }
889  if ((N[n1] == 0) || (N[n2] == 0)) continue;
890 
891  Direction =
892  (fXDir + fPitchDir[n1] * (1. / ViewGrad[n1]) +
893  fWireDir[n1] *
894  (((1. / ViewGrad[n2]) - fPitchDir[n1].Dot(fPitchDir[n2]) * (1. / ViewGrad[n1])) /
895  fWireDir[n1].Dot(fPitchDir[n2])))
896  .Unit();
897 
898  /*
899  Center2D[n] =
900  fXDir * 0.5 * (MeanTimeCoord[n1]+MeanTimeCoord[n2])
901  + fPitchDir[n1] * (MeanWireCoord[n1] + fWireZeroOffset[n1])
902  + fWireDir[n1] * ( ((MeanWireCoord[n2] + fWireZeroOffset[n2]) - ( MeanWireCoord[n1] + fWireZeroOffset[n1] )*fPitchDir[n1].Dot(fPitchDir[n2]))/(fPitchDir[n2].Dot(fWireDir[n1])) );
903  */
904 
905  double TimeCoord = 0.5 * (MeanTimeCoord[n1] + MeanTimeCoord[n2]);
906  double WireCoordIn1 = (TimeCoord - ViewOffset[n1]) / ViewGrad[n1] + fWireZeroOffset[n1];
907  double WireCoordIn2 = (TimeCoord - ViewOffset[n2]) / ViewGrad[n2] + fWireZeroOffset[n2];
908 
909  Center = fXDir * TimeCoord + fPitchDir[n1] * WireCoordIn1 +
910  fWireDir[n1] * ((WireCoordIn2 - WireCoordIn1 * fPitchDir[n1].Dot(fPitchDir[n2])) /
911  (fPitchDir[n2].Dot(fWireDir[n1])));
912 
913  ViewRMS[n] = -fabs(ViewRMS[n]);
914  ViewRMS[n1] = fabs(ViewRMS[n1]);
915  ViewRMS[n2] = fabs(ViewRMS[n2]);
916 
917  break;
918  }
919  }
920  }
921 
922  //-----------------------------------------------
924  {
926  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
927 
928  // Total number of channels in the detector
929  fNChannels = wireReadoutGeom.Nchannels();
930 
931  // Find pitch of each wireplane
932  fPitches.resize(3);
933  fPitches.at(0) = fabs(wireReadoutGeom.Plane({0, 0, geo::kU}).WirePitch());
934  fPitches.at(1) = fabs(wireReadoutGeom.Plane({0, 0, geo::kV}).WirePitch());
935  fPitches.at(2) = fabs(wireReadoutGeom.Plane({0, 0, geo::kW}).WirePitch());
936 
937  // Setup basis vectors
938  fXDir = TVector3(1, 0, 0);
939  fYDir = TVector3(0, 1, 0);
940  fZDir = TVector3(0, 0, 1);
941 
942  fWireDir.resize(3);
943  fPitchDir.resize(3);
944  fWireZeroOffset.resize(3);
945 
946  double xyzStart1[3], xyzStart2[3];
947  double xyzEnd1[3], xyzEnd2[3];
948 
949  // Calculate wire coordinate systems
950  for (unsigned int n = 0; n != 3u; ++n) {
951  wireReadoutGeom.WireEndPoints(geo::WireID{0, 0, n, 0}, xyzStart1, xyzEnd1);
952  wireReadoutGeom.WireEndPoints(geo::WireID{0, 0, n, 1}, xyzStart2, xyzEnd2);
953  fWireDir[n] =
954  TVector3(xyzEnd1[0] - xyzStart1[0], xyzEnd1[1] - xyzStart1[1], xyzEnd1[2] - xyzStart1[2])
955  .Unit();
956  fPitchDir[n] = fWireDir[n].Cross(fXDir).Unit();
957  if (fPitchDir[n].Dot(TVector3(
958  xyzEnd2[0] - xyzEnd1[0], xyzEnd2[1] - xyzEnd1[1], xyzEnd2[2] - xyzEnd1[2])) < 0)
959  fPitchDir[n] = -fPitchDir[n];
960 
961  fWireZeroOffset[n] =
962  xyzEnd1[0] * fPitchDir[n][0] + xyzEnd1[1] * fPitchDir[n][1] + xyzEnd1[2] * fPitchDir[n][2];
963  }
964  }
965 
966  //-----------------------------------------------
967 
969  detinfo::DetectorClocksData const& clockData,
970  detinfo::DetectorPropertiesData const& detProp,
973  unsigned int StopAfter) const
974  {
975  return FindSeeds(clockData, detProp, Hits, HitCatalogue, StopAfter);
976  }
977 
978  //---------------------------------------------
979 
980  std::vector<std::vector<recob::Seed>> SeedFinderAlgorithm::GetSeedsFromSortedHits(
981  detinfo::DetectorClocksData const& clockData,
982  detinfo::DetectorPropertiesData const& detProp,
985  unsigned int StopAfter) const
986  {
987 
988  std::vector<std::vector<recob::Seed>> ReturnVec;
989 
990  // This piece of code looks detector specific, but its not -
991  // it also works for 2 planes, but one vector is empty.
992 
993  if (!(fSptalg->enableU() && fSptalg->enableV() && fSptalg->enableW()))
994  mf::LogWarning("BezierTrackerModule")
995  << "Warning: SpacePointAlg is does not have three views enabled. This may cause unexpected "
996  "behaviour in the bezier tracker.";
997 
998  try {
1000  SortedHits.at(geo::kU).begin();
1001  itU != SortedHits.at(geo::kU).end();
1002  ++itU)
1004  SortedHits.at(geo::kV).begin();
1005  itV != SortedHits.at(geo::kV).end();
1006  ++itV)
1008  SortedHits.at(geo::kW).begin();
1009  itW != SortedHits.at(geo::kW).end();
1010  ++itW) {
1011  art::PtrVector<recob::Hit> HitsThisComboFlat;
1012 
1013  if (fSptalg->enableU())
1014  for (size_t i = 0; i != itU->size(); ++i)
1015  HitsThisComboFlat.push_back(itU->at(i));
1016 
1017  if (fSptalg->enableV())
1018  for (size_t i = 0; i != itV->size(); ++i)
1019  HitsThisComboFlat.push_back(itV->at(i));
1020 
1021  if (fSptalg->enableW())
1022  for (size_t i = 0; i != itW->size(); ++i)
1023  HitsThisComboFlat.push_back(itW->at(i));
1024 
1025  std::vector<art::PtrVector<recob::Hit>> CataloguedHits;
1026 
1027  std::vector<recob::Seed> Seeds =
1028  FindSeeds(clockData, detProp, HitsThisComboFlat, CataloguedHits, StopAfter);
1029 
1030  // Add this harvest to return vectors
1031  HitsPerSeed.push_back(CataloguedHits);
1032  ReturnVec.push_back(Seeds);
1033 
1034  // Tidy up
1035  CataloguedHits.clear();
1036  Seeds.clear();
1037  }
1038  }
1039  catch (...) {
1040  mf::LogWarning("SeedFinderTrackerModule")
1041  << " bailed during hit map lookup - have you enabled all 3 planes?";
1042  ReturnVec.push_back(std::vector<recob::Seed>());
1043  }
1044 
1045  return ReturnVec;
1046  }
1047 
1048 }
Float_t x
Definition: compare.C:6
TRandom r
Definition: spectrum.C:23
std::vector< double > fPitches
float Length(const PFPStruct &pfp)
Definition: PFPUtils.cxx:3267
SeedFinderAlgorithm(const fhicl::ParameterSet &pset)
bool enableW() const noexcept
Definition: SpacePointAlg.h:93
bool enableV() const noexcept
Definition: SpacePointAlg.h:92
Double_t xx
Definition: macro.C:12
static std::string ViewName(View_t view)
Returns the name of the specified view.
Definition: PlaneGeo.cxx:584
bool IsValid() const
Definition: Seed.cxx:58
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Planes which measure V.
Definition: geo_types.h:132
iterator begin()
Definition: PtrVector.h:217
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
void GetPoint(double *Pt, double *Err) const
Definition: Seed.cxx:89
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
void GetHitDistAndProj(detinfo::DetectorPropertiesData const &detProp, recob::Seed const &ASeed, art::Ptr< recob::Hit > const &AHit, double &disp, double &s) const
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:290
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
std::map< int, art::Ptr< recob::Hit > > HitMap
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
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:131
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:292
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
T get(std::string const &key) const
Definition: ParameterSet.h:314
iterator end()
Definition: PtrVector.h:231
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
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
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
float PeakTimeMinusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:300
size_type size() const
Definition: PtrVector.h:302
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
double ConvertTicksToX(double ticks, int p, int t, int c) const
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
void clearHitMap() const
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:226
Encapsulate the construction of a single detector plane .
double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Contains all timing reference information for the detector.
TDirectory * dir
Definition: macro.C:5
bool enableU() const noexcept
Definition: SpacePointAlg.h:91
std::vector< TVector3 > fWireDir
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
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:133
float PeakTimePlusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:295
Algorithm for generating space points from hits.
std::set< art::Ptr< recob::Hit > > HitList
void clear()
Definition: PtrVector.h:533
void SetValidity(bool Validity)
Definition: Seed.cxx:74
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
void GetDirection(double *Dir, double *Err) const
Definition: Seed.cxx:80
std::vector< double > fWireZeroOffset
art framework interface to geometry description
double GetLength() const
Definition: Seed.cxx:132
void SetDirection(double *Dir, double *Err)
Definition: Seed.cxx:112
std::vector< double > fMaxViewRMS