diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/AnnulusBounds.h b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/AnnulusBounds.h
index 9e017a7d3de3c3c72641ef10e75b0f2a823d992b..26bb21aedeca83ac965d36454cbeb142c00f003c 100644
--- a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/AnnulusBounds.h
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/AnnulusBounds.h
@@ -203,65 +203,7 @@ private:
   std::vector<TDD_real_t> m_solution_R_max;
 };
 
-inline AnnulusBounds*
-AnnulusBounds::clone() const
-{
-  return new AnnulusBounds(*this);
-}
-
-inline double
-AnnulusBounds::minR() const
-{
-  return m_boundValues[AnnulusBounds::bv_minR];
-}
-
-inline double
-AnnulusBounds::maxR() const
-{
-  return m_boundValues[AnnulusBounds::bv_maxR];
-}
-
-inline double
-AnnulusBounds::R() const
-{
-  return m_boundValues[AnnulusBounds::bv_R];
-}
-
-inline double
-AnnulusBounds::phi() const
-{
-  return m_boundValues[AnnulusBounds::bv_phi];
-}
-
-inline double
-AnnulusBounds::phiS() const
-{
-  return m_boundValues[AnnulusBounds::bv_phiS];
-}
-
-inline double
-AnnulusBounds::r() const
-
-{
-  return AnnulusBounds::bv_maxR;
-} // MW to be fixed!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-inline bool
-AnnulusBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
-{
-  return (locpo[locX] > std::min(m_solution_L_min[0], m_solution_L_max[0]) - tol1 &&
-          locpo[locX] < std::max(m_solution_R_min[0], m_solution_R_max[0]) + tol1);
-}
-// MW Fix it
-
-inline bool
-AnnulusBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
-{
-  return (locpo[locY] > std::min(m_solution_L_min[1], m_solution_L_max[1]) - tol2 &&
-          locpo[locY] < std::max(m_solution_R_min[1], m_solution_R_max[1]) + tol2);
-}
-// MW Fix it
-
 } // end of namespace
 
+#include "TrkSurfaces/AnnulusBounds.icc"
 #endif // TRKSURFACES_ANNULUSBOUNDS_H
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/AnnulusBounds.icc b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/AnnulusBounds.icc
new file mode 100644
index 0000000000000000000000000000000000000000..e0d4e411108c0dc574ca109a05c29208367227f3
--- /dev/null
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/AnnulusBounds.icc
@@ -0,0 +1,68 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+namespace Trk {
+inline AnnulusBounds*
+AnnulusBounds::clone() const
+{
+  return new AnnulusBounds(*this);
+}
+
+inline double
+AnnulusBounds::minR() const
+{
+  return m_boundValues[AnnulusBounds::bv_minR];
+}
+
+inline double
+AnnulusBounds::maxR() const
+{
+  return m_boundValues[AnnulusBounds::bv_maxR];
+}
+
+inline double
+AnnulusBounds::R() const
+{
+  return m_boundValues[AnnulusBounds::bv_R];
+}
+
+inline double
+AnnulusBounds::phi() const
+{
+  return m_boundValues[AnnulusBounds::bv_phi];
+}
+
+inline double
+AnnulusBounds::phiS() const
+{
+  return m_boundValues[AnnulusBounds::bv_phiS];
+}
+
+inline double
+AnnulusBounds::r() const
+
+{
+  return AnnulusBounds::bv_maxR;
+} // MW to be fixed!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+inline bool
+AnnulusBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
+{
+  return (
+    locpo[locX] > std::min(m_solution_L_min[0], m_solution_L_max[0]) - tol1 &&
+    locpo[locX] < std::max(m_solution_R_min[0], m_solution_R_max[0]) + tol1);
+}
+// MW Fix it
+
+inline bool
+AnnulusBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
+{
+  return (
+    locpo[locY] > std::min(m_solution_L_min[1], m_solution_L_max[1]) - tol2 &&
+    locpo[locY] < std::max(m_solution_R_min[1], m_solution_R_max[1]) + tol2);
+}
+// MW Fix it
+
+}
+
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/BoundaryCheck.h b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/BoundaryCheck.h
index 84459bae0138f3ae5c026d5715be89bec6e356d2..525dbb8b6a48d35b823aa4c7cccd55e121e9b0fa 100644
--- a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/BoundaryCheck.h
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/BoundaryCheck.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 ///////////////////////////////////////////////////////////////////
@@ -110,7 +110,8 @@ public:
 
   /** Each Bounds has a method inside, which checks if a LocalPosition is inside the bounds.
             Inside can be called without/with boundary check */
-  void ComputeKDOP(const std::vector<Amg::Vector2D> &v, const std::vector<Amg::Vector2D> &KDOPAxes, std::vector<KDOP>& kdop) const;
+  void ComputeKDOP(const std::vector<Amg::Vector2D> &v, const std::vector<Amg::Vector2D> &KDOPAxes, 
+                   std::vector<KDOP>& kdop) const;
 
   std::vector<Amg::Vector2D> EllipseToPoly(int resolution = 3) const;
 
@@ -121,178 +122,7 @@ public:
   sincosCache FastSinCos(double x) const;
 };
 
-// should have maximum (average) error of 0.0015 (0.00089) radians or 0.0859 (0.0509) degrees, fine for us and much
-// faster (>4 times)
-inline double
-BoundaryCheck::FastArcTan(double x) const
-{
-  double y;
-  bool complement = false; // true if arg was >1
-  bool sign = false;       // true if arg was < 0
-  if (x < 0.) {
-    x = -x;
-    sign = true; // arctan(-x)=-arctan(x)
-  }
-  if (x > 1.) {
-    x = 1. / x; // keep arg between 0 and 1
-    complement = true;
-  }
-  y = M_PI_4 * x - x * (fabs(x) - 1) * (0.2447 + 0.0663 * fabs(x));
-  if (complement)
-    y = M_PI_2 - y; // correct for 1/x if we did that
-  if (sign)
-    y = -y; // correct for negative arg
-  return y;
-}
-
-// should have maximum (average) error of 0.001 (0.0005) radians or 0.0573 (0.029) degrees, fine for us and much faster
-// (>8 times)
-inline sincosCache
-BoundaryCheck::FastSinCos(double x) const
-{
-  sincosCache tmp;
-  // always wrap input angle to -PI..PI
-  if (x < -M_PI)
-    x += 2. * M_PI;
-  else if (x > M_PI)
-    x -= 2. * M_PI;
-
-  // compute sine
-  if (x < 0.) {
-    tmp.sinC = 1.27323954 * x + .405284735 * x * x;
-
-    if (tmp.sinC < 0.)
-      tmp.sinC = .225 * (tmp.sinC * -tmp.sinC - tmp.sinC) + tmp.sinC;
-    else
-      tmp.sinC = .225 * (tmp.sinC * tmp.sinC - tmp.sinC) + tmp.sinC;
-  } else {
-    tmp.sinC = 1.27323954 * x - 0.405284735 * x * x;
-
-    if (tmp.sinC < 0.)
-      tmp.sinC = .225 * (tmp.sinC * -tmp.sinC - tmp.sinC) + tmp.sinC;
-    else
-      tmp.sinC = .225 * (tmp.sinC * tmp.sinC - tmp.sinC) + tmp.sinC;
-  }
-
-  // compute cosine: sin(x + PI/2) = cos(x)
-  x += M_PI_2;
-  if (x > M_PI)
-    x -= 2. * M_PI;
-
-  if (x < 0.) {
-    tmp.cosC = 1.27323954 * x + 0.405284735 * x * x;
-
-    if (tmp.cosC < 0.)
-      tmp.cosC = .225 * (tmp.cosC * -tmp.cosC - tmp.cosC) + tmp.cosC;
-    else
-      tmp.cosC = .225 * (tmp.cosC * tmp.cosC - tmp.cosC) + tmp.cosC;
-  } else {
-    tmp.cosC = 1.27323954 * x - 0.405284735 * x * x;
-
-    if (tmp.cosC < 0.)
-      tmp.cosC = .225 * (tmp.cosC * -tmp.cosC - tmp.cosC) + tmp.cosC;
-    else
-      tmp.cosC = .225 * (tmp.cosC * tmp.cosC - tmp.cosC) + tmp.cosC;
-  }
-  return tmp;
-}
-
-// does the conversion of an ellipse of height h and width w to an polygon with 4 + 4*resolution points
-inline std::vector<Amg::Vector2D>
-BoundaryCheck::EllipseToPoly(int resolution) const
-{
-  const double h = nSigmas * sqrt(lCovariance(1, 1));
-  const double w = nSigmas * sqrt(lCovariance(0, 0));
-
-  // first add the four vertices
-  std::vector<Amg::Vector2D> v((1 + resolution) * 4);
-  Amg::Vector2D p;
-  p << w, 0;
-  v[0] = p;
-  p << -w, 0;
-  v[1] = p;
-  p << 0, h;
-  v[2] = p;
-  p << 0, -h;
-  v[3] = p;
-
-  // now add a number, equal to the resolution, of equidistant points  in each quadrant
-  // resolution == 3 seems to be a solid working point, but possibility open to change to any number in the future
-  Amg::Vector2D t(1, 1);
-  AmgSymMatrix(2) t1;
-  t1 << 1, 0, 0, -1;
-  AmgSymMatrix(2) t2;
-  t2 << -1, 0, 0, -1;
-  AmgSymMatrix(2) t3;
-  t3 << -1, 0, 0, 1;
-  if (resolution != 3) {
-    sincosCache scResult;
-    for (int i = 1; i <= resolution; i++) {
-      scResult = FastSinCos(M_PI_2 * i / (resolution + 1));
-      t << w * scResult.sinC, h * scResult.cosC;
-      v[i * 4 + 0] = t;
-      v[i * 4 + 1] = t1 * t;
-      v[i * 4 + 2] = t2 * t;
-      v[i * 4 + 3] = t3 * t;
-    }
-  } else {
-    t << w * s_cos22, h * s_cos67;
-    v[4] = t;
-    v[5] = t1 * t;
-    v[6] = t2 * t;
-    v[7] = t3 * t;
-    t << w * s_cos45, h * s_cos45;
-    v[8] = t;
-    v[9] = t1 * t;
-    v[10] = t2 * t;
-    v[11] = t3 * t;
-    t << w * s_cos67, h * s_cos22;
-    v[12] = t;
-    v[13] = t1 * t;
-    v[14] = t2 * t;
-    v[15] = t3 * t;
-  }
-  return v;
-}
-
-// calculates KDOP object from given polygon and set of axes
-inline void
-BoundaryCheck::ComputeKDOP(const std::vector<Amg::Vector2D> &v,
-                           const std::vector<Amg::Vector2D> &KDOPAxes,
-                           std::vector<KDOP>& kdop) const
-{
-  // initialize KDOP to first point
-  unsigned int k = KDOPAxes.size();
-  for (unsigned int i = 0; i < k; i++) {
-    kdop[i].max = KDOPAxes[i](0, 0) * v[0](0, 0) + KDOPAxes[i](1, 0) * v[0](1, 0);
-    kdop[i].min = KDOPAxes[i](0, 0) * v[0](0, 0) + KDOPAxes[i](1, 0) * v[0](1, 0);
-  }
-  // now for each additional point, update KDOP bounds if necessary
-  float value;
-  for (unsigned int i = 1; i < v.size(); i++) {
-    for (unsigned int j = 0; j < k; j++) {
-      value = KDOPAxes[j](0, 0) * v[i](0, 0) + KDOPAxes[j](1, 0) * v[i](1, 0);
-      if (value < kdop[j].min)
-        kdop[j].min = value;
-      else if (value > kdop[j].max)
-        kdop[j].max = value;
-    }
-  }
-}
-
-// this is the method to actually check if two KDOPs overlap
-inline bool
-BoundaryCheck::TestKDOPKDOP(const std::vector<KDOP>& a, const std::vector<KDOP>& b) const
-{
-  int k = a.size();
-  // check if any intervals are non-overlapping, return if so
-  for (int i = 0; i < k; i++)
-    if (a[i].min > b[i].max || a[i].max < b[i].min)
-      return false;
-  // all intervals are overlapping, so KDOPs must intersect
-  return true;
-}
-
 }
 
+#include "TrkSurfaces/BoundaryCheck.icc"
 #endif // TRKSURFACES_BOUNDARYCHECK_H
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/BoundaryCheck.icc b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/BoundaryCheck.icc
new file mode 100644
index 0000000000000000000000000000000000000000..89d18b4645e436244e9850a077b7136ae0e06665
--- /dev/null
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/BoundaryCheck.icc
@@ -0,0 +1,185 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+namespace Trk {
+// should have maximum (average) error of 0.0015 (0.00089) radians or 0.0859
+// (0.0509) degrees, fine for us and much faster (>4 times)
+inline double
+BoundaryCheck::FastArcTan(double x) const
+{
+  double y;
+  bool complement = false; // true if arg was >1
+  bool sign = false;       // true if arg was < 0
+  if (x < 0.) {
+    x = -x;
+    sign = true; // arctan(-x)=-arctan(x)
+  }
+  if (x > 1.) {
+    x = 1. / x; // keep arg between 0 and 1
+    complement = true;
+  }
+  y = M_PI_4 * x - x * (fabs(x) - 1) * (0.2447 + 0.0663 * fabs(x));
+  if (complement)
+    y = M_PI_2 - y; // correct for 1/x if we did that
+  if (sign)
+    y = -y; // correct for negative arg
+  return y;
+}
+
+// should have maximum (average) error of 0.001 (0.0005) radians or 0.0573
+// (0.029) degrees, fine for us and much faster
+// (>8 times)
+inline sincosCache
+BoundaryCheck::FastSinCos(double x) const
+{
+  sincosCache tmp;
+  // always wrap input angle to -PI..PI
+  if (x < -M_PI)
+    x += 2. * M_PI;
+  else if (x > M_PI)
+    x -= 2. * M_PI;
+
+  // compute sine
+  if (x < 0.) {
+    tmp.sinC = 1.27323954 * x + .405284735 * x * x;
+
+    if (tmp.sinC < 0.)
+      tmp.sinC = .225 * (tmp.sinC * -tmp.sinC - tmp.sinC) + tmp.sinC;
+    else
+      tmp.sinC = .225 * (tmp.sinC * tmp.sinC - tmp.sinC) + tmp.sinC;
+  } else {
+    tmp.sinC = 1.27323954 * x - 0.405284735 * x * x;
+
+    if (tmp.sinC < 0.)
+      tmp.sinC = .225 * (tmp.sinC * -tmp.sinC - tmp.sinC) + tmp.sinC;
+    else
+      tmp.sinC = .225 * (tmp.sinC * tmp.sinC - tmp.sinC) + tmp.sinC;
+  }
+
+  // compute cosine: sin(x + PI/2) = cos(x)
+  x += M_PI_2;
+  if (x > M_PI)
+    x -= 2. * M_PI;
+
+  if (x < 0.) {
+    tmp.cosC = 1.27323954 * x + 0.405284735 * x * x;
+
+    if (tmp.cosC < 0.)
+      tmp.cosC = .225 * (tmp.cosC * -tmp.cosC - tmp.cosC) + tmp.cosC;
+    else
+      tmp.cosC = .225 * (tmp.cosC * tmp.cosC - tmp.cosC) + tmp.cosC;
+  } else {
+    tmp.cosC = 1.27323954 * x - 0.405284735 * x * x;
+
+    if (tmp.cosC < 0.)
+      tmp.cosC = .225 * (tmp.cosC * -tmp.cosC - tmp.cosC) + tmp.cosC;
+    else
+      tmp.cosC = .225 * (tmp.cosC * tmp.cosC - tmp.cosC) + tmp.cosC;
+  }
+  return tmp;
+}
+
+// does the conversion of an ellipse of height h and width w to an polygon with
+// 4 + 4*resolution points
+inline std::vector<Amg::Vector2D>
+BoundaryCheck::EllipseToPoly(int resolution) const
+{
+  const double h = nSigmas * sqrt(lCovariance(1, 1));
+  const double w = nSigmas * sqrt(lCovariance(0, 0));
+
+  // first add the four vertices
+  std::vector<Amg::Vector2D> v((1 + resolution) * 4);
+  Amg::Vector2D p;
+  p << w, 0;
+  v[0] = p;
+  p << -w, 0;
+  v[1] = p;
+  p << 0, h;
+  v[2] = p;
+  p << 0, -h;
+  v[3] = p;
+
+  // now add a number, equal to the resolution, of equidistant points  in each
+  // quadrant resolution == 3 seems to be a solid working point, but possibility
+  // open to change to any number in the future
+  Amg::Vector2D t(1, 1);
+  AmgSymMatrix(2) t1;
+  t1 << 1, 0, 0, -1;
+  AmgSymMatrix(2) t2;
+  t2 << -1, 0, 0, -1;
+  AmgSymMatrix(2) t3;
+  t3 << -1, 0, 0, 1;
+  if (resolution != 3) {
+    sincosCache scResult;
+    for (int i = 1; i <= resolution; i++) {
+      scResult = FastSinCos(M_PI_2 * i / (resolution + 1));
+      t << w * scResult.sinC, h * scResult.cosC;
+      v[i * 4 + 0] = t;
+      v[i * 4 + 1] = t1 * t;
+      v[i * 4 + 2] = t2 * t;
+      v[i * 4 + 3] = t3 * t;
+    }
+  } else {
+    t << w * s_cos22, h * s_cos67;
+    v[4] = t;
+    v[5] = t1 * t;
+    v[6] = t2 * t;
+    v[7] = t3 * t;
+    t << w * s_cos45, h * s_cos45;
+    v[8] = t;
+    v[9] = t1 * t;
+    v[10] = t2 * t;
+    v[11] = t3 * t;
+    t << w * s_cos67, h * s_cos22;
+    v[12] = t;
+    v[13] = t1 * t;
+    v[14] = t2 * t;
+    v[15] = t3 * t;
+  }
+  return v;
+}
+
+// calculates KDOP object from given polygon and set of axes
+inline void
+BoundaryCheck::ComputeKDOP(const std::vector<Amg::Vector2D>& v,
+                           const std::vector<Amg::Vector2D>& KDOPAxes,
+                           std::vector<KDOP>& kdop) const
+{
+  // initialize KDOP to first point
+  unsigned int k = KDOPAxes.size();
+  for (unsigned int i = 0; i < k; i++) {
+    kdop[i].max =
+      KDOPAxes[i](0, 0) * v[0](0, 0) + KDOPAxes[i](1, 0) * v[0](1, 0);
+    kdop[i].min =
+      KDOPAxes[i](0, 0) * v[0](0, 0) + KDOPAxes[i](1, 0) * v[0](1, 0);
+  }
+  // now for each additional point, update KDOP bounds if necessary
+  float value;
+  for (unsigned int i = 1; i < v.size(); i++) {
+    for (unsigned int j = 0; j < k; j++) {
+      value = KDOPAxes[j](0, 0) * v[i](0, 0) + KDOPAxes[j](1, 0) * v[i](1, 0);
+      if (value < kdop[j].min)
+        kdop[j].min = value;
+      else if (value > kdop[j].max)
+        kdop[j].max = value;
+    }
+  }
+}
+
+// this is the method to actually check if two KDOPs overlap
+inline bool
+BoundaryCheck::TestKDOPKDOP(const std::vector<KDOP>& a,
+                            const std::vector<KDOP>& b) const
+{
+  int k = a.size();
+  // check if any intervals are non-overlapping, return if so
+  for (int i = 0; i < k; i++)
+    if (a[i].min > b[i].max || a[i].max < b[i].min)
+      return false;
+  // all intervals are overlapping, so KDOPs must intersect
+  return true;
+}
+
+}
+
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/ConeBounds.h b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/ConeBounds.h
index df5fd340dae08e2f8a0312359a495459d422dcbe..b83ea854ec34e1310a55c2196bb717b5a86abdcf 100644
--- a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/ConeBounds.h
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/ConeBounds.h
@@ -158,141 +158,7 @@ private:
   }
 };
 
-inline ConeBounds*
-ConeBounds::clone() const
-{
-  return new ConeBounds(*this);
-}
-
-inline bool
-ConeBounds::inside(const Amg::Vector2D& locpo, double tol1, double tol2) const
-{
-  double z = locpo[locZ];
-  bool insideZ = z > (m_boundValues[ConeBounds::bv_minZ] - tol2) && z < (m_boundValues[ConeBounds::bv_maxZ] + tol2);
-  if (!insideZ)
-    return false;
-  // TODO: Do we need some sort of "R" tolerance also here (take
-  // it off the z tol2 in that case?) or does the rphi tol1 cover
-  // this? (Could argue either way)
-  double coneR = z * m_tanAlpha;
-  double minRPhi = coneR * minPhi() - tol1;
-  double maxRPhi = coneR * maxPhi() + tol1;
-  return minRPhi < locpo[locRPhi] && locpo[locRPhi] < maxRPhi;
-}
-
-inline bool
-ConeBounds::inside(const Amg::Vector3D& glopo, double tol1, double tol2) const
-{
-  // coords are (rphi,z)
-  return inside(Amg::Vector2D(glopo.perp() * glopo.phi(), glopo.z()), tol1, tol2);
-}
-
-inline bool
-ConeBounds::inside(const Amg::Vector3D& glopo, const BoundaryCheck& bchk) const
-{
-  // coords are (rphi,z)
-  return inside(Amg::Vector2D(glopo.perp() * glopo.phi(), glopo.z()), bchk.toleranceLoc1, bchk.toleranceLoc2);
-}
-
-inline bool
-ConeBounds::inside(const Amg::Vector2D& locpo, const BoundaryCheck& bchk) const
-{
-  return ConeBounds::inside(locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
-}
-
-inline bool
-ConeBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
-{
-  double z = locpo[locZ];
-  double coneR = z * m_tanAlpha;
-  double minRPhi = coneR * minPhi() - tol1;
-  double maxRPhi = coneR * maxPhi() + tol1;
-  return minRPhi < locpo[locRPhi] && locpo[locRPhi] < maxRPhi;
-}
-
-inline bool
-ConeBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
-{
-  double z = locpo[locZ];
-  return (z > (m_boundValues[ConeBounds::bv_minZ] - tol2) && z < (m_boundValues[ConeBounds::bv_maxZ] + tol2));
-}
-
-inline double
-ConeBounds::r() const
-{
-  double z = fabs(m_boundValues[ConeBounds::bv_maxZ]);
-  double mz = fabs(m_boundValues[ConeBounds::bv_minZ]);
-  if (mz > z)
-    z = mz;
-  return fabs(z * m_tanAlpha);
-}
-
-inline double
-ConeBounds::r(double z) const
-{
-  return fabs(z * m_tanAlpha);
-}
-
-inline double
-ConeBounds::tanAlpha() const
-{
-  return m_tanAlpha;
-}
-
-inline double
-ConeBounds::sinAlpha() const
-{
-  return m_sinAlpha;
-}
-
-inline double
-ConeBounds::cosAlpha() const
-{
-  return m_cosAlpha;
-}
-
-inline double
-ConeBounds::alpha() const
-{
-  return m_boundValues[ConeBounds::bv_alpha];
-}
-
-inline double
-ConeBounds::minZ() const
-{
-  return m_boundValues[ConeBounds::bv_minZ];
-}
-
-inline double
-ConeBounds::maxZ() const
-{
-  return m_boundValues[ConeBounds::bv_maxZ];
-}
-
-inline double
-ConeBounds::averagePhi() const
-{
-  return m_boundValues[ConeBounds::bv_averagePhi];
-}
-
-inline double
-ConeBounds::halfPhiSector() const
-{
-  return m_boundValues[ConeBounds::bv_halfPhiSector];
-}
-
-inline void
-ConeBounds::initCache()
-{
-  m_tanAlpha = tan(m_boundValues[ConeBounds::bv_alpha]);
-  m_sinAlpha = sin(m_boundValues[ConeBounds::bv_alpha]);
-  m_cosAlpha = cos(m_boundValues[ConeBounds::bv_alpha]);
-  // validate the halfphi
-  if (m_boundValues[ConeBounds::bv_halfPhiSector] < 0.)
-    m_boundValues[ConeBounds::bv_halfPhiSector] = -m_boundValues[ConeBounds::bv_halfPhiSector];
-  if (m_boundValues[ConeBounds::bv_halfPhiSector] > M_PI)
-    m_boundValues[ConeBounds::bv_halfPhiSector] = M_PI;
-}
 }
 
+#include "TrkSurfaces/ConeBounds.icc"
 #endif // TRKSURFACES_CONEBOUNDS_H
diff --git a/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/ConeBounds.icc b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/ConeBounds.icc
new file mode 100644
index 0000000000000000000000000000000000000000..907affe4e8ee9d8131788873f221b9caac8026b8
--- /dev/null
+++ b/Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/ConeBounds.icc
@@ -0,0 +1,149 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+namespace Trk {
+inline ConeBounds*
+ConeBounds::clone() const
+{
+  return new ConeBounds(*this);
+}
+
+inline bool
+ConeBounds::inside(const Amg::Vector2D& locpo, double tol1, double tol2) const
+{
+  double z = locpo[locZ];
+  bool insideZ = z > (m_boundValues[ConeBounds::bv_minZ] - tol2) &&
+                 z < (m_boundValues[ConeBounds::bv_maxZ] + tol2);
+  if (!insideZ)
+    return false;
+  // TODO: Do we need some sort of "R" tolerance also here (take
+  // it off the z tol2 in that case?) or does the rphi tol1 cover
+  // this? (Could argue either way)
+  double coneR = z * m_tanAlpha;
+  double minRPhi = coneR * minPhi() - tol1;
+  double maxRPhi = coneR * maxPhi() + tol1;
+  return minRPhi < locpo[locRPhi] && locpo[locRPhi] < maxRPhi;
+}
+
+inline bool
+ConeBounds::inside(const Amg::Vector3D& glopo, double tol1, double tol2) const
+{
+  // coords are (rphi,z)
+  return inside(
+    Amg::Vector2D(glopo.perp() * glopo.phi(), glopo.z()), tol1, tol2);
+}
+
+inline bool
+ConeBounds::inside(const Amg::Vector3D& glopo, const BoundaryCheck& bchk) const
+{
+  // coords are (rphi,z)
+  return inside(Amg::Vector2D(glopo.perp() * glopo.phi(), glopo.z()),
+                bchk.toleranceLoc1,
+                bchk.toleranceLoc2);
+}
+
+inline bool
+ConeBounds::inside(const Amg::Vector2D& locpo, const BoundaryCheck& bchk) const
+{
+  return ConeBounds::inside(locpo, bchk.toleranceLoc1, bchk.toleranceLoc2);
+}
+
+inline bool
+ConeBounds::insideLoc1(const Amg::Vector2D& locpo, double tol1) const
+{
+  double z = locpo[locZ];
+  double coneR = z * m_tanAlpha;
+  double minRPhi = coneR * minPhi() - tol1;
+  double maxRPhi = coneR * maxPhi() + tol1;
+  return minRPhi < locpo[locRPhi] && locpo[locRPhi] < maxRPhi;
+}
+
+inline bool
+ConeBounds::insideLoc2(const Amg::Vector2D& locpo, double tol2) const
+{
+  double z = locpo[locZ];
+  return (z > (m_boundValues[ConeBounds::bv_minZ] - tol2) &&
+          z < (m_boundValues[ConeBounds::bv_maxZ] + tol2));
+}
+
+inline double
+ConeBounds::r() const
+{
+  double z = fabs(m_boundValues[ConeBounds::bv_maxZ]);
+  double mz = fabs(m_boundValues[ConeBounds::bv_minZ]);
+  if (mz > z)
+    z = mz;
+  return fabs(z * m_tanAlpha);
+}
+
+inline double
+ConeBounds::r(double z) const
+{
+  return fabs(z * m_tanAlpha);
+}
+
+inline double
+ConeBounds::tanAlpha() const
+{
+  return m_tanAlpha;
+}
+
+inline double
+ConeBounds::sinAlpha() const
+{
+  return m_sinAlpha;
+}
+
+inline double
+ConeBounds::cosAlpha() const
+{
+  return m_cosAlpha;
+}
+
+inline double
+ConeBounds::alpha() const
+{
+  return m_boundValues[ConeBounds::bv_alpha];
+}
+
+inline double
+ConeBounds::minZ() const
+{
+  return m_boundValues[ConeBounds::bv_minZ];
+}
+
+inline double
+ConeBounds::maxZ() const
+{
+  return m_boundValues[ConeBounds::bv_maxZ];
+}
+
+inline double
+ConeBounds::averagePhi() const
+{
+  return m_boundValues[ConeBounds::bv_averagePhi];
+}
+
+inline double
+ConeBounds::halfPhiSector() const
+{
+  return m_boundValues[ConeBounds::bv_halfPhiSector];
+}
+
+inline void
+ConeBounds::initCache()
+{
+  m_tanAlpha = tan(m_boundValues[ConeBounds::bv_alpha]);
+  m_sinAlpha = sin(m_boundValues[ConeBounds::bv_alpha]);
+  m_cosAlpha = cos(m_boundValues[ConeBounds::bv_alpha]);
+  // validate the halfphi
+  if (m_boundValues[ConeBounds::bv_halfPhiSector] < 0.)
+    m_boundValues[ConeBounds::bv_halfPhiSector] =
+      -m_boundValues[ConeBounds::bv_halfPhiSector];
+  if (m_boundValues[ConeBounds::bv_halfPhiSector] > M_PI)
+    m_boundValues[ConeBounds::bv_halfPhiSector] = M_PI;
+}
+
+}
+