29 : fTpcGeo(
art::ServiceHandle<
geo::Geometry const>()->TPC())
61 unsigned int lastPlane =
geo::kZ;
85 if (d < dmin) dmin =
d;
88 if (d < dmin) dmin =
d;
90 if (d < dmin) dmin =
d;
93 if (d < dmin) dmin =
d;
95 if (d < dmin) dmin =
d;
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)))
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)))
122 bool trimmed =
false;
186 if (h->IsEnabled()) {
187 unsigned int view = h->View2D();
190 h->GetSigmaFactor() *
205 throw cet::exception(
"Node3D") <<
"Isolated vertex." << std::endl;
210 throw cet::exception(
"Node3D") <<
"Corrupted vertex." << std::endl;
246 double v0Norm = v0.Mod();
247 double v1Norm = v1.Mod();
248 double mag = v0Norm * v1Norm;
250 if (mag != 0.0) cosine = v0 * v1 / mag;
259 mag = v0Norm * vN.Mod();
260 double cosineN = 0.0;
261 if (mag != 0.0) cosineN = v0 * vN / mag;
272 float b = (float)(v0Norm * cosine / v1Norm);
309 double mag = sqrt(v1.Mag2() * v2.Mag2());
311 if (mag != 0.0) cosine = v1.Dot(v2) / mag;
315 mf::LogError(
"pma::Node3D") <<
"pma::Node3D::SegmentCos(): neighbours not initialized.";
327 double mag = sqrt(v1.Mag2() * v2.Mag2());
329 if (mag != 0.0) cosine = v1.Dot(v2) / mag;
333 mf::LogError(
"pma::Node3D") <<
"pma::Node3D::SegmentCosZX(): neighbours not initialized.";
345 double mag = sqrt(v1.Mag2() * v2.Mag2());
347 if (mag != 0.0) cosine = v1.Dot(v2) / mag;
351 mf::LogError(
"pma::Node3D") <<
"pma::Node3D::SegmentCosZY(): neighbours not initialized.";
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;
399 unsigned int nseg = 1;
421 if (nnext > 1)
return true;
438 if ((segPrev->
TPC() < 0) || (segNext->
TPC() < 0))
return true;
445 std::vector<pma::Track3D*> branches;
448 for (
size_t i = 0; i <
NextCount(); ++i) {
450 branches.push_back(seg->
Parent());
465 if ((segPrev->
TPC() < 0) || (segNext->
TPC() < 0))
470 double lAsymmFactor = 0.0;
472 double lPrev = segPrev->
Length();
473 double lNext = segNext->
Length();
474 double lSum = lPrev + lNext;
476 double lAsymm = (1.0 - segCos) * (lPrev - lNext) / lSum;
477 lAsymmFactor = 0.05 * lAsymm * lAsymm;
484 return scale * (1.0 + segCos + lAsymmFactor) *
Length2();
487 double pi_result = 0.0;
488 unsigned int nSeg = 0;
494 if (prevVtx->
Prev()) nSeg++;
504 mf::LogWarning(
"pma::Node3D") <<
"pma::Node3D::Pi(): an isolated vertex?";
507 if (nSeg == 1) pi_result = endSegWeight * seg->
Length2();
514 unsigned int nseg = 1;
515 double penalty =
Pi(endSegWeight,
true);
518 for (
unsigned int i = 0; i <
NextCount(); i++) {
520 penalty += v->
Pi(endSegWeight,
false);
525 penalty += v->
Pi(endSegWeight,
false);
528 return penalty / nseg;
537 for (
unsigned int i = 0; i <
NextCount(); i++) {
560 double l1 = 0.0, l2 = 0.0, minLength2 = 0.0;
572 if ((l1 > 0.01) && (l1 < l2))
574 else if ((l2 > 0.01) && (l2 < l1))
579 double dxi = 0.001 * sqrt(minLength2);
581 if (dxi < 6.0
E-37)
return 0.0;
590 gpoint[0] =
tmp[0] + dxi;
595 gpoint[0] =
tmp[0] - dxi;
605 gpoint[1] =
tmp[1] + dxi;
610 gpoint[1] =
tmp[1] - dxi;
620 gpoint[2] =
tmp[2] + dxi;
626 gpoint[2] =
tmp[2] - dxi;
636 if (
fGradient.Mag2() < 6.0E-37)
return 0.0;
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;
649 if (g < zero_tol)
return 0.0;
673 return (g0 - g3) / g3;
682 if (g3 < zero_tol)
return 0.0;
684 if (++steps > 1000) {
699 t2 = (t1 * g3 + t3 * g1) / (g1 + g3);
702 t2 = 0.05 * t3 + 0.95 *
t2;
708 if (fabs(t2 - t1) < tol) {
720 return (g0 - g2) / g2;
724 return (g0 - g1) / g1;
729 return (g0 - g3) / g3;
738 if (g2 < zero_tol)
return 0.0;
745 while (fabs(t1 - t3) > tol) {
747 if ((fabs(t2 - t1) < eps) || (fabs(t2 - t3) < eps))
break;
748 if ((fabs(g2 - g1) < eps) && (fabs(g2 - g3) < eps))
break;
750 p1 = (t2 -
t1) * (g2 - g3);
751 p2 = (t2 - t3) * (g2 - g1);
752 if (fabs(p1 - p2) < eps)
break;
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);
763 if ((g < g0) && (g < g1) && (g < g3))
765 else if ((g1 < g0) && (g1 < g3)) {
768 return (g0 - g1) / g1;
773 return (g0 - g3) / g3;
782 if (g < zero_tol)
return 0.0;
787 if (fabs(t - t2) < 0.2 * tol)
break;
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);
833 std::vector<pma::Track3D*> to_check;
837 if (seg->
Parent() != trk) to_check.push_back(seg->
Parent());
839 for (
unsigned int i = 0; i <
NextCount(); i++) {
841 if (seg->
Parent() != trk) to_check.push_back(seg->
Parent());
847 for (
size_t t = 0; t < to_check.size(); t++)
864 for (
size_t t = 0; (t < to_check.size()) && !found; t++)
865 for (
size_t i = 0; i < to_check[t]->size(); i++) {
size_t NPrecalcEnabledHits(void) const
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
TVector3 const & Point3D() const
TVector2 const & Projection2D(unsigned int view) const
double Length2() const override
std::vector< TVector3 * > fAssignedPoints
double Penalty(float endSegWeight) const
Implementation of the Projection Matching Algorithm.
TVector2 const & Point2D() const noexcept
double Dist2(const TVector2 &v1, const TVector2 &v2)
unsigned int Nplanes() const
Number of planes in this tpc.
void ClearAssigned(pma::Track3D *trk=0) override
Implementation of the Projection Matching Algorithm.
pma::SortedObjectBase * prev
double GetDistance2To(const TVector3 &p3d) const override
Distance [cm] from the 3D point to the point 3D.
pma::Vector3D GetDirection3D() const override
bool LimitPoint3D()
Returns true if node position was trimmed to its TPC volume + fMargin.
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
virtual unsigned int NextCount(void) const
std::vector< pma::Track3D * > GetBranches() const
void Optimize(float penaltyValue, float endSegWeight)
Planes which measure Z direction.
virtual pma::SortedObjectBase * Next(unsigned int=0) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double Length2(void) const override
bool IsTPCEdge() const
Is the first/last in this TPC?
recob::tracking::Vector_t Vector3D
static bool fGradFixed[3]
virtual pma::Vector3D GetDirection3D(void) const =0
Get 3D direction cosines corresponding to this element.
double SegmentCosWirePlane() const
int TPC(void) const
TPC index or -1 if out of any TPC.
double PenaltyInWirePlane() const
double GetDistToWall() const
virtual unsigned int NextCount(void) const
double PiInWirePlane() const
double SegmentCos() const
Cosine of 3D angle between connected segments.
bool SetPoint3D(const TVector3 &p3d)
double MinZ() const
Returns the world z coordinate of the start of the box.
double StepWithGradient(float alfa, float tol, float penalty, float weight)
double SegmentCosTransverse() const
unsigned int NumberTimeSamples() const
double EndPtCos2Transverse() const
bool SameTPC(const TVector3 &p3d, float margin=0.0F) const
Check if p3d is in the same TPC as the node.
void SetProjection(const TVector2 &p, float b)
double SumDist2(void) const
Detector simulation of raw signals on wires.
double MaxY() const
Returns the world y coordinate of the end of the box.
std::vector< pma::Hit3D * > fAssignedHits
double ConvertTicksToX(double ticks, int p, int t, int c) const
double MakeGradient(float penaltyValue, float endSegWeight)
unsigned int View2D() const noexcept
Implementation of the Projection Matching Algorithm.
virtual pma::SortedObjectBase * Next(unsigned int index=0) const
Encapsulate the construction of a single detector plane.
double GetObjFunction(float penaltyValue, float endSegWeight) const
Objective function minimized during oprimization.
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)
ROOT::Math::DisplacementVector2D< ROOT::Math::Cartesian2D< double > > Vector2D
bool IsBranching() const
Belongs to more than one track?
pma::SortedObjectBase * next
double SumDist2Hits() const override
geo::TPCGeo const & fTpcGeo
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.
virtual pma::SortedObjectBase * Prev(void) const
Namespace collecting geometry-related classes utilities.
void SetProjection(pma::Hit3D &h) const override
Set hit 3D position and its 2D projection to the vertex.
double MinY() const
Returns the world y coordinate of the start of the box.
pma::Track3D * Parent(void) const
double Pi(float endSegWeight, bool doAsymm) const
double PlaneCoordinate(geo::Point_t const &point) const
Returns the coordinate of the point on the plane.
double Length(void) const
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
static float OptFactor(unsigned int view)
Encapsulate the construction of a single detector plane.