From b79b863d16e57d0a5d32584ddddff0d090670d66 Mon Sep 17 00:00:00 2001
From: Christos Anastopoulos <christos.anastopoulos@cern.ch>
Date: Tue, 3 Nov 2020 16:14:41 +0000
Subject: [PATCH] Surfaces : Move inline methods to separate .icc files

---
 .../TrkSurfaces/TrkSurfaces/TrapezoidBounds.h | 136 +--------------
 .../TrkSurfaces/TrapezoidBounds.icc           | 159 ++++++++++++++++++
 .../TrkSurfaces/TrkSurfaces/TriangleBounds.h  | 122 +-------------
 .../TrkSurfaces/TriangleBounds.icc            | 149 ++++++++++++++++
 4 files changed, 310 insertions(+), 256 deletions(-)
 create mode 100644 Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TrapezoidBounds.icc
 create mode 100644 Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TriangleBounds.icc

diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TrapezoidBounds.h b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TrapezoidBounds.h
index e6b530c1439..cf95ee03b19 100644
--- a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TrapezoidBounds.h
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TrapezoidBounds.h
@@ -172,140 +172,6 @@ private:
   TDD_real_t m_beta;
 };
 
-inline TrapezoidBounds*
-TrapezoidBounds::clone() const
-{
-  return new TrapezoidBounds(*this);
-}
-
-inline double
-TrapezoidBounds::minHalflengthX() const
-{
-  return m_boundValues[TrapezoidBounds::bv_minHalfX];
-}
-
-inline double
-TrapezoidBounds::maxHalflengthX() const
-{
-  return m_boundValues[TrapezoidBounds::bv_maxHalfX];
-}
-
-inline double
-TrapezoidBounds::halflengthY() const
-{
-  return m_boundValues[TrapezoidBounds::bv_halfY];
-}
-
-inline double
-TrapezoidBounds::minHalflengthPhi() const
-{
-  return minHalflengthX();
-}
-
-inline double
-TrapezoidBounds::maxHalflengthPhi() const
-{
-  return maxHalflengthX();
-}
-
-inline double
-TrapezoidBounds::halflengthEta() const
-{
-  return halflengthY();
-}
-
-inline double
-TrapezoidBounds::alpha() const
-{
-  return m_alpha;
-}
-
-inline double
-TrapezoidBounds::beta() const
-{
-  return m_beta;
-}
-
-inline double
-TrapezoidBounds::r() const
-{
-  return sqrt(m_boundValues[TrapezoidBounds::bv_maxHalfX] * m_boundValues[TrapezoidBounds::bv_maxHalfX] +
-              m_boundValues[TrapezoidBounds::bv_halfY] * m_boundValues[TrapezoidBounds::bv_halfY]);
-}
-
-inline bool
-TrapezoidBounds::inside(const Amg::Vector2D& locpo, const BoundaryCheck& bchk) const
-{
-  if (bchk.bcType == 0)
-    return TrapezoidBounds::inside(locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
-
-  // a fast FALSE
-  double fabsY = fabs(locpo[Trk::locY]);
-  double max_ell = bchk.lCovariance(0, 0) > bchk.lCovariance(1, 1) ? bchk.lCovariance(0, 0) : bchk.lCovariance(1, 1);
-  double limit = bchk.nSigmas * sqrt(max_ell);
-  if (fabsY > (m_boundValues[TrapezoidBounds::bv_halfY] + limit))
-    return false;
-  // a fast FALSE
-  double fabsX = fabs(locpo[Trk::locX]);
-  if (fabsX > (m_boundValues[TrapezoidBounds::bv_maxHalfX] + limit))
-    return false;
-  // a fast TRUE
-  double min_ell = bchk.lCovariance(0, 0) < bchk.lCovariance(1, 1) ? bchk.lCovariance(0, 0) : bchk.lCovariance(1, 1);
-  limit = bchk.nSigmas * sqrt(min_ell);
-  if (fabsX < (m_boundValues[TrapezoidBounds::bv_minHalfX] + limit) &&
-      fabsY < (m_boundValues[TrapezoidBounds::bv_halfY] + limit))
-    return true;
-
-  // compute KDOP and axes for surface polygon
-  std::vector<KDOP> elementKDOP(3);
-  std::vector<Amg::Vector2D> elementP(4);
-  float theta = (bchk.lCovariance(1, 0) != 0 && (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)) != 0)
-                  ? .5 * bchk.FastArcTan(2 * bchk.lCovariance(1, 0) / (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)))
-                  : 0.;
-  sincosCache scResult = bchk.FastSinCos(theta);
-  AmgMatrix(2, 2) rotMatrix;
-  rotMatrix << scResult.cosC, scResult.sinC, -scResult.sinC, scResult.cosC;
-  AmgMatrix(2, 2) normal;
-  normal << 0, -1, 1, 0;
-  // ellipse is always at (0,0), surface is moved to ellipse position and then rotated
-  Amg::Vector2D p;
-  p << m_boundValues[TrapezoidBounds::bv_minHalfX], -m_boundValues[TrapezoidBounds::bv_halfY];
-  elementP[0] = (rotMatrix * (p - locpo));
-  p << -m_boundValues[TrapezoidBounds::bv_minHalfX], -m_boundValues[TrapezoidBounds::bv_halfY];
-  elementP[1] = (rotMatrix * (p - locpo));
-  scResult = bchk.FastSinCos(m_beta);
-  p << m_boundValues[TrapezoidBounds::bv_minHalfX] +
-         (2. * m_boundValues[TrapezoidBounds::bv_halfY]) * (scResult.sinC / scResult.cosC),
-    m_boundValues[TrapezoidBounds::bv_halfY];
-  elementP[2] = (rotMatrix * (p - locpo));
-  scResult = bchk.FastSinCos(m_alpha);
-  p << -(m_boundValues[TrapezoidBounds::bv_minHalfX] +
-         (2. * m_boundValues[TrapezoidBounds::bv_halfY]) * (scResult.sinC / scResult.cosC)),
-    m_boundValues[TrapezoidBounds::bv_halfY];
-  elementP[3] = (rotMatrix * (p - locpo));
-  std::vector<Amg::Vector2D> axis = { normal * (elementP[1] - elementP[0]),
-                                      normal * (elementP[3] - elementP[1]),
-                                      normal * (elementP[2] - elementP[0]) };
-  bchk.ComputeKDOP(elementP, axis, elementKDOP);
-  // compute KDOP for error ellipse
-  std::vector<KDOP> errelipseKDOP(3);
-  bchk.ComputeKDOP(bchk.EllipseToPoly(3), axis, errelipseKDOP);
-  // check if KDOPs overlap and return result
-  return bchk.TestKDOPKDOP(elementKDOP, errelipseKDOP);
-}
-
-inline bool
-TrapezoidBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
-{
-  return (fabs(locpo[locX]) < m_boundValues[TrapezoidBounds::bv_maxHalfX] + tol1);
-}
-
-inline bool
-TrapezoidBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
-{
-  return (fabs(locpo[locY]) < m_boundValues[TrapezoidBounds::bv_halfY] + tol2);
-}
-
 } // end of namespace
-
+#include "TrkSurfaces/TrapezoidBounds.icc"
 #endif // TRKSURFACES_TRAPEZOIDBOUNDS_H
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TrapezoidBounds.icc b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TrapezoidBounds.icc
new file mode 100644
index 00000000000..5b800e5f4e8
--- /dev/null
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TrapezoidBounds.icc
@@ -0,0 +1,159 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+namespace Trk {
+
+inline TrapezoidBounds*
+TrapezoidBounds::clone() const
+{
+  return new TrapezoidBounds(*this);
+}
+
+inline double
+TrapezoidBounds::minHalflengthX() const
+{
+  return m_boundValues[TrapezoidBounds::bv_minHalfX];
+}
+
+inline double
+TrapezoidBounds::maxHalflengthX() const
+{
+  return m_boundValues[TrapezoidBounds::bv_maxHalfX];
+}
+
+inline double
+TrapezoidBounds::halflengthY() const
+{
+  return m_boundValues[TrapezoidBounds::bv_halfY];
+}
+
+inline double
+TrapezoidBounds::minHalflengthPhi() const
+{
+  return minHalflengthX();
+}
+
+inline double
+TrapezoidBounds::maxHalflengthPhi() const
+{
+  return maxHalflengthX();
+}
+
+inline double
+TrapezoidBounds::halflengthEta() const
+{
+  return halflengthY();
+}
+
+inline double
+TrapezoidBounds::alpha() const
+{
+  return m_alpha;
+}
+
+inline double
+TrapezoidBounds::beta() const
+{
+  return m_beta;
+}
+
+inline double
+TrapezoidBounds::r() const
+{
+  return sqrt(m_boundValues[TrapezoidBounds::bv_maxHalfX] *
+                m_boundValues[TrapezoidBounds::bv_maxHalfX] +
+              m_boundValues[TrapezoidBounds::bv_halfY] *
+                m_boundValues[TrapezoidBounds::bv_halfY]);
+}
+
+inline bool
+TrapezoidBounds::inside(const Amg::Vector2D& locpo,
+                        const BoundaryCheck& bchk) const
+{
+  if (bchk.bcType == 0)
+    return TrapezoidBounds::inside(
+      locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
+
+  // a fast FALSE
+  double fabsY = fabs(locpo[Trk::locY]);
+  double max_ell = bchk.lCovariance(0, 0) > bchk.lCovariance(1, 1)
+                     ? bchk.lCovariance(0, 0)
+                     : bchk.lCovariance(1, 1);
+  double limit = bchk.nSigmas * sqrt(max_ell);
+  if (fabsY > (m_boundValues[TrapezoidBounds::bv_halfY] + limit))
+    return false;
+  // a fast FALSE
+  double fabsX = fabs(locpo[Trk::locX]);
+  if (fabsX > (m_boundValues[TrapezoidBounds::bv_maxHalfX] + limit))
+    return false;
+  // a fast TRUE
+  double min_ell = bchk.lCovariance(0, 0) < bchk.lCovariance(1, 1)
+                     ? bchk.lCovariance(0, 0)
+                     : bchk.lCovariance(1, 1);
+  limit = bchk.nSigmas * sqrt(min_ell);
+  if (fabsX < (m_boundValues[TrapezoidBounds::bv_minHalfX] + limit) &&
+      fabsY < (m_boundValues[TrapezoidBounds::bv_halfY] + limit))
+    return true;
+
+  // compute KDOP and axes for surface polygon
+  std::vector<KDOP> elementKDOP(3);
+  std::vector<Amg::Vector2D> elementP(4);
+  float theta =
+    (bchk.lCovariance(1, 0) != 0 &&
+     (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)) != 0)
+      ? .5 * bchk.FastArcTan(2 * bchk.lCovariance(1, 0) /
+                             (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)))
+      : 0.;
+  sincosCache scResult = bchk.FastSinCos(theta);
+  AmgMatrix(2, 2) rotMatrix;
+  rotMatrix << scResult.cosC, scResult.sinC, -scResult.sinC, scResult.cosC;
+  AmgMatrix(2, 2) normal;
+  normal << 0, -1, 1, 0;
+  // ellipse is always at (0,0), surface is moved to ellipse position and then
+  // rotated
+  Amg::Vector2D p;
+  p << m_boundValues[TrapezoidBounds::bv_minHalfX],
+    -m_boundValues[TrapezoidBounds::bv_halfY];
+  elementP[0] = (rotMatrix * (p - locpo));
+  p << -m_boundValues[TrapezoidBounds::bv_minHalfX],
+    -m_boundValues[TrapezoidBounds::bv_halfY];
+  elementP[1] = (rotMatrix * (p - locpo));
+  scResult = bchk.FastSinCos(m_beta);
+  p << m_boundValues[TrapezoidBounds::bv_minHalfX] +
+         (2. * m_boundValues[TrapezoidBounds::bv_halfY]) *
+           (scResult.sinC / scResult.cosC),
+    m_boundValues[TrapezoidBounds::bv_halfY];
+  elementP[2] = (rotMatrix * (p - locpo));
+  scResult = bchk.FastSinCos(m_alpha);
+  p << -(m_boundValues[TrapezoidBounds::bv_minHalfX] +
+         (2. * m_boundValues[TrapezoidBounds::bv_halfY]) *
+           (scResult.sinC / scResult.cosC)),
+    m_boundValues[TrapezoidBounds::bv_halfY];
+  elementP[3] = (rotMatrix * (p - locpo));
+  std::vector<Amg::Vector2D> axis = { normal * (elementP[1] - elementP[0]),
+                                      normal * (elementP[3] - elementP[1]),
+                                      normal * (elementP[2] - elementP[0]) };
+  bchk.ComputeKDOP(elementP, axis, elementKDOP);
+  // compute KDOP for error ellipse
+  std::vector<KDOP> errelipseKDOP(3);
+  bchk.ComputeKDOP(bchk.EllipseToPoly(3), axis, errelipseKDOP);
+  // check if KDOPs overlap and return result
+  return bchk.TestKDOPKDOP(elementKDOP, errelipseKDOP);
+}
+
+inline bool
+TrapezoidBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
+{
+  return (fabs(locpo[locX]) <
+          m_boundValues[TrapezoidBounds::bv_maxHalfX] + tol1);
+}
+
+inline bool
+TrapezoidBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
+{
+  return (fabs(locpo[locY]) < m_boundValues[TrapezoidBounds::bv_halfY] + tol2);
+}
+
+} // end of namespace
+
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TriangleBounds.h b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TriangleBounds.h
index 9117cfc67b8..c7e0858df56 100644
--- a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TriangleBounds.h
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TriangleBounds.h
@@ -111,127 +111,7 @@ public:
 private:
   std::vector<TDD_real_t> m_boundValues;
 };
-
-inline TriangleBounds*
-TriangleBounds::clone() const
-{
-  return new TriangleBounds(*this);
-}
-
-inline bool
-TriangleBounds::inside(const Amg::Vector2D& locpo, double tol1, double tol2) const
-{
-  std::pair<double, double> locB(m_boundValues[TriangleBounds::bv_x2] - m_boundValues[TriangleBounds::bv_x1],
-                                 m_boundValues[TriangleBounds::bv_y2] - m_boundValues[TriangleBounds::bv_y1]);
-  std::pair<double, double> locT(m_boundValues[TriangleBounds::bv_x3] - locpo[0],
-                                 m_boundValues[TriangleBounds::bv_y3] - locpo[1]);
-  std::pair<double, double> locV(m_boundValues[TriangleBounds::bv_x1] - locpo[0],
-                                 m_boundValues[TriangleBounds::bv_y1] - locpo[1]);
-
-  // special case :: third vertex ?
-  if (locT.first * locT.first + locT.second * locT.second < tol1 * tol1)
-    return true;
-
-  // special case : lies on base ?
-  double db = locB.first * locV.second - locB.second * locV.first;
-  if (fabs(db) < tol1) {
-    double a = (locB.first != 0) ? -locV.first / locB.first : -locV.second / locB.second;
-    return a > -tol2 && a - 1. < tol2;
-  }
-
-  double dn = locB.first * locT.second - locB.second * locT.first;
-
-  if (fabs(dn) > fabs(tol1)) {
-    double t = (locB.first * locV.second - locB.second * locV.first) / dn;
-    if (t > 0.)
-      return false;
-
-    double a =
-      (locB.first != 0.) ? (t * locT.first - locV.first) / locB.first : (t * locT.second - locV.second) / locB.second;
-    if (a < -tol2 || a - 1. > tol2)
-      return false;
-  } else {
-    return false;
-  }
-  return true;
-}
-
-inline bool
-TriangleBounds::inside(const Amg::Vector2D& locpo, const BoundaryCheck& bchk) const
-{
-  if (bchk.bcType == 0)
-    return TriangleBounds::inside(locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
-
-  // a fast FALSE
-  double fabsR = sqrt(locpo[Trk::locX] * locpo[Trk::locX] + locpo[Trk::locY] * locpo[Trk::locY]);
-  double max_ell = bchk.lCovariance(0, 0) > bchk.lCovariance(1, 1) ? bchk.lCovariance(0, 0) : bchk.lCovariance(1, 1);
-  double limit = bchk.nSigmas * sqrt(max_ell);
-  double r_max = TriangleBounds::r();
-  if (fabsR > (r_max + limit))
-    return false;
-
-  // compute KDOP and axes for surface polygon
-  std::vector<KDOP> elementKDOP(3);
-  std::vector<Amg::Vector2D> elementP(3);
-  float theta = (bchk.lCovariance(1, 0) != 0 && (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)) != 0)
-                  ? .5 * bchk.FastArcTan(2 * bchk.lCovariance(1, 0) / (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)))
-                  : 0.;
-  sincosCache scResult = bchk.FastSinCos(theta);
-  AmgMatrix(2, 2) rotMatrix;
-  rotMatrix << scResult.cosC, scResult.sinC, -scResult.sinC, scResult.cosC;
-  AmgMatrix(2, 2) normal;
-  normal << 0, -1, 1, 0;
-  // ellipse is always at (0,0), surface is moved to ellipse position and then rotated
-  Amg::Vector2D p;
-  p << m_boundValues[TriangleBounds::bv_x1], m_boundValues[TriangleBounds::bv_y1];
-  elementP[0] = (rotMatrix * (p - locpo));
-  p << m_boundValues[TriangleBounds::bv_x2], m_boundValues[TriangleBounds::bv_y2];
-  elementP[1] = (rotMatrix * (p - locpo));
-  p << m_boundValues[TriangleBounds::bv_x3], m_boundValues[TriangleBounds::bv_y3];
-  elementP[2] = (rotMatrix * (p - locpo));
-  std::vector<Amg::Vector2D> axis = { normal * (elementP[1] - elementP[0]),
-                                      normal * (elementP[2] - elementP[1]),
-                                      normal * (elementP[2] - elementP[0]) };
-  bchk.ComputeKDOP(elementP, axis, elementKDOP);
-  // compute KDOP for error ellipse
-  std::vector<KDOP> errelipseKDOP(3);
-  bchk.ComputeKDOP(bchk.EllipseToPoly(3), axis, errelipseKDOP);
-  // check if KDOPs overlap and return result
-  return bchk.TestKDOPKDOP(elementKDOP, errelipseKDOP);
-}
-
-inline bool
-TriangleBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
-{
-  return inside(locpo, tol1, tol1);
-}
-
-inline bool
-TriangleBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
-{
-  return inside(locpo, tol2, tol2);
-}
-
-inline std::vector<std::pair<TDD_real_t, TDD_real_t>>
-TriangleBounds::vertices() const
-{
-  std::vector<std::pair<TDD_real_t, TDD_real_t>> vertices;
-  vertices.resize(3);
-  for (size_t iv = 0; iv < 3; iv++)
-    vertices.emplace_back(m_boundValues[2 * iv], m_boundValues[2 * iv + 1]);
-  return vertices;
-}
-
-inline double
-TriangleBounds::r() const
-{
-  double rmax = 0.;
-  for (size_t iv = 0; iv < 3; iv++)
-    rmax =
-      fmax(rmax, m_boundValues[2 * iv] * m_boundValues[2 * iv] + m_boundValues[2 * iv + 1] * m_boundValues[2 * iv + 1]);
-  return sqrt(rmax);
-}
-
 } // end of namespace
 
+#include "TrkSurfaces/TriangleBounds.icc"
 #endif // TRKSURFACES_RECTANGLEBOUNDS_H
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TriangleBounds.icc b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TriangleBounds.icc
new file mode 100644
index 00000000000..acb149216c8
--- /dev/null
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/TriangleBounds.icc
@@ -0,0 +1,149 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+namespace Trk {
+
+inline TriangleBounds*
+TriangleBounds::clone() const
+{
+  return new TriangleBounds(*this);
+}
+
+inline bool
+TriangleBounds::inside(const Amg::Vector2D& locpo,
+                       double tol1,
+                       double tol2) const
+{
+  std::pair<double, double> locB(m_boundValues[TriangleBounds::bv_x2] -
+                                   m_boundValues[TriangleBounds::bv_x1],
+                                 m_boundValues[TriangleBounds::bv_y2] -
+                                   m_boundValues[TriangleBounds::bv_y1]);
+  std::pair<double, double> locT(
+    m_boundValues[TriangleBounds::bv_x3] - locpo[0],
+    m_boundValues[TriangleBounds::bv_y3] - locpo[1]);
+  std::pair<double, double> locV(
+    m_boundValues[TriangleBounds::bv_x1] - locpo[0],
+    m_boundValues[TriangleBounds::bv_y1] - locpo[1]);
+
+  // special case :: third vertex ?
+  if (locT.first * locT.first + locT.second * locT.second < tol1 * tol1)
+    return true;
+
+  // special case : lies on base ?
+  double db = locB.first * locV.second - locB.second * locV.first;
+  if (fabs(db) < tol1) {
+    double a =
+      (locB.first != 0) ? -locV.first / locB.first : -locV.second / locB.second;
+    return a > -tol2 && a - 1. < tol2;
+  }
+
+  double dn = locB.first * locT.second - locB.second * locT.first;
+
+  if (fabs(dn) > fabs(tol1)) {
+    double t = (locB.first * locV.second - locB.second * locV.first) / dn;
+    if (t > 0.)
+      return false;
+
+    double a = (locB.first != 0.)
+                 ? (t * locT.first - locV.first) / locB.first
+                 : (t * locT.second - locV.second) / locB.second;
+    if (a < -tol2 || a - 1. > tol2)
+      return false;
+  } else {
+    return false;
+  }
+  return true;
+}
+
+inline bool
+TriangleBounds::inside(const Amg::Vector2D& locpo,
+                       const BoundaryCheck& bchk) const
+{
+  if (bchk.bcType == 0)
+    return TriangleBounds::inside(
+      locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
+
+  // a fast FALSE
+  double fabsR = sqrt(locpo[Trk::locX] * locpo[Trk::locX] +
+                      locpo[Trk::locY] * locpo[Trk::locY]);
+  double max_ell = bchk.lCovariance(0, 0) > bchk.lCovariance(1, 1)
+                     ? bchk.lCovariance(0, 0)
+                     : bchk.lCovariance(1, 1);
+  double limit = bchk.nSigmas * sqrt(max_ell);
+  double r_max = TriangleBounds::r();
+  if (fabsR > (r_max + limit))
+    return false;
+
+  // compute KDOP and axes for surface polygon
+  std::vector<KDOP> elementKDOP(3);
+  std::vector<Amg::Vector2D> elementP(3);
+  float theta =
+    (bchk.lCovariance(1, 0) != 0 &&
+     (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)) != 0)
+      ? .5 * bchk.FastArcTan(2 * bchk.lCovariance(1, 0) /
+                             (bchk.lCovariance(1, 1) - bchk.lCovariance(0, 0)))
+      : 0.;
+  sincosCache scResult = bchk.FastSinCos(theta);
+  AmgMatrix(2, 2) rotMatrix;
+  rotMatrix << scResult.cosC, scResult.sinC, -scResult.sinC, scResult.cosC;
+  AmgMatrix(2, 2) normal;
+  normal << 0, -1, 1, 0;
+  // ellipse is always at (0,0), surface is moved to ellipse position and then
+  // rotated
+  Amg::Vector2D p;
+  p << m_boundValues[TriangleBounds::bv_x1],
+    m_boundValues[TriangleBounds::bv_y1];
+  elementP[0] = (rotMatrix * (p - locpo));
+  p << m_boundValues[TriangleBounds::bv_x2],
+    m_boundValues[TriangleBounds::bv_y2];
+  elementP[1] = (rotMatrix * (p - locpo));
+  p << m_boundValues[TriangleBounds::bv_x3],
+    m_boundValues[TriangleBounds::bv_y3];
+  elementP[2] = (rotMatrix * (p - locpo));
+  std::vector<Amg::Vector2D> axis = { normal * (elementP[1] - elementP[0]),
+                                      normal * (elementP[2] - elementP[1]),
+                                      normal * (elementP[2] - elementP[0]) };
+  bchk.ComputeKDOP(elementP, axis, elementKDOP);
+  // compute KDOP for error ellipse
+  std::vector<KDOP> errelipseKDOP(3);
+  bchk.ComputeKDOP(bchk.EllipseToPoly(3), axis, errelipseKDOP);
+  // check if KDOPs overlap and return result
+  return bchk.TestKDOPKDOP(elementKDOP, errelipseKDOP);
+}
+
+inline bool
+TriangleBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
+{
+  return inside(locpo, tol1, tol1);
+}
+
+inline bool
+TriangleBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
+{
+  return inside(locpo, tol2, tol2);
+}
+
+inline std::vector<std::pair<TDD_real_t, TDD_real_t>>
+TriangleBounds::vertices() const
+{
+  std::vector<std::pair<TDD_real_t, TDD_real_t>> vertices;
+  vertices.resize(3);
+  for (size_t iv = 0; iv < 3; iv++)
+    vertices.emplace_back(m_boundValues[2 * iv], m_boundValues[2 * iv + 1]);
+  return vertices;
+}
+
+inline double
+TriangleBounds::r() const
+{
+  double rmax = 0.;
+  for (size_t iv = 0; iv < 3; iv++)
+    rmax = fmax(rmax,
+                m_boundValues[2 * iv] * m_boundValues[2 * iv] +
+                  m_boundValues[2 * iv + 1] * m_boundValues[2 * iv + 1]);
+  return sqrt(rmax);
+}
+
+} // end of namespace
+
-- 
GitLab