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