LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PmaNode3D.cxx
Go to the documentation of this file.
1 
14 
19 
22 
23 // Fixed optimization directions: X Y Z
24 bool pma::Node3D::fGradFixed[3] = {false, false, false};
25 
26 double pma::Node3D::fMargin = 3.0;
27 
29  : fTpcGeo(art::ServiceHandle<geo::Geometry const>()->TPC())
30  , fMinX(0)
31  , fMaxX(0)
32  , fMinY(0)
33  , fMaxY(0)
34  , fMinZ(0)
35  , fMaxZ(0)
36  , fPoint3D(0, 0, 0)
37  , fDriftOffset(0)
38  , fIsVertex(false)
39 {
40  fTPC = 0;
41  fCryo = 0;
42 
43  fProj2D[0].Set(0);
44  fProj2D[1].Set(0);
45  fProj2D[2].Set(0);
46 }
47 
49  const TVector3& p3d,
50  unsigned int tpc,
51  unsigned int cryo,
52  bool vtx,
53  double xshift)
54  : fTpcGeo(art::ServiceHandle<geo::Geometry const>()->TPC(geo::TPCID(cryo, tpc)))
55  , fDriftOffset(xshift)
56  , fIsVertex(vtx)
57 {
58  fTPC = tpc;
59  fCryo = cryo;
60 
61  unsigned int lastPlane = geo::kZ;
62  while ((lastPlane > 0) && !fTpcGeo.HasPlane(lastPlane))
63  lastPlane--;
64 
65  fMinX = detProp.ConvertTicksToX(0, lastPlane, tpc, cryo);
66  fMaxX = detProp.ConvertTicksToX(detProp.NumberTimeSamples() - 1, lastPlane, tpc, cryo);
67  if (fMaxX < fMinX) {
68  double tmp = fMaxX;
69  fMaxX = fMinX;
70  fMinX = tmp;
71  }
72 
73  fMinY = fTpcGeo.MinY();
74  fMaxY = fTpcGeo.MaxY();
75  fMinZ = fTpcGeo.MinZ();
76  fMaxZ = fTpcGeo.MaxZ();
77 
78  SetPoint3D(p3d);
79 }
80 
82 {
83  double d, dmin = fPoint3D.X() - fMinX;
84  d = fMaxX - fPoint3D.X();
85  if (d < dmin) dmin = d;
86 
87  d = fPoint3D.Y() - fMinY;
88  if (d < dmin) dmin = d;
89  d = fMaxY - fPoint3D.Y();
90  if (d < dmin) dmin = d;
91 
92  d = fPoint3D.Z() - fMinZ;
93  if (d < dmin) dmin = d;
94  d = fMaxZ - fPoint3D.Z();
95  if (d < dmin) dmin = d;
96 
97  return dmin;
98 }
99 
100 bool pma::Node3D::SameTPC(const TVector3& p3d, float margin) const
101 {
102  if (((fMinX - margin) <= p3d.X()) && (p3d.X() <= (fMaxX + margin)) &&
103  ((fMinY - margin) <= p3d.Y()) && (p3d.Y() <= (fMaxY + margin)) &&
104  ((fMinZ - margin) <= p3d.Z()) && (p3d.Z() <= (fMaxZ + margin)))
105  return true;
106  else
107  return false;
108 }
109 
110 bool pma::Node3D::SameTPC(const pma::Vector3D& p3d, float margin) const
111 {
112  if (((fMinX - margin) <= p3d.X()) && (p3d.X() <= (fMaxX + margin)) &&
113  ((fMinY - margin) <= p3d.Y()) && (p3d.Y() <= (fMaxY + margin)) &&
114  ((fMinZ - margin) <= p3d.Z()) && (p3d.Z() <= (fMaxZ + margin)))
115  return true;
116  else
117  return false;
118 }
119 
121 {
122  bool trimmed = false;
123 
124  if (fPoint3D.X() < fMinX - fMargin) {
125  fPoint3D.SetX(fMinX - fMargin);
126  trimmed = true;
127  }
128  if (fPoint3D.X() > fMaxX + fMargin) {
129  fPoint3D.SetX(fMaxX + fMargin);
130  trimmed = true;
131  }
132 
133  if (fPoint3D.Y() < fMinY - fMargin) {
134  fPoint3D.SetY(fMinY - fMargin);
135  trimmed = true;
136  }
137  if (fPoint3D.Y() > fMaxY + fMargin) {
138  fPoint3D.SetY(fMaxY + fMargin);
139  trimmed = true;
140  }
141 
142  if (fPoint3D.Z() < fMinZ - fMargin) {
143  fPoint3D.SetZ(fMinZ - fMargin);
144  trimmed = true;
145  }
146  if (fPoint3D.Z() > fMaxZ + fMargin) {
147  fPoint3D.SetZ(fMaxZ + fMargin);
148  trimmed = true;
149  }
150 
151  return trimmed;
152 }
153 
155 {
156  for (size_t i = 0; i < fTpcGeo.Nplanes(); ++i) {
158  fPoint3D.X() - fDriftOffset);
159  }
160 }
161 
162 bool pma::Node3D::SetPoint3D(const TVector3& p3d)
163 {
164  fPoint3D = p3d;
165 
166  bool accepted = !LimitPoint3D();
167  UpdateProj2D();
168 
169  return accepted;
170 }
171 
172 double pma::Node3D::GetDistance2To(const TVector3& p3d) const
173 {
174  return pma::Dist2(fPoint3D, p3d);
175 }
176 
177 double pma::Node3D::GetDistance2To(const TVector2& p2d, unsigned int view) const
178 {
179  return pma::Dist2(fProj2D[view], p2d);
180 }
181 
183 {
184  double sum = 0.0F;
185  for (auto h : fAssignedHits) {
186  if (h->IsEnabled()) {
187  unsigned int view = h->View2D();
188 
189  sum += OptFactor(view) * // alpha_i
190  h->GetSigmaFactor() * // hit_amp / hit_max_amp
191  pma::Dist2(h->Point2D(), fProj2D[view]);
192  }
193  }
194  return sum;
195 }
196 
198 {
199  pma::Element3D* seg = 0;
200  if (next) { seg = dynamic_cast<pma::Element3D*>(next); }
201  else if (prev) {
202  seg = dynamic_cast<pma::Element3D*>(prev);
203  }
204  else {
205  throw cet::exception("Node3D") << "Isolated vertex." << std::endl;
206  }
207 
208  if (seg) { return seg->GetDirection3D(); }
209  else {
210  throw cet::exception("Node3D") << "Corrupted vertex." << std::endl;
211  }
212 }
213 
215 {
216  TVector2 gstart;
217  TVector3 g3d;
218  if (prev) {
219  pma::Node3D* vtx = static_cast<pma::Node3D*>(prev->Prev());
220  gstart = vtx->Projection2D(h.View2D());
221  if (!next) g3d = vtx->Point3D();
222  }
223  else if (next) {
224  pma::Node3D* vtx = static_cast<pma::Node3D*>(next->Next());
225  gstart = Projection2D(h.View2D());
226  gstart -= vtx->Projection2D(h.View2D()) - Projection2D(h.View2D());
227  if (!prev) {
228  g3d = fPoint3D;
229  g3d -= vtx->Point3D() - fPoint3D;
230  }
231  }
232  else {
233  mf::LogError("pma::Node3D") << "Isolated vertex.";
234  TVector2 p(Projection2D(h.View2D()));
235  h.SetProjection(p, 0.0F);
236  h.SetPoint3D(fPoint3D);
237  return;
238  }
239 
240  TVector2 v0(h.Point2D());
241  v0 -= Projection2D(h.View2D());
242 
243  TVector2 v1(gstart);
244  v1 -= Projection2D(h.View2D());
245 
246  double v0Norm = v0.Mod();
247  double v1Norm = v1.Mod();
248  double mag = v0Norm * v1Norm;
249  double cosine = 0.0;
250  if (mag != 0.0) cosine = v0 * v1 / mag;
251 
252  TVector2 p(Projection2D(h.View2D()));
253 
254  if (prev && next) {
255  pma::Node3D* vNext = static_cast<pma::Node3D*>(next->Next());
256  TVector2 vN(vNext->Projection2D(h.View2D()));
257  vN -= Projection2D(h.View2D());
258 
259  mag = v0Norm * vN.Mod();
260  double cosineN = 0.0;
261  if (mag != 0.0) cosineN = v0 * vN / mag;
262 
263  // hit on the previous segment side, sorting on the -cosine(prev_seg, point) /max.val. = 1/
264  if (cosineN <= cosine) h.SetProjection(p, -(float)cosine);
265  // hit on the next segment side, sorting on the 1+cosine(next_seg, point) /min.val. = 1/
266  else
267  h.SetProjection(p, 2.0F + (float)cosineN);
268 
269  h.SetPoint3D(fPoint3D);
270  }
271  else {
272  float b = (float)(v0Norm * cosine / v1Norm);
273  if (fFrozen) // limit 3D positions to outermose node if frozen
274  {
275  h.SetPoint3D(fPoint3D);
276  }
277  else // or set 3D positions along the line of outermost segment
278  {
279  g3d -= fPoint3D;
280  h.SetPoint3D(fPoint3D + (g3d * b));
281 
282  p += (v1 * b);
283  }
284  h.SetProjection(p, -b);
285  }
286 }
287 
288 double pma::Node3D::Length2() const
289 {
290  double l = 0.0;
291  if (next) l += (static_cast<pma::Segment3D*>(next))->Length();
292  if (prev) l += (static_cast<pma::Segment3D*>(prev))->Length();
293 
294  if (next && prev)
295  return 0.25 * l * l;
296  else
297  return l * l;
298 }
299 
301 {
302  if (prev && next) {
303  auto const& vStop1 = static_cast<pma::Node3D*>(prev->Prev())->fPoint3D;
304  auto const& vStop2 = static_cast<pma::Node3D*>(next->Next())->fPoint3D;
305  pma::Vector3D v1(
306  vStop1.X() - fPoint3D.X(), vStop1.Y() - fPoint3D.Y(), vStop1.Z() - fPoint3D.Z());
307  pma::Vector3D v2(
308  vStop2.X() - fPoint3D.X(), vStop2.Y() - fPoint3D.Y(), vStop2.Z() - fPoint3D.Z());
309  double mag = sqrt(v1.Mag2() * v2.Mag2());
310  double cosine = 0.0;
311  if (mag != 0.0) cosine = v1.Dot(v2) / mag;
312  return cosine;
313  }
314  else {
315  mf::LogError("pma::Node3D") << "pma::Node3D::SegmentCos(): neighbours not initialized.";
316  return -1.0;
317  }
318 }
319 
321 {
322  if (prev && next) {
323  auto const& vStop1 = static_cast<pma::Node3D*>(prev->Prev())->fPoint3D;
324  auto const& vStop2 = static_cast<pma::Node3D*>(next->Next())->fPoint3D;
325  pma::Vector2D v1(vStop1.Y() - fPoint3D.Y(), vStop1.Z() - fPoint3D.Z());
326  pma::Vector2D v2(vStop2.Y() - fPoint3D.Y(), vStop2.Z() - fPoint3D.Z());
327  double mag = sqrt(v1.Mag2() * v2.Mag2());
328  double cosine = 0.0;
329  if (mag != 0.0) cosine = v1.Dot(v2) / mag;
330  return cosine;
331  }
332  else {
333  mf::LogError("pma::Node3D") << "pma::Node3D::SegmentCosZX(): neighbours not initialized.";
334  return -1.0;
335  }
336 }
337 
339 {
340  if (prev && next) {
341  auto const& vStop1 = static_cast<pma::Node3D*>(prev->Prev())->fPoint3D;
342  auto const& vStop2 = static_cast<pma::Node3D*>(next->Next())->fPoint3D;
343  pma::Vector2D v1(vStop1.X() - fPoint3D.X(), vStop1.Z() - fPoint3D.Z());
344  pma::Vector2D v2(vStop2.X() - fPoint3D.X(), vStop2.Z() - fPoint3D.Z());
345  double mag = sqrt(v1.Mag2() * v2.Mag2());
346  double cosine = 0.0;
347  if (mag != 0.0) cosine = v1.Dot(v2) / mag;
348  return cosine;
349  }
350  else {
351  mf::LogError("pma::Node3D") << "pma::Node3D::SegmentCosZY(): neighbours not initialized.";
352  return -1.0;
353  }
354 }
355 
356 // *** Note: should be changed / generalized for horizontal wire planes (e.g. 2-phase LAr). ***
358 {
359  if (prev && next) {
360  auto const& vStart = static_cast<pma::Node3D*>(prev->Prev())->fPoint3D;
361  auto const& vStop = static_cast<pma::Node3D*>(next->Next())->fPoint3D;
362 
363  double dy = vStop.X() - vStart.X();
364  double dz = vStop.Z() - vStart.Z();
365  double len2 = dy * dy + dz * dz;
366  double cosine2 = 0.0;
367  if (len2 > 0.0) cosine2 = dz * dz / len2;
368  return cosine2;
369  }
370  else
371  return 0.0;
372 }
373 
375 {
376  if (prev && NextCount()) {
377  pma::Segment3D* seg0 = dynamic_cast<pma::Segment3D*>(prev);
378  pma::Segment3D* seg1 = dynamic_cast<pma::Segment3D*>(Next(0));
379  unsigned int nInd1 = NHits(geo::kU) + seg0->NHits(geo::kU) + seg1->NHits(geo::kU);
380 
381  if (fHitsRadius > 0.0F)
382  return (1.0 + SegmentCosWirePlane()) * fHitsRadius * fHitsRadius / (4 * nInd1 + 1.0);
383  else
384  return (1.0 + SegmentCosWirePlane()) * Length2() / (4 * nInd1 + 1.0);
385  }
386  else
387  return 0.0;
388 }
389 
390 // Constraint on two segments angle in projection to plane parallel to wire plane, suppressed by
391 // the orientation in the plane transverse to wire plane (only sections with low variation of
392 // drift time are penalized with this constraint); PiInWirePlane() components are reduced if
393 // there are Ind1 / geo::kU hits which add information to the object shape.
394 // *** Note: should be changed / generalized for horizontal wire planes (e.g. 2-phase LAr). ***
396 {
397  if (fIsVertex) return 0.0;
398 
399  unsigned int nseg = 1;
400  double penalty = PiInWirePlane();
401  pma::Node3D* v;
402  if (next) {
403  v = static_cast<pma::Node3D*>(next->Next());
404  penalty += v->PiInWirePlane();
405  nseg++;
406  }
407  if (prev) {
408  v = static_cast<pma::Node3D*>(prev->Prev());
409  penalty += v->PiInWirePlane();
410  nseg++;
411  }
412  if (penalty > 0.0)
413  return pow(EndPtCos2Transverse(), 10) * penalty / nseg;
414  else
415  return 0.0;
416 }
417 
419 {
420  size_t nnext = NextCount();
421  if (nnext > 1) return true; // 1 trk -> vtx -> n*trk
422 
423  if (prev && nnext) {
424  pma::Segment3D* segPrev = static_cast<pma::Segment3D*>(prev);
425  pma::Segment3D* segNext = static_cast<pma::Segment3D*>(Next(0));
426  if (segNext->Parent() != segPrev->Parent()) // 1 trk -> vtx -> 1 trk
427  return true;
428  }
429  return false;
430 }
431 
433 {
434  if (prev && (NextCount() == 1)) {
435  pma::Segment3D* segPrev = static_cast<pma::Segment3D*>(prev);
436  pma::Segment3D* segNext = static_cast<pma::Segment3D*>(Next(0));
437 
438  if ((segPrev->TPC() < 0) || (segNext->TPC() < 0)) return true;
439  }
440  return false;
441 }
442 
443 std::vector<pma::Track3D*> pma::Node3D::GetBranches() const
444 {
445  std::vector<pma::Track3D*> branches;
446  if (NextCount()) {
447  branches.reserve(NextCount());
448  for (size_t i = 0; i < NextCount(); ++i) {
449  pma::Segment3D* seg = static_cast<pma::Segment3D*>(Next(i));
450  branches.push_back(seg->Parent());
451  }
452  }
453  return branches;
454 }
455 
456 double pma::Node3D::Pi(float endSegWeight, bool doAsymm) const
457 {
458  if (fIsVertex) return 0.0;
459 
460  if (prev && NextCount()) {
461  pma::Segment3D* segPrev = static_cast<pma::Segment3D*>(prev);
462  pma::Segment3D* segNext = static_cast<pma::Segment3D*>(Next(0));
463 
464  double scale = 1.0;
465  if ((segPrev->TPC() < 0) || (segNext->TPC() < 0))
466  scale = 0.5; // lower penalty on segments between tpc's
467 
468  double segCos = SegmentCos();
469 
470  double lAsymmFactor = 0.0;
471  if (doAsymm) {
472  double lPrev = segPrev->Length();
473  double lNext = segNext->Length();
474  double lSum = lPrev + lNext;
475  if (lSum > 0.1) {
476  double lAsymm = (1.0 - segCos) * (lPrev - lNext) / lSum;
477  lAsymmFactor = 0.05 * lAsymm * lAsymm;
478  }
479  }
480 
481  if (fHitsRadius > 0.0F)
482  return scale * (1.0 + segCos + lAsymmFactor) * fHitsRadius * fHitsRadius;
483  else
484  return scale * (1.0 + segCos + lAsymmFactor) * Length2();
485  }
486  else {
487  double pi_result = 0.0;
488  unsigned int nSeg = 0;
489  pma::Segment3D* seg = 0;
490  if (prev) {
491  seg = static_cast<pma::Segment3D*>(prev);
492 
493  SortedObjectBase* prevVtx = seg->Prev();
494  if (prevVtx->Prev()) nSeg++;
495  nSeg += prevVtx->NextCount();
496  }
497  else if (next) {
498  seg = static_cast<pma::Segment3D*>(next);
499 
500  SortedObjectBase* nextVtx = seg->Next();
501  nSeg += nextVtx->NextCount() + 1;
502  }
503  else {
504  mf::LogWarning("pma::Node3D") << "pma::Node3D::Pi(): an isolated vertex?";
505  return 0.0;
506  }
507  if (nSeg == 1) pi_result = endSegWeight * seg->Length2();
508  return pi_result;
509  }
510 }
511 
512 double pma::Node3D::Penalty(float endSegWeight) const
513 {
514  unsigned int nseg = 1;
515  double penalty = Pi(endSegWeight, true);
516 
517  pma::Node3D* v;
518  for (unsigned int i = 0; i < NextCount(); i++) {
519  v = static_cast<pma::Node3D*>(Next(i)->Next());
520  penalty += v->Pi(endSegWeight, false);
521  nseg++;
522  }
523  if (prev) {
524  v = static_cast<pma::Node3D*>(prev->Prev());
525  penalty += v->Pi(endSegWeight, false);
526  nseg++;
527  }
528  return penalty / nseg;
529 }
530 
531 double pma::Node3D::Mse() const
532 {
533  unsigned int nhits = NPrecalcEnabledHits(); //NEnabledHits();
534  double mse = SumDist2();
535 
536  pma::Segment3D* seg;
537  for (unsigned int i = 0; i < NextCount(); i++) {
538  seg = static_cast<pma::Segment3D*>(Next(i));
539  nhits += seg->NPrecalcEnabledHits(); //NEnabledHits();
540  mse += seg->SumDist2();
541  }
542  if (prev) {
543  seg = static_cast<pma::Segment3D*>(prev);
544  nhits += seg->NPrecalcEnabledHits(); //NEnabledHits();
545  mse += seg->SumDist2();
546  }
547  if (!nhits)
548  return 0.0;
549  else
550  return mse / nhits;
551 }
552 
553 double pma::Node3D::GetObjFunction(float penaltyValue, float endSegWeight) const
554 {
555  return Mse() + penaltyValue * (Penalty(endSegWeight) + PenaltyInWirePlane());
556 }
557 
558 double pma::Node3D::MakeGradient(float penaltyValue, float endSegWeight)
559 {
560  double l1 = 0.0, l2 = 0.0, minLength2 = 0.0;
561  TVector3 tmp(fPoint3D), gpoint(fPoint3D);
562 
563  pma::Segment3D* seg;
564  if (prev) {
565  seg = static_cast<pma::Segment3D*>(prev);
566  l1 = seg->Length2();
567  }
568  if (next) {
569  seg = static_cast<pma::Segment3D*>(next);
570  l2 = seg->Length2();
571  }
572  if ((l1 > 0.01) && (l1 < l2))
573  minLength2 = l1;
574  else if ((l2 > 0.01) && (l2 < l1))
575  minLength2 = l2;
576  else
577  minLength2 = 0.01;
578 
579  double dxi = 0.001 * sqrt(minLength2);
580 
581  if (dxi < 6.0E-37) return 0.0;
582 
583  double gi, g0, gz;
584  gz = g0 = GetObjFunction(penaltyValue, endSegWeight);
585 
586  //if (fQPenaltyFactor > 0.0F) gz += fQPenaltyFactor * QPenalty(); <----------------------- maybe later..
587 
588  if (!fGradFixed[0]) // gradX
589  {
590  gpoint[0] = tmp[0] + dxi;
591  SetPoint3D(gpoint);
592  gi = GetObjFunction(penaltyValue, endSegWeight);
593  fGradient[0] = (g0 - gi) / dxi;
594 
595  gpoint[0] = tmp[0] - dxi;
596  SetPoint3D(gpoint);
597  gi = GetObjFunction(penaltyValue, endSegWeight);
598  fGradient[0] = 0.5 * (fGradient[0] + (gi - g0) / dxi);
599 
600  gpoint[0] = tmp[0];
601  }
602 
603  if (!fGradFixed[1]) // gradY
604  {
605  gpoint[1] = tmp[1] + dxi;
606  SetPoint3D(gpoint);
607  gi = GetObjFunction(penaltyValue, endSegWeight);
608  fGradient[1] = (g0 - gi) / dxi;
609 
610  gpoint[1] = tmp[1] - dxi;
611  SetPoint3D(gpoint);
612  gi = GetObjFunction(penaltyValue, endSegWeight);
613  fGradient[1] = 0.5 * (fGradient[1] + (gi - g0) / dxi);
614 
615  gpoint[1] = tmp[1];
616  }
617 
618  if (!fGradFixed[2]) // gradZ
619  {
620  gpoint[2] = tmp[2] + dxi;
621  SetPoint3D(gpoint);
622  gi = GetObjFunction(penaltyValue, endSegWeight);
623  //if (fQPenaltyFactor > 0.0F) gi += fQPenaltyFactor * QPenalty();
624  fGradient[2] = (gz - gi) / dxi;
625 
626  gpoint[2] = tmp[2] - dxi;
627  SetPoint3D(gpoint);
628  gi = GetObjFunction(penaltyValue, endSegWeight);
629  //if (fQPenaltyFactor > 0.0F) gi += fQPenaltyFactor * QPenalty();
630  fGradient[2] = 0.5 * (fGradient[2] + (gi - gz) / dxi);
631 
632  gpoint[2] = tmp[2];
633  }
634 
635  SetPoint3D(tmp);
636  if (fGradient.Mag2() < 6.0E-37) return 0.0;
637 
638  return g0;
639 }
640 
641 double pma::Node3D::StepWithGradient(float alfa, float tol, float penalty, float weight)
642 {
643  unsigned int steps = 0;
644  double t, t1, t2, t3, g, g0, g1, g2, g3, p1, p2;
645  double eps = 6.0E-37, zero_tol = 1.0E-15;
646  TVector3 tmp(fPoint3D), gpoint(fPoint3D);
647 
648  g = MakeGradient(penalty, weight);
649  if (g < zero_tol) return 0.0;
650  g0 = g;
651 
652  //**** first three points ****//
653  alfa *= 0.8F;
654  t2 = 0.0;
655  g2 = g;
656  t3 = 0.0;
657  g3 = g;
658  do {
659  t1 = t2;
660  g1 = g2;
661  t2 = t3;
662  g2 = g3;
663 
664  alfa *= 1.25F;
665  t3 += alfa;
666  gpoint = tmp;
667  gpoint += (fGradient * t3);
668  if (!SetPoint3D(gpoint)) // stepped out of allowed volume
669  {
670  //std::cout << "**** SetPoint trimmed 1 ****" << std::endl;
671  g3 = GetObjFunction(penalty, weight);
672  if (g3 < g2)
673  return (g0 - g3) / g3; // exit with the node at the border
674  else {
675  SetPoint3D(tmp);
676  return 0.0;
677  } // exit with restored original position
678  }
679 
680  g3 = GetObjFunction(penalty, weight);
681 
682  if (g3 < zero_tol) return 0.0;
683 
684  if (++steps > 1000) {
685  SetPoint3D(tmp);
686  return 0.0;
687  }
688 
689  } while (g3 < g2);
690  //****************************//
691 
692  //*** first step overshoot ***//
693  if (steps == 1) {
694  t2 = t3;
695  g2 = g3;
696  do {
697  t3 = t2;
698  g3 = g2;
699  t2 = (t1 * g3 + t3 * g1) / (g1 + g3);
700 
701  // small shift...
702  t2 = 0.05 * t3 + 0.95 * t2;
703 
704  // break: starting point is at the minimum
705  //if (t2 == t1) { SetPoint3D(tmp); return 0.0F; }
706 
707  // break: starting point is very close to the minimum
708  if (fabs(t2 - t1) < tol) {
709  SetPoint3D(tmp);
710  return 0.0;
711  }
712 
713  gpoint = tmp;
714  gpoint += (fGradient * t2);
715  if (!SetPoint3D(gpoint)) // select the best point to exit
716  {
717  //std::cout << "**** SetPoint trimmed 2 ****" << std::endl;
718  g2 = GetObjFunction(penalty, weight);
719  if (g2 < g0)
720  return (g0 - g2) / g2; // exit with the node at the border
721  else if (g1 < g0) {
722  gpoint = tmp;
723  gpoint += (fGradient * t1);
724  return (g0 - g1) / g1;
725  }
726  else if (g3 < g0) {
727  gpoint = tmp;
728  gpoint += (fGradient * t3);
729  return (g0 - g3) / g3;
730  }
731  else {
732  SetPoint3D(tmp);
733  return 0.0;
734  }
735  }
736  g2 = GetObjFunction(penalty, weight);
737 
738  if (g2 < zero_tol) return 0.0;
739  steps++;
740 
741  } while (g2 >= g1);
742  }
743  //****************************//
744 
745  while (fabs(t1 - t3) > tol) {
746  //*** 3-point minimization ***//
747  if ((fabs(t2 - t1) < eps) || (fabs(t2 - t3) < eps)) break; // minimum on the edge
748  if ((fabs(g2 - g1) < eps) && (fabs(g2 - g3) < eps)) break; // ~singularity
749 
750  p1 = (t2 - t1) * (g2 - g3);
751  p2 = (t2 - t3) * (g2 - g1);
752  if (fabs(p1 - p2) < eps) break; // ~linearity
753 
754  t = t2 + ((t2 - t1) * p1 - (t2 - t3) * p2) / (2 * (p2 - p1));
755  if ((t <= t1) || (t >= t3)) t = (t1 * g3 + t3 * g1) / (g1 + g3);
756 
757  gpoint = tmp;
758  gpoint += (fGradient * t);
759  if (!SetPoint3D(gpoint)) // select the best point to exit
760  {
761  //std::cout << "**** SetPoint trimmed 3 ****" << std::endl;
762  g = GetObjFunction(penalty, weight);
763  if ((g < g0) && (g < g1) && (g < g3))
764  return (g0 - g) / g; // exit with the node at the border
765  else if ((g1 < g0) && (g1 < g3)) {
766  gpoint = tmp;
767  gpoint += (fGradient * t1);
768  return (g0 - g1) / g1;
769  }
770  else if (g3 < g0) {
771  gpoint = tmp;
772  gpoint += (fGradient * t3);
773  return (g0 - g3) / g3;
774  }
775  else {
776  SetPoint3D(tmp);
777  return 0.0;
778  }
779  }
780 
781  g = GetObjFunction(penalty, weight);
782  if (g < zero_tol) return 0.0;
783  steps++;
784  //****************************//
785 
786  //*** select next 3 points ***//
787  if (fabs(t - t2) < 0.2 * tol) break; // start in a new direction
788  if (g < g2) {
789  if (t < t2) {
790  t3 = t2;
791  g3 = g2;
792  }
793  else {
794  t1 = t2;
795  g1 = g2;
796  }
797  t2 = t;
798  g2 = g;
799  }
800  else {
801  if (t < t2) {
802  t1 = t;
803  g1 = g;
804  }
805  else {
806  t3 = t;
807  g3 = g;
808  }
809  }
810  //****************************//
811  }
812 
813  return (g0 - g) / g;
814 }
815 
816 void pma::Node3D::Optimize(float penaltyValue, float endSegWeight)
817 {
818  if (!fFrozen) {
819  double dg = StepWithGradient(0.1F, 0.002F, penaltyValue, endSegWeight);
820  if (dg > 0.01) dg = StepWithGradient(0.03F, 0.0001F, penaltyValue, endSegWeight);
821  if (dg > 0.0) dg = StepWithGradient(0.03F, 0.0001F, penaltyValue, endSegWeight);
822  }
823 }
824 
826 {
827  if (!trk) {
828  // like in the base class:
829  fAssignedPoints.clear();
830  fAssignedHits.clear();
831  }
832  else {
833  std::vector<pma::Track3D*> to_check;
834  pma::Segment3D* seg;
835  if (Prev()) {
836  seg = static_cast<pma::Segment3D*>(Prev());
837  if (seg->Parent() != trk) to_check.push_back(seg->Parent());
838  }
839  for (unsigned int i = 0; i < NextCount(); i++) {
840  seg = static_cast<pma::Segment3D*>(Next(i));
841  if (seg->Parent() != trk) to_check.push_back(seg->Parent());
842  }
843 
844  unsigned int p = 0;
845  while (p < fAssignedPoints.size()) {
846  bool found = false;
847  for (size_t t = 0; t < to_check.size(); t++)
848  if (to_check[t]->HasRefPoint(fAssignedPoints[p])) {
849  found = true;
850  break;
851  }
852 
853  if (!found)
854  fAssignedPoints.erase(fAssignedPoints.begin() + p);
855  else
856  p++;
857  }
858 
859  unsigned int h = 0;
860  while (h < fAssignedHits.size()) {
861  bool found = false;
863 
864  for (size_t t = 0; (t < to_check.size()) && !found; t++)
865  for (size_t i = 0; i < to_check[t]->size(); i++) {
866  pma::Hit3D* pmaHit = static_cast<pma::Hit3D*>((*(to_check[t]))[i]);
867  if (hit == pmaHit) {
868  found = true;
869  break;
870  }
871  }
872 
873  if (!found)
874  fAssignedHits.erase(fAssignedHits.begin() + h);
875  else
876  h++;
877  }
878  }
879 
880  fHitsRadius = 0.0F;
881 }
size_t NPrecalcEnabledHits(void) const
Definition: PmaElement3D.h:78
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:147
TVector3 const & Point3D() const
Definition: PmaNode3D.h:44
TVector2 const & Projection2D(unsigned int view) const
Definition: PmaNode3D.h:50
double Length2() const override
Definition: PmaNode3D.cxx:288
TTree * t1
Definition: plottest35.C:26
std::vector< TVector3 * > fAssignedPoints
Definition: PmaElement3D.h:123
double Penalty(float endSegWeight) const
Definition: PmaNode3D.cxx:512
Implementation of the Projection Matching Algorithm.
TVector2 const & Point2D() const noexcept
Definition: PmaHit3D.h:55
double fMinZ
Definition: PmaNode3D.h:149
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:137
void ClearAssigned(pma::Track3D *trk=0) override
Definition: PmaNode3D.cxx:825
Implementation of the Projection Matching Algorithm.
double fDriftOffset
Definition: PmaNode3D.h:155
pma::SortedObjectBase * prev
Definition: SortedObjects.h:54
double GetDistance2To(const TVector3 &p3d) const override
Distance [cm] from the 3D point to the point 3D.
Definition: PmaNode3D.cxx:172
pma::Vector3D GetDirection3D() const override
Definition: PmaNode3D.cxx:197
bool LimitPoint3D()
Returns true if node position was trimmed to its TPC volume + fMargin.
Definition: PmaNode3D.cxx:120
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
virtual unsigned int NextCount(void) const
Definition: SortedObjects.h:84
std::vector< pma::Track3D * > GetBranches() const
Definition: PmaNode3D.cxx:443
void Optimize(float penaltyValue, float endSegWeight)
Definition: PmaNode3D.cxx:816
Planes which measure Z direction.
Definition: geo_types.h:138
Float_t tmp
Definition: plot.C:35
virtual pma::SortedObjectBase * Next(unsigned int=0) const
Definition: SortedObjects.h:43
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double Length2(void) const override
bool IsTPCEdge() const
Is the first/last in this TPC?
Definition: PmaNode3D.cxx:432
recob::tracking::Vector_t Vector3D
Definition: Utilities.h:34
static bool fGradFixed[3]
Definition: PmaNode3D.h:160
virtual pma::Vector3D GetDirection3D(void) const =0
Get 3D direction cosines corresponding to this element.
double SegmentCosWirePlane() const
Definition: PmaNode3D.cxx:320
int TPC(void) const
TPC index or -1 if out of any TPC.
Definition: PmaElement3D.h:35
double PenaltyInWirePlane() const
Definition: PmaNode3D.cxx:395
TVector2 fProj2D[3]
Definition: PmaNode3D.h:153
double Mse() const
Definition: PmaNode3D.cxx:531
Planes which measure U.
Definition: geo_types.h:135
Double_t scale
Definition: plot.C:24
double GetDistToWall() const
Definition: PmaNode3D.cxx:81
virtual unsigned int NextCount(void) const
Definition: SortedObjects.h:44
double fMaxX
Definition: PmaNode3D.h:149
double PiInWirePlane() const
Definition: PmaNode3D.cxx:374
Float_t E
Definition: plot.C:20
void UpdateProj2D()
Definition: PmaNode3D.cxx:154
double SegmentCos() const
Cosine of 3D angle between connected segments.
Definition: PmaNode3D.cxx:300
bool SetPoint3D(const TVector3 &p3d)
Definition: PmaNode3D.cxx:162
double MinZ() const
Returns the world z coordinate of the start of the box.
double StepWithGradient(float alfa, float tol, float penalty, float weight)
Definition: PmaNode3D.cxx:641
Float_t d
Definition: plot.C:235
double SegmentCosTransverse() const
Definition: PmaNode3D.cxx:338
double EndPtCos2Transverse() const
Definition: PmaNode3D.cxx:357
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
Definition: PmaNode3D.cxx:100
double fMaxY
Definition: PmaNode3D.h:149
void SetProjection(const TVector2 &p, float b)
Definition: PmaHit3D.h:75
double SumDist2(void) const
Detector simulation of raw signals on wires.
double fMaxZ
Definition: PmaNode3D.h:149
TTree * t2
Definition: plottest35.C:36
double MaxY() const
Returns the world y coordinate of the end of the box.
std::vector< pma::Hit3D * > fAssignedHits
Definition: PmaElement3D.h:122
double ConvertTicksToX(double ticks, int p, int t, int c) const
double MakeGradient(float penaltyValue, float endSegWeight)
Definition: PmaNode3D.cxx:558
double fMinY
Definition: PmaNode3D.h:149
unsigned int View2D() const noexcept
Definition: PmaHit3D.h:60
double weight
Definition: plottest35.C:25
Implementation of the Projection Matching Algorithm.
virtual pma::SortedObjectBase * Next(unsigned int index=0) const
Definition: SortedObjects.h:77
Encapsulate the construction of a single detector plane.
double GetObjFunction(float penaltyValue, float endSegWeight) const
Objective function minimized during oprimization.
Definition: PmaNode3D.cxx:553
double MaxZ() const
Returns the world z coordinate of the end of the box.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void SetPoint3D(const TVector3 &p3d)
Definition: PmaHit3D.h:52
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double > > Vector2D
Definition: Utilities.h:33
double fHitsRadius
Definition: PmaElement3D.h:128
bool IsBranching() const
Belongs to more than one track?
Definition: PmaNode3D.cxx:418
pma::SortedObjectBase * next
Definition: SortedObjects.h:53
double SumDist2Hits() const override
Definition: PmaNode3D.cxx:182
Definition: MVAAlg.h:12
geo::TPCGeo const & fTpcGeo
Definition: PmaNode3D.h:147
TVector3 fPoint3D
Definition: PmaNode3D.h:152
IDparameter< geo::TPCID > TPCID
Member type of validated geo::TPCID parameter.
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:252
virtual pma::SortedObjectBase * Prev(void) const
Definition: SortedObjects.h:42
Double_t sum
Definition: plot.C:31
static double fMargin
Definition: PmaNode3D.h:161
Namespace collecting geometry-related classes utilities.
void SetProjection(pma::Hit3D &h) const override
Set hit 3D position and its 2D projection to the vertex.
Definition: PmaNode3D.cxx:214
double MinY() const
Returns the world y coordinate of the start of the box.
pma::Track3D * Parent(void) const
Definition: PmaSegment3D.h:66
size_t NHits(void) const
Definition: PmaElement3D.h:76
double Pi(float endSegWeight, bool doAsymm) const
Definition: PmaNode3D.cxx:456
double PlaneCoordinate(geo::Point_t const &point) const
Returns the coordinate of the point on the plane.
Definition: PlaneGeo.h:777
bool fIsVertex
Definition: PmaNode3D.h:158
double Length(void) const
Definition: PmaElement3D.h:53
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
static float OptFactor(unsigned int view)
Definition: PmaElement3D.h:112
double fMinX
Definition: PmaNode3D.h:149
Encapsulate the construction of a single detector plane.
TVector3 fGradient
Definition: PmaNode3D.h:157