LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
PmaTrack3D.cxx
Go to the documentation of this file.
1 
12 #include "larreco/RecoAlg/PMAlg/PmaTrack3D.h"
13 #include "larreco/RecoAlg/PMAlg/Utilities.h"
14 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
15 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
16 #include "larcore/Geometry/Geometry.h"
17 #include "messagefacility/MessageLogger/MessageLogger.h"
18 
19 #include "TMath.h"
20 
21 pma::Track3D::Track3D(void) :
22  fMaxHitsPerSeg(70),
23  fPenaltyFactor(1.0F),
24  fMaxSegStopFactor(8.0F),
25 
26  fSegStopValue(2), fMinSegStop(2), fMaxSegStop(2),
27  fSegStopFactor(0.2F),
28  fPenaltyValue(0.1F),
29  fEndSegWeight(0.05F),
30  fHitsRadius(1.0F),
31 
32  fT0(0.0),
33  fT0Flag(false),
34 
35  fTag(pma::Track3D::kNotTagged)
36 {
37 }
38 
39 pma::Track3D::Track3D(const Track3D& src) :
40  fMaxHitsPerSeg(src.fMaxHitsPerSeg),
41  fPenaltyFactor(src.fPenaltyFactor),
42  fMaxSegStopFactor(src.fMaxSegStopFactor),
43 
44  fSegStopValue(src.fSegStopValue),
45  fMinSegStop(src.fMinSegStop),
46  fMaxSegStop(src.fMaxSegStop),
47  fSegStopFactor(src.fSegStopFactor),
48  fPenaltyValue(src.fPenaltyValue),
49  fEndSegWeight(src.fEndSegWeight),
50  fHitsRadius(src.fHitsRadius),
51 
52  fT0(src.fT0),
53  fT0Flag(src.fT0Flag),
54 
55  fTag(src.fTag)
56 {
57  fHits.reserve(src.fHits.size());
58  for (auto const& hit : src.fHits)
59  {
60  pma::Hit3D* h3d = new pma::Hit3D(*hit);
61  h3d->fParent = this;
62  fHits.push_back(h3d);
63  }
64 
65  fNodes.reserve(src.fNodes.size());
66  for (auto const& node : src.fNodes) fNodes.push_back(new pma::Node3D(node->Point3D(), node->TPC(), node->Cryo(), node->IsVertex(), node->GetDriftShift()));
67 
68  for (auto const& point : src.fAssignedPoints) fAssignedPoints.push_back(new TVector3(*point));
69 
70  RebuildSegments();
71  MakeProjection();
72 }
73 
74 pma::Track3D::~Track3D(void)
75 {
76  for (size_t i = 0; i < fHits.size(); i++) delete fHits[i];
77  for (size_t i = 0; i < fAssignedPoints.size(); i++) delete fAssignedPoints[i];
78 
79  for (size_t i = 0; i < fSegments.size(); i++) delete fSegments[i];
80  for (size_t i = 0; i < fNodes.size(); i++)
81  if (!fNodes[i]->NextCount() && !fNodes[i]->Prev()) delete fNodes[i];
82 }
83 
84 bool pma::Track3D::Initialize(float initEndSegW)
85 {
86  if (!HasTwoViews(2))
87  {
88  mf::LogError("pma::Track3D") << "Need min. 2 hits per view, at least two views.";
89  return false;
90  }
91 
92  auto cryos = Cryos();
93  if (cryos.size() > 1)
94  {
95  mf::LogError("pma::Track3D") << "Only one cryostat for now, please.";
96  return false;
97  }
98  int cryo = cryos.front();
99 
100  auto tpcs = TPCs();
101  if (tpcs.size() > 1)
102  {
103  mf::LogError("pma::Track3D") << "Only one TPC, please.";
104  return false;
105  }
106  // single tpc, many tpc's are ok, but need to be handled from ProjectionMatchingAlg::buildMultiTPCTrack()
107  int tpc = tpcs.front();
108 
109  if (InitFromRefPoints(tpc, cryo)) mf::LogVerbatim("pma::Track3D") << "Track initialized with 3D reference points.";
110  else
111  {
112  if (InitFromHits(tpc, cryo, initEndSegW)) mf::LogVerbatim("pma::Track3D") << "Track initialized with hit positions.";
113  else { InitFromMiddle(tpc, cryo); mf::LogVerbatim("pma::Track3D") << "Track initialized in the module center."; }
114  }
115 
116  UpdateHitsRadius();
117  return true;
118 }
119 
120 void pma::Track3D::ClearNodes(void)
121 {
122  for (size_t i = 0; i < fNodes.size(); i++) delete fNodes[i];
123  fNodes.clear();
124 }
125 
126 bool pma::Track3D::InitFromHits(int tpc, int cryo, float initEndSegW)
127 {
128  auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
129  art::ServiceHandle<geo::Geometry> geom;
130 
131  float wtmp = fEndSegWeight;
132  fEndSegWeight = initEndSegW;
133 
134  // endpoints for the first combination:
135  TVector3 v3d_1(0., 0., 0.), v3d_2(0., 0., 0.);
136  double x, y, z;
137 
138  pma::Hit3D* hit0_a = front();
139  pma::Hit3D* hit0_b = 0;
140 
141  pma::Hit3D* hit1_a = front();
142  pma::Hit3D* hit1_b = 0;
143 
144  pma::Hit3D* hit = 0;
145  float minX = detprop->ConvertTicksToX(hit0_a->PeakTime(), hit0_a->View2D(), tpc, cryo);
146  float maxX = detprop->ConvertTicksToX(hit1_a->PeakTime(), hit1_a->View2D(), tpc, cryo);
147  for (size_t i = 1; i < size(); i++)
148  {
149  hit = fHits[i];
150  x = detprop->ConvertTicksToX(hit->PeakTime(), hit->View2D(), tpc, cryo);
151  if (x < minX) { minX = x; hit0_a = hit; }
152  if (x > maxX) { maxX = x; hit1_a = hit; }
153  }
154 
155  float diff, minDiff0 = 5000, minDiff1 = 5000;
156  for (size_t i = 0; i < size(); i++)
157  {
158  hit = fHits[i];
159  x = detprop->ConvertTicksToX(hit->PeakTime(), hit->View2D(), tpc, cryo);
160  diff = fabs(x - minX);
161  if ((diff < minDiff0) && (hit->View2D() != hit0_a->View2D()))
162  {
163  minDiff0 = diff; hit0_b = hit;
164  }
165  diff = fabs(x - maxX);
166  if ((diff < minDiff1) && (hit->View2D() != hit1_a->View2D()))
167  {
168  minDiff1 = diff; hit1_b = hit;
169  }
170  }
171 
172  if (hit0_a && hit0_b && hit1_a && hit1_b)
173  {
174  x = 0.5 * (
175  detprop->ConvertTicksToX(hit0_a->PeakTime(), hit0_a->View2D(), tpc, cryo) +
176  detprop->ConvertTicksToX(hit0_b->PeakTime(), hit0_b->View2D(), tpc, cryo));
177  geom->IntersectionPoint(hit0_a->Wire(), hit0_b->Wire(),
178  hit0_a->View2D(), hit0_b->View2D(), cryo, tpc, y, z);
179  v3d_1.SetXYZ(x, y, z);
180 
181  x = 0.5 * (
182  detprop->ConvertTicksToX(hit1_a->PeakTime(), hit1_a->View2D(), tpc, cryo) +
183  detprop->ConvertTicksToX(hit1_b->PeakTime(), hit1_b->View2D(), tpc, cryo));
184  geom->IntersectionPoint(hit1_a->Wire(), hit1_b->Wire(),
185  hit1_a->View2D(), hit1_b->View2D(), cryo, tpc, y, z);
186  v3d_2.SetXYZ(x, y, z);
187 
188  ClearNodes();
189  AddNode(v3d_1, tpc, cryo);
190  AddNode(v3d_2, tpc, cryo);
191 
192  MakeProjection();
193  UpdateHitsRadius();
194  Optimize(0, 0.01F, false, true, 100);
195  SelectAllHits();
196  }
197  else
198  {
199  mf::LogVerbatim("pma::Track3D") << "Good hits not found.";
200  fEndSegWeight = wtmp;
201  return false;
202  }
203 
204  if (pma::Dist2(fNodes.front()->Point3D(), fNodes.back()->Point3D()) < (0.3*0.3))
205  {
206  mf::LogVerbatim("pma::Track3D") << "Short initial segment.";
207  fEndSegWeight = wtmp;
208  return false;
209  }
210 
211  fEndSegWeight = wtmp;
212  return true;
213 }
214 
215 bool pma::Track3D::InitFromRefPoints(int tpc, int cryo)
216 {
217  if (fAssignedPoints.size() < 2) return false;
218 
219  ClearNodes();
220 
221  TVector3 mean(0., 0., 0.), stdev(0., 0., 0.), p(0., 0., 0.);
222  for (size_t i = 0; i < fAssignedPoints.size(); i++)
223  {
224  p = *(fAssignedPoints[i]);
225  mean += p;
226  p.SetXYZ( p.X()*p.X(), p.Y()*p.Y(), p.Z()*p.Z() );
227  stdev += p;
228  }
229  stdev *= 1.0 / fAssignedPoints.size();
230  mean *= 1.0 / fAssignedPoints.size();
231  p = mean;
232  p.SetXYZ( p.X()*p.X(), p.Y()*p.Y(), p.Z()*p.Z() );
233  stdev -= p;
234 
235  double sx = stdev.X(), sy = stdev.Y(), sz = stdev.Z();
236  if (sx >= 0.0) sx = sqrt(sx);
237  else sx = 0.0;
238  if (sy >= 0.0) sy = sqrt(sy);
239  else sy = 0.0;
240  if (sz >= 0.0) sz = sqrt(sz);
241  else sz = 0.0;
242  stdev.SetXYZ(sx, sy, sz);
243 
244  double scale = 2.0 * stdev.Mag();
245  double iscale = 1.0 / scale;
246 
247  size_t max_index = 0;
248  double norm2, max_norm2 = 0.0;
249  std::vector< TVector3 > data;
250  for (size_t i = 0; i < fAssignedPoints.size(); i++)
251  {
252  p = *(fAssignedPoints[i]);
253  p -= mean;
254  p *= iscale;
255  norm2 = p.Mag2();
256  if (norm2 > max_norm2)
257  {
258  max_norm2 = norm2;
259  max_index = i;
260  }
261  data.push_back(p);
262  }
263 
264  double y = 0.0, kappa = 1.0, prev_kappa, kchg = 1.0;
265  TVector3 w(data[max_index]);
266 
267  while (kchg > 0.0001)
268  for (size_t i = 0; i < data.size(); i++)
269  {
270  y = (data[i] * w);
271  w += (y/kappa) * (data[i] - y*w);
272 
273  prev_kappa = kappa;
274  kappa += y*y;
275  kchg = fabs((kappa - prev_kappa) / prev_kappa);
276  }
277  w *= 1.0 / w.Mag();
278 
279  TVector3 v1(w), v2(w);
280  v1 *= scale; v1 += mean;
281  v2 *= -scale; v2 += mean;
282  std::sort(fAssignedPoints.begin(), fAssignedPoints.end(), pma::bSegmentProjLess(v1, v2));
283  for (size_t i = 0; i < fAssignedPoints.size(); i++)
284  {
285  AddNode(*(fAssignedPoints[i]), tpc, cryo);
286  }
287 
288  RebuildSegments();
289  MakeProjection();
290 
291  if (size()) UpdateHitsRadius();
292 
293  Optimize(0, 0.01F, false, true, 100);
294  SelectAllHits();
295 
296  return true;
297 }
298 
299 void pma::Track3D::InitFromMiddle(int tpc, int cryo)
300 {
301  art::ServiceHandle<geo::Geometry> geom;
302 
303  const auto& tpcGeo = geom->TPC(tpc, cryo);
304 
305  double minX = tpcGeo.MinX(), maxX = tpcGeo.MaxX();
306  double minY = tpcGeo.MinY(), maxY = tpcGeo.MaxY();
307  double minZ = tpcGeo.MinZ(), maxZ = tpcGeo.MaxZ();
308 
309  TVector3 v3d_1(0.5 * (minX + maxX), 0.5 * (minY + maxY), 0.5 * (minZ + maxZ));
310  TVector3 v3d_2(v3d_1);
311 
312  TVector3 shift(5.0, 5.0, 5.0);
313  v3d_1 += shift;
314  v3d_2 -= shift;
315 
316  ClearNodes();
317  AddNode(v3d_1, tpc, cryo);
318  AddNode(v3d_2, tpc, cryo);
319 
320  MakeProjection();
321  UpdateHitsRadius();
322 
323  Optimize(0, 0.01F);
324 }
325 
326 int pma::Track3D::index_of(const pma::Node3D* n) const
327 {
328  for (size_t i = 0; i < fNodes.size(); ++i)
329  if (fNodes[i] == n) return (int)i;
330  return -1;
331 }
332 
333 int pma::Track3D::index_of(const pma::Hit3D* hit) const
334 {
335  for (size_t i = 0; i < size(); i++)
336  if (fHits[i] == hit) return (int)i;
337  return -1;
338 }
339 
340 pma::Hit3D* pma::Track3D::release_at(size_t index)
341 {
342  pma::Hit3D* h3d = fHits[index];
343  fHits.erase(fHits.begin() + index);
344  return h3d;
345 }
346 
347 bool pma::Track3D::push_back(const art::Ptr< recob::Hit > & hit)
348 {
349  for (auto const& trk_hit : fHits)
350  {
351  if (trk_hit->fHit == hit) return false;
352  }
353  pma::Hit3D* h3d = new pma::Hit3D(hit);
354  h3d->fParent = this;
355  fHits.push_back(h3d);
356  return true;
357 }
358 
359 bool pma::Track3D::erase(const art::Ptr< recob::Hit > & hit)
360 {
361  for (size_t i = 0; i < size(); i++)
362  {
363  if (hit == fHits[i]->fHit)
364  {
365  pma::Hit3D* h3d = fHits[i];
366  fHits.erase(fHits.begin() + i);
367  delete h3d; return true;
368  }
369  }
370  return false;
371 }
372 
373 pma::Vector3D pma::Track3D::GetDirection3D(size_t index) const
374 {
375  pma::Hit3D* h = fHits[index];
376 
377  for (auto s : fSegments) { if (s->HasHit(h)) return s->GetDirection3D(); }
378  for (auto n : fNodes) { if (n->HasHit(h)) return n->GetDirection3D(); }
379 
380  auto pe = GetNearestElement(h->Point2D(), h->View2D(), h->TPC());
381  if (pe)
382  {
383  mf::LogWarning("pma::Track3D") << "GetDirection3D(): had to update hit assignment to segment/node.";
384  pe->AddHit(h); return pe-> GetDirection3D();
385  }
386  else
387  {
388  throw cet::exception("pma::Track3D") << "GetDirection3D(): direction of a not assigned hit "
389  << index << " (size: " << fHits.size() << ")" << std::endl;
390  }
391 }
392 
393 void pma::Track3D::AddHits(const std::vector< art::Ptr<recob::Hit> >& hits)
394 {
395  fHits.reserve(fHits.size() + hits.size());
396  for (auto const& hit : hits) push_back(hit);
397 }
398 
399 void pma::Track3D::RemoveHits(const std::vector< art::Ptr<recob::Hit> >& hits)
400 {
401  unsigned int n = 0;
402  for (auto const& hit : hits) { if (erase(hit)) n++; }
403  if (n) MakeProjection();
404 }
405 
406 unsigned int pma::Track3D::NHits(unsigned int view) const
407 {
408  unsigned int n = 0;
409  for (size_t i = 0; i < fHits.size(); i++)
410  {
411  if (fHits[i]->View2D() == view) n++;
412  }
413  return n;
414 }
415 
416 unsigned int pma::Track3D::NEnabledHits(unsigned int view) const
417 {
418  unsigned int n = 0;
419  for (size_t i = 0; i < size(); i++)
420  {
421  pma::Hit3D* hit = fHits[i];
422  if (hit->IsEnabled() &&
423  ((view == geo::kUnknown) || (view == hit->View2D()))) n++;
424  }
425  return n;
426 }
427 
428 bool pma::Track3D::HasTwoViews(size_t nmin) const
429 {
430  unsigned int nviews = 0;
431  if (NHits(geo::kU) >= nmin) nviews++;
432  if (NHits(geo::kV) >= nmin) nviews++;
433  if (NHits(geo::kZ) >= nmin) nviews++;
434  return (nviews > 1);
435 }
436 
437 std::vector< unsigned int > pma::Track3D::TPCs(void) const
438 {
439  std::vector< unsigned int > tpc_idxs;
440  for (size_t i = 0; i < size(); i++)
441  {
442  unsigned int tpc = fHits[i]->TPC();
443 
444  bool found = false;
445  for (size_t j = 0; j < tpc_idxs.size(); j++)
446  if (tpc_idxs[j] == tpc) { found = true; break; }
447 
448  if (!found) tpc_idxs.push_back(tpc);
449  }
450  return tpc_idxs;
451 }
452 
453 std::vector< unsigned int > pma::Track3D::Cryos(void) const
454 {
455  std::vector< unsigned int > cryo_idxs;
456  for (size_t i = 0; i < size(); i++)
457  {
458  unsigned int cryo = fHits[i]->Cryo();
459 
460  bool found = false;
461  for (size_t j = 0; j < cryo_idxs.size(); j++)
462  if (cryo_idxs[j] == cryo) { found = true; break; }
463 
464  if (!found) cryo_idxs.push_back(cryo);
465  }
466  return cryo_idxs;
467 }
468 
469 std::pair< TVector2, TVector2 > pma::Track3D::WireDriftRange(unsigned int view, unsigned int tpc, unsigned int cryo) const
470 {
471  std::pair< TVector2, TVector2 > range(TVector2(0., 0.), TVector2(0., 0.));
472 
473  size_t n0 = 0;
474  while ((n0 < fNodes.size()) && (fNodes[n0]->TPC() != (int)tpc)) n0++;
475 
476  if (n0 < fNodes.size())
477  {
478  size_t n1 = n0;
479  while ((n1 < fNodes.size()) && (fNodes[n1]->TPC() == (int)tpc)) n1++;
480 
481  if (n0 > 0) n0--;
482  if (n1 == fNodes.size()) n1--;
483 
484  TVector2 p0 = fNodes[n0]->Projection2D(view);
485  p0 = pma::CmToWireDrift(p0.X(), p0.Y(), view, tpc, cryo);
486 
487  TVector2 p1 = fNodes[n1]->Projection2D(view);
488  p1 = pma::CmToWireDrift(p1.X(), p1.Y(), view, tpc, cryo);
489 
490  if (p0.X() > p1.X()) { double tmp = p0.X(); p0.Set(p1.X(), p0.Y()); p1.Set(tmp, p1.Y()); }
491  if (p0.Y() > p1.Y()) { double tmp = p0.Y(); p0.Set(p0.X(), p1.Y()); p1.Set(p1.X(), tmp); }
492 
493  range.first = p0;
494  range.second = p1;
495  }
496  return range;
497 }
498 
499 bool pma::Track3D::Flip(std::vector< pma::Track3D* >& allTracks)
500 {
501  if (fNodes.size() < 2) { return true; }
502 
503  std::vector< pma::Track3D* > toSort;
504 
505  pma::Node3D* n = fNodes.front();
506  if (n->Prev())
507  {
508  pma::Segment3D* s = static_cast< pma::Segment3D* >(n->Prev());
509  pma::Track3D* t = s->Parent();
510 
511  if (t->NextSegment(n)) // starts from middle of another track: need to break that one first
512  {
513  int idx = t->index_of(n);
514  if (idx >= 0)
515  {
516  pma::Track3D* u = t->Split(idx, false); // u is in front of t
517  if (u)
518  {
519  allTracks.push_back(u);
520  if (u->Flip(allTracks))
521  {
522  InternalFlip(toSort);
523  toSort.push_back(this);
524  }
525  else
526  {
527  mf::LogWarning("pma::Track3D") << "Flip(): Could not flip after split (but new track is preserved!!).";
528  return false;
529  }
530  }
531  else { return false; } // could not flip/break associated tracks, so give up on this one
532  }
533  else { throw cet::exception("pma::Track3D") << "Node not found." << std::endl; }
534  }
535  else // flip root
536  {
537  if (t->Flip(allTracks)) // all ok, so can flip this track
538  {
539  InternalFlip(toSort);
540  toSort.push_back(this);
541  }
542  else { return false; } // could not flip/break associated tracks, so give up on this one
543  }
544  }
545  else // simply flip
546  {
547  InternalFlip(toSort);
548  toSort.push_back(this);
549  }
550 
551  for (size_t t = 0; t < toSort.size(); t++)
552  {
553  bool sorted = false;
554  for (size_t u = 0; u < t; u++)
555  if (toSort[u] == toSort[t]) { sorted = true; break; }
556  if (!sorted)
557  {
558  toSort[t]->MakeProjection();
559  toSort[t]->SortHits();
560  }
561  }
562  return true;
563 }
564 
565 void pma::Track3D::InternalFlip(std::vector< pma::Track3D* >& toSort)
566 {
567  for (size_t i = 0; i < fNodes.size() - 1; i++)
568  if (fNodes[i]->NextCount() > 1)
569  {
570  for (size_t j = 0; j < fNodes[i]->NextCount(); j++)
571  {
572  pma::Segment3D* s = static_cast< pma::Segment3D* >(fNodes[i]->Next(j));
573  if (s->Parent() != this) toSort.push_back(s->Parent());
574  }
575  }
576 
577  if (fNodes.back()->NextCount())
578  for (size_t j = 0; j < fNodes.back()->NextCount(); j++)
579  {
580  pma::Segment3D* s = static_cast< pma::Segment3D* >(fNodes.back()->Next(j));
581  toSort.push_back(s->Parent());
582  }
583 
584  if (fNodes.front()->Prev())
585  {
586  pma::Segment3D* s = static_cast< pma::Segment3D* >(fNodes.front()->Prev());
587  toSort.push_back(s->Parent());
588  s->Parent()->InternalFlip(toSort);
589  }
590 
591  std::reverse(fNodes.begin(), fNodes.end());
592  toSort.push_back(this);
593  RebuildSegments();
594 }
595 
596 void pma::Track3D::Flip(void)
597 {
598  std::vector< pma::Track3D* > toSort;
599  InternalFlip(toSort);
600  toSort.push_back(this);
601 
602  for (size_t t = 0; t < toSort.size(); t++)
603  {
604  bool sorted = false;
605  for (size_t u = 0; u < t; u++)
606  if (toSort[u] == toSort[t]) { sorted = true; break; }
607  if (!sorted)
608  {
609  toSort[t]->MakeProjection();
610  toSort[t]->SortHits();
611  }
612  }
613 }
614 
615 bool pma::Track3D::CanFlip(void) const
616 {
617  if (!fNodes.size()) { return false; }
618 
619  pma::Node3D* n = fNodes.front();
620  if (n->Prev())
621  {
622  pma::Segment3D* s = static_cast< pma::Segment3D* >(n->Prev());
623  pma::Track3D* t = s->Parent();
624  if (t->NextSegment(n)) { return false; } // cannot flip if starts from middle of another track
625  else { return t->CanFlip(); } // check root
626  }
627  else { return true; }
628 }
629 
630 void pma::Track3D::AutoFlip(pma::Track3D::EDirection dir, double thr, unsigned int n)
631 {
632  unsigned int nViews = 3;
633  pma::dedx_map dedxInViews[3];
634  for (unsigned int i = 0; i < nViews; i++)
635  {
636  GetRawdEdxSequence(dedxInViews[i], i, 1);
637  }
638  unsigned int bestView = 2;
639  if (dedxInViews[0].size() > 2 * dedxInViews[2].size()) bestView = 0;
640  if (dedxInViews[1].size() > 2 * dedxInViews[2].size()) bestView = 1;
641 
642  std::vector< std::vector<double> > dedx;
643  for (size_t i = 0; i < size(); i++)
644  {
645  auto it = dedxInViews[bestView].find(i);
646  if (it != dedxInViews[bestView].end())
647  {
648  dedx.push_back(it->second);
649  }
650  }
651  if (!dedx.empty()) dedx.pop_back();
652 
653  float dEdxStart = 0.0F, dEdxStop = 0.0F;
654  float dEStart = 0.0F, dxStart = 0.0F;
655  float dEStop = 0.0F, dxStop = 0.0F;
656  if (dedx.size() > 4)
657  {
658  if (!n) // use default options
659  {
660  if (dedx.size() > 30) n = 12;
661  else if (dedx.size() > 20) n = 8;
662  else if (dedx.size() > 10) n = 4;
663  else n = 3;
664  }
665 
666  size_t k = (dedx.size() - 2) >> 1;
667  if (n > k) n = k;
668 
669  for (size_t i = 1, j = 0; j < n; i++, j++)
670  {
671  dEStart += dedx[i][5]; dxStart += dedx[i][6];
672  }
673  if (dxStart > 0.0F) dEdxStart = dEStart / dxStart;
674 
675  for (size_t i = dedx.size() - 2, j = 0; j < n; i--, j++)
676  {
677  dEStop += dedx[i][5]; dxStop += dedx[i][6];
678  }
679  if (dxStop > 0.0F) dEdxStop = dEStop / dxStop;
680  }
681  else if (dedx.size() == 4)
682  {
683  dEStart = dedx[0][5] + dedx[1][5]; dxStart = dedx[0][6] + dedx[1][6];
684  dEStop = dedx[2][5] + dedx[3][5]; dxStop = dedx[2][6] + dedx[3][6];
685  if (dxStart > 0.0F) dEdxStart = dEStart / dxStart;
686  if (dxStop > 0.0F) dEdxStop = dEStop / dxStop;
687 
688  }
689  else if (dedx.size() > 1)
690  {
691  if (dedx.front()[2] > 0.0F) dEdxStart = dedx.front()[5] / dedx.front()[6];
692  if (dedx.back()[2] > 0.0F) dEdxStop = dedx.back()[5] / dedx.back()[6];
693  }
694  else return;
695 
696  if ((dir == pma::Track3D::kForward) && ((1.0 + thr) * dEdxStop < dEdxStart))
697  {
698  mf::LogVerbatim("pma::Track3D") << "Auto-flip fired (1), thr: " << (1.0+thr) << ", value: " << dEdxStart/dEdxStop;
699  Flip(); // particle stops at the end of the track
700  }
701  if ((dir == pma::Track3D::kBackward) && (dEdxStop > (1.0 + thr) * dEdxStart))
702  {
703  mf::LogVerbatim("pma::Track3D") << "Auto-flip fired (2), thr: " << (1.0+thr) << ", value: " << dEdxStop/dEdxStart;
704  Flip(); // particle stops at the front of the track
705  }
706 }
707 
708 bool pma::Track3D::AutoFlip(std::vector< pma::Track3D* >& allTracks, pma::Track3D::EDirection dir, double thr, unsigned int n)
709 {
710  unsigned int nViews = 3;
711  pma::dedx_map dedxInViews[3];
712  for (unsigned int i = 0; i < nViews; i++)
713  {
714  GetRawdEdxSequence(dedxInViews[i], i, 1);
715  }
716  unsigned int bestView = 2;
717  if (dedxInViews[0].size() > 2 * dedxInViews[2].size()) bestView = 0;
718  if (dedxInViews[1].size() > 2 * dedxInViews[2].size()) bestView = 1;
719 
720  std::vector< std::vector<double> > dedx;
721  for (size_t i = 0; i < size(); i++)
722  {
723  auto it = dedxInViews[bestView].find(i);
724  if (it != dedxInViews[bestView].end())
725  {
726  dedx.push_back(it->second);
727  }
728  }
729  if (!dedx.empty()) dedx.pop_back();
730 
731  float dEdxStart = 0.0F, dEdxStop = 0.0F;
732  float dEStart = 0.0F, dxStart = 0.0F;
733  float dEStop = 0.0F, dxStop = 0.0F;
734  if (dedx.size() > 4)
735  {
736  if (!n) // use default options
737  {
738  if (dedx.size() > 30) n = 12;
739  else if (dedx.size() > 20) n = 8;
740  else if (dedx.size() > 10) n = 4;
741  else n = 3;
742  }
743 
744  size_t k = (dedx.size() - 2) >> 1;
745  if (n > k) n = k;
746 
747  for (size_t i = 1, j = 0; j < n; i++, j++)
748  {
749  dEStart += dedx[i][5]; dxStart += dedx[i][6];
750  }
751  if (dxStart > 0.0F) dEdxStart = dEStart / dxStart;
752 
753  for (size_t i = dedx.size() - 2, j = 0; j < n; i--, j++)
754  {
755  dEStop += dedx[i][5]; dxStop += dedx[i][6];
756  }
757  if (dxStop > 0.0F) dEdxStop = dEStop / dxStop;
758  }
759  else if (dedx.size() == 4)
760  {
761  dEStart = dedx[0][5] + dedx[1][5]; dxStart = dedx[0][6] + dedx[1][6];
762  dEStop = dedx[2][5] + dedx[3][5]; dxStop = dedx[2][6] + dedx[3][6];
763  if (dxStart > 0.0F) dEdxStart = dEStart / dxStart;
764  if (dxStop > 0.0F) dEdxStop = dEStop / dxStop;
765 
766  }
767  else if (dedx.size() > 1)
768  {
769  if (dedx.front()[2] > 0.0F) dEdxStart = dedx.front()[5] / dedx.front()[6];
770  if (dedx.back()[2] > 0.0F) dEdxStop = dedx.back()[5] / dedx.back()[6];
771  }
772  else return false;
773 
774  bool done = false;
775  if ((dir == pma::Track3D::kForward) && ((1.0 + thr) * dEdxStop < dEdxStart))
776  {
777  mf::LogVerbatim("pma::Track3D") << "Auto-flip fired (1), thr: " << (1.0+thr) << ", value: " << dEdxStart/dEdxStop;
778  done = Flip(allTracks); // particle stops at the end of the track
779  }
780  if ((dir == pma::Track3D::kBackward) && (dEdxStop > (1.0 + thr) * dEdxStart))
781  {
782  mf::LogVerbatim("pma::Track3D") << "Auto-flip fired (2), thr: " << (1.0+thr) << ", value: " << dEdxStop/dEdxStart;
783  done = Flip(allTracks); // particle stops at the front of the track
784  }
785  return done;
786 }
787 
788 double pma::Track3D::TestHitsMse(const std::vector< art::Ptr<recob::Hit> >& hits, bool normalized) const
789 {
790  if (hits.empty())
791  {
792  mf::LogWarning("pma::Track3D") << "TestHitsMse(): Empty cluster.";
793  return -1.0;
794  }
795 
796  double mse = 0.0;
797  for (auto const & h : hits)
798  {
799  unsigned int cryo = h->WireID().Cryostat;
800  unsigned int tpc = h->WireID().TPC;
801  unsigned int view = h->WireID().Plane;
802  unsigned int wire = h->WireID().Wire;
803  float drift = h->PeakTime();
804 
805  mse += Dist2(pma::WireDriftToCm(wire, drift, view, tpc, cryo), view, tpc, cryo);
806  }
807  if (normalized) return mse / hits.size();
808  else return mse;
809 }
810 
811 unsigned int pma::Track3D::TestHits(const std::vector< art::Ptr<recob::Hit> >& hits, double dist) const
812 {
813  if (hits.empty())
814  {
815  mf::LogWarning("pma::Track3D") << "TestHits(): Empty cluster.";
816  return 0;
817  }
818 
819  TVector3 p3d;
820  double tst, d2 = dist * dist;
821  unsigned int nhits = 0;
822  for (auto const & h : hits)
823  if (GetUnconstrainedProj3D(h, p3d, tst))
824  {
825  if (tst < d2) nhits++;
826  }
827 
828  return nhits;
829 }
830 
831 double pma::Track3D::Length(size_t start, size_t stop, size_t step) const
832 {
833  if (size() < 2) return 0.0;
834 
835  if (start > stop)
836  {
837  size_t tmp = stop; stop = start; start = tmp;
838  }
839  if (start >= size() - 1) return 0.0;
840  if (stop >= size()) stop = size() - 1;
841 
842  double result = 0.0;
843  pma::Hit3D *h1 = 0, *h0 = fHits[start];
844 
845  if (!step) step = 1;
846  size_t i = start + step;
847  while (i <= stop)
848  {
849  h1 = fHits[i];
850  result += sqrt(pma::Dist2(h1->Point3D(), h0->Point3D()));
851  h0 = h1;
852  i += step;
853  }
854  if (i - step < stop) // last step jumped beyond the stop point
855  { // so need to add the last short segment
856  result += sqrt(pma::Dist2(h0->Point3D(), back()->Point3D()));
857  }
858  return result;
859 }
860 
861 int pma::Track3D::NextHit(int index, unsigned int view, bool inclDisabled) const
862 {
863  pma::Hit3D* hit = 0;
864  if (index < -1) index = -1;
865  while (++index < (int)size()) // look for the next index of hit from the view
866  {
867  hit = fHits[index];
868  if (hit->View2D() == view)
869  {
870  if (inclDisabled) break;
871  else if (hit->IsEnabled()) break;
872  }
873  }
874  return index;
875 }
876 
877 int pma::Track3D::PrevHit(int index, unsigned int view, bool inclDisabled) const
878 {
879  pma::Hit3D* hit = 0;
880  if (index > (int)size()) index = (int)size();
881  while (--index >= 0) // look for the prev index of hit from the view
882  {
883  hit = fHits[index];
884  if (hit->View2D() == view)
885  {
886  if (inclDisabled) break;
887  else if (hit->IsEnabled()) break;
888  }
889  }
890  return index;
891 }
892 
893 double pma::Track3D::HitDxByView(size_t index, unsigned int view,
894  pma::Track3D::EDirection dir, bool secondDir) const
895 {
896  pma::Hit3D* nexthit = 0;
897  pma::Hit3D* hit = fHits[index];
898 
899  if (hit->View2D() != view)
900  {
901  mf::LogWarning("pma::Track3D") << "Function used with the hit not matching specified view.";
902  }
903 
904  double dx = 0.0; // [cm]
905  bool hitFound = false;
906  int i = index;
907  switch (dir)
908  {
909  case pma::Track3D::kForward:
910  while (!hitFound && (++i < (int)size()))
911  {
912  nexthit = fHits[i];
913  dx += sqrt(pma::Dist2(hit->Point3D(), nexthit->Point3D()));
914 
915  if (nexthit->View2D() == view) hitFound = true;
916  else hitFound = false;
917 
918  hit = nexthit;
919  }
920  if (!hitFound)
921  {
922  if (!secondDir) dx = 0.5 * HitDxByView(index, view, pma::Track3D::kBackward, true);
923  else { dx = Length(); mf::LogWarning("pma::Track3D") << "Single hit in this view."; }
924  }
925  break;
926 
927  case pma::Track3D::kBackward:
928  while (!hitFound && (--i >= 0))
929  {
930  nexthit = fHits[i];
931  dx += sqrt(pma::Dist2(hit->Point3D(), nexthit->Point3D()));
932 
933  if (nexthit->View2D() == view) hitFound = true;
934  else hitFound = false;
935 
936  hit = nexthit;
937  }
938  if (!hitFound)
939  {
940  if (!secondDir) dx = 0.5 * HitDxByView(index, view, pma::Track3D::kForward, true);
941  else { dx = Length(); mf::LogWarning("pma::Track3D") << "Single hit in this view."; }
942  }
943  break;
944 
945  default:
946  mf::LogError("pma::Track3D") << "Direction undefined.";
947  break;
948  }
949  return dx;
950 }
951 
952 double pma::Track3D::HitDxByView(size_t index, unsigned int view) const
953 {
954  if (index < size())
955  {
956  return 0.5 * (HitDxByView(index, view, pma::Track3D::kForward)
957  + HitDxByView(index, view, pma::Track3D::kBackward));
958  }
959  else
960  {
961  mf::LogError("pma::Track3D") << "Hit index out of range.";
962  return 0.0;
963  }
964 }
965 
966 pma::Segment3D* pma::Track3D::NextSegment(pma::Node3D* vtx) const
967 {
968  pma::Segment3D* seg = 0;
969  unsigned int nCount = vtx->NextCount();
970  unsigned int k = 0;
971  while (k < nCount)
972  {
973  seg = static_cast< pma::Segment3D* >(vtx->Next(k));
974  if (seg && (seg->Parent() == this)) return seg;
975  k++;
976  }
977  return 0;
978 }
979 
980 pma::Segment3D* pma::Track3D::PrevSegment(pma::Node3D* vtx) const
981 {
982  if (vtx->Prev())
983  {
984  pma::Segment3D* seg = static_cast< pma::Segment3D* >(vtx->Prev());
985  if (seg->Parent() == this) return seg;
986  }
987  return 0;
988 }
989 
990 double pma::Track3D::GetRawdEdxSequence(
991  std::map< size_t, std::vector<double> > & dedx,
992  unsigned int view, unsigned int skip,
993  bool inclDisabled) const
994 {
995  dedx.clear();
996 
997  if (!size()) return 0.0;
998 
999  size_t step = 1;
1000 
1001  pma::Hit3D* hit = 0;
1002 
1003  double rv, dr, dR, dq, dEq, qSkipped = 0.0;
1004 
1005  size_t j = NextHit(-1, view, inclDisabled), s = skip;
1006  if (j >= size()) return 0.0F; // no charged hits at all
1007  while (j < size()) // look for the first hit index
1008  {
1009  hit = fHits[j];
1010  dq = hit->SummedADC();
1011  if (s) { qSkipped += dq; s--; }
1012  else break;
1013 
1014  j = NextHit(j, view, inclDisabled);
1015  }
1016 
1017  size_t jmax = PrevHit(size(), view, inclDisabled);
1018 
1019  std::vector< size_t > indexes;
1020  TVector3 p0(0., 0., 0.), p1(0., 0., 0.);
1021  TVector2 c0(0., 0.), c1(0., 0.);
1022  while (j <= jmax)
1023  {
1024  indexes.clear(); // prepare to collect hit indexes used for this dE/dx entry
1025 
1026  indexes.push_back(j);
1027  hit = fHits[j];
1028 
1029  p0 = hit->Point3D();
1030  p1 = hit->Point3D();
1031 
1032  c0.Set(hit->Wire(), hit->PeakTime());
1033  c1.Set(hit->Wire(), hit->PeakTime());
1034 
1035  dEq = hit->SummedADC(); // [now it is ADC sum]
1036 
1037  dr = HitDxByView(j, view, pma::Track3D::kForward); // protection against hits on the same position
1038  rv = HitDxByView(j, view); // dx seen by j-th hit
1039  fHits[j]->fDx = rv;
1040  dR = rv;
1041 
1042  size_t m = 1; // number of hits with charge > 0
1043  while (((m < step) || (dR < 0.1) || (dr == 0.0)) && (j <= jmax))
1044  {
1045  j = NextHit(j, view); // just next, even if tagged as outlier
1046  if (j > jmax) break; // no more hits in this view
1047 
1048  hit = fHits[j];
1049  if (!inclDisabled && !hit->IsEnabled())
1050  {
1051  if (dr == 0.0) continue;
1052  else break;
1053  }
1054  indexes.push_back(j);
1055 
1056  p1 = hit->Point3D();
1057 
1058  c1.Set(hit->Wire(), hit->PeakTime());
1059 
1060  dq = hit->SummedADC();
1061 
1062  dEq += dq;
1063 
1064  dr = HitDxByView(j, view, pma::Track3D::kForward);
1065  rv = HitDxByView(j, view);
1066  fHits[j]->fDx = rv;
1067  dR += rv;
1068  m++;
1069  }
1070  p0 += p1; p0 *= 0.5;
1071  c0 += c1; c0 *= 0.5;
1072 
1073  double range = Length(0, j);
1074 
1075  std::vector<double> trk_section;
1076  trk_section.push_back(c0.X());
1077  trk_section.push_back(c0.Y());
1078  trk_section.push_back(p0.X());
1079  trk_section.push_back(p0.Y());
1080  trk_section.push_back(p0.Z());
1081  trk_section.push_back(dEq);
1082  trk_section.push_back(dR);
1083  trk_section.push_back(range);
1084 
1085  for (auto const idx : indexes) dedx[idx] = trk_section;
1086 
1087  j = NextHit(j, view, inclDisabled);
1088  }
1089 
1090  return qSkipped;
1091 }
1092 
1093 std::vector<float> pma::Track3D::DriftsOfWireIntersection(unsigned int wire, unsigned int view) const
1094 {
1095  std::vector<float> drifts;
1096  for (size_t i = 0; i < fNodes.size() - 1; i++)
1097  {
1098  int tpc = fNodes[i]->TPC(), cryo = fNodes[i]->Cryo();
1099  if ((tpc != fNodes[i + 1]->TPC()) || (cryo != fNodes[i + 1]->Cryo())) continue;
1100 
1101  TVector2 p0 = pma::CmToWireDrift(
1102  fNodes[i]->Projection2D(view).X(), fNodes[i]->Projection2D(view).Y(),
1103  view, fNodes[i]->TPC(), fNodes[i]->Cryo());
1104  TVector2 p1 = pma::CmToWireDrift(
1105  fNodes[i + 1]->Projection2D(view).X(), fNodes[i + 1]->Projection2D(view).Y(),
1106  view, fNodes[i + 1]->TPC(), fNodes[i + 1]->Cryo());
1107 
1108  if ((p0.X() - wire) * (p1.X() - wire) <= 0.0)
1109  {
1110  double frac_w = (wire - p0.X()) / (p1.X() - p0.X());
1111  double d = p0.Y() + frac_w * (p1.Y() - p0.Y());
1112  drifts.push_back((float)d);
1113  }
1114  }
1115  return drifts;
1116 }
1117 
1118 size_t pma::Track3D::CompleteMissingWires(unsigned int view)
1119 {
1120  int dPrev, dw, w, wx, wPrev, i = NextHit(-1, view);
1121  //TVector2 projPoint;
1122 
1123  pma::Hit3D* hitPrev = 0;
1124  pma::Hit3D* hit = 0;
1125 
1126  std::vector< pma::Hit3D* > missHits;
1127 
1128  while (i < (int)size())
1129  {
1130  hitPrev = hit;
1131  hit = fHits[i];
1132 
1133  //projPoint = pma::CmToWireDrift(
1134  // hit->Projection2D().X(), hit->Projection2D().Y(),
1135  // hit->View2D(), hit->TPC(), hit->Cryo());
1136 
1137  if (hit->View2D() == view) // &&
1138  //(projPoint.Y() >= hit->Hit2DPtr()->StartTick()-1) && // cannot use Hit2DPtr here
1139  //(projPoint.Y() <= hit->Hit2DPtr()->EndTick()+1))
1140  {
1141  w = hit->Wire();
1142  if (hitPrev)
1143  {
1144  wPrev = hitPrev->Wire();
1145  dPrev = (int)hitPrev->PeakTime();
1146  if (abs(w - wPrev) > 1)
1147  {
1148  if (w > wPrev) dw = 1;
1149  else dw = -1;
1150  wx = wPrev + dw;
1151  int k = 1;
1152  while (wx != w)
1153  {
1154  int peakTime = 0;
1155  unsigned int iWire = wx;
1156  std::vector<float> drifts = DriftsOfWireIntersection(wx, view);
1157  if (drifts.size())
1158  {
1159  peakTime = drifts[0];
1160  for (size_t d = 1; d < drifts.size(); d++)
1161  if (fabs(drifts[d] - dPrev) < fabs(peakTime - dPrev))
1162  peakTime = drifts[d];
1163  }
1164  else
1165  {
1166  mf::LogVerbatim("pma::Track3D")
1167  << "Track does not intersect with the wire " << wx;
1168  break;
1169  }
1170  pma::Hit3D* hmiss = new pma::Hit3D(
1171  iWire, view, hit->TPC(), hit->Cryo(), peakTime, 1.0, 1.0);
1172  missHits.push_back(hmiss);
1173  wx += dw; k++;
1174  }
1175  }
1176  }
1177  }
1178  else hit = hitPrev;
1179 
1180  i = NextHit(i, view, true);
1181  while ((i < (int)size()) && (hit->TPC() != fHits[i]->TPC()))
1182  {
1183  hitPrev = hit; hit = fHits[i];
1184  i = NextHit(i, view, true);
1185  }
1186  }
1187 
1188  if (missHits.size())
1189  {
1190  for (size_t hi = 0; hi < missHits.size(); hi++)
1191  {
1192  //mf::LogVerbatim("pma::Track3D") << "Fill missing hit on wire " << missHits[hi]->Wire();
1193  push_back(missHits[hi]);
1194  }
1195 
1196  MakeProjection();
1197  SortHits();
1198  }
1199 
1200  return missHits.size();
1201 }
1202 
1203 void pma::Track3D::AddNode(pma::Node3D* node)
1204 {
1205  fNodes.push_back(node);
1206  if (fNodes.size() > 1) RebuildSegments();
1207 }
1208 
1209 bool pma::Track3D::AddNode(void)
1210 {
1211  pma::Segment3D* seg;
1212  pma::Segment3D* maxSeg = 0;
1213 
1214  size_t si = 0;
1215  while (si < fSegments.size())
1216  {
1217  if (!fSegments[si]->IsFrozen())
1218  {
1219  maxSeg = fSegments[si];
1220  break;
1221  }
1222  else si++;
1223  }
1224  if (!maxSeg) return false;
1225 
1226  unsigned int nHitsByView[3];
1227  unsigned int nHits, maxHits = 0;
1228  unsigned int vIndex = 0, segHits, maxSegHits = 0;
1229  float segLength, maxLength = maxSeg->Length();
1230  for (unsigned int i = si + 1; i < fNodes.size(); i++)
1231  {
1232  seg = static_cast< pma::Segment3D* >(fNodes[i]->Prev());
1233  if (seg->IsFrozen()) continue;
1234 
1235  nHitsByView[0] = seg->NEnabledHits(geo::kU);
1236  nHitsByView[1] = seg->NEnabledHits(geo::kV);
1237  nHitsByView[2] = seg->NEnabledHits(geo::kZ);
1238  segHits = nHitsByView[0] + nHitsByView[1] + nHitsByView[2];
1239 
1240  if (segHits < 15)
1241  {
1242  if ((nHitsByView[0] == 0) && ((nHitsByView[1] < 4) || (nHitsByView[2] < 4))) continue;
1243  if ((nHitsByView[1] == 0) && ((nHitsByView[0] < 4) || (nHitsByView[2] < 4))) continue;
1244  if ((nHitsByView[2] == 0) && ((nHitsByView[0] < 4) || (nHitsByView[1] < 4))) continue;
1245  }
1246 
1247  nHits = fNodes[i]->NEnabledHits() + seg->NEnabledHits() + fNodes[i-1]->NEnabledHits();
1248 
1249  if (nHits > maxHits)
1250  {
1251  maxHits = nHits;
1252  maxLength = seg->Length();
1253  maxSegHits = segHits;
1254  maxSeg = seg;
1255  vIndex = i;
1256  }
1257  else if (nHits == maxHits)
1258  {
1259  segLength = seg->Length();
1260  if (segLength > maxLength)
1261  {
1262  maxLength = segLength;
1263  maxSegHits = segHits;
1264  maxSeg = seg;
1265  vIndex = i;
1266  }
1267  }
1268  }
1269 
1270  if (maxSegHits > 1)
1271  {
1272  maxSeg->SortHits();
1273 
1274  nHitsByView[0] = maxSeg->NEnabledHits(geo::kU);
1275  nHitsByView[1] = maxSeg->NEnabledHits(geo::kV);
1276  nHitsByView[2] = maxSeg->NEnabledHits(geo::kZ);
1277 
1278 
1279  unsigned int maxViewIdx = 2, midViewIdx = 2;
1280  if ((nHitsByView[2] >= nHitsByView[1]) && (nHitsByView[1] >= nHitsByView[0])) { maxViewIdx = 2; midViewIdx = 1; }
1281  else if ((nHitsByView[1] >= nHitsByView[2]) && (nHitsByView[2] >= nHitsByView[0])) { maxViewIdx = 1; midViewIdx = 2; }
1282  else if ((nHitsByView[0] >= nHitsByView[2]) && (nHitsByView[2] >= nHitsByView[1])) { maxViewIdx = 0; midViewIdx = 2; }
1283  else if ((nHitsByView[2] >= nHitsByView[0]) && (nHitsByView[0] >= nHitsByView[1])) { maxViewIdx = 2; midViewIdx = 0; }
1284  else if ((nHitsByView[0] >= nHitsByView[1]) && (nHitsByView[1] >= nHitsByView[2])) { maxViewIdx = 0; midViewIdx = 1; }
1285  else if ((nHitsByView[1] >= nHitsByView[0]) && (nHitsByView[0] >= nHitsByView[2])) { maxViewIdx = 1; midViewIdx = 0; }
1286  if (nHitsByView[midViewIdx] < 2) midViewIdx = maxViewIdx;
1287 
1288  if (nHitsByView[midViewIdx] < 2) { mf::LogVerbatim("pma::Track3D") << "AddNode(): too few hits."; return false; }
1289 
1290  unsigned int mHits[3] = { 0, 0, 0 };
1291  unsigned int halfIndex = (nHitsByView[midViewIdx] >> 1) - 1;
1292  unsigned int n = 0, i = 0, i0 = 0, i1 = 0;
1293  while (i < maxSeg->NHits() - 1)
1294  {
1295  if (maxSeg->Hit(i).IsEnabled())
1296  {
1297  mHits[maxSeg->Hit(i).View2D()]++;
1298  if (maxSeg->Hit(i).View2D() == midViewIdx)
1299  {
1300  if (n == halfIndex) break;
1301  n++;
1302  }
1303  }
1304  i++;
1305  }
1306 
1307  i0 = i; i++;
1308  while ((i < maxSeg->NHits()) && !((maxSeg->Hit(i).View2D() == midViewIdx) && maxSeg->Hit(i).IsEnabled()))
1309  {
1310  i++;
1311  }
1312  i1 = i;
1313 
1314  if (!nHitsByView[0])
1315  {
1316  if (nHitsByView[1] && (mHits[1] < 2)) { mf::LogVerbatim("pma::Track3D") << "AddNode(): low Ind2 hits."; return false; }
1317  if (nHitsByView[2] && (mHits[2] < 2)) { mf::LogVerbatim("pma::Track3D") << "AddNode(): low Coll hits."; return false; }
1318  }
1319 
1320  maxSeg->SetProjection(maxSeg->Hit(i0));
1321  maxSeg->SetProjection(maxSeg->Hit(i1));
1322 
1323  unsigned int tpc = maxSeg->Hit(i0).TPC();
1324  unsigned int cryo = maxSeg->Hit(i0).Cryo();
1325 
1326  pma::Node3D* p = new pma::Node3D((maxSeg->Hit(i0).Point3D() + maxSeg->Hit(i1).Point3D()) * 0.5, tpc, cryo, false, fNodes[vIndex]->GetDriftShift());
1327 
1328  //mf::LogVerbatim("pma::Track3D") << "add node x:" << p->Point3D().X()
1329  // << " y:" << p->Point3D().Y() << " z:" << p->Point3D().Z();
1330  fNodes.insert(fNodes.begin() + vIndex, p);
1331 
1332  maxSeg->AddNext(fNodes[vIndex]);
1333 
1334  seg = new pma::Segment3D(this, fNodes[vIndex], fNodes[vIndex + 1]);
1335  fSegments.insert(fSegments.begin() + vIndex, seg);
1336 
1337  return true;
1338  }
1339  else return false;
1340 }
1341 
1342 void pma::Track3D::InsertNode(
1343  TVector3 const & p3d, size_t at_idx,
1344  unsigned int tpc, unsigned int cryo)
1345 {
1346  //std::cout << " insert after " << at_idx << " in " << fNodes.size() << std::endl;
1347  pma::Node3D* vtx = new pma::Node3D(p3d, tpc, cryo, false, fNodes[at_idx]->GetDriftShift());
1348  fNodes.insert(fNodes.begin() + at_idx, vtx);
1349  //std::cout << " inserted " << std::endl;
1350 
1351  if (fNodes.size() > 1) RebuildSegments();
1352  //std::cout << " rebuild ok " << std::endl;
1353 }
1354 
1355 bool pma::Track3D::RemoveNode(size_t idx)
1356 {
1357  if ((fNodes.size() > 1) && (idx < fNodes.size()))
1358  {
1359  pma::Node3D* vtx = fNodes[idx];
1360  fNodes.erase(fNodes.begin() + idx);
1361  RebuildSegments();
1362 
1363  if (!vtx->NextCount() && !vtx->Prev()) { delete vtx; }
1364 
1365  return true;
1366  }
1367  else return false;
1368 }
1369 
1370 pma::Track3D* pma::Track3D::Split(size_t idx, bool try_start_at_idx)
1371 {
1372  if (!idx || (idx + 1 >= fNodes.size())) return 0;
1373 
1374  pma::Node3D* n = 0;
1375  pma::Track3D* t0 = new pma::Track3D();
1376  t0->fT0 = fT0; t0->fT0Flag = fT0Flag; t0->fTag = fTag;
1377 
1378  for (size_t i = 0; i < idx; ++i)
1379  {
1380  n = fNodes.front();
1381  n->ClearAssigned();
1382 
1383  pma::Segment3D* s = static_cast< pma::Segment3D* >(n->Prev());
1384  if (s && (s->Parent() == this)) s->RemoveNext(n);
1385 
1386  size_t k = 0;
1387  while (k < n->NextCount())
1388  {
1389  s = static_cast< pma::Segment3D* >(n->Next(k));
1390  if (s->Parent() == this) n->RemoveNext(s);
1391  else k++;
1392  }
1393 
1394  fNodes.erase(fNodes.begin());
1395  t0->fNodes.push_back(n);
1396  }
1397 
1398  n = fNodes.front();
1399  t0->fNodes.push_back(new pma::Node3D(n->Point3D(), n->TPC(), n->Cryo(), false, n->GetDriftShift()));
1400  t0->RebuildSegments();
1401  RebuildSegments();
1402 
1403  size_t h = 0;
1404  while (h < size())
1405  {
1406  pma::Hit3D* h3d = fHits[h];
1407  unsigned int view = h3d->View2D(), tpc = h3d->TPC(), cryo = h3d->Cryo();
1408  double dist2D_old = Dist2(h3d->Point2D(), view, tpc, cryo);
1409  double dist2D_new = t0->Dist2(h3d->Point2D(), view, tpc, cryo);
1410 
1411  if ((dist2D_new < dist2D_old) && t0->HasTPC(tpc)) t0->push_back(release_at(h));
1412  else if (!HasTPC(tpc) && t0->HasTPC(tpc)) t0->push_back(release_at(h));
1413  else h++;
1414  }
1415 
1416  bool passed = true;
1417  if (HasTwoViews() && t0->HasTwoViews())
1418  {
1419  if (try_start_at_idx && t0->CanFlip())
1420  {
1421  mf::LogVerbatim("pma::Track3D") << " attach new t0 and this trk to a common start node";
1422  t0->Flip();
1423  passed = t0->AttachTo(fNodes.front());
1424  }
1425  else
1426  {
1427  mf::LogVerbatim("pma::Track3D") << " attach this trk to the new t0 end node";
1428  passed = AttachTo(t0->fNodes.back());
1429  }
1430  }
1431  else
1432  {
1433  mf::LogVerbatim("pma::Track3D") << " single-view track";
1434  passed = false;
1435  }
1436 
1437  if (!passed)
1438  {
1439  mf::LogVerbatim("pma::Track3D") << " undo split";
1440  while (t0->size()) push_back(t0->release_at(0));
1441 
1442  for (size_t i = 0; i < idx; ++i)
1443  {
1444  fNodes.insert(fNodes.begin() + i, t0->fNodes.front());
1445  t0->fNodes.erase(t0->fNodes.begin());
1446  }
1447 
1448  RebuildSegments();
1449  delete t0;
1450  t0 = 0;
1451  }
1452 
1453  return t0;
1454 }
1455 
1456 bool pma::Track3D::AttachTo(pma::Node3D* vStart, bool noFlip)
1457 {
1458  pma::Node3D* vtx = fNodes.front();
1459 
1460  if (vtx == vStart) return true; // already connected!
1461 
1462  pma::Track3D* rootThis = GetRoot();
1463  if (vStart->Prev())
1464  {
1465  pma::Track3D* rootPrev = static_cast< pma::Segment3D* >(vStart->Prev())->Parent()->GetRoot();
1466  if (rootPrev->IsAttachedTo(rootThis))
1467  {
1468  //std::cout << "*** loop with prev ***" << std::endl;
1469  return false;
1470  }
1471  }
1472  else if (vStart->NextCount())
1473  {
1474  pma::Track3D* rootNext = static_cast< pma::Segment3D* >(vStart->Next(0))->Parent()->GetRoot();
1475  if (rootNext->IsAttachedTo(rootThis))
1476  {
1477  //std::cout << "*** loop with next ***" << std::endl;
1478  return false;
1479  }
1480  }
1481  else
1482  {
1483  //std::cout << "*** isolated vertex ***" << std::endl;
1484  return false;
1485  }
1486 
1487  for (auto n : fNodes) if (n == vStart) { mf::LogError("pma::Track3D") << "Don't create loops!"; return false; }
1488 
1489  if ( !noFlip && CanFlip() && (vStart->TPC() == fNodes.back()->TPC()) &&
1490  (pma::Dist2(vtx->Point3D(), vStart->Point3D()) > pma::Dist2(fNodes.back()->Point3D(), vStart->Point3D())) &&
1491  (fNodes.back()->NextCount() == 0) )
1492  {
1493  mf::LogError("pma::Track3D") << "Flip, endpoint closer to vStart.";
1494  //std::cout << "Flip, endpoint closer to vStart." << std::endl;
1495  Flip();
1496  }
1497 
1498  if (vStart->TPC() == vtx->TPC()) return AttachToSameTPC(vStart);
1499  else return AttachToOtherTPC(vStart);
1500 }
1501 
1502 bool pma::Track3D::AttachToOtherTPC(pma::Node3D* vStart)
1503 {
1504  if (fNodes.front()->Prev()) return false;
1505 
1506  mf::LogVerbatim("pma::Track3D") << "Attach to track in another TPC.";
1507 
1508  fNodes.insert(fNodes.begin(), vStart);
1509 
1510  fSegments.insert(fSegments.begin(), new pma::Segment3D(this, fNodes[0], fNodes[1]));
1511 
1512  return true;
1513 }
1514 
1515 bool pma::Track3D::AttachToSameTPC(pma::Node3D* vStart)
1516 {
1517  pma::Node3D* vtx = fNodes.front();
1518 
1519  if (vtx->Prev())
1520  {
1521  pma::Segment3D* segPrev = static_cast< pma::Segment3D* >(vtx->Prev());
1522  pma::Track3D* tpPrev = segPrev->Parent();
1523  if (tpPrev->NextSegment(vtx))
1524  {
1525  //std::cout << "Do not reattach from vtx inner in another track." << std::endl;
1526  return false;
1527  }
1528  else if (tpPrev->CanFlip())
1529  {
1530  //std::cout << "Flip local." << std::endl;
1531  tpPrev->Flip();
1532  } // flip in local vtx, no problem
1533  else
1534  {
1535  if (vStart->Prev())
1536  {
1537  pma::Segment3D* segNew = static_cast< pma::Segment3D* >(vStart->Prev());
1538  pma::Track3D* tpNew = segNew->Parent();
1539  if (tpNew->NextSegment(vStart))
1540  {
1541  //std::cout << "Do not reattach to vStart inner in another track when vtx also has Prev()." << std::endl;
1542  return false;
1543  }
1544  else if (tpNew->CanFlip()) // flip in remote vStart, no problem
1545  {
1546  tpNew->Flip();
1547 
1548  mf::LogVerbatim("pma::Track3D") << "Reconnect prev to vStart.";
1549  //std::cout << "Reconnect prev to vStart." << std::endl;
1550  tpPrev->fNodes[tpPrev->fNodes.size() - 1] = vStart;
1551  segPrev->AddNext(vStart);
1552  }
1553  else
1554  {
1555  mf::LogError("pma::Track3D") << "Flips in vtx and vStart not possible, cannot attach.";
1556  return false;
1557  }
1558  }
1559  else
1560  {
1561  mf::LogVerbatim("pma::Track3D") << "Reconnect prev to vStart.";
1562  //std::cout << "Reconnect prev to vStart." << std::endl;
1563  tpPrev->fNodes[tpPrev->fNodes.size() - 1] = vStart;
1564  segPrev->AddNext(vStart);
1565  }
1566  }
1567  }
1568 
1569  while (vtx->NextCount()) // reconnect nexts to vStart
1570  {
1571  pma::Segment3D* seg = static_cast< pma::Segment3D* >(vtx->Next(0));
1572  pma::Track3D* trk = seg->Parent();
1573 
1574  vtx->RemoveNext(seg);
1575  trk->fNodes[0] = vStart;
1576  vStart->AddNext(seg);
1577  }
1578 
1579  if (vtx->NextCount() || vtx->Prev()) // better throw here
1580  {
1581  throw cet::exception("pma::Track3D") << "Something is still using disconnected vertex.";
1582  }
1583  else delete vtx; // ok
1584  return true;
1585 }
1586 
1587 bool pma::Track3D::AttachBackTo(pma::Node3D* vStart)
1588 {
1589  pma::Node3D* vtx = fNodes.back();
1590 
1591  if (vtx == vStart) return true; // already connected!
1592 
1593  pma::Track3D* rootThis = GetRoot();
1594  if (vStart->Prev())
1595  {
1596  pma::Track3D* rootPrev = static_cast< pma::Segment3D* >(vStart->Prev())->Parent()->GetRoot();
1597  if (rootPrev->IsAttachedTo(rootThis))
1598  {
1599  //std::cout << "*** loop with prev ***" << std::endl;
1600  return false;
1601  }
1602  }
1603  else if (vStart->NextCount())
1604  {
1605  pma::Track3D* rootNext = static_cast< pma::Segment3D* >(vStart->Next(0))->Parent()->GetRoot();
1606  if (rootNext->IsAttachedTo(rootThis))
1607  {
1608  //std::cout << "*** loop with next ***" << std::endl;
1609  return false;
1610  }
1611  }
1612  else
1613  {
1614  //std::cout << "*** isolated vertex ***" << std::endl;
1615  return false;
1616  }
1617 
1618  for (auto n : fNodes) if (n == vStart) { mf::LogError("pma::Track3D") << "Don't create loops!"; return false; }
1619 
1620  if (vStart->TPC() == vtx->TPC()) return AttachBackToSameTPC(vStart);
1621  else return AttachBackToOtherTPC(vStart);
1622 }
1623 
1624 bool pma::Track3D::AttachBackToOtherTPC(pma::Node3D* vStart)
1625 {
1626  if (vStart->Prev()) return false;
1627 
1628  mf::LogVerbatim("pma::Track3D") << "Attach to track in another TPC.";
1629 
1630  fNodes.push_back(vStart);
1631 
1632  size_t idx = fNodes.size() - 1;
1633  fSegments.push_back(new pma::Segment3D(this, fNodes[idx - 1], fNodes[idx]));
1634 
1635  return true;
1636 }
1637 
1638 bool pma::Track3D::AttachBackToSameTPC(pma::Node3D* vStart)
1639 {
1640  pma::Node3D* vtx = fNodes.back();
1641 
1642  if (vStart->Prev())
1643  {
1644  pma::Segment3D* seg = static_cast< pma::Segment3D* >(vStart->Prev());
1645  pma::Track3D* tp = seg->Parent();
1646  if (tp->NextSegment(vStart))
1647  {
1648  mf::LogError("pma::Track3D") << "Cannot attach back to inner node of other track.";
1649  //std::cout << "Cannot attach back to inner node of other track." << std::endl;
1650  return false;
1651  }
1652 
1653  if (tp->CanFlip()) tp->Flip(); // flip in remote vStart, no problem
1654  else
1655  {
1656  mf::LogError("pma::Track3D") << "Flip not possible, cannot attach.";
1657  //std::cout << "Flip not possible, cannot attach." << std::endl;
1658  return false;
1659  }
1660  }
1661  fNodes[fNodes.size() - 1] = vStart;
1662  //RebuildSegments();
1663  fSegments[fSegments.size() - 1]->AddNext(vStart);
1664 
1665  while (vtx->NextCount()) // reconnect nexts to vStart
1666  {
1667  pma::Segment3D* seg = static_cast< pma::Segment3D* >(vtx->Next(0));
1668  pma::Track3D* trk = seg->Parent();
1669 
1670  vtx->RemoveNext(seg);
1671  trk->fNodes[0] = vStart;
1672  vStart->AddNext(seg);
1673  }
1674 
1675  if (vtx->NextCount() || vtx->Prev())
1676  {
1677  throw cet::exception("pma::Track3D") << "Something is still using disconnected vertex." << std::endl;
1678  }
1679  else delete vtx; // ok
1680 
1681  return true;
1682 }
1683 
1684 void pma::Track3D::ExtendWith(pma::Track3D* src)
1685 {
1686  if (src->fNodes.front()->Prev())
1687  {
1688  throw cet::exception("pma::Track3D") << "Cant extend with a track being a daughter of another track." << std::endl;
1689  }
1690 
1691  src->DeleteSegments(); // disconnect from vertices and delete
1692  for (size_t i = 0; i < src->fNodes.size(); ++i)
1693  {
1694  fNodes.push_back(src->fNodes[i]);
1695 
1696  size_t idx = fNodes.size() - 1;
1697  fSegments.push_back(new pma::Segment3D(this, fNodes[idx - 1], fNodes[idx]));
1698  }
1699  src->fNodes.clear(); // just empty the node collection
1700 
1701  while (src->size()) { push_back(src->release_at(0)); } // move hits
1702 
1703  for (auto p : src->fAssignedPoints) { fAssignedPoints.push_back(p); } // move 3D reference points
1704  src->fAssignedPoints.clear();
1705 
1706  MakeProjection();
1707  SortHits();
1708 
1709  delete src;
1710 }
1711 
1712 pma::Track3D* pma::Track3D::GetRoot(void)
1713 {
1714  pma::Track3D* trk = 0;
1715 
1716  if (fNodes.size())
1717  {
1718  pma::Segment3D* seg = static_cast< pma::Segment3D* >(fNodes.front()->Prev());
1719  if (seg) trk = seg->Parent()->GetRoot();
1720  if (!trk) trk = this;
1721  }
1722  else throw cet::exception("pma::Track3D") << "Broken track, no nodes.";
1723 
1724  return trk;
1725 }
1726 
1727 bool pma::Track3D::GetBranches(std::vector< pma::Track3D const * >& branches, bool skipFirst) const
1728 {
1729  for (auto trk : branches) if (trk == this)
1730  {
1731  //std::cout << "HAS LOOP" << std::endl;
1732  throw cet::exception("pma::Track3D") << "Track tree has loop.";
1733  //mf::LogError("pma::Track3D") << "Track tree has loop.";
1734  //return false;
1735  }
1736 
1737  branches.push_back(this);
1738 
1739  size_t i0 = 0;
1740  if (skipFirst) i0 = 1;
1741 
1742  for (size_t i = i0; i < fNodes.size(); ++i)
1743  for (size_t n = 0; n < fNodes[i]->NextCount(); n++)
1744  {
1745  pma::Segment3D* seg = static_cast< pma::Segment3D* >(fNodes[i]->Next(n));
1746  if (seg && (seg->Parent() != this)) seg->Parent()->GetBranches(branches, true);
1747  }
1748 
1749  return true;
1750 }
1751 
1752 bool pma::Track3D::IsAttachedTo(pma::Track3D const * trk) const
1753 {
1754  if (trk == this) return true;
1755 
1756  std::vector< pma::Track3D const * > branchesThis, branchesTrk;
1757 
1758  if (!trk->GetBranches(branchesTrk)) return true; // has loop - better check the reason
1759  if (!GetBranches(branchesThis)) return true; // has loop - better check the reason
1760 
1761  for (auto bThis : branchesThis)
1762  for (auto bTrk : branchesTrk)
1763  if (bThis == bTrk) return true;
1764  return false;
1765 }
1766 
1767 bool pma::Track3D::HasRefPoint(TVector3* p) const
1768 {
1769  for (auto point : fAssignedPoints)
1770  if (point == p) return true;
1771  return false;
1772 }
1773 
1774 double pma::Track3D::GetMse(unsigned int view) const
1775 {
1776  double sumMse = 0.0;
1777  unsigned int nEnabledHits = 0;
1778  for (auto n : fNodes)
1779  {
1780  sumMse += n->SumDist2(view);
1781  nEnabledHits += n->NEnabledHits(view);
1782  }
1783  for (auto s : fSegments)
1784  {
1785  sumMse += s->SumDist2(view);
1786  nEnabledHits += s->NEnabledHits(view);
1787  }
1788 
1789  if (nEnabledHits) return sumMse / nEnabledHits;
1790  else return 0.0;
1791 }
1792 
1793 double pma::Track3D::GetObjFunction(float penaltyFactor) const
1794 {
1795  double sum = 0.0;
1796  float p = penaltyFactor * fPenaltyValue;
1797  for (size_t i = 0; i < fNodes.size(); i++)
1798  {
1799  sum += fNodes[i]->GetObjFunction(p, fEndSegWeight);
1800  }
1801  return sum / fNodes.size();
1802 }
1803 
1804 double pma::Track3D::Optimize(int nNodes, double eps, bool selAllHits, bool setAllNodes, size_t selSegHits, size_t selVtxHits)
1805 {
1806  if (!fNodes.size()) { mf::LogError("pma::Track3D") << "Track not initialized."; return 0.0; }
1807 
1808  if (!UpdateParams()) { mf::LogError("pma::Track3D") << "Track empty."; return 1.0e10; }
1809  double g0 = GetObjFunction(), g1 = 0.0;
1810  if (g0 == 0.0) return g0;
1811 
1812  //mf::LogVerbatim("pma::Track3D") << "objective function at opt start: " << g0;
1813 
1814  // set branching flag only at the beginning, optimization is not changing
1815  // that and new nodes are not branching
1816  for (auto n : fNodes) n->SetVertexToBranching(setAllNodes);
1817 
1818  bool stop = false;
1819  fMinSegStop = fSegments.size();
1820  fMaxSegStop = (int)(size() / fMaxSegStopFactor) + 1;
1821  do
1822  {
1823  if (selSegHits || selVtxHits) SelectRndHits(selSegHits, selVtxHits);
1824 
1825  bool stepDone = true;
1826  unsigned int stepIter = 0;
1827  do
1828  {
1829  double gstep = 1.0;
1830  unsigned int iter = 0;
1831  while ((gstep > eps) && (iter < 1000))
1832  {
1833  if ((fNodes.size() < 4) || (iter % 10 == 0)) MakeProjection();
1834  else MakeFastProjection();
1835 
1836  if (!UpdateParams()) { mf::LogError("pma::Track3D") << "Track empty."; return 0.0; }
1837 
1838  for (auto n : fNodes) n->Optimize(fPenaltyValue, fEndSegWeight);
1839 
1840  g1 = g0; g0 = GetObjFunction();
1841 
1842  //mf::LogVerbatim("pma::Track3D") << "obj fn: " << g0;
1843  if (g0 == 0.0F) { MakeProjection(); break; }
1844  gstep = fabs(g0 - g1) / g0;
1845  iter++;
1846  }
1847 
1848  stepIter++;
1849  if (fNodes.size() > 2)
1850  {
1851  stepDone = !(CheckEndSegment(pma::Track3D::kEnd) ||
1852  CheckEndSegment(pma::Track3D::kBegin));
1853  }
1854  } while (!stepDone && (stepIter < 5));
1855 
1856  if (selAllHits)
1857  {
1858  selAllHits = false;
1859  selSegHits = 0;
1860  selVtxHits = 0;
1861  if (SelectAllHits()) continue;
1862  }
1863 
1864  switch (nNodes)
1865  {
1866  case 0: stop = true; break; // just optimize existing vertices
1867 
1868  case -1: // grow and optimize until automatic stop condition
1869  mf::LogVerbatim("pma::Track3D") << "optimized segments: " << fSegments.size();
1870  if ((fSegments.size() >= fSegStopValue) ||
1871  (fSegments.size() >= fMaxSegStop))
1872  {
1873  stop = true;
1874  }
1875  else
1876  {
1877  if (!AddNode()) stop = true;
1878  }
1879  break;
1880 
1881  default: // grow and optimize until fixed number of vertices is added
1882  if (nNodes > 12)
1883  {
1884  if (AddNode()) { MakeProjection(); nNodes--; }
1885  else { mf::LogVerbatim("pma::Track3D") << "stop (3)"; stop = true; break; }
1886 
1887  if (AddNode())
1888  {
1889  MakeProjection(); nNodes--;
1890  if (AddNode()) nNodes--;
1891  }
1892  }
1893  else if (nNodes > 4)
1894  {
1895  if (AddNode()) { MakeProjection(); nNodes--; }
1896  else { mf::LogVerbatim("pma::Track3D") << "stop (2)"; stop = true; break; }
1897 
1898  if (AddNode()) nNodes--;
1899  }
1900  else
1901  {
1902  if (AddNode()) { nNodes--; }
1903  else { mf::LogVerbatim("pma::Track3D") << "stop (1)"; stop = true; break; }
1904  }
1905  break;
1906  }
1907 
1908  } while (!stop);
1909  //mf::LogVerbatim("pma::Track3D") << "Done (optimized segments: " << fSegments.size() << ").";
1910 
1911  MakeProjection();
1912  return GetObjFunction();
1913 }
1914 
1915 bool pma::Track3D::UpdateParamsInTree(bool skipFirst)
1916 {
1917  const size_t maxTreeDepth = 100; // really big tree...
1918  static size_t depth;
1919 
1920  pma::Node3D* vtx = fNodes.front();
1921  pma::Segment3D* segThis = 0;
1922  pma::Segment3D* seg = 0;
1923 
1924  bool isOK = true;
1925 
1926  if (skipFirst)
1927  {
1928  segThis = NextSegment(vtx);
1929  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
1930 
1931  if (++depth > maxTreeDepth)
1932  {
1933  mf::LogError("pma::Track3D") << "Tree depth above allowed max.";
1934  }
1935  }
1936  else { depth = 1; }
1937 
1938  while (vtx)
1939  {
1940  segThis = NextSegment(vtx);
1941 
1942  for (size_t i = 0; i < vtx->NextCount(); i++)
1943  {
1944  seg = static_cast< pma::Segment3D* >(vtx->Next(i));
1945  if (seg != segThis) { isOK &= seg->Parent()->UpdateParamsInTree(true); }
1946  }
1947 
1948  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
1949  else break;
1950  }
1951 
1952  if (!UpdateParams()) { mf::LogError("pma::Track3D") << "Track empty."; isOK = false; }
1953 
1954  depth--;
1955 
1956  return isOK;
1957 }
1958 
1959 double pma::Track3D::TuneSinglePass(bool skipFirst)
1960 {
1961  pma::Node3D* vtx = fNodes.front();
1962  pma::Segment3D* segThis = 0;
1963  pma::Segment3D* seg = 0;
1964 
1965  if (skipFirst)
1966  {
1967  segThis = NextSegment(vtx);
1968  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
1969  }
1970 
1971  double g = 0.0;
1972  while (vtx)
1973  {
1974  vtx->Optimize(fPenaltyValue, fEndSegWeight);
1975  segThis = NextSegment(vtx);
1976 
1977  for (size_t i = 0; i < vtx->NextCount(); i++)
1978  {
1979  seg = static_cast< pma::Segment3D* >(vtx->Next(i));
1980  if (seg != segThis) g += seg->Parent()->TuneSinglePass(true);
1981  }
1982 
1983  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
1984  else break;
1985  }
1986 
1987  return g + GetObjFunction();
1988 }
1989 
1990 pma::Track3D* pma::Track3D::GetNearestTrkInTree(
1991  const TVector3& p3d_cm, double& dist, bool skipFirst)
1992 {
1993  pma::Node3D* vtx = fNodes.front();
1994  pma::Segment3D* segThis = 0;
1995  pma::Segment3D* seg = 0;
1996 
1997  if (skipFirst)
1998  {
1999  segThis = NextSegment(vtx);
2000  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2001  }
2002 
2003  pma::Track3D* result = this;
2004  dist = sqrt( Dist2(p3d_cm) );
2005 
2006  pma::Track3D* candidate = 0;
2007  while (vtx)
2008  {
2009  segThis = NextSegment(vtx);
2010 
2011  double d;
2012  for (size_t i = 0; i < vtx->NextCount(); i++)
2013  {
2014  seg = static_cast< pma::Segment3D* >(vtx->Next(i));
2015  if (seg != segThis)
2016  {
2017  candidate = seg->Parent()->GetNearestTrkInTree(p3d_cm, d, true);
2018  if (d < dist) { dist = d; result = candidate; }
2019  }
2020  }
2021 
2022  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2023  else break;
2024  }
2025 
2026  return result;
2027 }
2028 
2029 pma::Track3D* pma::Track3D::GetNearestTrkInTree(
2030  const TVector2& p2d_cm, unsigned view, unsigned int tpc, unsigned int cryo,
2031  double& dist, bool skipFirst)
2032 {
2033  pma::Node3D* vtx = fNodes.front();
2034  pma::Segment3D* segThis = 0;
2035  pma::Segment3D* seg = 0;
2036 
2037  if (skipFirst)
2038  {
2039  segThis = NextSegment(vtx);
2040  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2041  }
2042 
2043  pma::Track3D* result = this;
2044  dist = Dist2(p2d_cm, view, tpc, cryo);
2045 
2046  pma::Track3D* candidate = 0;
2047  while (vtx)
2048  {
2049  segThis = NextSegment(vtx);
2050 
2051  double d;
2052  for (size_t i = 0; i < vtx->NextCount(); i++)
2053  {
2054  seg = static_cast< pma::Segment3D* >(vtx->Next(i));
2055  if (seg != segThis)
2056  {
2057  candidate = seg->Parent()->GetNearestTrkInTree(p2d_cm, view, tpc, cryo, d, true);
2058  if (d < dist) { dist = d; result = candidate; }
2059  }
2060  }
2061 
2062  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2063  else break;
2064  }
2065 
2066  dist = sqrt(dist);
2067  return result;
2068 }
2069 
2070 void pma::Track3D::ReassignHitsInTree(pma::Track3D* trkRoot)
2071 {
2072  bool skipFirst;
2073  if (trkRoot) skipFirst = true;
2074  else { trkRoot = this; skipFirst = false; }
2075 
2076  pma::Node3D* vtx = fNodes.front();
2077  pma::Segment3D* segThis = 0;
2078  pma::Segment3D* seg = 0;
2079 
2080  if (skipFirst)
2081  {
2082  segThis = NextSegment(vtx);
2083  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2084  }
2085 
2086  while (vtx)
2087  {
2088  segThis = NextSegment(vtx);
2089 
2090  for (size_t i = 0; i < vtx->NextCount(); i++)
2091  {
2092  seg = static_cast< pma::Segment3D* >(vtx->Next(i));
2093  if (seg != segThis) seg->Parent()->ReassignHitsInTree(trkRoot);
2094  }
2095 
2096  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2097  else break;
2098  }
2099 
2100  double d0, dmin;
2101  pma::Hit3D* hit = 0;
2102  pma::Track3D* nearestTrk = 0;
2103  size_t i = 0;
2104  while (HasTwoViews(2) && (i < size()))
2105  {
2106  hit = fHits[i];
2107  d0 = hit->GetDistToProj();
2108 
2109  nearestTrk = GetNearestTrkInTree(hit->Point2D(), hit->View2D(), hit->TPC(), hit->Cryo(), dmin);
2110  if ((nearestTrk != this) && (dmin < 0.5 * d0))
2111  {
2112  nearestTrk->push_back(release_at(i));
2113  mf::LogVerbatim("pma::Track3D") << "*** hit moved to another track ***";
2114  }
2115  else i++;
2116  }
2117  if (!size())
2118  {
2119  mf::LogError("pma::Track3D") << "ALL hits moved to other tracks.";
2120  }
2121 }
2122 
2123 void pma::Track3D::MakeProjectionInTree(bool skipFirst)
2124 {
2125  pma::Node3D* vtx = fNodes.front();
2126  pma::Segment3D* segThis = 0;
2127  pma::Segment3D* seg = 0;
2128 
2129  if (skipFirst)
2130  {
2131  segThis = NextSegment(vtx);
2132  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2133  }
2134 
2135  while (vtx)
2136  {
2137  segThis = NextSegment(vtx);
2138 
2139  for (size_t i = 0; i < vtx->NextCount(); i++)
2140  {
2141  seg = static_cast< pma::Segment3D* >(vtx->Next(i));
2142  if (seg != segThis) seg->Parent()->MakeProjectionInTree(true);
2143  }
2144 
2145  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2146  else break;
2147  }
2148 
2149  MakeProjection();
2150 }
2151 
2152 void pma::Track3D::SortHitsInTree(bool skipFirst)
2153 {
2154  pma::Node3D* vtx = fNodes.front();
2155  pma::Segment3D* segThis = 0;
2156  pma::Segment3D* seg = 0;
2157 
2158  if (skipFirst)
2159  {
2160  segThis = NextSegment(vtx);
2161  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2162  }
2163 
2164  while (vtx)
2165  {
2166  segThis = NextSegment(vtx);
2167 
2168  for (size_t i = 0; i < vtx->NextCount(); i++)
2169  {
2170  seg = static_cast< pma::Segment3D* >(vtx->Next(i));
2171  if (seg != segThis) seg->Parent()->SortHitsInTree(true);
2172  }
2173 
2174  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2175  else break;
2176  }
2177 
2178  SortHits();
2179 }
2180 
2181 double pma::Track3D::GetObjFnInTree(bool skipFirst)
2182 {
2183  pma::Node3D* vtx = fNodes.front();
2184  pma::Segment3D* segThis = 0;
2185  pma::Segment3D* seg = 0;
2186 
2187  if (skipFirst)
2188  {
2189  segThis = NextSegment(vtx);
2190  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2191  }
2192 
2193  double g = 0.0;
2194  while (vtx)
2195  {
2196  segThis = NextSegment(vtx);
2197 
2198  for (size_t i = 0; i < vtx->NextCount(); i++)
2199  {
2200  seg = static_cast< pma::Segment3D* >(vtx->Next(i));
2201  if (seg != segThis) g += seg->Parent()->GetObjFnInTree(true);
2202  }
2203 
2204  if (segThis) vtx = static_cast< pma::Node3D* >(segThis->Next());
2205  else break;
2206  }
2207 
2208  return g + GetObjFunction();
2209 }
2210 
2211 double pma::Track3D::TuneFullTree(double eps, double gmax)
2212 {
2213  if (!UpdateParamsInTree())
2214  {
2215  mf::LogError("pma::Track3D") << "TuneFullTree failed.";
2216  return -2; // negetive to tag destroyed tree
2217  }
2218 
2219  double g0 = GetObjFnInTree(), g1 = 0.0;
2220  if (!std::isfinite(g0))
2221  {
2222  mf::LogError("pma::Track3D") << "Infinifte value of g.";
2223  //std::cout << "Infinifte value of g." << std::endl;
2224  return -2; // negetive to tag destroyed tree
2225  }
2226  if (g0 > gmax)
2227  {
2228  mf::LogWarning("pma::Track3D") << "Too high value of g: " << g0;
2229  //std::cout << "Too high value of g: " << g0 << std::endl;
2230  return -1; // negetive to tag bad tree
2231  }
2232  if (g0 == 0.0) return g0;
2233 
2234  mf::LogVerbatim("pma::Track3D") << "Tune tree, g = " << g0;
2235  //std::cout << "Tune tree, g = " << g0 << std::endl;
2236  unsigned int stepIter = 0;
2237  do
2238  {
2239  double gstep = 1.0;
2240  unsigned int iter = 0;
2241  while ((gstep > eps) && (iter < 60))
2242  {
2243  g1 = g0;
2244  g0 = TuneSinglePass();
2245 
2246  MakeProjectionInTree();
2247  //ReassignHitsInTree();
2248 
2249  if (!UpdateParamsInTree()) { g0 = -2; break; } // negetive to tag destroyed tree
2250 
2251  if (g0 == 0.0F) break;
2252 
2253  gstep = fabs(g0 - g1) / g0;
2254  iter++;
2255  }
2256 
2257  stepIter++;
2258 
2259  } while (stepIter < 5);
2260 
2261  //ReassignHitsInTree();
2262  //TuneSinglePass();
2263  MakeProjectionInTree();
2264  SortHitsInTree();
2265 
2266  if (g0 >= 0) { mf::LogVerbatim("pma::Track3D") << " done, g = " << g0; }
2267  else { mf::LogError("pma::Track3D") << "TuneFullTree failed."; }
2268  return g0;
2269 }
2270 
2271 void pma::Track3D::ApplyDriftShiftInTree(double dx, bool skipFirst)
2272 {
2273  pma::Node3D* node = fNodes.front();
2274  pma::Segment3D* segThis = 0;
2275  pma::Segment3D* seg = 0;
2276 
2277  if (skipFirst)
2278  {
2279  segThis = NextSegment(node);
2280  if (segThis) node = static_cast< pma::Node3D* >(segThis->Next());
2281  }
2282 
2283  while (node)
2284  {
2285  node->ApplyDriftShift(dx);
2286 
2287  segThis = NextSegment(node);
2288  for (size_t i = 0; i < node->NextCount(); i++)
2289  {
2290  seg = static_cast< pma::Segment3D* >(node->Next(i));
2291  if (seg != segThis) seg->Parent()->ApplyDriftShiftInTree(dx, true);
2292  }
2293 
2294  if (segThis) node = static_cast< pma::Node3D* >(segThis->Next());
2295  else break;
2296  }
2297 
2298  for (auto h : fHits) { h->fPoint3D[0] += dx; }
2299 
2300  for (auto p : fAssignedPoints) { (*p)[0] += dx; }
2301 
2302  // For T0 we need to make sure we use the total shift, not just this current one in case of multiple shifts
2303  double newdx = fNodes.front()->GetDriftShift();
2304 
2305  // Now convert this newdx into T0 and store in fT0
2306  SetT0FromDx(newdx);
2307 }
2308 
2309 void pma::Track3D::SetT0FromDx(double dx){
2310 
2311  // Get all of the services we need
2312  auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
2313  auto const* detclock = lar::providerFrom<detinfo::DetectorClocksService>();
2314  auto const* geom = lar::providerFrom<geo::Geometry>();
2315  const geo::TPCGeo &tpcGeo = geom->TPC(fNodes.front()->TPC(), fNodes.front()->Cryo());
2316 
2317  // GetXTicksCoefficient has a sign that we don't care about. We need to decide
2318  // the sign for ourselves based on the drift direction of the TPC
2319  double correctedSign = 0;
2320  if(tpcGeo.DetectDriftDirection() > 0){
2321  if(dx > 0){
2322  correctedSign = 1.0;
2323  }
2324  else{
2325  correctedSign = -1.0;
2326  }
2327  }
2328  else{
2329  if(dx > 0){
2330  correctedSign = -1.0;
2331  }
2332  else{
2333  correctedSign = 1.0;
2334  }
2335  }
2336 
2337  // The magnitude of x in ticks is fine
2338  double dxInTicks = fabs(dx / detprop->GetXTicksCoefficient());
2339  // Now correct the sign
2340  dxInTicks = dxInTicks * correctedSign;
2341  // At this stage, we have xInTicks relative to the trigger time
2342  double dxInTime = dxInTicks * detclock->TPCClock().TickPeriod();
2343  // Reconstructed times are relative to the trigger (t=0), so this is our T0
2344  fT0 = dxInTime;
2345 
2346  mf::LogDebug("pma::Track3D") << dx << ", " << dxInTicks << ", " << correctedSign
2347  << ", " << fT0 << ", " << tpcGeo.DetectDriftDirection() << " :: "
2348  << detclock->TriggerTime() << ", " << detclock->TriggerOffsetTPC() << std::endl;
2349 
2350  // Leigh test
2351 // mf::LogDebug("pma::Track3D") << " Check x position for T = Trigger Offset " << detprop->ConvertTicksToX(detprop->TriggerOffset(),2,tpcGeo.ID().TPC,tpcGeo.ID().Cryostat) << std::endl;
2352 // float apaTime = detprop->GetXTicksOffset(2,tpcGeo.ID().TPC,tpcGeo.ID().Cryostat) - detprop->TriggerOffset();
2353 // mf::LogDebug("pma::Track3D") << " Check x position for T = Trigger Offset + DriftTime " << detprop->ConvertTicksToX(detprop->TriggerOffset() + apaTime,2,tpcGeo.ID().TPC,tpcGeo.ID().Cryostat) << std::endl;
2354 
2355  // Mark this track as having a measured T0
2356  fT0Flag = true;
2357 }
2358 
2359 void pma::Track3D::DeleteSegments(void)
2360 {
2361  for (size_t i = 0; i < fSegments.size(); i++) delete fSegments[i];
2362  fSegments.clear();
2363 }
2364 
2365 void pma::Track3D::RebuildSegments(void)
2366 {
2367  DeleteSegments();
2368 
2369  fSegments.reserve(fNodes.size() - 1);
2370  for (size_t i = 1; i < fNodes.size(); i++)
2371  {
2372  fSegments.push_back(new pma::Segment3D(this, fNodes[i - 1], fNodes[i]));
2373  }
2374 }
2375 
2376 void pma::Track3D::CleanupTails(void)
2377 {
2378  unsigned int nhits = 0;
2379  while (!nhits && (fNodes.size() > 2) && !fNodes.front()->IsBranching())
2380  {
2381  pma::Node3D* vtx = fNodes.front();
2382  nhits = vtx->NHits();
2383 
2384  pma::Segment3D* seg = NextSegment(vtx);
2385  nhits += seg->NHits();
2386 
2387  if (!nhits)
2388  {
2389  if (vtx->NextCount() < 2) delete vtx;
2390  fNodes.erase(fNodes.begin());
2391 
2392  delete seg;
2393  fSegments.erase(fSegments.begin());
2394 
2395  MakeProjection();
2396  }
2397  }
2398 
2399  nhits = 0;
2400  while (!nhits && (fNodes.size() > 2) && !fNodes.back()->IsBranching())
2401  {
2402  pma::Node3D* vtx = fNodes.back();
2403  nhits = vtx->NHits();
2404 
2405  pma::Segment3D* seg = PrevSegment(vtx);
2406  nhits += seg->NHits();
2407 
2408  if (!nhits)
2409  {
2410  if (vtx->NextCount() < 1) delete fNodes.back();
2411  fNodes.pop_back();
2412 
2413  delete seg;
2414  fSegments.pop_back();
2415 
2416  MakeProjection();
2417  }
2418  }
2419 }
2420 
2421 
2422 bool pma::Track3D::ShiftEndsToHits(void)
2423 {
2424  pma::Element3D* el;
2425  pma::Node3D* vtx;
2426 
2427  if (!(fNodes.front()->Prev()))
2428  {
2429  el = GetNearestElement(front()->Point3D());
2430  vtx = dynamic_cast< pma::Node3D* >(el);
2431  if (vtx)
2432  {
2433  if (vtx == fNodes.front()) fNodes.front()->SetPoint3D(front()->Point3D());
2434  else
2435  {
2436  mf::LogWarning("pma::Track3D") << "First hit is projected to inner node.";
2437  return false;
2438  }
2439  }
2440  else
2441  {
2442  pma::Segment3D* seg = dynamic_cast< pma::Segment3D* >(el);
2443  if (seg)
2444  {
2445  if (seg->Prev() == fNodes.front())
2446  {
2447  double l0 = seg->Length();
2448  fNodes.front()->SetPoint3D(front()->Point3D());
2449  if ((seg->Length() < 0.2 * l0) && (fNodes.size() > 2))
2450  {
2451  pma::Node3D* toRemove = fNodes[1];
2452  if (toRemove->NextCount() == 1)
2453  {
2454  mf::LogWarning("pma::Track3D") << "ShiftEndsToHits(): Short start segment, node removed.";
2455  //std::cout << "ShiftEndsToHits(): Short start segment, node removed." << std::endl;
2456  fNodes.erase(fNodes.begin() + 1);
2457 
2458  RebuildSegments();
2459  MakeProjection();
2460 
2461  delete toRemove;
2462  }
2463  else
2464  {
2465  mf::LogWarning("pma::Track3D") << "Branching node, not removed.";
2466  return false;
2467  }
2468  }
2469  }
2470  else
2471  {
2472  mf::LogWarning("pma::Track3D") << "First hit is projected to inner segment.";
2473  return false;
2474  }
2475  }
2476  }
2477  }
2478 
2479  if (!(fNodes.back()->NextCount()))
2480  {
2481  el = GetNearestElement(back()->Point3D());
2482  vtx = dynamic_cast< pma::Node3D* >(el);
2483  if (vtx)
2484  {
2485  if (vtx == fNodes.back()) fNodes.back()->SetPoint3D(back()->Point3D());
2486  else
2487  {
2488  mf::LogWarning("pma::Track3D") << "First hit is projected to inner node.";
2489  return false;
2490  }
2491  }
2492  else
2493  {
2494  pma::Segment3D* seg = dynamic_cast< pma::Segment3D* >(el);
2495  if (seg)
2496  {
2497  if (seg->Next() == fNodes.back())
2498  {
2499  double l0 = seg->Length();
2500  fNodes.back()->SetPoint3D(back()->Point3D());
2501  if ((seg->Length() < 0.2 * l0) && (fNodes.size() > 2))
2502  {
2503  size_t idx = fNodes.size() - 2;
2504  pma::Node3D* toRemove = fNodes[idx];
2505  if (toRemove->NextCount() == 1)
2506  {
2507  mf::LogWarning("pma::Track3D") << "ShiftEndsToHits(): Short end segment, node removed.";
2508  //std::cout << "ShiftEndsToHits(): Short end segment, node removed." << std::endl;
2509  fNodes.erase(fNodes.begin() + idx);
2510 
2511  RebuildSegments();
2512  MakeProjection();
2513 
2514  delete toRemove;
2515  }
2516  else
2517  {
2518  mf::LogWarning("pma::Track3D") << "Branching node, not removed.";
2519  return false;
2520  }
2521  }
2522  }
2523  else
2524  {
2525  mf::LogWarning("pma::Track3D") << "First hit is projected to inner segment.";
2526  return false;
2527  }
2528  }
2529  }
2530  }
2531 
2532  return true;
2533 }
2534 
2535 double pma::Track3D::Dist2(const TVector2& p2d, unsigned int view, unsigned int tpc, unsigned int cryo) const
2536 {
2537  double dist, min_dist = 1.0e12;
2538 
2539  int t = (int)tpc, c = (int)cryo;
2540  auto n0 = fNodes.front();
2541  if ((n0->TPC() == t) && (n0->Cryo() == c))
2542  {
2543  dist = n0->GetDistance2To(p2d, view);
2544  if (dist < min_dist) min_dist = dist;
2545  }
2546  auto n1 = fNodes.back();
2547  if ((n1->TPC() == t) && (n1->Cryo() == c))
2548  {
2549  dist = n1->GetDistance2To(p2d, view);
2550  if (dist < min_dist) min_dist = dist;
2551  }
2552 
2553  for (auto s : fSegments)
2554  if ((s->TPC() == t) && (s->Cryo() == c))
2555  {
2556  dist = s->GetDistance2To(p2d, view);
2557  if (dist < min_dist) min_dist = dist;
2558  }
2559  return min_dist;
2560 }
2561 
2562 double pma::Track3D::Dist2(const TVector3& p3d) const
2563 {
2564  double dist, min_dist = fSegments.front()->GetDistance2To(p3d);
2565  for (size_t i = 1; i < fSegments.size(); i++)
2566  {
2567  dist = fSegments[i]->GetDistance2To(p3d);
2568  if (dist < min_dist) min_dist = dist;
2569  }
2570  return min_dist;
2571 }
2572 
2573 pma::Element3D* pma::Track3D::GetNearestElement(
2574  const TVector2& p2d, unsigned int view, int tpc,
2575  bool skipFrontVtx, bool skipBackVtx) const
2576 {
2577  if (fSegments.front()->TPC() < 0) skipFrontVtx = false;
2578  if (fSegments.back()->TPC() < 0) skipBackVtx = false;
2579 
2580  if (skipFrontVtx && skipBackVtx && (fSegments.size() == 1))
2581  return fSegments.front(); // no need for searching...
2582 
2583  size_t v0 = 0, v1 = fNodes.size();
2584  if (skipFrontVtx) v0 = 1;
2585  if (skipBackVtx) --v1;
2586 
2587  pma::Element3D* pe_min = 0;
2588  double dist, min_dist = 1.0e9;
2589  for (size_t i = v0; i < v1; i++)
2590  //if ((tpc == -1) || (fNodes[i]->TPC() == tpc))
2591  if (fNodes[i]->TPC() == tpc)
2592  {
2593  dist = fNodes[i]->GetDistance2To(p2d, view);
2594  if (dist < min_dist)
2595  {
2596  min_dist = dist; pe_min = fNodes[i];
2597  }
2598  }
2599  for (size_t i = 0; i < fSegments.size(); i++)
2600  //if ((tpc == -1) || (fSegments[i]->TPC() == tpc))
2601  if (fSegments[i]->TPC() == tpc)
2602  {
2603  if (fSegments[i]->TPC() < 0) continue; // segment between TPC's
2604 
2605  dist = fSegments[i]->GetDistance2To(p2d, view);
2606  if (dist < min_dist)
2607  {
2608  min_dist = dist; pe_min = fSegments[i];
2609  }
2610  }
2611  if (!pe_min) throw cet::exception("pma::Track3D") << "Nearest element not found." << std::endl;
2612  return pe_min;
2613 }
2614 
2615 pma::Element3D* pma::Track3D::GetNearestElement(const TVector3& p3d) const
2616 {
2617  pma::Element3D* pe_min = fNodes.front();
2618  double dist, min_dist = pe_min->GetDistance2To(p3d);
2619  for (size_t i = 1; i < fNodes.size(); i++)
2620  {
2621  dist = fNodes[i]->GetDistance2To(p3d);
2622  if (dist < min_dist)
2623  {
2624  min_dist = dist; pe_min = fNodes[i];
2625  }
2626  }
2627  for (size_t i = 0; i < fSegments.size(); i++)
2628  {
2629  dist = fSegments[i]->GetDistance2To(p3d);
2630  if (dist < min_dist)
2631  {
2632  min_dist = dist; pe_min = fSegments[i];
2633  }
2634  }
2635  return pe_min;
2636 }
2637 
2638 bool pma::Track3D::GetUnconstrainedProj3D(art::Ptr<recob::Hit> hit, TVector3& p3d, double& dist2) const
2639 {
2640  TVector2 p2d = pma::WireDriftToCm(
2641  hit->WireID().Wire, hit->PeakTime(),
2642  hit->WireID().Plane, hit->WireID().TPC, hit->WireID().Cryostat);
2643 
2644  pma::Segment3D* seg = 0;
2645  double d2, min_d2 = 1.0e100;
2646  int tpc = (int)hit->WireID().TPC;
2647  int view = hit->WireID().Plane;
2648  for (size_t i = 0; i < fSegments.size(); i++)
2649  {
2650  if (fSegments[i]->TPC() != tpc) continue;
2651 
2652  d2 = fSegments[i]->GetDistance2To(p2d, view);
2653  if (d2 < min_d2)
2654  {
2655  min_d2 = d2; seg = fSegments[i];
2656  }
2657  }
2658  if (seg)
2659  {
2660  p3d = seg->GetUnconstrainedProj3D(p2d, view);
2661  dist2 = min_d2;
2662 
2663  pma::Node3D* prev = static_cast< pma::Node3D* >(seg->Prev());
2664  return prev->SameTPC(p3d); // 3D can be beyond the segment endpoints => in other TPC
2665  }
2666 
2667  return false;
2668 }
2669 
2670 void pma::Track3D::SortHits(void)
2671 {
2672  std::vector< pma::Hit3D* > hits_tmp;
2673  hits_tmp.reserve(size());
2674 
2675  pma::Node3D* vtx = fNodes.front();
2676  pma::Segment3D* seg = NextSegment(vtx);
2677  while (vtx)
2678  {
2679  vtx->SortHits();
2680  for (size_t i = 0; i < vtx->NHits(); i++)
2681  {
2682  pma::Hit3D* h3d = &(vtx->Hit(i));
2683  if (h3d->fParent == this) hits_tmp.push_back(h3d);
2684  }
2685 
2686  if (seg)
2687  {
2688  seg->SortHits();
2689  for (size_t i = 0; i < seg->NHits(); i++)
2690  {
2691  pma::Hit3D* h3d = &(seg->Hit(i));
2692  if (h3d->fParent == this) hits_tmp.push_back(h3d);
2693  }
2694 
2695  vtx = static_cast< pma::Node3D* >(seg->Next());
2696  seg = NextSegment(vtx);
2697  }
2698  else break;
2699  }
2700 
2701  if (size() == hits_tmp.size())
2702  {
2703  for (size_t i = 0; i < size(); i++)
2704  {
2705  fHits[i] = hits_tmp[i];
2706  }
2707  }
2708  else mf::LogError("pma::Track3D") << "Hit sorting problem.";
2709 }
2710 
2711 unsigned int pma::Track3D::DisableSingleViewEnds(void)
2712 {
2713  SortHits();
2714 
2715  unsigned int nDisabled = 0;
2716 
2717  int hasHits[4];
2718 
2719  pma::Hit3D* nextHit = 0;
2720  int hitIndex = -1;
2721 
2722  bool rebuild = false, stop = false;
2723  int nViews = 0;
2724  hasHits[0] = hasHits[1] = hasHits[2] = hasHits[3] = 0;
2725  do
2726  {
2727  pma::Node3D* vtx = fNodes.front();
2728  pma::Segment3D* seg = NextSegment(vtx);
2729  if (!seg) break;
2730 
2731  if (vtx->NPoints() + seg->NPoints() > 0) hasHits[3] = 1;
2732 
2733  for (size_t i = 0; i < vtx->NHits(); i++)
2734  {
2735  hitIndex = index_of(&(vtx->Hit(i)));
2736  if ((hitIndex >= 0) && (hitIndex + 1 < (int)size())) nextHit = fHits[hitIndex + 1];
2737  else nextHit = 0;
2738 
2739  if (vtx->Hit(i).IsEnabled()) hasHits[vtx->Hit(i).View2D()] = 1;
2740  if (nextHit && nextHit->IsEnabled()) hasHits[nextHit->View2D()] = 1;
2741  nViews = hasHits[0] + hasHits[1] + hasHits[2] + hasHits[3];
2742  if (nViews < 2)
2743  {
2744  if (vtx->Hit(i).IsEnabled())
2745  {
2746  vtx->Hit(i).SetEnabled(false);
2747  nDisabled++;
2748  }
2749  }
2750  }
2751  for (size_t i = 0; i < seg->NHits(); i++)
2752  {
2753  hitIndex = index_of(&(seg->Hit(i)));
2754  if ((hitIndex >= 0) && (hitIndex + 1 < (int)size())) nextHit = fHits[hitIndex + 1];
2755  else nextHit = 0;
2756 
2757  if (seg->Hit(i).IsEnabled()) hasHits[seg->Hit(i).View2D()] = 1;
2758  if (nextHit && nextHit->IsEnabled()) hasHits[nextHit->View2D()] = 1;
2759  nViews = hasHits[0] + hasHits[1] + hasHits[2] + hasHits[3];
2760  if (nViews < 2)
2761  {
2762  if (seg->Hit(i).IsEnabled())
2763  {
2764  seg->Hit(i).SetEnabled(false);
2765  nDisabled++;
2766  }
2767  }
2768  }
2769 
2770  if (fNodes.size() < 3) break;
2771 
2772  nViews = hasHits[1] + hasHits[2] + hasHits[3];
2773  if (hasHits[0] || (nViews > 1)) stop = true;
2774  else if (!fNodes.front()->IsBranching())
2775  {
2776  pma::Node3D* vtx_front = fNodes.front();
2777  fNodes.erase(fNodes.begin());
2778  delete vtx_front;
2779  rebuild = true;
2780  }
2781 
2782  } while (!stop);
2783 
2784  stop = false;
2785  nViews = 0;
2786  hasHits[0] = hasHits[1] = hasHits[2] = hasHits[3] = 0;
2787  do
2788  {
2789  pma::Node3D* vtx = fNodes.back();
2790  pma::Segment3D* seg = PrevSegment(vtx);
2791  if (!seg) break;
2792 
2793  if (vtx->NPoints() || seg->NPoints()) hasHits[3] = 1;
2794 
2795  for (int i = vtx->NHits() - 1; i >= 0; i--)
2796  {
2797  hitIndex = index_of(&(vtx->Hit(i)));
2798  if ((hitIndex >= 0) && (hitIndex - 1 >= 0)) nextHit = fHits[hitIndex - 1];
2799  else nextHit = 0;
2800 
2801  if (vtx->Hit(i).IsEnabled()) hasHits[vtx->Hit(i).View2D()] = 1;
2802  if (nextHit && nextHit->IsEnabled()) hasHits[nextHit->View2D()] = 1;
2803  nViews = hasHits[0] + hasHits[1] + hasHits[2] + hasHits[3];
2804  if (nViews < 2)
2805  {
2806  if (vtx->Hit(i).IsEnabled())
2807  {
2808  vtx->Hit(i).SetEnabled(false);
2809  nDisabled++;
2810  }
2811  }
2812  }
2813  for (int i = seg->NHits() - 1; i >= 0; i--)
2814  {
2815  hitIndex = index_of(&(seg->Hit(i)));
2816  if ((hitIndex >= 0) && (hitIndex - 1 >= 0)) nextHit = fHits[hitIndex - 1];
2817  else nextHit = 0;
2818 
2819  if (seg->Hit(i).IsEnabled()) hasHits[seg->Hit(i).View2D()] = 1;
2820  if (nextHit && nextHit->IsEnabled()) hasHits[nextHit->View2D()] = 1;
2821  nViews = hasHits[0] + hasHits[1] + hasHits[2] + hasHits[3];
2822  if (nViews < 2)
2823  {
2824  if (seg->Hit(i).IsEnabled())
2825  {
2826  seg->Hit(i).SetEnabled(false);
2827  nDisabled++;
2828  }
2829  }
2830  }
2831 
2832  if (fNodes.size() < 3) break;
2833 
2834  nViews = hasHits[1] + hasHits[2] + hasHits[3];
2835  if (hasHits[0] || (nViews > 1)) stop = true;
2836  else if (!fNodes.back()->IsBranching())
2837  {
2838  pma::Node3D* vtx_back = fNodes.back();
2839  fNodes.pop_back();
2840  delete vtx_back;
2841  rebuild = true;
2842  }
2843 
2844  } while (!stop);
2845 
2846  if (rebuild)
2847  {
2848  RebuildSegments();
2849  MakeProjection();
2850  }
2851 
2852  return nDisabled;
2853 }
2854 
2855 bool pma::Track3D::SelectHits(float fraction)
2856 {
2857  if (fraction < 0.0F) fraction = 0.0F;
2858  if (fraction > 1.0F) fraction = 1.0F;
2859  if (fraction < 1.0F) std::sort(fHits.begin(), fHits.end(), pma::bTrajectory3DDistLess());
2860 
2861  size_t nHitsColl = (size_t)(fraction * NHits(geo::kZ));
2862  size_t nHitsInd2 = (size_t)(fraction * NHits(geo::kV));
2863  size_t nHitsInd1 = (size_t)(fraction * NHits(geo::kU));
2864  size_t coll = 0, ind2 = 0, ind1 = 0;
2865 
2866  bool b, changed = false;
2867  for (auto h : fHits)
2868  {
2869  b = h->IsEnabled();
2870  if (fraction < 1.0F)
2871  {
2872  h->SetEnabled(false);
2873  switch (h->View2D())
2874  {
2875  case geo::kZ: if (coll++ < nHitsColl) h->SetEnabled(true); break;
2876  case geo::kV: if (ind2++ < nHitsInd2) h->SetEnabled(true); break;
2877  case geo::kU: if (ind1++ < nHitsInd1) h->SetEnabled(true); break;
2878  }
2879  }
2880  else h->SetEnabled(true);
2881 
2882  changed |= (b != h->IsEnabled());
2883  }
2884 
2885  if (fraction < 1.0F)
2886  {
2887  FirstElement()->SelectAllHits();
2888  LastElement()->SelectAllHits();
2889  }
2890  return changed;
2891 }
2892 
2893 bool pma::Track3D::SelectRndHits(size_t segmax, size_t vtxmax)
2894 {
2895  bool changed = false;
2896  for (auto n : fNodes) changed |= n->SelectRndHits(vtxmax);
2897  for (auto s : fSegments) changed |= s->SelectRndHits(segmax);
2898  return changed;
2899 }
2900 
2901 bool pma::Track3D::SelectAllHits(void)
2902 {
2903  bool changed = false;
2904  for (auto h : fHits)
2905  {
2906  changed |= !(h->IsEnabled());
2907  h->SetEnabled(true);
2908  }
2909  return changed;
2910 }
2911 
2912 void pma::Track3D::MakeProjection(void)
2913 {
2914  for (auto n : fNodes) n->ClearAssigned(this);
2915  for (auto s : fSegments) s->ClearAssigned(this);
2916 
2917  pma::Element3D* pe = 0;
2918 
2919  bool skipFrontVtx = false, skipBackVtx = false;
2920 
2921  // skip outermost vertices if not branching
2922  if (!(fNodes.front()->IsFrozen()) &&
2923  !(fNodes.front()->Prev()) &&
2924  (fNodes.front()->NextCount() == 1))
2925  {
2926  skipFrontVtx = true;
2927  }
2928  if (!(fNodes.front()->IsFrozen()) &&
2929  (fNodes.back()->NextCount() == 0))
2930  {
2931  skipBackVtx = true;
2932  }
2933 
2934  for (auto h : fHits) // assign hits to nodes/segments
2935  {
2936  pe = GetNearestElement(h->Point2D(), h->View2D(), h->TPC(), skipFrontVtx, skipBackVtx);
2937  pe->AddHit(h);
2938  }
2939 
2940  for (auto p : fAssignedPoints) // assign ref points to nodes/segments
2941  {
2942  pe = GetNearestElement(*p);
2943  pe->AddPoint(p);
2944  }
2945 
2946  for (auto n : fNodes) n->UpdateHitParams();
2947  for (auto s : fSegments) s->UpdateHitParams();
2948 }
2949 
2950 void pma::Track3D::MakeFastProjection(void)
2951 {
2952  std::vector< std::pair< pma::Hit3D*, pma::Element3D* > > assignments;
2953  assignments.reserve(fHits.size());
2954 
2955  for (auto hi : fHits)
2956  {
2957  pma::Element3D* pe = 0;
2958 
2959  for (auto s : fSegments)
2960  {
2961  for (size_t j = 0; j < s->NHits(); ++j)
2962  if (hi == s->Hits()[j]) // look at next/prev vtx,seg,vtx
2963  {
2964  pe = s;
2965  double d2, min_d2 = s->GetDistance2To(hi->Point2D(), hi->View2D());
2966  int tpc = hi->TPC();
2967 
2968  pma::Node3D* nnext = static_cast< pma::Node3D* >(s->Next());
2969  if (nnext->TPC() == tpc)
2970  {
2971  d2 = nnext->GetDistance2To(hi->Point2D(), hi->View2D());
2972  if (d2 < min_d2) { min_d2 = d2; pe = nnext; }
2973 
2974  pma::Segment3D* snext = NextSegment(nnext);
2975  if (snext && (snext->TPC() == tpc))
2976  {
2977  d2 = snext->GetDistance2To(hi->Point2D(), hi->View2D());
2978  if (d2 < min_d2) { min_d2 = d2; pe = snext; }
2979 
2980  nnext = static_cast< pma::Node3D* >(snext->Next());
2981  if (nnext->TPC() == tpc)
2982  {
2983  d2 = nnext->GetDistance2To(hi->Point2D(), hi->View2D());
2984  if (d2 < min_d2) { min_d2 = d2; pe = nnext; }
2985  }
2986  }
2987  }
2988 
2989  pma::Node3D* nprev = static_cast< pma::Node3D* >(s->Prev());
2990  if (nprev->TPC() == tpc)
2991  {
2992  d2 = nprev->GetDistance2To(hi->Point2D(), hi->View2D());
2993  if (d2 < min_d2) { min_d2 = d2; pe = nprev; }
2994 
2995  pma::Segment3D* sprev = PrevSegment(nprev);
2996  if (sprev && (sprev->TPC() == tpc))
2997  {
2998  d2 = sprev->GetDistance2To(hi->Point2D(), hi->View2D());
2999  if (d2 < min_d2) { min_d2 = d2; pe = sprev; }
3000 
3001  nprev = static_cast< pma::Node3D* >(sprev->Prev());
3002  if (nprev->TPC() == tpc)
3003  {
3004  d2 = nprev->GetDistance2To(hi->Point2D(), hi->View2D());
3005  if (d2 < min_d2) { min_d2 = d2; pe = nprev; }
3006  }
3007  }
3008  }
3009 
3010  s->RemoveHitAt(j); break;
3011  }
3012  if (pe) break;
3013  }
3014 
3015  if (!pe) for (auto n : fNodes)
3016  {
3017  for (size_t j = 0; j < n->NHits(); ++j)
3018  if (hi == n->Hits()[j]) // look at next/prev seg,vtx,seg
3019  {
3020  pe = n;
3021  double d2, min_d2 = n->GetDistance2To(hi->Point2D(), hi->View2D());
3022  int tpc = hi->TPC();
3023 
3024  pma::Segment3D* snext = NextSegment(n);
3025  if (snext && (snext->TPC() == tpc))
3026  {
3027  d2 = snext->GetDistance2To(hi->Point2D(), hi->View2D());
3028  if (d2 < min_d2) { min_d2 = d2; pe = snext; }
3029 
3030  pma::Node3D* nnext = static_cast< pma::Node3D* >(snext->Next());
3031  if (nnext->TPC() == tpc)
3032  {
3033  d2 = nnext->GetDistance2To(hi->Point2D(), hi->View2D());
3034  if (d2 < min_d2) { min_d2 = d2; pe = nnext; }
3035 
3036  snext = NextSegment(nnext);
3037  if (snext && (snext->TPC() == tpc))
3038  {
3039  d2 = snext->GetDistance2To(hi->Point2D(), hi->View2D());
3040  if (d2 < min_d2) { min_d2 = d2; pe = snext; }
3041  }
3042  }
3043  }
3044 
3045  pma::Segment3D* sprev = PrevSegment(n);
3046  if (sprev && (sprev->TPC() == tpc))
3047  {
3048  d2 = sprev->GetDistance2To(hi->Point2D(), hi->View2D());
3049  if (d2 < min_d2) { min_d2 = d2; pe = sprev; }
3050 
3051  pma::Node3D* nprev = static_cast< pma::Node3D* >(sprev->Prev());
3052  if (nprev->TPC() == tpc)
3053  {
3054  d2 = nprev->GetDistance2To(hi->Point2D(), hi->View2D());
3055  if (d2 < min_d2) { min_d2 = d2; pe = nprev; }
3056 
3057  sprev = PrevSegment(nprev);
3058  if (sprev && (sprev->TPC() == tpc))
3059  {
3060  d2 = sprev->GetDistance2To(hi->Point2D(), hi->View2D());
3061  if (d2 < min_d2) { min_d2 = d2; pe = sprev; }
3062  }
3063  }
3064  }
3065 
3066  n->RemoveHitAt(j); break;
3067  }
3068  if (pe) break;
3069  }
3070 
3071  if (pe) assignments.emplace_back(hi, pe);
3072  else mf::LogWarning("pma::Track3D") << "Hit was not assigned to any element.";
3073  }
3074 
3075  for (auto const & a : assignments) a.second->AddHit(a.first);
3076 
3077  for (auto n : fNodes) n->UpdateHitParams();
3078  for (auto s : fSegments) s->UpdateHitParams();
3079 }
3080 
3081 void pma::Track3D::UpdateProjection(void)
3082 {
3083  for (auto n : fNodes) n->UpdateProjection();
3084  for (auto s : fSegments) s->UpdateProjection();
3085 }
3086 
3087 double pma::Track3D::AverageDist2(void) const
3088 {
3089  double sum = 0.0;
3090  unsigned int count = 0;
3091 
3092  for (auto n : fNodes)
3093  {
3094  sum += n->SumDist2();
3095  count += n->NEnabledHits();
3096  }
3097 
3098  for (auto s : fSegments)
3099  {
3100  sum += s->SumDist2();
3101  count += s->NEnabledHits();
3102  }
3103 
3104  if (count) { return sum / count; }
3105  else
3106  {
3107  mf::LogError("pma::Track3D") << "0 enabled hits in AverageDist2 calculation.";
3108  return 0;
3109  }
3110 }
3111 
3112 bool pma::Track3D::UpdateParams(void)
3113 {
3114  size_t n = size();
3115  if (n == 0)
3116  {
3117  fPenaltyValue = 1;
3118  fSegStopValue = 1;
3119  return false;
3120  }
3121 
3122  float nCubeRoot = pow((double)n, 1.0/3.0);
3123  float avgDist2Root = sqrt(AverageDist2());
3124  if (avgDist2Root > 0)
3125  {
3126  fPenaltyValue = fPenaltyFactor * pow((double)fSegments.size(), 1.8) * avgDist2Root / (fHitsRadius * nCubeRoot);
3127 
3128  fSegStopValue = (int)(fSegStopFactor * nCubeRoot * fHitsRadius / avgDist2Root);
3129  if (fSegStopValue < fMinSegStop) fSegStopValue = fMinSegStop;
3130  return true;
3131  }
3132  else
3133  {
3134  fPenaltyValue = 1;
3135  fSegStopValue = 1;
3136  return false;
3137  }
3138 }
3139 
3140 bool pma::Track3D::SwapVertices(size_t v0, size_t v1)
3141 {
3142  if (v0 == v1) return false;
3143 
3144  if (v0 > v1)
3145  {
3146  size_t vx = v0;
3147  v0 = v1; v1 = vx;
3148  }
3149 
3150  pma::Node3D* vtmp;
3151  if (v1 - v0 == 1)
3152  {
3153  pma::Segment3D* midSeg = NextSegment(fNodes[v0]);
3154  pma::Segment3D* prevSeg = PrevSegment(fNodes[v0]);
3155  pma::Segment3D* nextSeg = NextSegment(fNodes[v1]);
3156 
3157  fNodes[v1]->RemoveNext(nextSeg);
3158  midSeg->Disconnect();
3159 
3160  vtmp = fNodes[v0];
3161  fNodes[v0] = fNodes[v1];
3162  fNodes[v1] = vtmp;
3163 
3164  if (prevSeg) prevSeg->AddNext(fNodes[v0]);
3165  fNodes[v0]->AddNext(midSeg);
3166  midSeg->AddNext(fNodes[v1]);
3167  if (nextSeg) fNodes[v1]->AddNext(nextSeg);
3168 
3169  return false;
3170  }
3171  else
3172  {
3173  vtmp = fNodes[v0];
3174  fNodes[v0] = fNodes[v1];
3175  fNodes[v1] = vtmp;
3176  return true;
3177  }
3178 }
3179 
3180 bool pma::Track3D::CheckEndSegment(pma::Track3D::ETrackEnd endCode)
3181 {
3182  unsigned int v1, v2;
3183  switch (endCode)
3184  {
3185  case pma::Track3D::kBegin:
3186  if (fSegments.front()->IsFrozen()) return false;
3187  if (fNodes.front()->NextCount() > 1) return false;
3188  v1 = 0; v2 = 1; break;
3189  case pma::Track3D::kEnd:
3190  if (fSegments.back()->IsFrozen()) return false;
3191  if (fNodes.back()->NextCount() > 1) return false;
3192  v1 = fNodes.size() - 1;
3193  v2 = fNodes.size() - 2;
3194  break;
3195  default: return false;
3196  }
3197  if (fNodes[v1]->TPC() != fNodes[v2]->TPC()) return false;
3198 
3199  double g1, g0 = GetObjFunction();
3200 
3201  if (SwapVertices(v1, v2)) RebuildSegments();
3202  MakeProjection();
3203  g1 = GetObjFunction();
3204 
3205  if (g1 >= g0)
3206  {
3207  if (SwapVertices(v1, v2)) RebuildSegments();
3208  MakeProjection();
3209  return false;
3210  }
3211  else return true;
3212 }
3213 
3214 void pma::Track3D::UpdateHitsRadius(void)
3215 {
3216  std::vector< pma::Hit3D* > hitsColl, hitsInd1, hitsInd2;
3217  for (size_t i = 0; i < size(); i++)
3218  {
3219  pma::Hit3D* hit = (*this)[i];
3220  switch (hit->View2D())
3221  {
3222  case geo::kZ: hitsColl.push_back(hit); break;
3223  case geo::kV: hitsInd2.push_back(hit); break;
3224  case geo::kU: hitsInd1.push_back(hit); break;
3225  }
3226  }
3227  fHitsRadius = pma::GetHitsRadius2D(hitsColl, true);
3228  float r = pma::GetHitsRadius2D(hitsInd2, true);
3229  if (r > fHitsRadius) fHitsRadius = r;
3230  r = pma::GetHitsRadius2D(hitsInd1, true);
3231  if (r > fHitsRadius) fHitsRadius = r;
3232 }
3233 
x
Float_t x
Definition: compare.C:6
pma::Track3D::TestHits
unsigned int TestHits(const std::vector< art::Ptr< recob::Hit > > &hits, double dist=0.4) const
Count close 2D hits.
Definition: PmaTrack3D.cxx:811
pma::Track3D::SelectHits
bool SelectHits(float fraction=1.0F)
Definition: PmaTrack3D.cxx:2855
pma::Node3D::ApplyDriftShift
void ApplyDriftShift(double dx)
Definition: PmaNode3D.h:101
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:18
s
Float_t s
Definition: plot.C:23
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:219
pma::Track3D::fSegments
std::vector< pma::Segment3D * > fSegments
Definition: PmaTrack3D.h:331
pma::Hit3D::TPC
unsigned int TPC(void) const
Definition: PmaHit3D.h:57
pma::Track3D::LastElement
pma::Node3D * LastElement(void) const
Definition: PmaTrack3D.h:234
pma::Track3D::TPCs
std::vector< unsigned int > TPCs(void) const
Definition: PmaTrack3D.cxx:437
pma::Track3D::AutoFlip
void AutoFlip(pma::Track3D::EDirection dir, double thr=0.0, unsigned int n=0)
Definition: PmaTrack3D.cxx:630
pma::Track3D::CleanupTails
void CleanupTails(void)
Cut out tails with no hits assigned.
Definition: PmaTrack3D.cxx:2376
pma::Track3D::InitFromHits
bool InitFromHits(int tpc, int cryo, float initEndSegW=0.05F)
Definition: PmaTrack3D.cxx:126
pma::Track3D::HasTPC
bool HasTPC(int tpc) const
Definition: PmaTrack3D.h:110
pma::Hit3D::GetDistToProj
double GetDistToProj(void) const
Definition: PmaHit3D.h:69
pma::SortedBranchBase::AddNext
virtual bool AddNext(pma::SortedObjectBase *nextElement)
Definition: SortedObjects.cxx:71
minY
double minY
Definition: plot_hist.C:9
pma::Track3D::DriftsOfWireIntersection
std::vector< float > DriftsOfWireIntersection(unsigned int wire, unsigned int view) const
Definition: PmaTrack3D.cxx:1093
pma::Track3D::GetMse
double GetMse(unsigned int view=geo::kUnknown) const
MSE of hits weighted with hit amplidudes and wire plane coefficients.
Definition: PmaTrack3D.cxx:1774
pma::Track3D::SortHitsInTree
void SortHitsInTree(bool skipFirst=false)
Definition: PmaTrack3D.cxx:2152
pma::Track3D::AttachTo
bool AttachTo(pma::Node3D *vStart, bool noFlip=false)
Definition: PmaTrack3D.cxx:1456
pma::GetHitsRadius2D
double GetHitsRadius2D(const std::vector< pma::Hit3D * > &hits, bool exact=false)
Definition: Utilities.cxx:78
pma::Track3D::UpdateProjection
void UpdateProjection(void)
Definition: PmaTrack3D.cxx:3081
pma::Track3D::Dist2
double Dist2(const TVector2 &p2d, unsigned int view, unsigned int tpc, unsigned int cryo) const
Definition: PmaTrack3D.cxx:2535
pma::Hit3D::Cryo
unsigned int Cryo(void) const
Definition: PmaHit3D.h:56
PmaTrack3D.h
Implementation of the Projection Matching Algorithm.
pma::Track3D::InsertNode
void InsertNode(TVector3 const &p3d, size_t at_idx, unsigned int tpc, unsigned int cryo)
Definition: PmaTrack3D.cxx:1342
pma::Track3D::IsAttachedTo
bool IsAttachedTo(pma::Track3D const *trk) const
Definition: PmaTrack3D.cxx:1752
pma::Hit3D::View2D
unsigned int View2D(void) const
Definition: PmaHit3D.h:58
pma::Hit3D::fParent
pma::Track3D * fParent
Definition: PmaHit3D.h:106
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:19
pma::Track3D::PrevHit
int PrevHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:877
pma::Track3D::Flip
bool Flip(std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:499
geo::kV
Planes which measure V.
Definition: geo_types.h:77
pma::Node3D::Point3D
TVector3 const & Point3D(void) const
Definition: PmaNode3D.h:33
geo::kUnknown
Unknown view.
Definition: geo_types.h:83
pma::Hit3D::IsEnabled
bool IsEnabled(void) const
Definition: PmaHit3D.h:82
recob::Hit::WireID
geo::WireID WireID() const
Initial tdc tick for hit.
Definition: Hit.h:234
pma::Node3D::ClearAssigned
void ClearAssigned(pma::Track3D *trk=0) override
Definition: PmaNode3D.cxx:766
pma::Element3D::AddPoint
void AddPoint(TVector3 *p)
Definition: PmaElement3D.h:77
pma::SortedBranchBase::RemoveNext
virtual int RemoveNext(pma::SortedObjectBase *nextElement)
Definition: SortedObjects.cxx:107
y
Float_t y
Definition: compare.C:6
lar::fhicl::geo::Point3D
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
Definition: geo_vectors_fhicl.h:38
z
Double_t z
Definition: plot.C:279
Y
Float_t Y
Definition: plot.C:39
geo::BoxBoundedGeo::MinX
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:90
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:37
pma::Track3D::index_of
int index_of(const pma::Hit3D *hit) const
Definition: PmaTrack3D.cxx:333
pma::Node3D::GetDistance2To
double GetDistance2To(const TVector3 &p3d) const override
Distance [cm] from the 3D point to the point 3D.
Definition: PmaNode3D.cxx:133
pma::Track3D::AttachBackToSameTPC
bool AttachBackToSameTPC(pma::Node3D *vStart)
Definition: PmaTrack3D.cxx:1638
pma::Track3D::AddNode
bool AddNode(void)
Definition: PmaTrack3D.cxx:1209
pma::Track3D::GetBranches
bool GetBranches(std::vector< pma::Track3D const * > &branches, bool skipFirst=false) const
Definition: PmaTrack3D.cxx:1727
pma::SortedBranchBase::NextCount
virtual unsigned int NextCount(void) const
Definition: SortedObjects.h:84
pma::Track3D::AttachToOtherTPC
bool AttachToOtherTPC(pma::Node3D *vStart)
Definition: PmaTrack3D.cxx:1502
pma::Track3D::CheckEndSegment
bool CheckEndSegment(pma::Track3D::ETrackEnd endCode)
Definition: PmaTrack3D.cxx:3180
pma::Track3D::GetObjFunction
double GetObjFunction(float penaltyFactor=1.0F) const
Objective function optimized in track reconstruction.
Definition: PmaTrack3D.cxx:1793
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:130
pma::Node3D::Optimize
void Optimize(float penaltyValue, float endSegWeight)
Definition: PmaNode3D.cxx:756
pma::Track3D::GetUnconstrainedProj3D
bool GetUnconstrainedProj3D(art::Ptr< recob::Hit > hit, TVector3 &p3d, double &dist2) const
Definition: PmaTrack3D.cxx:2638
geo::kZ
Planes which measure Z direction.
Definition: geo_types.h:79
maxY
double maxY
Definition: plot_hist.C:10
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:313
pma::Track3D::DeleteSegments
void DeleteSegments(void)
Definition: PmaTrack3D.cxx:2359
tmp
Float_t tmp
Definition: plot.C:37
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:214
pma::Track3D::Split
pma::Track3D * Split(size_t idx, bool try_start_at_idx=true)
Definition: PmaTrack3D.cxx:1370
pma::Track3D::UpdateParams
bool UpdateParams(void)
Definition: PmaTrack3D.cxx:3112
pma::Track3D::AttachBackToOtherTPC
bool AttachBackToOtherTPC(pma::Node3D *vStart)
Definition: PmaTrack3D.cxx:1624
pma::Track3D::fPenaltyValue
float fPenaltyValue
Definition: PmaTrack3D.h:342
pma::Track3D::RemoveHits
void RemoveHits(const std::vector< art::Ptr< recob::Hit > > &hits)
Remove hits; removes also hit->node/seg assignments.
Definition: PmaTrack3D.cxx:399
pma::Hit3D::Wire
unsigned int Wire(void) const
Definition: PmaHit3D.h:59
pma::Vector3D
recob::tracking::Vector_t Vector3D
Definition: Utilities.h:29
pma::Track3D::TuneFullTree
double TuneFullTree(double eps=0.001, double gmax=50.0)
Definition: PmaTrack3D.cxx:2211
pma::Track3D::GetDirection3D
pma::Vector3D GetDirection3D(size_t index) const
Get trajectory direction at given hit index.
Definition: PmaTrack3D.cxx:373
pma::Track3D::TuneSinglePass
double TuneSinglePass(bool skipFirst=false)
Definition: PmaTrack3D.cxx:1959
pma::Track3D::fEndSegWeight
float fEndSegWeight
Definition: PmaTrack3D.h:343
pma::Track3D::AverageDist2
double AverageDist2(void) const
Definition: PmaTrack3D.cxx:3087
pma::Track3D::NEnabledHits
unsigned int NEnabledHits(unsigned int view=geo::kUnknown) const
Definition: PmaTrack3D.cxx:416
pma::Element3D::TPC
int TPC(void) const
TPC index or -1 if out of any TPC.
Definition: PmaElement3D.h:31
pma::Track3D::ShiftEndsToHits
bool ShiftEndsToHits(void)
Definition: PmaTrack3D.cxx:2422
pma::Track3D::~Track3D
~Track3D(void)
Definition: PmaTrack3D.cxx:74
pma::Track3D::fMaxHitsPerSeg
unsigned int fMaxHitsPerSeg
Definition: PmaTrack3D.h:333
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
pma::Track3D::WireDriftRange
std::pair< TVector2, TVector2 > WireDriftRange(unsigned int view, unsigned int tpc, unsigned int cryo) const
Definition: PmaTrack3D.cxx:469
geo::kU
Planes which measure U.
Definition: geo_types.h:76
scale
Double_t scale
Definition: plot.C:25
pma::Element3D::NPoints
size_t NPoints(void) const
Definition: PmaElement3D.h:76
pma::Track3D::fSegStopValue
unsigned int fSegStopValue
Definition: PmaTrack3D.h:337
pma::Track3D::ExtendWith
void ExtendWith(pma::Track3D *src)
Extend the track with everything from src, delete the src;.
Definition: PmaTrack3D.cxx:1684
hits
void hits()
Definition: readHits.C:15
pma::Track3D::release_at
pma::Hit3D * release_at(size_t index)
Definition: PmaTrack3D.cxx:340
geo::vect::isfinite
bool isfinite(Vector const &v)
Returns whether all components of the vector are finite.
Definition: geo_vectors_utils.h:1759
c1
TCanvas * c1
Definition: plotHisto.C:7
pma::Track3D::NHits
unsigned int NHits(unsigned int view) const
Definition: PmaTrack3D.cxx:406
pma::SortedObjectBase::RemoveNext
virtual int RemoveNext(pma::SortedObjectBase *nextElement)
Definition: SortedObjects.cxx:47
pma::Track3D::MakeFastProjection
void MakeFastProjection(void)
Definition: PmaTrack3D.cxx:2950
pma::Track3D::InitFromMiddle
void InitFromMiddle(int tpc, int cryo)
Definition: PmaTrack3D.cxx:299
pma::Track3D::fT0
double fT0
Definition: PmaTrack3D.h:346
pma::Track3D::HasRefPoint
bool HasRefPoint(TVector3 *p) const
Definition: PmaTrack3D.cxx:1767
pma::Track3D::CanFlip
bool CanFlip(void) const
Check if the track can be flipped without breaking any other track.
Definition: PmaTrack3D.cxx:615
pma::Track3D::CompleteMissingWires
size_t CompleteMissingWires(unsigned int view)
Definition: PmaTrack3D.cxx:1118
pma::Track3D::back
pma::Hit3D *& back()
Definition: PmaTrack3D.h:74
d
Float_t d
Definition: plot.C:237
pma::Element3D::Hit
pma::Hit3D & Hit(size_t index)
Definition: PmaElement3D.h:59
pma::Hit3D::SetEnabled
void SetEnabled(bool state)
Definition: PmaHit3D.h:83
pma::Track3D::Optimize
double Optimize(int nNodes=-1, double eps=0.01, bool selAllHits=true, bool setAllNodes=true, size_t selSegHits=0, size_t selVtxHits=0)
Main optimization method.
Definition: PmaTrack3D.cxx:1804
pma::Track3D::UpdateHitsRadius
void UpdateHitsRadius(void)
Definition: PmaTrack3D.cxx:3214
pma::Track3D::fPenaltyFactor
float fPenaltyFactor
Definition: PmaTrack3D.h:334
pma::Node3D::SameTPC
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
Definition: PmaNode3D.cxx:83
pma::Hit3D::SummedADC
float SummedADC(void) const
Definition: PmaHit3D.h:62
pma::Track3D::GetObjFnInTree
double GetObjFnInTree(bool skipFirst=false)
Definition: PmaTrack3D.cxx:2181
pma::Track3D::InternalFlip
void InternalFlip(std::vector< pma::Track3D * > &toSort)
Definition: PmaTrack3D.cxx:565
pma::Track3D::erase
bool erase(const art::Ptr< recob::Hit > &hit)
Definition: PmaTrack3D.cxx:359
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:258
pma::Hit3D::Point2D
TVector2 const & Point2D(void) const
Definition: PmaHit3D.h:53
pma::Track3D::AttachToSameTPC
bool AttachToSameTPC(pma::Node3D *vStart)
Definition: PmaTrack3D.cxx:1515
pma::Node3D::GetDriftShift
double GetDriftShift(void) const
Definition: PmaNode3D.h:102
pma::Track3D::Cryos
std::vector< unsigned int > Cryos(void) const
Definition: PmaTrack3D.cxx:453
pma::Track3D::RebuildSegments
void RebuildSegments(void)
Definition: PmaTrack3D.cxx:2365
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
pma::Track3D::SelectRndHits
bool SelectRndHits(size_t segmax, size_t vtxmax)
Definition: PmaTrack3D.cxx:2893
pma::Track3D::TestHitsMse
double TestHitsMse(const std::vector< art::Ptr< recob::Hit > > &hits, bool normalized=true) const
MSE of 2D hits.
Definition: PmaTrack3D.cxx:788
pma::Track3D::DisableSingleViewEnds
unsigned int DisableSingleViewEnds(void)
Definition: PmaTrack3D.cxx:2711
pma::Track3D::ClearNodes
void ClearNodes(void)
Definition: PmaTrack3D.cxx:120
pma::Track3D::push_back
void push_back(pma::Hit3D *hit)
Definition: PmaTrack3D.h:66
pma::Hit3D::Point3D
TVector3 const & Point3D(void) const
Definition: PmaHit3D.h:48
pma::Track3D::HasTwoViews
bool HasTwoViews(size_t nmin=1) const
Definition: PmaTrack3D.cxx:428
pma::Element3D::Cryo
int Cryo(void) const
Cryostat index or -1 if out of any cryostat.
Definition: PmaElement3D.h:33
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:15
Utilities.h
Implementation of the Projection Matching Algorithm.
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:219
pma::Track3D::MakeProjectionInTree
void MakeProjectionInTree(bool skipFirst=false)
Definition: PmaTrack3D.cxx:2123
pma::Track3D::GetRoot
pma::Track3D * GetRoot(void)
Definition: PmaTrack3D.cxx:1712
pma::Element3D::SortHits
void SortHits(void)
Definition: PmaElement3D.cxx:44
pma::SortedBranchBase::Next
virtual pma::SortedObjectBase * Next(unsigned int index=0) const
Definition: SortedObjects.h:79
pma::Track3D::Length
double Length(size_t step=1) const
Definition: PmaTrack3D.h:81
pma::Element3D::IsFrozen
bool IsFrozen(void) const
Check if the vertex 3D position is fixed.
Definition: PmaElement3D.h:96
pma::Track3D::UpdateParamsInTree
bool UpdateParamsInTree(bool skipFirst=false)
Definition: PmaTrack3D.cxx:1915
dir
TDirectory * dir
Definition: macro.C:5
pma::Track3D::PrevSegment
pma::Segment3D * PrevSegment(pma::Node3D *vtx) const
Definition: PmaTrack3D.cxx:980
pma::Track3D::RemoveNode
bool RemoveNode(size_t idx)
Definition: PmaTrack3D.cxx:1355
pma::Track3D::ReassignHitsInTree
void ReassignHitsInTree(pma::Track3D *plRoot=0)
Definition: PmaTrack3D.cxx:2070
pma::WireDriftToCm
TVector2 WireDriftToCm(unsigned int wire, float drift, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:296
mf::LogDebug
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Definition: MessageLogger.h:211
pma::Track3D::Initialize
bool Initialize(float initEndSegW=0.05F)
Definition: PmaTrack3D.cxx:84
pma::Track3D::fSegStopFactor
float fSegStopFactor
Definition: PmaTrack3D.h:341
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:213
pma::Track3D::NextHit
int NextHit(int index, unsigned int view=geo::kZ, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:861
pma::Element3D::GetDistance2To
virtual double GetDistance2To(const TVector3 &p3d) const =0
Distance [cm] from the 3D point to the object 3D.
pma::Track3D::MakeProjection
void MakeProjection(void)
Definition: PmaTrack3D.cxx:2912
pma::Track3D::SetT0FromDx
void SetT0FromDx(double dx)
Function to convert dx into dT0.
Definition: PmaTrack3D.cxx:2309
pma::Track3D::GetRawdEdxSequence
double GetRawdEdxSequence(std::map< size_t, std::vector< double > > &dedx, unsigned int view=geo::kZ, unsigned int skip=0, bool inclDisabled=false) const
Definition: PmaTrack3D.cxx:990
pma::Track3D::fNodes
std::vector< pma::Node3D * > fNodes
Definition: PmaTrack3D.h:330
h1
TH1F * h1
Definition: plot.C:43
pma::Track3D::HitDxByView
double HitDxByView(size_t index, unsigned int view) const
Definition: PmaTrack3D.cxx:952
pma::CmToWireDrift
TVector2 CmToWireDrift(float xw, float yd, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:307
pma::Segment3D::GetDistance2To
double GetDistance2To(const TVector3 &p3d) const override
Distance [cm] from the 3D segment to the point 3D.
Definition: PmaSegment3D.cxx:24
pma::Element3D::AddHit
void AddHit(pma::Hit3D *h)
Definition: PmaElement3D.h:65
n
Char_t n[5]
Definition: comparison_ascii.C:55
pma::Track3D::SwapVertices
bool SwapVertices(size_t v0, size_t v1)
Definition: PmaTrack3D.cxx:3140
pma::Hit3D::PeakTime
float PeakTime(void) const
Definition: PmaHit3D.h:60
pma::Node3D::fPoint3D
TVector3 fPoint3D
Definition: PmaNode3D.h:129
pma::SortedObjectBase::Disconnect
virtual void Disconnect(void)
Definition: SortedObjects.cxx:26
evd::details::end
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
Definition: RawDataDrawer.cxx:317
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2393
pma::dedx_map
std::map< size_t, std::vector< double > > dedx_map
Definition: Utilities.h:31
pma::Track3D::GetNearestTrkInTree
pma::Track3D * GetNearestTrkInTree(const TVector3 &p3d_cm, double &dist, bool skipFirst=false)
Definition: PmaTrack3D.cxx:1990
pma::SortedObjectBase::AddNext
virtual bool AddNext(pma::SortedObjectBase *nextElement)
Definition: SortedObjects.cxx:32
pma::Track3D::SortHits
void SortHits(void)
Definition: PmaTrack3D.cxx:2670
pma::Track3D::fMaxSegStopFactor
float fMaxSegStopFactor
Definition: PmaTrack3D.h:335
pma::Element3D::SelectAllHits
bool SelectAllHits(void)
Definition: PmaElement3D.cxx:296
pma::SortedObjectBase::Prev
virtual pma::SortedObjectBase * Prev(void) const
Definition: SortedObjects.h:44
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:203
pma::Track3D::fMinSegStop
unsigned int fMinSegStop
Definition: PmaTrack3D.h:338
pma::Track3D::fHitsRadius
float fHitsRadius
Definition: PmaTrack3D.h:344
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:76
X
Float_t X
Definition: plot.C:39
w
Float_t w
Definition: plot.C:23
pma::Segment3D::Parent
pma::Track3D * Parent(void) const
Definition: PmaSegment3D.h:64
pma::Element3D::NHits
size_t NHits(void) const
Definition: PmaElement3D.h:71
pma::Track3D::front
pma::Hit3D *& front()
Definition: PmaTrack3D.h:72
pma::Track3D::fHits
std::vector< pma::Hit3D * > fHits
Definition: PmaTrack3D.h:322
pma::SortedObjectBase::Next
virtual pma::SortedObjectBase * Next(unsigned int index=0) const
Definition: SortedObjects.h:45
pma::Track3D::AttachBackTo
bool AttachBackTo(pma::Node3D *vStart)
Definition: PmaTrack3D.cxx:1587
pma::Track3D::fMaxSegStop
unsigned int fMaxSegStop
Definition: PmaTrack3D.h:339
pma::Track3D::FirstElement
pma::Node3D * FirstElement(void) const
Definition: PmaTrack3D.h:233
pma::Segment3D::SetProjection
void SetProjection(pma::Hit3D &h) const override
Set hit 3D position and its 2D projection to the vertex.
Definition: PmaSegment3D.cxx:149
pma::Track3D::GetNearestElement
pma::Element3D * GetNearestElement(const TVector2 &p2d, unsigned int view, int tpc=-1, bool skipFrontVtx=false, bool skipBackVtx=false) const
Definition: PmaTrack3D.cxx:2573
pma::Track3D::NextSegment
pma::Segment3D * NextSegment(pma::Node3D *vtx) const
Definition: PmaTrack3D.cxx:966
pma::Element3D::Length
double Length(void) const
Definition: PmaElement3D.h:49
Geometry.h
art framework interface to geometry description
pma::Element3D::NEnabledHits
size_t NEnabledHits(unsigned int view=geo::kUnknown) const
Definition: PmaElement3D.cxx:35
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
pma::Track3D::ApplyDriftShiftInTree
void ApplyDriftShiftInTree(double dx, bool skipFirst=false)
Adjust track tree position in the drift direction (when T0 is being corrected).
Definition: PmaTrack3D.cxx:2271
pma::Track3D::SelectAllHits
bool SelectAllHits(void)
Definition: PmaTrack3D.cxx:2901
pma::Track3D::Flip
void Flip(void)
Definition: PmaTrack3D.cxx:596
pma::Track3D::fAssignedPoints
std::vector< TVector3 * > fAssignedPoints
Definition: PmaTrack3D.h:324
pma::Track3D::AddHits
void AddHits(const std::vector< art::Ptr< recob::Hit > > &hits)
Definition: PmaTrack3D.cxx:393
pma::Track3D::InitFromRefPoints
bool InitFromRefPoints(int tpc, int cryo)
Definition: PmaTrack3D.cxx:215