LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
BezierTrackerAlgorithm.cxx
Go to the documentation of this file.
1 //
2 // Name: BezierTrackerAlgorithm.cxx
3 //
4 // Purpose: Implementation file for module BezierTrackerAlgorithm.
5 //
6 // Ben Jones, MIT
7 //
8 
9 #include <vector>
10 #include "messagefacility/MessageLogger/MessageLogger.h"
11 #include "larreco/RecoAlg/BezierTrackerAlgorithm.h"
12 #include "lardataobj/RecoBase/Hit.h"
13 #include "lardataobj/RecoBase/Seed.h"
14 #include "lardataobj/RecoBase/Vertex.h"
15 #include "lardataobj/RecoBase/Track.h"
16 #include "lardataobj/RecoBase/SpacePoint.h"
17 #include "larreco/Deprecated/BezierTrack.h"
18 #include "larreco/Deprecated/BezierCurveHelper.h"
19 #include "lardata/Utilities/AssociationUtil.h"
20 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
21 
22 
23 #include "art/Framework/Principal/Event.h"
24 #include "art/Framework/Services/Registry/ServiceHandle.h"
25 
26 #include "TVector3.h"
27 
28 namespace trkf {
29 
30  bool BTrack_SeedCountComparator(std::vector<recob::Seed> const& s1, std::vector<recob::Seed> const& s2);
31 
32 
33  BezierTrackerAlgorithm::BezierTrackerAlgorithm(const fhicl::ParameterSet& pset)
34  {
35 
36  reconfigure(pset);
37  }
38 
39  //------------------------------------------------
40 
41  BezierTrackerAlgorithm::~BezierTrackerAlgorithm()
42  {
43  }
44 
45  //------------------------------------------------
46 
47 
48  void BezierTrackerAlgorithm::reconfigure(fhicl::ParameterSet const& pset)
49  {
50 
51  CalculateGeometricalElements();
52 
53  fTrackJoinAngle = pset.get<double>("TrackJoinAngle");
54  fDirectJoinDistance= pset.get<double>("DirectJoinDistance");
55 
56  fVertexImpactThreshold = pset.get<double>("VertexImpactThreshold");
57  fVertexExtrapDistance = pset.get<double>("VertexExtrapDistance");
58 
59  fOccupancyThresh = pset.get<double>("OccupancyThresh");
60  fTrackResolution = pset.get<double>("TrackResolution");
61  fOverlapCut = pset.get<double>("OverlapCut");
62 
63 
64 
65  fTheSeedFinder = new SeedFinderAlgorithm(pset.get<fhicl::ParameterSet>("SeedFinder"));
66 
67  }
68 
69 
70  //-----------------------------------------------
71  std::vector<trkf::BezierTrack> BezierTrackerAlgorithm::MakeTracks(std::vector< std::vector<art::PtrVector<recob::Hit> > >& SortedHits, std::vector<art::PtrVector<recob::Hit> >& HitAssocs, std::vector<std::vector<std::pair<int,int> > >& ClustersUsed)
72  {
73 
74  std::vector<trkf::BezierTrack> ReturnVector;
75 
76  size_t UEntries = SortedHits[0].size();
77  size_t VEntries = SortedHits[1].size();
78  size_t WEntries = SortedHits[2].size();
79 
80  // Get seeds from the hit collection
81 
82  std::vector< std::vector< std::vector<int> > > OrgHitsU(UEntries);
83  std::vector< std::vector< std::vector<int> > > OrgHitsV(VEntries);
84  std::vector< std::vector< std::vector<int> > > OrgHitsW(WEntries);
85 
86 
87  art::ServiceHandle<geo::Geometry> geo;
88  size_t NChannels = geo->Nchannels();
89 
90  std::vector<double> MinUTime(UEntries,10000);
91  std::vector<double> MaxUTime(UEntries,0);
92  std::vector<double> MinVTime(VEntries,10000);
93  std::vector<double> MaxVTime(VEntries,0);
94  std::vector<double> MinWTime(WEntries,10000);
95  std::vector<double> MaxWTime(WEntries,0);
96 
97  std::vector<int> MinUWire(UEntries,10000);
98  std::vector<int> MaxUWire(UEntries,0);
99  std::vector<int> MinVWire(VEntries,10000);
100  std::vector<int> MaxVWire(VEntries,0);
101  std::vector<int> MinWWire(WEntries,10000);
102  std::vector<int> MaxWWire(WEntries,0);
103 
104 
105  for(size_t nU=0; nU!=UEntries; ++nU)
106  {
107  OrgHitsU[nU].resize(NChannels);
108  for(size_t iH=0; iH!=SortedHits[0][nU].size(); ++iH)
109  {
110  OrgHitsU[nU][SortedHits[0][nU][iH]->Channel()].push_back(iH);
111  double Time = SortedHits[0][nU][iH]->PeakTime();
112  if(Time<MinUTime[nU]) MinUTime[nU]=Time;
113  if(Time>MaxUTime[nU]) MaxUTime[nU]=Time;
114  int Wire = SortedHits[0][nU][iH]->WireID().Wire;
115  if(Wire<MinUWire[nU]) MinUWire[nU]=Wire;
116  if(Wire>MaxUWire[nU]) MaxUWire[nU]=Wire;
117 
118  }
119  MaxUTime[nU] -= fPlaneTimeOffsets[0];
120  MinUTime[nU] -= fPlaneTimeOffsets[0];
121  }
122 
123  for(size_t nV=0; nV!=VEntries; ++nV)
124  {
125  OrgHitsV[nV].resize(NChannels);
126  for(size_t iH=0; iH!=SortedHits[1][nV].size(); ++iH)
127  {
128  OrgHitsV[nV][SortedHits[1][nV][iH]->Channel()].push_back(iH);
129  double Time = SortedHits[1][nV][iH]->PeakTime();
130  if(Time<MinVTime[nV]) MinVTime[nV]=Time;
131  if(Time>MaxVTime[nV]) MaxVTime[nV]=Time;
132  int Wire = SortedHits[1][nV][iH]->WireID().Wire;
133  if(Wire<MinVWire[nV]) MinVWire[nV]=Wire;
134  if(Wire>MaxVWire[nV]) MaxVWire[nV]=Wire;
135  }
136  MaxVTime[nV] -= fPlaneTimeOffsets[1];
137  MinVTime[nV] -= fPlaneTimeOffsets[1];
138  }
139 
140  for(size_t nW=0; nW!=WEntries; ++nW)
141  {
142  OrgHitsW[nW].resize(NChannels);
143  for(size_t iH=0; iH!=SortedHits[2][nW].size(); ++iH)
144  {
145  OrgHitsW[nW][SortedHits[2][nW][iH]->Channel()].push_back(iH);
146  double Time = SortedHits[2][nW][iH]->PeakTime();
147  if(Time<MinWTime[nW]) MinWTime[nW]=Time;
148  if(Time>MaxWTime[nW]) MaxWTime[nW]=Time;
149  int Wire = SortedHits[2][nW][iH]->WireID().Wire;
150  if(Wire<MinWWire[nW]) MinWWire[nW]=Wire;
151  if(Wire>MaxWWire[nW]) MaxWWire[nW]=Wire;
152  }
153  MaxWTime[nW] -= fPlaneTimeOffsets[2];
154  MinWTime[nW] -= fPlaneTimeOffsets[2];
155  }
156 
157  for(size_t nU =0; nU != UEntries; ++nU)
158  for(size_t nV =0; nV != VEntries; ++nV)
159  for(size_t nW =0; nW != WEntries; ++nW)
160 
161  {
162  // First check if there is some chance of overlap. If not, move on.
163 
164  // check time region
165 
166  if(MaxUTime[nU] < MinWTime[nW] ) continue;
167  if(MaxUTime[nU] < MinVTime[nV] ) continue;
168  if(MaxVTime[nV] < MinUTime[nU] ) continue;
169  if(MaxVTime[nV] < MinWTime[nW] ) continue;
170  if(MaxWTime[nW] < MinUTime[nU] ) continue;
171  if(MaxWTime[nW] < MinVTime[nV] ) continue;
172 
173  // We won't need to include hits outside this time range
174  double SpacePointRes = GetSeedFinderAlgorithm()->GetSpacePointAlg()->maxDT();
175  double MinTime = std::max(std::max(MinUTime[nU],MinVTime[nV]),MinWTime[nW])-SpacePointRes;
176  double MaxTime = std::min(std::min(MaxUTime[nU],MaxVTime[nV]),MaxWTime[nW])+SpacePointRes;
177 
178  // check wire region
179  double y=0,z=0;
180  geo->IntersectionPoint(MaxUWire[nU],MaxVWire[nV],0,1,0,0,y,z);
181  double MaxUVOnW = ( z - fWireZeroOffset[2] ) / fPitches[2];
182  if(MaxUVOnW < MinWWire[nW]) continue;
183  geo->IntersectionPoint(MinUWire[nU],MinVWire[nV],0,1,0,0,y,z);
184  double MinUVOnW = ( z - fWireZeroOffset[2] ) / fPitches[2];
185  if(MinUVOnW > MaxWWire[nW]) continue;
186 
187 
188 
189  art::PtrVector<recob::Hit> HitsFlat;
190 
191  double EffMinTime = MinTime + fPlaneTimeOffsets[0];
192  double EffMaxTime = MaxTime + fPlaneTimeOffsets[0];
193 
194  for(size_t i=0; i!=SortedHits[0][nU].size(); ++i)
195  if((SortedHits[0][nU][i]->PeakTimePlusRMS(+1.) > EffMinTime)
196  &&(SortedHits[0][nU][i]->PeakTimePlusRMS(-1.)<EffMaxTime))
197  HitsFlat.push_back(SortedHits[0][nU][i]);
198 
199  EffMinTime = MinTime + fPlaneTimeOffsets[1];
200  EffMaxTime = MaxTime + fPlaneTimeOffsets[1];
201 
202  for(size_t i=0; i!=SortedHits[1][nV].size(); ++i)
203  if((SortedHits[1][nV][i]->PeakTimePlusRMS(+1.)>EffMinTime)
204  &&(SortedHits[1][nV][i]->PeakTimePlusRMS(-1.)<EffMaxTime))
205  HitsFlat.push_back(SortedHits[1][nV][i]);
206 
207  EffMinTime = MinTime + fPlaneTimeOffsets[2];
208  EffMaxTime = MaxTime + fPlaneTimeOffsets[2];
209 
210  for(size_t i=0; i!=SortedHits[2][nW].size(); ++i)
211  if((SortedHits[2][nW][i]->PeakTimePlusRMS(+1.)>EffMinTime)
212  &&(SortedHits[2][nW][i]->PeakTimePlusRMS(-1.)<EffMaxTime))
213  HitsFlat.push_back(SortedHits[2][nW][i]);
214 
215  std::vector<art::PtrVector<recob::Hit> > HitsPerSeed;
216 
217  std::vector<recob::Seed> SeedsThisCombo =
218  GetSeedFinderAlgorithm()->GetSeedsFromUnSortedHits(HitsFlat, HitsPerSeed);
219 
220 
221  if(SeedsThisCombo.size()>0)
222  {
223 
224 
225  std::vector<art::PtrVector<recob::Hit>* > HitStruct(3);
226  HitStruct[0] = & SortedHits[0].at(nU);
227  HitStruct[1] = & SortedHits[1].at(nV);
228  HitStruct[2] = & SortedHits[2].at(nW);
229 
230  std::vector<std::vector< std::vector<int> >* > OrgHits(3);
231  OrgHits[0] = & OrgHitsU[nU];
232  OrgHits[1] = & OrgHitsV[nV];
233  OrgHits[2] = & OrgHitsW[nW];
234 
235 
236 
237  std::vector<std::vector<std::vector<int > > > HitsPerTrack;
238 
239  std::vector<std::vector<recob::Seed> > SeedsForTracks = OrganizeSeedsIntoTracks(SeedsThisCombo, HitStruct, HitsPerSeed, OrgHits, HitsPerTrack);
240  for(size_t iTrack=0; iTrack!=SeedsForTracks.size(); ++iTrack)
241  {
242  std::vector<std::pair<int, int> > ClustersThisTrack;
243  ClustersThisTrack.push_back(std::pair<int,int>(0,nU));
244  ClustersThisTrack.push_back(std::pair<int,int>(1,nV));
245  ClustersThisTrack.push_back(std::pair<int,int>(2,nW));
246 
247  ClustersUsed.push_back(ClustersThisTrack);
248 
249  ReturnVector.push_back(trkf::BezierTrack(SeedsForTracks.at(iTrack)));
250  // Fill up the hit vector
251  art::PtrVector<recob::Hit> HitsForThisTrack;
252 
253  for(size_t iH=0; iH!=HitsPerTrack.at(iTrack)[0].size(); ++iH)
254  HitsForThisTrack.push_back(SortedHits[geo::kU][nU][HitsPerTrack.at(iTrack)[0].at(iH)]);
255 
256  for(size_t iH=0; iH!=HitsPerTrack.at(iTrack)[1].size(); ++iH)
257  HitsForThisTrack.push_back(SortedHits[geo::kV][nV][HitsPerTrack.at(iTrack)[1].at(iH)]);
258 
259  for(size_t iH=0; iH!=HitsPerTrack.at(iTrack)[2].size(); ++iH)
260  HitsForThisTrack.push_back(SortedHits[geo::kW][nW][HitsPerTrack.at(iTrack)[2].at(iH)]);
261 
262  HitAssocs.push_back(HitsForThisTrack);
263 
264  }
265 
266  }
267 
268  }
269 
270  return ReturnVector;
271  }
272 
273 
274  //------------------------------------------------
275 
276  void BezierTrackerAlgorithm::FilterOverlapTracks(std::vector<trkf::BezierTrack>& BTracks, std::vector<art::PtrVector<recob::Hit> > & HitVecs, std::vector<std::vector<std::pair<int, int> > >& ClustersUsed)
277  {
278 
279  std::map<int, bool> ToErase;
280 
281  std::vector<double> Lengths;
282  std::vector<TVector3> End1Points, End2Points;
283 
284 
285  for(size_t i=0; i!=BTracks.size(); ++i)
286  {
287  End1Points.push_back(BTracks.at(i).GetTrackPointV(0));
288  End2Points.push_back(BTracks.at(i).GetTrackPointV(1));
289  Lengths.push_back(BTracks.at(i).GetLength());
290  }
291 
292  for(size_t t1=0; t1!=BTracks.size(); ++t1)
293  for(size_t t2=0; t2!=BTracks.size(); ++t2)
294  {
295  if(t1!=t2)
296  {
297  if(Lengths.at(t1)>Lengths.at(t2))
298  {
299  double s1,d1,s2,d2;
300  BTracks.at(t1).GetClosestApproach(End1Points.at(t2),s1,d1);
301  BTracks.at(t1).GetClosestApproach(End2Points.at(t2),s2,d2);
302  if((s1>0.01)&&(s1<0.99)&&(s2>0.01)&&(s2<0.99)&&
303  (d1<fTrackResolution)&&(d2<fTrackResolution))
304  {
305  // Track t2 is a fraud - toss it
306  ToErase[t2]=true;
307  }
308  }
309  else
310  {
311  double s1,d1,s2,d2;
312  BTracks.at(t2).GetClosestApproach(End1Points.at(t1),s1,d1);
313  BTracks.at(t2).GetClosestApproach(End2Points.at(t1),s2,d2);
314  if((s1>0.01)&&(s1<0.99)&&(s2>0.01)&&(s2<0.99)&&
315  (d1<fTrackResolution)&&(d2<fTrackResolution))
316  {
317  // Track t2 is a fraud - toss it
318  ToErase[t1]=true;
319  }
320  }
321  }
322  }
323 
324  // Remove the tracks we marked for deletion
325  for(int i=BTracks.size()-1; i >= 0; --i)
326  {
327  if(ToErase[i])
328  {
329  BTracks.erase(BTracks.begin()+i);
330  HitVecs.erase(HitVecs.begin()+i);
331  ClustersUsed.erase(ClustersUsed.begin()+i);
332  }
333  }
334  }
335 
336 
337 
338  //------------------------------------------------
339 
340  void BezierTrackerAlgorithm::MakeOverlapJoins(std::vector<trkf::BezierTrack>& BTracks, std::vector<art::PtrVector<recob::Hit> > & HitVecs, std::vector<std::vector<std::pair<int, int> > >& ClustersUsed)
341  {
342 
343  std::vector<double> Lengths;
344  std::vector<TVector3> End1Directions, End1Points, End2Directions, End2Points;
345  for(size_t i=0; i!=BTracks.size(); ++i)
346  {
347  Lengths.push_back(BTracks.at(i).GetLength());
348  End1Points.push_back(BTracks.at(i).GetTrackPointV(0));
349  End2Points.push_back(BTracks.at(i).GetTrackPointV(1));
350  End1Directions.push_back(BTracks.at(i).GetTrackDirectionV(0).Unit());
351  End2Directions.push_back(BTracks.at(i).GetTrackDirectionV(1).Unit());
352  }
353 
354  int LoopStatus=1;
355  while(LoopStatus!=0)
356  {
357  LoopStatus=0;
358  for(size_t t1=0; t1!=BTracks.size(); ++t1)
359  {
360  for(size_t t2=0; t2!=BTracks.size(); ++t2)
361  {
362  if(t1!=t2)
363  {
364  double s12,d12;
365  double s21,d21;
366  BTracks.at(t1).GetClosestApproach(End1Points.at(t2),s12,d12);
367  BTracks.at(t2).GetClosestApproach(End2Points.at(t1),s21,d21);
368 
369  if((s12>0.0001)&&(s12<0.999)
370  &&(s21>0.0001)&&(s21<0.999)
371  &&(d12<fTrackResolution)
372  &&(d21<fTrackResolution))
373  {
374  // In this situation we want to join track 1 onto the end of track 2
375  trkf::BezierTrack Part1 = BTracks.at(t1).GetPartialTrack(0,s21);
376  trkf::BezierTrack Part2 = BTracks.at(t2).GetPartialTrack(s12,1);
377  BTracks.at(t1) = JoinTracks(Part1, Part2);
378  BTracks.erase(BTracks.begin()+t2);
379 
380  // Add the pointervectors into 1 and remove 2
381  AddPtrVectors(HitVecs.at(t1), HitVecs.at(t2));
382  ClustersUsed.at(t1).insert(ClustersUsed.at(t1).begin(),
383  ClustersUsed.at(t2).begin(),
384  ClustersUsed.at(t2).end());
385 
386  HitVecs.erase(HitVecs.begin()+t2);
387  ClustersUsed.erase(ClustersUsed.begin()+t2);
388 
389  LoopStatus=3;
390  }
391  }
392  if(LoopStatus==3) break;
393  }
394  if(LoopStatus==3) break;
395  }
396  }
397  }
398 
399 
400  //--------------------------------------------
401  void BezierTrackerAlgorithm::SortTracksByLength(std::vector<trkf::BezierTrack>& BTracks, std::vector<art::PtrVector<recob::Hit> > & HitVecs, std::vector<std::vector<std::pair<int, int> > >& ClustersUsed)
402  {
403 
404  std::vector<double> Lengths;
405  std::vector<int> SortedOrder;
406  std::vector<bool> Counted(BTracks.size());
407 
408  for(size_t i=0; i!=BTracks.size(); ++i)
409  {
410  Lengths.push_back(BTracks.at(i).GetLength());
411  if(BTracks.at(i).GetTrackDirectionV(0.0).Dot(fZDir)<0.0) BTracks.at(i)=BTracks.at(i).Reverse();
412  }
413  while(SortedOrder.size()<BTracks.size())
414  {
415  int Longest=-1;
416  for(size_t i=0; i!=BTracks.size(); ++i)
417  {
418  if(!Counted[i])
419  {
420  if(Longest==-1) Longest=i;
421  else
422  if(Lengths.at(i)>Lengths.at(Longest)) Longest=i;
423  }
424  }
425  Counted[Longest]=true;
426  SortedOrder.push_back(Longest);
427  }
428 
429  std::vector<trkf::BezierTrack> BTracksNew;
430  std::vector<art::PtrVector<recob::Hit> > HitVecsNew;
431  std::vector<std::vector<std::pair<int, int> > > ClustersUsedNew;
432  for(size_t i=0; i!=SortedOrder.size(); ++i)
433  {
434  BTracksNew.push_back(BTracks.at(SortedOrder.at(i)));
435  HitVecsNew.push_back(HitVecs.at(SortedOrder.at(i)));
436  ClustersUsedNew.push_back(ClustersUsed.at(SortedOrder.at(i)));
437  }
438 
439  BTracks=BTracksNew;
440  HitVecs=HitVecsNew;
441  ClustersUsed=ClustersUsedNew;
442  }
443 
444 
445  //--------------------------------------------
446 
447 
448  void BezierTrackerAlgorithm::MakeDirectJoins(std::vector<trkf::BezierTrack>& BTracks, std::vector<art::PtrVector<recob::Hit> > & HitVecs, std::vector<std::vector<std::pair<int, int> > >& ClustersUsed)
449  {
450 
451  int LoopStatus=1;
452 
453  // 0 : Quit loop
454  // 1 : Keep looping
455  // 2 : Deleted a seed : start again from the beginning
456 
457  while(LoopStatus>0)
458  {
459  LoopStatus = 0;
460  for(size_t t1=0; t1!=BTracks.size(); ++t1)
461  {
462  for(size_t t2=0; t2!=BTracks.size(); ++t2)
463  {
464  if(t1==t2) continue;
465 
466  TVector3 End1Dir_i = BTracks.at(t1).GetTrackDirectionV(1.0).Unit();
467  TVector3 End0Dir_j = BTracks.at(t2).GetTrackDirectionV(0.0).Unit();
468 
469  TVector3 End0Point_i = BTracks.at(t1).GetTrackPointV(0.0);
470  TVector3 End0Point_j = BTracks.at(t2).GetTrackPointV(0.0);
471  TVector3 End1Point_i = BTracks.at(t1).GetTrackPointV(1.0);
472  TVector3 End1Point_j = BTracks.at(t2).GetTrackPointV(1.0);
473 
474 
475 
476  double JoinAngle = End1Dir_i.Angle(End0Dir_j);
477  double JoinSign = (End0Point_j-End1Point_i).Dot(End1Dir_i);
478  double Colinearity = ((End0Point_j-End1Point_i) - End1Dir_i*((End0Point_j-End1Point_i).Dot(End1Dir_i))).Mag();
479 
480 
481  if((JoinAngle<fTrackJoinAngle)&&(JoinSign>0)&&(Colinearity<fTrackResolution))
482  {
483  double EndSep10 = (End1Point_i-End0Point_j).Mag();
484  double EndSep01 = (End0Point_i-End1Point_j).Mag();
485  double EndSep11 = (End1Point_i-End1Point_j).Mag();
486  double EndSep00 = (End0Point_i-End0Point_j).Mag();
487 
488 
489  if( (EndSep10 < fDirectJoinDistance)
490  &&(EndSep10 < EndSep11)
491  &&(EndSep10 < EndSep00)
492  &&(EndSep10 < EndSep01) )
493  {
494 
495  BTracks.at(t1) = JoinTracks(BTracks.at(t1), BTracks.at(t2));
496  BTracks.erase(BTracks.begin()+t2);
497 
498  // Add the pointervectors into 1 and remove 2
499  AddPtrVectors(HitVecs.at(t1), HitVecs.at(t2));
500  ClustersUsed.at(t1).insert(ClustersUsed.at(t1).begin(),
501  ClustersUsed.at(t2).begin(),
502  ClustersUsed.at(t2).end());
503 
504  HitVecs.erase(HitVecs.begin()+t2);
505  ClustersUsed.erase(ClustersUsed.begin()+t2);
506  LoopStatus=3;
507  }
508  }
509  if(LoopStatus==3) break;
510  }
511  if(LoopStatus==3) break;
512  }
513  }
514  }
515 
516 
517  //-----------------------------------
518 
519  trkf::BezierTrack BezierTrackerAlgorithm::JoinTracks(trkf::BezierTrack& BT1, trkf::BezierTrack& BT2)
520  {
521 
522  std::vector<recob::Seed> SeedCol1 = BT1.GetSeedVector();
523  std::vector<recob::Seed> SeedCol2 = BT2.GetSeedVector();
524 
525  SeedCol1.insert(SeedCol1.end(),SeedCol2.begin(),SeedCol2.end());
526 
527  return trkf::BezierTrack(SeedCol1);
528  }
529 
530 
531 
532  //-----------------------------------------------
533  void BezierTrackerAlgorithm::AddPtrVectors(art::PtrVector<recob::Hit>& Receiver, art::PtrVector<recob::Hit> const & ToAdd)
534  {
535  for(size_t i=0; i!=ToAdd.size(); ++i)
536  {
537  Receiver.push_back(ToAdd.at(i));
538  }
539  }
540 
541 
542 
543  //-----------------------------------------------------
544 
545  std::vector<std::vector< recob::Seed > > BezierTrackerAlgorithm::OrganizeSeedsIntoTracks(std::vector<recob::Seed > & AllSeeds, std::vector<art::PtrVector<recob::Hit>*>& AllHits, std::vector<art::PtrVector<recob::Hit> >& WhichHitsPerSeed, std::vector<std::vector< std::vector<int> >* >& OrgHits, std::vector<std::vector<std::vector<int> > >& HitsWithTracks )
546  {
547  mf::LogVerbatim("BezierTrackerAlgorithm")<<"Organizing seed collection into track collections ";
548 
549  std::vector<std::vector<uint32_t> > MaxSeedChanPerView(3);
550  std::vector<std::vector<uint32_t> > MinSeedChanPerView(3);
551 
552  // This vector will keep track of which hits we are still allowed
553  // to use. Key :
554  // 0 - hit available
555  // 1 - hit taken by an extrapolated bezier curve
556  // 2 - hit taken by a seed on the track
557 
558  std::vector<std::vector<int> > HitStatus(3);
559  for(size_t n=0; n!=3; ++n) HitStatus.at(n).resize(AllHits[n]->size());
560 
561 
562  // This map keeps track of the indices of the hits (in the AllHits vector)
563  // which go with each seed, in each view
564  // SeedToHitMap[per_view][per_seed][1_to_n_with_this_seed] = index
565 
566  std::vector<std::vector<std::vector<int> > > SeedToHitMap(3);
567 
568  // Fill Seed <--> Hit lookups
569  for(size_t iSeed=0; iSeed!=AllSeeds.size(); ++iSeed)
570  {
571  // first make sure we have at least one hit in each view.
572 
573  bool CompleteSet = false;
574  std::vector<bool> AtLeastOneHitThisView(3,false);
575  for(size_t iHit=0; iHit!=WhichHitsPerSeed.at(iSeed).size(); ++iHit)
576  {
577  if(!AtLeastOneHitThisView[WhichHitsPerSeed.at(iSeed).at(iHit)->View()])
578  {
579  AtLeastOneHitThisView[WhichHitsPerSeed.at(iSeed).at(iHit)->View()] = true;
580  if(AtLeastOneHitThisView[0]&&AtLeastOneHitThisView[1]&&AtLeastOneHitThisView[2])
581  {
582  CompleteSet = true;
583  break;
584  }
585  }
586  }
587  if(!CompleteSet)
588  {
589  // If we don't have at least one hit in each view from this
590  // seed, we're not going to be able to use it.
591  // ditch it (and its corresponding hits) from the game.
592  AllSeeds.erase(AllSeeds.begin() + iSeed);
593  WhichHitsPerSeed.erase(WhichHitsPerSeed.begin() + iSeed);
594  --iSeed;
595  continue;
596  }
597  else
598  {
599  // if we didn't throw out that seed, push back a new entry in
600  // the seed <-> hit map vector, and
601  // Find the min and max channels for each seed
602  for(size_t n=0; n!=3; ++n)
603  {
604  SeedToHitMap[n].push_back(std::vector<int>());
605  MinSeedChanPerView[n].push_back(10000);
606  MaxSeedChanPerView[n].push_back(0);
607  }
608  for(size_t iHit=0; iHit!=WhichHitsPerSeed.at(iSeed).size(); ++iHit)
609  {
610  uint32_t Channel = WhichHitsPerSeed.at(iSeed).at(iHit)->Channel();
611  geo::View_t View = WhichHitsPerSeed.at(iSeed).at(iHit)->View();
612 
613  int ViewID=0;
614  if(View==geo::kU) ViewID=0;
615  else if(View==geo::kV) ViewID=1;
616  else if(View==geo::kW) ViewID=2;
617 
618  double eta = 0.01;
619 
620  for(size_t iH=0; iH!=OrgHits[ViewID]->at(Channel).size();++iH)
621  {
622  double PeakTime1 = AllHits[ViewID]->at(OrgHits[ViewID]->at(Channel).at(iH))->PeakTime();
623  double PeakTime2 = WhichHitsPerSeed.at(iSeed).at(iHit)->PeakTime();
624 
625  if( fabs(PeakTime1-PeakTime2)<eta)
626 
627  SeedToHitMap[ViewID][iSeed].push_back(OrgHits[ViewID]->at(Channel).at(iH));
628  }
629 
630 
631  if(Channel>MaxSeedChanPerView[ViewID][iSeed]) MaxSeedChanPerView[ViewID][iSeed] = Channel;
632  if(Channel<MinSeedChanPerView[ViewID][iSeed]) MinSeedChanPerView[ViewID][iSeed] = Channel;
633 
634  }
635  }
636  }
637 
638  // Now we've set up the book keeping, we get into the real
639  // work of track finding.
640 
641  // This will store our proposed track vectors
642  std::vector<std::vector<recob::Seed > > OrganizedByTrack;
643 
644  TVector3 AverageDir;
645 
646  std::vector<TVector3> VSeedDirs;
647  std::vector<TVector3> VSeedPoss;
648  std::vector<std::vector<int> > SeedSDirs;
649  std::vector<std::vector<uint32_t> > SeedHighSChans;
650  std::vector<std::vector<uint32_t> > SeedLowSChans;
651 
652 
653 
654  for(size_t i=0; i!=AllSeeds.size(); ++i)
655  {
656  double SeedDir[3], SeedPos[3], Err[3];
657 
658  AllSeeds.at(i).GetDirection(SeedDir,Err);
659  AllSeeds.at(i).GetPoint(SeedPos,Err);
660  TVector3 VSeedDir = TVector3(SeedDir[0],SeedDir[1],SeedDir[2]);
661  VSeedDirs.push_back(VSeedDir);
662  TVector3 VSeedPos = TVector3(SeedPos[0],SeedPos[1],SeedPos[2]);
663  VSeedPoss.push_back(VSeedPos);
664 
665  // This is a two directional degeneracy
666  // This is one way to resolve it (maybe not the best)
667 
668  if(VSeedDir.Dot(fZDir)<0)
669  AverageDir -= VSeedDir;
670  else
671  AverageDir += VSeedDir;
672 
673  // Work out which direction the seed points out in each plane
674  std::vector<double> TimeDir, WireDir;
675 
676  GetSeedDirProjected(AllSeeds.at(i), WireDir, TimeDir);
677 
678  SeedSDirs.push_back(std::vector<int>(3));
679  SeedHighSChans.push_back(std::vector<uint32_t>(3));
680  SeedLowSChans.push_back(std::vector<uint32_t>(3));
681 
682 
683  for(size_t n=0; n!=3; ++n)
684  {
685  if(WireDir[n]>0)
686  {
687  SeedSDirs[i][n]=1;
688  SeedHighSChans[i][n] = MaxSeedChanPerView[n][i];
689  SeedLowSChans[i][n] = MinSeedChanPerView[n][i];
690  }
691  else
692  {
693  SeedSDirs[i][n]=-1;
694  SeedHighSChans[i][n] = MinSeedChanPerView[n][i];
695  SeedLowSChans[i][n] = MaxSeedChanPerView[n][i];
696  }
697  }
698  }
699 
700  AverageDir *= (1.0/AverageDir.Mag());
701 
702  std::vector<double> AngleToAve;
703  std::vector<double> ProjAlongAve;
704 
705  for(size_t i=0; i!=AllSeeds.size(); ++i)
706  {
707  ProjAlongAve.push_back(VSeedPoss.at(i).Dot(AverageDir));
708  AngleToAve.push_back(VSeedDirs.at(i).Dot(AverageDir)/VSeedDirs.at(i).Mag());
709  }
710 
711 
712  // For now we'll only look for 1 long track per cluster combo.
713  // There is room for improvement here if we need it.
714 
715  std::vector<size_t> TrackParts;
716  std::vector<size_t> DoesNotFit;
717  std::vector<size_t> AlreadyUsed;
718  std::vector<size_t> ThrownByHits;
719  std::map<int, bool> SeedNotAvailable;
720 
721  // This loop keeps running whilst there are still tracks to make
722  while((AlreadyUsed.size()+ThrownByHits.size())<AllSeeds.size())
723  {
724 
725  // This loop keeps running whilst there are still seeds
726  // to check for this track
727 
728 
729  while((ThrownByHits.size()+TrackParts.size()+DoesNotFit.size()+AlreadyUsed.size())<AllSeeds.size())
730  {
731 
732  // First find the lowest track proj remaining
733  int LowestLeft=-1;
734  for(size_t i=0; i!=AllSeeds.size(); ++i)
735  {
736  if(!SeedNotAvailable[i])
737  {
738  if(LowestLeft < 0) LowestLeft = i;
739  if(ProjAlongAve.at(i)<ProjAlongAve.at(LowestLeft))
740  LowestLeft = i;
741  }
742  }
743 
744  std::vector<std::vector<int> > InterpHits(3);
745 
746  if(TrackParts.size()==0)
747  {
748  // If there are no seeds on the track, accept this seed
749  // as the first one to start running with.
750  TrackParts.push_back(LowestLeft);
751  SeedNotAvailable[LowestLeft]=true;
752 
753  // Mark the hits that go with this seed as taken
754  for(size_t View=0; View!=3; ++View)
755  for(size_t iH=0; iH!=SeedToHitMap[View][LowestLeft].size();++iH)
756  HitStatus[View][SeedToHitMap[View][LowestLeft][iH]] = 2;
757 
758  }
759  else if(TrackParts.size()>0)
760  {
761  int a = TrackParts.at(TrackParts.size()-1);
762  int b = LowestLeft;
763 
764  double Angle = AllSeeds.at(a).GetAngle(AllSeeds.at(b));
765  double ProjDisc = AllSeeds.at(a).GetProjAngleDiscrepancy(AllSeeds.at(b));
766  int PointingSign_ab = AllSeeds.at(a).GetPointingSign(AllSeeds.at(b));
767  int PointingSign_ba = AllSeeds.at(b).GetPointingSign(AllSeeds.at(a));
768  int PointingSign_ac = 0;
769 
770  if(TrackParts.size()>1) PointingSign_ac = AllSeeds.at(a).GetPointingSign(AllSeeds.at(TrackParts.size()-2));
771 
772  bool SeedFits=false;
773 
774  // Check seed join conditions (no hits yet)
775  bool JoinAngleOK = fabs(Angle)<fTrackJoinAngle;
776  bool PointingSignOK = (PointingSign_ac * PointingSign_ab)<0.5;
777  bool ProjDiscOK = false;
778  if(PointingSign_ab > 0)
779  {
780  ProjDiscOK = fabs(ProjDisc) < fTrackJoinAngle;
781  }
782  else
783  {
784  ProjDiscOK = fabs(ProjDisc - 3.142) < fTrackJoinAngle;
785  }
786 
787  if(JoinAngleOK && PointingSignOK && ProjDiscOK)
788  {
789  std::vector<uint32_t> EndChan_a(3), EndChan_b(3), LowChan(3), HighChan(3);
790  std::vector<std::vector<int> > HitsUnderConsideration(3);
791 
792  for(size_t n=0; n!=3; ++n)
793  {
794  if(PointingSign_ab > 0)
795  EndChan_a[n] = SeedHighSChans[a][n];
796  else
797  EndChan_a[n] = SeedLowSChans[a][n];
798 
799  if(PointingSign_ba > 0)
800  EndChan_b[n] = SeedHighSChans[b][n];
801  else
802  EndChan_b[n] = SeedLowSChans[b][n];
803 
804  LowChan[n] = std::min(EndChan_a[n], EndChan_b[n]) ;
805  HighChan[n] = std::max(EndChan_a[n], EndChan_b[n]) ;
806  }
807 
808  std::vector<std::vector<int> > TheseHits(3);
809 
810 
811  SeedFits = EvaluateOccupancy(AllSeeds.at(a), AllSeeds.at(b), fTrackResolution, AllHits, OrgHits, LowChan, HighChan, HitStatus, TheseHits);
812 
813  if(SeedFits)
814  {
815  // This seed is good! Store it on the track
816  TrackParts.push_back(LowestLeft);
817  SeedNotAvailable[LowestLeft]=true;
818 
819  // Remove the hits that went with this seed from consideration
820  for(size_t View=0; View!=3; ++View)
821  for(size_t iH=0; iH!=SeedToHitMap[View][LowestLeft].size();++iH)
822  HitStatus[View][SeedToHitMap[View][LowestLeft][iH]] = 2;
823  // Remove the hits that were interpolated through from consideration
824  for(size_t n=0; n!=3; ++n)
825  for(size_t i=0; i!=InterpHits[n].size(); ++i)
826  {
827  HitStatus[n][InterpHits[n][i]] = 1;
828  }
829  // Kick out any seeds which are not invalidated by their
830  // hits being stolen
831  for(size_t j=0; j!=AllSeeds.size(); ++j)
832  {
833  if(!SeedNotAvailable[j])
834  {
835  double OverlapFrac[3];
836  for(size_t n=0; n!=3; ++n)
837  {
838  int TotalInSeed=0;
839  int TotalRejected=0;
840  for(size_t iH=0; iH!=SeedToHitMap[n][j].size(); ++iH)
841  {
842  TotalInSeed++;
843  if(HitStatus[n][SeedToHitMap[n][j][iH]]!=0)
844  TotalRejected++;
845  }
846  OverlapFrac[n] = float(TotalRejected)/float(TotalInSeed);
847  }
848  for(size_t n=0; n!=3; ++n)
849  {
850  size_t n1=(n+1)%3; size_t n2=(n+2)%3;
851  if( (OverlapFrac[n]<=OverlapFrac[n1]) &&
852  (OverlapFrac[n]<=OverlapFrac[n2]) )
853  {
854  if(OverlapFrac[n] > fOverlapCut)
855  {
856  SeedNotAvailable[j]=true;
857  ThrownByHits.push_back(j);
858  break;
859  }
860  }
861 
862  }// End view loop (seed tossing)
863  } // End checking seed availability (seed tossing)
864  } // End seed loop (seed tossing)
865  } // End if positive seed fit
866  }
867  if(!SeedFits)
868  {
869  // If seed doesn't fit, mark it as such, and make it
870  // temporarily unavailable.
871  DoesNotFit.push_back(LowestLeft);
872  SeedNotAvailable[LowestLeft]=true;
873  }
874  } // End if(not first seed)
875  } // End loop over all remaining seeds
876  if(TrackParts.size()>0)
877  {
878  mf::LogVerbatim("BezierTrackerAlgorithm")<<" Prototrack runs through seeds: " ;
879  std::vector<recob::Seed> TheTrackSeeds;
880 
881  TVector3 AveTrackDir =
882  VSeedPoss.at(TrackParts.at(TrackParts.size()-1))
883  - VSeedPoss.at(TrackParts.at(0));
884 
885  for(size_t i=0; i!=TrackParts.size(); ++i)
886  {
887  mf::LogVerbatim("BezierTrackerAlgorithm")<<" " << TrackParts.at(i)<<", ";
888  if(VSeedDirs.at(TrackParts.at(i)).Dot(AveTrackDir)>0)
889  TheTrackSeeds.push_back(AllSeeds.at(TrackParts.at(i)));
890  else
891  TheTrackSeeds.push_back(AllSeeds.at(TrackParts.at(i)).Reverse());
892  AlreadyUsed.push_back(TrackParts.at(i));
893  }
894  OrganizedByTrack.push_back(TheTrackSeeds);
895 
896  HitsWithTracks.push_back(std::vector<std::vector<int> >(3));
897  for(size_t n=0; n!=3; ++n)
898  for(size_t iH=0; iH!=HitStatus[n].size(); ++iH)
899  {
900  if((HitStatus[n][iH]==1)||(HitStatus[n][iH]==2))
901  {
902  HitStatus[n][iH]=3;
903  HitsWithTracks[HitsWithTracks.size()-1][n].push_back(iH);
904  }
905  }
906  }
907 
908 
909  // Reset the count of unavailable seeds
910  SeedNotAvailable.clear();
911 
912  // Anything which is either part of a track,
913  // or which has had its hits used, is not fair game.
914  for(size_t j=0; j!=AlreadyUsed.size(); ++j)
915  {
916  SeedNotAvailable[AlreadyUsed[j]]=true;
917  }
918 
919  // For all the seeds which didn't fit, see how many
920  // are now ruled out by hit overlaps
921  for(size_t i=0; i!=DoesNotFit.size(); ++i)
922  {
923  size_t j=DoesNotFit[i];
924  double OverlapFrac[3];
925  for(size_t n=0; n!=3; ++n)
926  {
927  int TotalInSeed=0;
928  int TotalRejected=0;
929  for(size_t iH=0; iH!=SeedToHitMap[n][j].size(); ++iH)
930  {
931  TotalInSeed++;
932  if(HitStatus[n][SeedToHitMap[n][j][iH]]!=0)
933  TotalRejected++;
934  }
935  OverlapFrac[n] = float(TotalRejected)/float(TotalInSeed);
936  }
937  for(size_t n=0; n!=3; ++n)
938  {
939  size_t n1=(n+1)%3; size_t n2=(n+2)%3;
940  if( (OverlapFrac[n]<=OverlapFrac[n1]) &&
941  (OverlapFrac[n]<=OverlapFrac[n2]) )
942  {
943  if(OverlapFrac[n] > fOverlapCut)
944  {
945  SeedNotAvailable[j]=true;
946  ThrownByHits.push_back(j);
947  break;
948  }
949  }
950  }// End view loop (seed tossing)
951  } // Loop over DoesNotFit
952 
953 
954  for(size_t i=0; i!=ThrownByHits.size(); ++i)
955  SeedNotAvailable[ThrownByHits[i]]=true;
956 
957  // Ready to go around again
958 
959  TrackParts.clear();
960  DoesNotFit.clear();
961  }
962 
963  return OrganizedByTrack;
964  }
965 
966 
967  //-------------------------------------------------
968 
969  bool BezierTrackerAlgorithm::EvaluateOccupancy(recob::Seed& Seed1, recob::Seed& Seed2, double dThresh, std::vector<art::PtrVector<recob::Hit>*>& AllHits, std::vector<std::vector< std::vector<int> >* >& OrgHits, std::vector<uint32_t>& LowChan, std::vector<uint32_t>& HighChan, std::vector<std::vector<int> >& HitStatus, std::vector<std::vector<int> >& TheseHits)
970  {
971  const detinfo::DetectorProperties* det = lar::providerFrom<detinfo::DetectorPropertiesService>();
972 
973  int NSteps = 2 * Seed1.GetDistance(Seed2) / dThresh;
974  if(NSteps<2) NSteps=2;
975  BezierCurveHelper bhlp(NSteps);
976 
977  std::vector<TVector3> Pts = bhlp.GetBezierPoints(Seed2, Seed1, NSteps);
978 
979 
980  for(size_t n=0; n!=3; ++n)
981  {
982  if((HighChan[n]-LowChan[n])<=2) continue;
983 
984  size_t EmptyChanLimit = int((1.-fOccupancyThresh)*(HighChan[n]-LowChan[n]-2));
985  size_t NEmptyChannels=0;
986 
987  for(uint32_t iChan = LowChan[n]+1; iChan < (HighChan[n]-1); ++iChan)
988  {
989  bool GotThisChannel=false;
990 
991  for(size_t iH = 0; iH!=OrgHits[n]->at(iChan).size(); ++iH)
992  {
993  if(HitStatus[n][OrgHits[n]->at(iChan).at(iH)]==0)
994  {
995  art::Ptr<recob::Hit> ThisHit = AllHits[n]->at(OrgHits[n]->at(iChan)[iH]);
996 
997  TVector3 WirePoint1 = ( fPitchDir[n]*(fWireZeroOffset[n] + fPitches[n] * ThisHit->WireID().Wire) )+ fXDir * det->ConvertTicksToX(ThisHit->PeakTime(), n, 0, 0);
998 
999  for(size_t ipt=0; ipt!=Pts.size(); ++ipt)
1000  {
1001  float d = ((WirePoint1-Pts.at(ipt))-fWireDir[n]*((WirePoint1-Pts.at(ipt)).Dot(fWireDir[n]))).Mag();
1002  if(d < dThresh)
1003  {
1004  TheseHits[n].push_back(OrgHits[n]->at(iChan).at(iH));
1005  GotThisChannel=true;
1006  break;
1007  }
1008  }
1009  }
1010  }
1011  if(!GotThisChannel) NEmptyChannels++;
1012  if(NEmptyChannels > EmptyChanLimit) return false;
1013  }
1014  }
1015  return true;
1016 
1017  }
1018 
1019 
1020  //-------------------------------------------------
1021 
1022  void BezierTrackerAlgorithm::GetSeedDirProjected(recob::Seed const& TheSeed, std::vector<double>& WireCoord, std::vector<double>& TimeCoord)
1023  {
1024  const detinfo::DetectorProperties* det = lar::providerFrom<detinfo::DetectorPropertiesService>();
1025 
1026  WireCoord.clear();
1027  WireCoord.resize(3);
1028  TimeCoord.clear();
1029  TimeCoord.resize(3);
1030 
1031  double dir[3], err[3];
1032 
1033  TheSeed.GetDirection(dir,err);
1034  TVector3 DirVec(dir[0],dir[1],dir[2]);
1035 
1036  for(size_t n=0; n!=3; ++n)
1037  {
1038  WireCoord[n] = DirVec.Dot(fPitchDir[n]) / fPitches[n];
1039  TimeCoord[n] = det->ConvertXToTicks(DirVec.Dot(fXDir),n,0,0);
1040  }
1041 
1042  }
1043 
1044 
1045 
1046 
1047  //-------------------------------------------------
1048 
1049 
1050  void BezierTrackerAlgorithm::CalculateGeometricalElements()
1051  {
1052  art::ServiceHandle<geo::Geometry> geom;
1053  const detinfo::DetectorProperties* det = lar::providerFrom<detinfo::DetectorPropertiesService>();
1054 
1055 
1056  // Find pitch of each wireplane
1057  fPitches.resize(3);
1058  fPitches.at(0) = fabs(geom->WirePitch(geo::kU));
1059  fPitches.at(1) = fabs(geom->WirePitch(geo::kV));
1060  fPitches.at(2) = fabs(geom->WirePitch(geo::kW));
1061 
1062  fPlaneTimeOffsets.resize(3);
1063  for(size_t n=0; n!=3; ++n)
1064  fPlaneTimeOffsets.at(n) = det->GetXTicksOffset(n,0,0);
1065 
1066  // Setup basis vectors
1067  fXDir = TVector3(1,0,0);
1068  fYDir = TVector3(0,1,0);
1069  fZDir = TVector3(0,0,1);
1070 
1071  fWireDir.resize(3);
1072  fPitchDir.resize(3);
1073  fWireZeroOffset.resize(3);
1074 
1075  double xyzStart1[3], xyzStart2[3];
1076  double xyzEnd1[3], xyzEnd2[3];
1077 
1078  // Calculate wire coordinate systems
1079  for(size_t n=0; n!=3; ++n)
1080  {
1081  geom->WireEndPoints(0,0,n,0,xyzStart1,xyzEnd1);
1082  geom->WireEndPoints(0,0,n,1,xyzStart2,xyzEnd2);
1083  fWireDir[n] = TVector3(xyzEnd1[0] - xyzStart1[0],
1084  xyzEnd1[1] - xyzStart1[1],
1085  xyzEnd1[2] - xyzStart1[2]).Unit();
1086  fPitchDir[n] = fWireDir[n].Cross(fXDir).Unit();
1087  if(fPitchDir[n].Dot(TVector3(xyzEnd2[0] - xyzEnd1[0],
1088  xyzEnd2[1] - xyzEnd1[1],
1089  xyzEnd2[2] - xyzEnd1[2]))<0) fPitchDir[n] = -fPitchDir[n];
1090 
1091  fWireZeroOffset[n] =
1092  xyzEnd1[0]*fPitchDir[n][0] +
1093  xyzEnd1[1]*fPitchDir[n][1] +
1094  xyzEnd1[2]*fPitchDir[n][2];
1095 
1096  }
1097 
1098 
1099  }
1100 
1101  //-----------------------------------------------------
1102 
1103  void BezierTrackerAlgorithm::MakeVertexJoins(std::vector<trkf::BezierTrack>& BTracks, std::vector<recob::Vertex>& Vertices, std::vector<std::vector<int> >& Mapping)
1104  {
1105 
1106  std::vector<TVector3> TrackEnd1s;
1107  std::vector<TVector3> TrackEnd2s;
1108  std::vector<TVector3> TrackDir1s;
1109  std::vector<TVector3> TrackDir2s;
1110 
1111 
1112  for(size_t i=0; i!=BTracks.size(); ++i)
1113  {
1114  TrackEnd1s.push_back( BTracks.at(i).GetTrackPointV(0));
1115  TrackEnd2s.push_back( BTracks.at(i).GetTrackPointV(1));
1116  TrackDir1s.push_back( BTracks.at(i).GetTrackDirectionV(0));
1117  TrackDir2s.push_back( BTracks.at(i).GetTrackDirectionV(1));
1118  }
1119 
1120  std::vector<std::map<int,int> > TracksMeetAtPoints;
1121  std::vector<std::vector<TVector3> > MeetPoints;
1122 
1123  // The connection[i][j] gives impact parameters for track ends
1124  // The sign tells us which end of track j
1125  for(size_t i=0; i!=BTracks.size(); ++i)
1126  {
1127  for(size_t j=0; j!=BTracks.size(); ++j)
1128  {
1129  if(i!=j)
1130  {
1131  double impact, disti, distj;
1132  GetImpact(TrackEnd1s[i],TrackDir1s[i],TrackEnd1s[j],TrackDir1s[j], impact, disti, distj);
1133 
1134  if((impact < fVertexImpactThreshold)
1135  && (fabs(disti) < fVertexExtrapDistance )
1136  && (fabs(distj) < fVertexExtrapDistance )
1137  && (disti < 0 )
1138  && (distj < 0 ))
1139  {
1140  bool ThisTrackAlreadyMet =false;
1141  TVector3 ExtrapPoint1 = TrackEnd1s[i] + TrackDir1s[j].Unit()*disti;
1142  TVector3 ExtrapPoint2 = TrackEnd1s[j] + TrackDir1s[j].Unit()*distj;
1143 
1144  for(size_t pt=0; pt!=TracksMeetAtPoints.size(); ++pt)
1145  {
1146  if(TracksMeetAtPoints.at(pt)[i]==-1)
1147  {
1148  TracksMeetAtPoints.at(pt)[j]=-1;
1149  MeetPoints.at(pt).push_back(ExtrapPoint1);
1150  MeetPoints.at(pt).push_back(ExtrapPoint2);
1151  ThisTrackAlreadyMet=true;
1152 
1153  }
1154  else if(TracksMeetAtPoints.at(pt)[j]==-1)
1155  {
1156  TracksMeetAtPoints.at(pt)[i]=-1;
1157  MeetPoints.at(pt).push_back(ExtrapPoint1);
1158  MeetPoints.at(pt).push_back(ExtrapPoint2);
1159  ThisTrackAlreadyMet=true;
1160  }
1161  }
1162  if(!ThisTrackAlreadyMet)
1163  {
1164  size_t pt = TracksMeetAtPoints.size();
1165  TracksMeetAtPoints.push_back(std::map<int,int>());
1166  TracksMeetAtPoints.at(pt)[i]=-1;
1167  TracksMeetAtPoints.at(pt)[j]=-1;
1168  MeetPoints.push_back(std::vector<TVector3>());
1169  MeetPoints.at(pt).push_back(ExtrapPoint1);
1170  MeetPoints.at(pt).push_back(ExtrapPoint2);
1171  }
1172 
1173  }
1174 
1175  GetImpact(TrackEnd2s[i],TrackDir2s[i],TrackEnd2s[j],TrackDir2s[j], impact, disti, distj);
1176  if((impact < fVertexImpactThreshold)
1177  && (fabs(disti) < fVertexExtrapDistance )
1178  && (fabs(distj) < fVertexExtrapDistance )
1179  && (disti > 0 )
1180  && (distj > 0 ))
1181  {
1182  bool ThisTrackAlreadyMet =false;
1183  TVector3 ExtrapPoint1 = TrackEnd1s[i] + TrackDir2s[j].Unit()*disti;
1184  TVector3 ExtrapPoint2 = TrackEnd2s[j] + TrackDir2s[j].Unit()*distj;
1185 
1186  for(size_t pt=0; pt!=TracksMeetAtPoints.size(); ++pt)
1187  {
1188  if(TracksMeetAtPoints.at(pt)[i]==1)
1189  {
1190  TracksMeetAtPoints.at(pt)[j]=1;
1191  MeetPoints.at(pt).push_back(ExtrapPoint1);
1192  MeetPoints.at(pt).push_back(ExtrapPoint2);
1193  ThisTrackAlreadyMet=true;
1194 
1195  }
1196  else if(TracksMeetAtPoints.at(pt)[j]==1)
1197  {
1198  TracksMeetAtPoints.at(pt)[i]=1;
1199  MeetPoints.at(pt).push_back(ExtrapPoint1);
1200  MeetPoints.at(pt).push_back(ExtrapPoint2);
1201  ThisTrackAlreadyMet=true;
1202  }
1203  }
1204  if(!ThisTrackAlreadyMet)
1205  {
1206  size_t pt = TracksMeetAtPoints.size();
1207  TracksMeetAtPoints.push_back(std::map<int,int>());
1208  TracksMeetAtPoints.at(pt)[i]=1;
1209  TracksMeetAtPoints.at(pt)[j]=1;
1210  MeetPoints.push_back(std::vector<TVector3>());
1211  MeetPoints.at(pt).push_back(ExtrapPoint1);
1212  MeetPoints.at(pt).push_back(ExtrapPoint2);
1213  }
1214 
1215 
1216  }
1217 
1218  GetImpact(TrackEnd1s[i],TrackDir1s[i],TrackEnd2s[j],TrackDir2s[j], impact, disti, distj);
1219  if((impact < fVertexImpactThreshold)
1220  && (fabs(disti) < fVertexExtrapDistance )
1221  && (fabs(distj) < fVertexExtrapDistance )
1222  && (disti < 0 )
1223  && (distj > 0 ))
1224  {
1225  bool ThisTrackAlreadyMet =false;
1226  TVector3 ExtrapPoint1 = TrackEnd1s[i] + TrackDir2s[j].Unit()*disti;
1227  TVector3 ExtrapPoint2 = TrackEnd2s[j] + TrackDir2s[j].Unit()*distj;
1228 
1229  for(size_t pt=0; pt!=TracksMeetAtPoints.size(); ++pt)
1230  {
1231  if(TracksMeetAtPoints.at(pt)[i]==-1)
1232  {
1233  TracksMeetAtPoints.at(pt)[j]=1;
1234  MeetPoints.at(pt).push_back(ExtrapPoint1);
1235  MeetPoints.at(pt).push_back(ExtrapPoint2);
1236  ThisTrackAlreadyMet=true;
1237 
1238  }
1239  else if(TracksMeetAtPoints.at(pt)[j]==1)
1240  {
1241  TracksMeetAtPoints.at(pt)[i]=-1;
1242  MeetPoints.at(pt).push_back(ExtrapPoint1);
1243  MeetPoints.at(pt).push_back(ExtrapPoint2);
1244  ThisTrackAlreadyMet=true;
1245  }
1246  }
1247  if(!ThisTrackAlreadyMet)
1248  {
1249  size_t pt = TracksMeetAtPoints.size();
1250  TracksMeetAtPoints.push_back(std::map<int,int>());
1251  TracksMeetAtPoints.at(pt)[i]=-1;
1252  TracksMeetAtPoints.at(pt)[j]=1;
1253  MeetPoints.push_back(std::vector<TVector3>());
1254  MeetPoints.at(pt).push_back(ExtrapPoint1);
1255  MeetPoints.at(pt).push_back(ExtrapPoint2);
1256  }
1257 
1258 
1259  }
1260 
1261  }
1262  }
1263  }
1264 
1265 
1266  for(size_t pt=0; pt!=TracksMeetAtPoints.size(); ++pt)
1267  {
1268  std::vector<int> TracksThisVertex;
1269  mf::LogVerbatim("BezierTrackerAlgorithm")<<" Making track adjustments for vertex " <<pt<<std::endl;
1270  TVector3 FinalMeetPt(0,0,0);
1271  for(size_t i=0; i!=MeetPoints.at(pt).size(); ++i)
1272  {
1273  FinalMeetPt += MeetPoints.at(pt).at(i);
1274  }
1275  FinalMeetPt *= (1./float(MeetPoints.at(pt).size()));
1276 
1277  for(size_t i=0; i!=BTracks.size(); ++i)
1278  {
1279  if(TracksMeetAtPoints.at(pt)[i]==1)
1280  {
1281  TracksThisVertex.push_back(i);
1282  TVector3 NewDirection = (FinalMeetPt - TrackEnd2s[i])*0.5;
1283  double NewPos[3];
1284  double NewDir[3];
1285  double Err[3];
1286 
1287  for(size_t n=0; n!=3; ++n)
1288  {
1289  NewPos[n]=FinalMeetPt[n]-NewDirection[n];
1290  NewDir[n]=NewDirection[n];
1291  }
1292 
1293  recob::Seed NewSeed(NewPos, NewDir, Err, Err);
1294  std::vector<recob::Seed> SeedCol = BTracks.at(i).GetSeedVector();
1295  SeedCol.at(SeedCol.size()-1) = NewSeed;
1296  BTracks.at(i) = trkf::BezierTrack(SeedCol);
1297  }
1298  else if(TracksMeetAtPoints.at(pt)[i]==-1)
1299  {
1300  TracksThisVertex.push_back(i);
1301  TVector3 NewDirection = -(FinalMeetPt - TrackEnd1s[i])*0.5;
1302  double NewPos[3];
1303  double NewDir[3];
1304  double Err[3];
1305 
1306  for(size_t n=0; n!=3; ++n)
1307  {
1308  NewPos[n]=FinalMeetPt[n] + NewDirection[n];
1309  NewDir[n]=NewDirection[n];
1310  }
1311 
1312  recob::Seed NewSeed(NewPos, NewDir, Err, Err);
1313  std::vector<recob::Seed> SeedCol = BTracks.at(i).GetSeedVector();
1314  SeedCol.at(0) = NewSeed;
1315  BTracks.at(i) = trkf::BezierTrack(SeedCol);
1316  }
1317 
1318 
1319  }
1320  Mapping.push_back(TracksThisVertex);
1321  double VtxPoint[3];
1322  for(size_t n=0; n!=3; ++n)
1323  VtxPoint[n]=FinalMeetPt[n];
1324  Vertices.push_back(recob::Vertex(VtxPoint, pt));
1325  }
1326 
1327  }
1328 
1329 
1330  //------------------------------------------------
1331 
1332  void BezierTrackerAlgorithm::GetImpact(TVector3 t1pt, TVector3 t1dir, TVector3 t2pt, TVector3 t2dir, double& ImpactParam, double& Dist1, double& Dist2)
1333  {
1334  double lambda1 = ((t2pt-t1pt).Cross(t2dir)).Dot(t1dir.Cross(t2dir)) / (t1dir.Cross(t2dir)).Mag2();
1335  double lambda2 = ((t1pt-t2pt).Cross(t1dir)).Dot(t2dir.Cross(t1dir)) / (t2dir.Cross(t1dir)).Mag2();
1336  TVector3 c = t1pt + t1dir*lambda1 - t2pt - t2dir*lambda2;
1337 
1338  ImpactParam = c.Mag();
1339  Dist1 = lambda1 * t1dir.Mag();
1340  Dist2 = lambda2 * t2dir.Mag();
1341  }
1342 
1343 
1344  //-------------------------------------------
1345  bool BTrack_SeedCountComparator(std::vector<recob::Seed> const& s1, std::vector<recob::Seed> const& s2)
1346  {
1347  return s1.size()>s2.size();
1348  }
1349 }
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:219
trkf::BezierTrackerAlgorithm::MakeOverlapJoins
void MakeOverlapJoins(std::vector< trkf::BezierTrack > &BTracks, std::vector< art::PtrVector< recob::Hit > > &HitVecs, std::vector< std::vector< std::pair< int, int > > > &ClustersUsed)
Definition: BezierTrackerAlgorithm.cxx:340
t1
TTree * t1
Definition: plottest35.C:26
trkf::BezierTrackerAlgorithm::GetSeedDirProjected
void GetSeedDirProjected(recob::Seed const &TheSeed, std::vector< double > &WireCoord, std::vector< double > &TimeCoord)
Definition: BezierTrackerAlgorithm.cxx:1022
trkf::BTrack_SeedCountComparator
bool BTrack_SeedCountComparator(std::vector< recob::Seed > const &s1, std::vector< recob::Seed > const &s2)
Definition: BezierTrackerAlgorithm.cxx:1345
trkf::BezierTrack::GetSeedVector
std::vector< recob::Seed > GetSeedVector()
Definition: BezierTrack.cxx:108
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:19
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
geo::kV
Planes which measure V.
Definition: geo_types.h:77
recob::Hit::WireID
geo::WireID WireID() const
Initial tdc tick for hit.
Definition: Hit.h:234
Hit.h
Declaration of signal hit object.
y
Float_t y
Definition: compare.C:6
z
Double_t z
Definition: plot.C:279
trkf::BezierTrackerAlgorithm::SortTracksByLength
void SortTracksByLength(std::vector< trkf::BezierTrack > &BTracks, std::vector< art::PtrVector< recob::Hit > > &HitVecs, std::vector< std::vector< std::pair< int, int > > > &ClustersUsed)
Definition: BezierTrackerAlgorithm.cxx:401
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:313
recob::Vertex
Definition of vertex object for LArSoft.
Definition: Vertex.h:35
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:783
detinfo::DetectorProperties::ConvertXToTicks
virtual double ConvertXToTicks(double X, int p, int t, int c) const =0
geo::GeometryCore::Nchannels
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
Definition: GeometryCore.cxx:191
trkf::BezierTrackerAlgorithm::MakeTracks
std::vector< trkf::BezierTrack > MakeTracks(std::vector< std::vector< art::PtrVector< recob::Hit > > > &SortedHits, std::vector< art::PtrVector< recob::Hit > > &HitAssocs, std::vector< std::vector< std::pair< int, int > > > &ClustersUsed)
Definition: BezierTrackerAlgorithm.cxx:71
trkf::SeedFinderAlgorithm::GetSpacePointAlg
SpacePointAlg * GetSpacePointAlg() const
Definition: SeedFinderAlgorithm.h:65
Event.h
trkf::BezierTrackerAlgorithm::fPitchDir
std::vector< TVector3 > fPitchDir
Definition: BezierTrackerAlgorithm.h:103
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
geo::kU
Planes which measure U.
Definition: geo_types.h:76
trkf::BezierTrackerAlgorithm::fPlaneTimeOffsets
std::vector< double > fPlaneTimeOffsets
Definition: BezierTrackerAlgorithm.h:102
trkf::BezierTrackerAlgorithm::fTheSeedFinder
SeedFinderAlgorithm * fTheSeedFinder
Definition: BezierTrackerAlgorithm.h:114
trkf::BezierTrackerAlgorithm::JoinTracks
trkf::BezierTrack JoinTracks(trkf::BezierTrack &BT1, trkf::BezierTrack &BT2)
Definition: BezierTrackerAlgorithm.cxx:519
max
Int_t max
Definition: plot.C:27
trkf::BezierTrackerAlgorithm::MakeVertexJoins
void MakeVertexJoins(std::vector< trkf::BezierTrack > &BTracks, std::vector< recob::Vertex > &Vertices, std::vector< std::vector< int > > &Mapping)
Definition: BezierTrackerAlgorithm.cxx:1103
trkf::BezierTrackerAlgorithm::reconfigure
void reconfigure(fhicl::ParameterSet const &pset)
Definition: BezierTrackerAlgorithm.cxx:48
Track.h
Provides recob::Track data product.
art::PtrVector::push_back
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
BezierCurveHelper.h
Service for finding 3D track seeds.
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:231
trkf::BezierTrackerAlgorithm::OrganizeSeedsIntoTracks
std::vector< std::vector< recob::Seed > > OrganizeSeedsIntoTracks(std::vector< recob::Seed > &AllSeeds, std::vector< art::PtrVector< recob::Hit > * > &AllHits, std::vector< art::PtrVector< recob::Hit > > &WhichHitsPerSeed, std::vector< std::vector< std::vector< int > > * > &OrgHits, std::vector< std::vector< std::vector< int > > > &WhichHitsPerTrack)
Definition: BezierTrackerAlgorithm.cxx:545
pt
TMarker * pt
Definition: egs.C:25
d
Float_t d
Definition: plot.C:237
trkf::SpacePointAlg::maxDT
double maxDT() const
Definition: SpacePointAlg.h:90
recob::Seed::GetDistance
double GetDistance(Seed const &AnotherSeed) const
Definition: Seed.cxx:243
art::PtrVector::at
reference at(size_type n)
Definition: PtrVector.h:365
geo::GeometryCore::IntersectionPoint
bool IntersectionPoint(geo::WireID const &wid1, geo::WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
Definition: GeometryCore.cxx:1749
trkf::BezierTrackerAlgorithm::AddPtrVectors
void AddPtrVectors(art::PtrVector< recob::Hit > &Receiever, art::PtrVector< recob::Hit > const &ToAdd)
Definition: BezierTrackerAlgorithm.cxx:533
art::PtrVector::size
size_type size() const
Definition: PtrVector.h:308
t2
TTree * t2
Definition: plottest35.C:36
trkf::BezierTrackerAlgorithm::MakeDirectJoins
void MakeDirectJoins(std::vector< trkf::BezierTrack > &BTracks, std::vector< art::PtrVector< recob::Hit > > &HitVecs, std::vector< std::vector< std::pair< int, int > > > &ClustersUsed)
Definition: BezierTrackerAlgorithm.cxx:448
detinfo::DetectorProperties::ConvertTicksToX
virtual double ConvertTicksToX(double ticks, int p, int t, int c) const =0
trkf::BezierTrackerAlgorithm::fWireZeroOffset
std::vector< double > fWireZeroOffset
Definition: BezierTrackerAlgorithm.h:105
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:219
AssociationUtil.h
Utility object to perform functions of association.
dir
TDirectory * dir
Definition: macro.C:5
min
Int_t min
Definition: plot.C:26
trkf::SeedFinderAlgorithm::GetSeedsFromUnSortedHits
std::vector< recob::Seed > GetSeedsFromUnSortedHits(art::PtrVector< recob::Hit > const &, std::vector< art::PtrVector< recob::Hit > > &, unsigned int StopAfter=0)
Definition: SeedFinderAlgorithm.cxx:1134
trkf::BezierTrackerAlgorithm::FilterOverlapTracks
void FilterOverlapTracks(std::vector< trkf::BezierTrack > &BTracks, std::vector< art::PtrVector< recob::Hit > > &HitVecs, std::vector< std::vector< std::pair< int, int > > > &ClustersUsed)
Definition: BezierTrackerAlgorithm.cxx:276
trkf::BezierTrackerAlgorithm::GetSeedFinderAlgorithm
trkf::SeedFinderAlgorithm * GetSeedFinderAlgorithm()
Definition: BezierTrackerAlgorithm.h:53
trkf::BezierTrackerAlgorithm::fWireDir
std::vector< TVector3 > fWireDir
Definition: BezierTrackerAlgorithm.h:104
geo::GeometryCore::WireEndPoints
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
Definition: GeometryCore.cxx:1307
n
Char_t n[5]
Definition: comparison_ascii.C:55
geo::kW
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:78
trkf::BezierCurveHelper::GetBezierPoints
std::vector< TVector3 > GetBezierPoints(recob::Seed const &s1, recob::Seed const &s2, int N=100)
Definition: BezierCurveHelper.cxx:225
geo
Namespace collecting geometry-related classes utilities.
Definition: TrackContainmentAlg.hh:20
recob::Seed::GetDirection
void GetDirection(double *Dir, double *Err) const
Definition: Seed.cxx:102
trkf::BezierTrackerAlgorithm::BezierTrackerAlgorithm
BezierTrackerAlgorithm(fhicl::ParameterSet const &pset)
Definition: BezierTrackerAlgorithm.cxx:33
trkf::BezierTrackerAlgorithm::EvaluateOccupancy
bool EvaluateOccupancy(recob::Seed &Seed1, recob::Seed &Seed2, double dThresh, std::vector< art::PtrVector< recob::Hit > * > &AllHits, std::vector< std::vector< std::vector< int > > * > &OrgHits, std::vector< uint32_t > &LowChan, std::vector< uint32_t > &HighChan, std::vector< std::vector< int > > &HitStatus, std::vector< std::vector< int > > &TheseHits)
Definition: BezierTrackerAlgorithm.cxx:969
trkf::BezierTrackerAlgorithm::GetImpact
void GetImpact(TVector3 t1pt, TVector3 t1dir, TVector3 t2pt, TVector3 t2dir, double &ImpactParam, double &Dist1, double &Dist2)
Definition: BezierTrackerAlgorithm.cxx:1332
detinfo::DetectorProperties::GetXTicksOffset
virtual double GetXTicksOffset(int p, int t, int c) const =0