diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/CylinderBounds.h b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/CylinderBounds.h
index 124e8171a3dbb6a221633e13d65002bab043b8ea..517ef5816aa532c4122b93964001049be0235815 100644
--- a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/CylinderBounds.h
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/CylinderBounds.h
@@ -140,133 +140,6 @@ private:
   std::vector<TDD_real_t> m_boundValues;
   bool m_checkPhi;
 };
-
-inline CylinderBounds*
-CylinderBounds::clone() const
-{
-  return new CylinderBounds(*this);
-}
-
-inline bool
-CylinderBounds::inside(const Amg::Vector2D& locpo, double tol1, double tol2) const
-{
-  double z = locpo[locZ];
-  bool insideZ = insideLocZ(z, tol2);
-  if (!insideZ)
-    return false;
-  // no check on Phi neccesary
-  if (!m_checkPhi)
-    return true;
-  // now check insidePhi
-  double localPhi =
-    (locpo[locRPhi] / m_boundValues[CylinderBounds::bv_radius]) - m_boundValues[CylinderBounds::bv_averagePhi];
-  localPhi -= (localPhi > M_PI) ? 2. * M_PI : 0.;
-  return (localPhi * localPhi < (m_boundValues[CylinderBounds::bv_halfPhiSector] + tol1) *
-                                  (m_boundValues[CylinderBounds::bv_halfPhiSector] + tol1));
-}
-
-inline bool
-CylinderBounds::inside(const Amg::Vector2D& locpo, const BoundaryCheck& bchk) const
-{
-  if (bchk.bcType == 0 || bchk.nSigmas == 0 || m_boundValues[CylinderBounds::bv_halfPhiSector] != M_PI)
-    return CylinderBounds::inside(locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
-
-  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);
-  double dphi = scResult.sinC * scResult.sinC * bchk.lCovariance(0, 0);
-  double dz = scResult.cosC * scResult.cosC * bchk.lCovariance(0, 1);
-  double max_ell = dphi > dz ? dphi : dz;
-  double limit = bchk.nSigmas * sqrt(max_ell);
-  return insideLocZ(locpo[locZ], limit);
-}
-
-inline bool
-CylinderBounds::inside3D(const Amg::Vector3D& glopo, double tol1, double tol2) const
-{
-  return inside(glopo.perp(), glopo.phi(), glopo.z(), tol1, tol2);
-}
-
-inline bool
-CylinderBounds::inside(double r, double phi, double z, double tol1, double tol2) const
-{
-
-  bool insideZ = insideLocZ(z, tol2);
-  if (!insideZ)
-    return false;
-  double diffR = (m_boundValues[CylinderBounds::bv_radius] - r);
-  bool insideR = diffR * diffR < tol1 * tol1;
-  if (!insideR)
-    return false;
-  // now check insidePhi if needed
-  if (!m_checkPhi)
-    return true;
-  // phi needs to be checked
-  double localPhi = phi - m_boundValues[CylinderBounds::bv_averagePhi];
-  localPhi -= (localPhi > M_PI) ? 2. * M_PI : 0.;
-  return (localPhi * localPhi <
-          m_boundValues[CylinderBounds::bv_halfPhiSector] * m_boundValues[CylinderBounds::bv_halfPhiSector]);
-}
-
-inline bool
-CylinderBounds::insideLocZ(double z, double tol2) const
-{
-  return (m_boundValues[CylinderBounds::bv_halfZ] + tol2) - fabs(z) > 0.;
-}
-
-inline bool
-CylinderBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
-{
-  bool insideRphi = false;
-  if (fabs(m_boundValues[CylinderBounds::bv_averagePhi]) < 10e-7)
-    insideRphi = (fabs(locpo[locRPhi] / m_boundValues[CylinderBounds::bv_radius]) <
-                  (m_boundValues[CylinderBounds::bv_halfPhiSector] + tol1));
-  else {
-    double localPhi =
-      (locpo[locRPhi] / m_boundValues[CylinderBounds::bv_radius]) - m_boundValues[CylinderBounds::bv_averagePhi];
-    localPhi -= (localPhi > M_PI) ? 2. * M_PI : 0.;
-    insideRphi = (localPhi < (m_boundValues[CylinderBounds::bv_halfPhiSector] + tol1));
-  }
-  return (insideRphi);
 }
-
-inline bool
-CylinderBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
-{
-  return insideLocZ(locpo[locZ], tol2);
-}
-
-inline bool
-CylinderBounds::insideRadius(const Amg::Vector2D& locpo, double tol) const
-{
-  return (this->inside(locpo, tol, 0) && fabs(locpo[locR]) < m_boundValues[CylinderBounds::bv_radius] + tol);
-}
-
-inline double
-CylinderBounds::r() const
-{
-  return m_boundValues[CylinderBounds::bv_radius];
-}
-
-inline double
-CylinderBounds::averagePhi() const
-{
-  return m_boundValues[CylinderBounds::bv_averagePhi];
-}
-
-inline double
-CylinderBounds::halfPhiSector() const
-{
-  return m_boundValues[CylinderBounds::bv_halfPhiSector];
-}
-
-inline double
-CylinderBounds::halflengthZ() const
-{
-  return m_boundValues[CylinderBounds::bv_halfZ];
-}
-
-}
-
+#include "TrkSurfaces/CylinderBounds.icc"
 #endif // TRKSURFACES_CYLINDERBOUNDS_H
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/CylinderBounds.icc b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/CylinderBounds.icc
new file mode 100644
index 0000000000000000000000000000000000000000..4a94a2b71a341779783b6376983cd42cf82f8f12
--- /dev/null
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/CylinderBounds.icc
@@ -0,0 +1,151 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+namespace Trk {
+inline CylinderBounds*
+CylinderBounds::clone() const
+{
+  return new CylinderBounds(*this);
+}
+
+inline bool
+CylinderBounds::inside(const Amg::Vector2D& locpo,
+                       double tol1,
+                       double tol2) const
+{
+  double z = locpo[locZ];
+  bool insideZ = insideLocZ(z, tol2);
+  if (!insideZ)
+    return false;
+  // no check on Phi neccesary
+  if (!m_checkPhi)
+    return true;
+  // now check insidePhi
+  double localPhi =
+    (locpo[locRPhi] / m_boundValues[CylinderBounds::bv_radius]) -
+    m_boundValues[CylinderBounds::bv_averagePhi];
+  localPhi -= (localPhi > M_PI) ? 2. * M_PI : 0.;
+  return (localPhi * localPhi <
+          (m_boundValues[CylinderBounds::bv_halfPhiSector] + tol1) *
+            (m_boundValues[CylinderBounds::bv_halfPhiSector] + tol1));
+}
+
+inline bool
+CylinderBounds::inside(const Amg::Vector2D& locpo,
+                       const BoundaryCheck& bchk) const
+{
+  if (bchk.bcType == 0 || bchk.nSigmas == 0 ||
+      m_boundValues[CylinderBounds::bv_halfPhiSector] != M_PI)
+    return CylinderBounds::inside(
+      locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
+
+  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);
+  double dphi = scResult.sinC * scResult.sinC * bchk.lCovariance(0, 0);
+  double dz = scResult.cosC * scResult.cosC * bchk.lCovariance(0, 1);
+  double max_ell = dphi > dz ? dphi : dz;
+  double limit = bchk.nSigmas * sqrt(max_ell);
+  return insideLocZ(locpo[locZ], limit);
+}
+
+inline bool
+CylinderBounds::inside3D(const Amg::Vector3D& glopo,
+                         double tol1,
+                         double tol2) const
+{
+  return inside(glopo.perp(), glopo.phi(), glopo.z(), tol1, tol2);
+}
+
+inline bool
+CylinderBounds::inside(double r, double phi, double z, double tol1, double tol2)
+  const
+{
+
+  bool insideZ = insideLocZ(z, tol2);
+  if (!insideZ)
+    return false;
+  double diffR = (m_boundValues[CylinderBounds::bv_radius] - r);
+  bool insideR = diffR * diffR < tol1 * tol1;
+  if (!insideR)
+    return false;
+  // now check insidePhi if needed
+  if (!m_checkPhi)
+    return true;
+  // phi needs to be checked
+  double localPhi = phi - m_boundValues[CylinderBounds::bv_averagePhi];
+  localPhi -= (localPhi > M_PI) ? 2. * M_PI : 0.;
+  return (localPhi * localPhi <
+          m_boundValues[CylinderBounds::bv_halfPhiSector] *
+            m_boundValues[CylinderBounds::bv_halfPhiSector]);
+}
+
+inline bool
+CylinderBounds::insideLocZ(double z, double tol2) const
+{
+  return (m_boundValues[CylinderBounds::bv_halfZ] + tol2) - fabs(z) > 0.;
+}
+
+inline bool
+CylinderBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
+{
+  bool insideRphi = false;
+  if (fabs(m_boundValues[CylinderBounds::bv_averagePhi]) < 10e-7)
+    insideRphi =
+      (fabs(locpo[locRPhi] / m_boundValues[CylinderBounds::bv_radius]) <
+       (m_boundValues[CylinderBounds::bv_halfPhiSector] + tol1));
+  else {
+    double localPhi =
+      (locpo[locRPhi] / m_boundValues[CylinderBounds::bv_radius]) -
+      m_boundValues[CylinderBounds::bv_averagePhi];
+    localPhi -= (localPhi > M_PI) ? 2. * M_PI : 0.;
+    insideRphi =
+      (localPhi < (m_boundValues[CylinderBounds::bv_halfPhiSector] + tol1));
+  }
+  return (insideRphi);
+}
+
+inline bool
+CylinderBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
+{
+  return insideLocZ(locpo[locZ], tol2);
+}
+
+inline bool
+CylinderBounds::insideRadius(const Amg::Vector2D& locpo, double tol) const
+{
+  return (this->inside(locpo, tol, 0) &&
+          fabs(locpo[locR]) < m_boundValues[CylinderBounds::bv_radius] + tol);
+}
+
+inline double
+CylinderBounds::r() const
+{
+  return m_boundValues[CylinderBounds::bv_radius];
+}
+
+inline double
+CylinderBounds::averagePhi() const
+{
+  return m_boundValues[CylinderBounds::bv_averagePhi];
+}
+
+inline double
+CylinderBounds::halfPhiSector() const
+{
+  return m_boundValues[CylinderBounds::bv_halfPhiSector];
+}
+
+inline double
+CylinderBounds::halflengthZ() const
+{
+  return m_boundValues[CylinderBounds::bv_halfZ];
+}
+
+}
+
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/DiamondBounds.h b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/DiamondBounds.h
index 2c375d619997ce689cf9d546f36b2e4525fb5af6..45017b57137d14c3e1979bdebc8d32ed99becc30 100644
--- a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/DiamondBounds.h
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/DiamondBounds.h
@@ -136,136 +136,7 @@ private:
   TDD_real_t m_alpha2;
 };
 
-inline DiamondBounds*
-DiamondBounds::clone() const
-{
-  return new DiamondBounds(*this);
-}
-
-inline double
-DiamondBounds::minHalflengthX() const
-{
-  return m_boundValues[DiamondBounds::bv_minHalfX];
-}
-
-inline double
-DiamondBounds::medHalflengthX() const
-{
-  return m_boundValues[DiamondBounds::bv_medHalfX];
-}
-
-inline double
-DiamondBounds::maxHalflengthX() const
-{
-  return m_boundValues[DiamondBounds::bv_maxHalfX];
-}
-
-inline double
-DiamondBounds::halflengthY1() const
-{
-  return m_boundValues[DiamondBounds::bv_halfY1];
-}
-
-inline double
-DiamondBounds::halflengthY2() const
-{
-  return m_boundValues[DiamondBounds::bv_halfY2];
-}
-
-inline double
-DiamondBounds::r() const
-{
-  return sqrt(m_boundValues[DiamondBounds::bv_medHalfX] * m_boundValues[DiamondBounds::bv_medHalfX] +
-              m_boundValues[DiamondBounds::bv_halfY1] * m_boundValues[DiamondBounds::bv_halfY1]);
-}
-
-inline bool
-DiamondBounds::inside(const Amg::Vector2D& locpo, const BoundaryCheck& bchk) const
-{
-  if (bchk.bcType == 0)
-    return DiamondBounds::inside(locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
-
-  // a fast FALSE
-  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 (locpo[Trk::locY] < -2 * m_boundValues[DiamondBounds::bv_halfY1] - limit)
-    return false;
-  if (locpo[Trk::locY] > 2 * m_boundValues[DiamondBounds::bv_halfY2] + limit)
-    return false;
-  // a fast FALSE
-  double fabsX = fabs(locpo[Trk::locX]);
-  if (fabsX > (m_boundValues[DiamondBounds::bv_medHalfX] + 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 < (fmin(m_boundValues[DiamondBounds::bv_minHalfX], m_boundValues[DiamondBounds::bv_maxHalfX]) - limit))
-    return true;
-  // a fast TRUE
-  if (fabs(locpo[Trk::locY]) <
-      (fmin(m_boundValues[DiamondBounds::bv_halfY1], m_boundValues[DiamondBounds::bv_halfY2]) - limit))
-    return true;
-
-  // compute KDOP and axes for surface polygon
-  std::vector<KDOP> elementKDOP(5);
-  std::vector<Amg::Vector2D> elementP(6);
-  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[DiamondBounds::bv_minHalfX], -2. * m_boundValues[DiamondBounds::bv_halfY1];
-  elementP[0] = (rotMatrix * (p - locpo));
-  p << -m_boundValues[DiamondBounds::bv_medHalfX], 0.;
-  elementP[1] = (rotMatrix * (p - locpo));
-  p << -m_boundValues[DiamondBounds::bv_maxHalfX], 2. * m_boundValues[DiamondBounds::bv_halfY2];
-  elementP[2] = (rotMatrix * (p - locpo));
-  p << m_boundValues[DiamondBounds::bv_maxHalfX], 2. * m_boundValues[DiamondBounds::bv_halfY2];
-  elementP[3] = (rotMatrix * (p - locpo));
-  p << m_boundValues[DiamondBounds::bv_medHalfX], 0.;
-  elementP[4] = (rotMatrix * (p - locpo));
-  p << m_boundValues[DiamondBounds::bv_minHalfX], -2. * m_boundValues[DiamondBounds::bv_halfY1];
-  elementP[5] = (rotMatrix * (p - locpo));
-  std::vector<Amg::Vector2D> axis = { normal * (elementP[1] - elementP[0]),
-                                      normal * (elementP[2] - elementP[1]),
-                                      normal * (elementP[3] - elementP[2]),
-                                      normal * (elementP[4] - elementP[3]),
-                                      normal * (elementP[5] - elementP[4]) };
-  bchk.ComputeKDOP(elementP, axis, elementKDOP);
-  // compute KDOP for error ellipse
-  std::vector<KDOP> errelipseKDOP(5);
-  bchk.ComputeKDOP(bchk.EllipseToPoly(3), axis, errelipseKDOP);
-  // check if KDOPs overlap and return result
-  return bchk.TestKDOPKDOP(elementKDOP, errelipseKDOP);
-}
-
-inline bool
-DiamondBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
-{
-  return (fabs(locpo[locX]) < m_boundValues[DiamondBounds::bv_medHalfX] + tol1);
-}
-
-inline bool
-DiamondBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
-{
-  return ((locpo[locY] > -2. * m_boundValues[DiamondBounds::bv_halfY1] - tol2) &&
-          (locpo[locY] < 2. * m_boundValues[DiamondBounds::bv_halfY2] + tol2));
-}
-
-inline void
-DiamondBounds::initCache()
-{
-  m_alpha1 = atan2(m_boundValues[DiamondBounds::bv_medHalfX] - m_boundValues[DiamondBounds::bv_minHalfX],
-                   2. * m_boundValues[DiamondBounds::bv_halfY1]);
-  m_alpha2 = atan2(m_boundValues[DiamondBounds::bv_medHalfX] - m_boundValues[DiamondBounds::bv_maxHalfX],
-                   2. * m_boundValues[DiamondBounds::bv_halfY2]);
-}
-
 } // end of namespace
 
+#include "TrkSurfaces/DiamondBounds.icc"
 #endif // TRKSURFACES_DIAMONDBOUNDS_H
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/DiamondBounds.icc b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/DiamondBounds.icc
new file mode 100644
index 0000000000000000000000000000000000000000..d2a631cc4b0111aed4bafcff7bb5f3baeb639d09
--- /dev/null
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/DiamondBounds.icc
@@ -0,0 +1,159 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+namespace Trk {
+
+inline DiamondBounds*
+DiamondBounds::clone() const
+{
+  return new DiamondBounds(*this);
+}
+
+inline double
+DiamondBounds::minHalflengthX() const
+{
+  return m_boundValues[DiamondBounds::bv_minHalfX];
+}
+
+inline double
+DiamondBounds::medHalflengthX() const
+{
+  return m_boundValues[DiamondBounds::bv_medHalfX];
+}
+
+inline double
+DiamondBounds::maxHalflengthX() const
+{
+  return m_boundValues[DiamondBounds::bv_maxHalfX];
+}
+
+inline double
+DiamondBounds::halflengthY1() const
+{
+  return m_boundValues[DiamondBounds::bv_halfY1];
+}
+
+inline double
+DiamondBounds::halflengthY2() const
+{
+  return m_boundValues[DiamondBounds::bv_halfY2];
+}
+
+inline double
+DiamondBounds::r() const
+{
+  return sqrt(m_boundValues[DiamondBounds::bv_medHalfX] *
+                m_boundValues[DiamondBounds::bv_medHalfX] +
+              m_boundValues[DiamondBounds::bv_halfY1] *
+                m_boundValues[DiamondBounds::bv_halfY1]);
+}
+
+inline bool
+DiamondBounds::inside(const Amg::Vector2D& locpo,
+                      const BoundaryCheck& bchk) const
+{
+  if (bchk.bcType == 0)
+    return DiamondBounds::inside(locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
+
+  // a fast FALSE
+  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 (locpo[Trk::locY] < -2 * m_boundValues[DiamondBounds::bv_halfY1] - limit)
+    return false;
+  if (locpo[Trk::locY] > 2 * m_boundValues[DiamondBounds::bv_halfY2] + limit)
+    return false;
+  // a fast FALSE
+  double fabsX = fabs(locpo[Trk::locX]);
+  if (fabsX > (m_boundValues[DiamondBounds::bv_medHalfX] + 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 < (fmin(m_boundValues[DiamondBounds::bv_minHalfX],
+                    m_boundValues[DiamondBounds::bv_maxHalfX]) -
+               limit))
+    return true;
+  // a fast TRUE
+  if (fabs(locpo[Trk::locY]) < (fmin(m_boundValues[DiamondBounds::bv_halfY1],
+                                     m_boundValues[DiamondBounds::bv_halfY2]) -
+                                limit))
+    return true;
+
+  // compute KDOP and axes for surface polygon
+  std::vector<KDOP> elementKDOP(5);
+  std::vector<Amg::Vector2D> elementP(6);
+  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[DiamondBounds::bv_minHalfX],
+    -2. * m_boundValues[DiamondBounds::bv_halfY1];
+  elementP[0] = (rotMatrix * (p - locpo));
+  p << -m_boundValues[DiamondBounds::bv_medHalfX], 0.;
+  elementP[1] = (rotMatrix * (p - locpo));
+  p << -m_boundValues[DiamondBounds::bv_maxHalfX],
+    2. * m_boundValues[DiamondBounds::bv_halfY2];
+  elementP[2] = (rotMatrix * (p - locpo));
+  p << m_boundValues[DiamondBounds::bv_maxHalfX],
+    2. * m_boundValues[DiamondBounds::bv_halfY2];
+  elementP[3] = (rotMatrix * (p - locpo));
+  p << m_boundValues[DiamondBounds::bv_medHalfX], 0.;
+  elementP[4] = (rotMatrix * (p - locpo));
+  p << m_boundValues[DiamondBounds::bv_minHalfX],
+    -2. * m_boundValues[DiamondBounds::bv_halfY1];
+  elementP[5] = (rotMatrix * (p - locpo));
+  std::vector<Amg::Vector2D> axis = { normal * (elementP[1] - elementP[0]),
+                                      normal * (elementP[2] - elementP[1]),
+                                      normal * (elementP[3] - elementP[2]),
+                                      normal * (elementP[4] - elementP[3]),
+                                      normal * (elementP[5] - elementP[4]) };
+  bchk.ComputeKDOP(elementP, axis, elementKDOP);
+  // compute KDOP for error ellipse
+  std::vector<KDOP> errelipseKDOP(5);
+  bchk.ComputeKDOP(bchk.EllipseToPoly(3), axis, errelipseKDOP);
+  // check if KDOPs overlap and return result
+  return bchk.TestKDOPKDOP(elementKDOP, errelipseKDOP);
+}
+
+inline bool
+DiamondBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
+{
+  return (fabs(locpo[locX]) < m_boundValues[DiamondBounds::bv_medHalfX] + tol1);
+}
+
+inline bool
+DiamondBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
+{
+  return (
+    (locpo[locY] > -2. * m_boundValues[DiamondBounds::bv_halfY1] - tol2) &&
+    (locpo[locY] < 2. * m_boundValues[DiamondBounds::bv_halfY2] + tol2));
+}
+
+inline void
+DiamondBounds::initCache()
+{
+  m_alpha1 = atan2(m_boundValues[DiamondBounds::bv_medHalfX] -
+                     m_boundValues[DiamondBounds::bv_minHalfX],
+                   2. * m_boundValues[DiamondBounds::bv_halfY1]);
+  m_alpha2 = atan2(m_boundValues[DiamondBounds::bv_medHalfX] -
+                     m_boundValues[DiamondBounds::bv_maxHalfX],
+                   2. * m_boundValues[DiamondBounds::bv_halfY2]);
+}
+
+} // end of namespace
+