diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 498eb1c..ba5b812 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -5,7 +5,7 @@ ci:
 
 repos:
   - repo: https://github.com/adamchainz/blacken-docs
-    rev: "1.18.0"
+    rev: "1.19.1"
     hooks:
       - id: blacken-docs
         additional_dependencies: [black==23.*]
@@ -42,14 +42,14 @@ repos:
         args: [--prose-wrap=always]
 
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: "v0.6.9"
+    rev: "v0.7.2"
     hooks:
       - id: ruff
         args: ["--fix", "--show-fixes"]
       - id: ruff-format
 
   - repo: https://github.com/pre-commit/mirrors-mypy
-    rev: "v1.11.2"
+    rev: "v1.13.0"
     hooks:
       - id: mypy
         files: src|tests
@@ -76,12 +76,12 @@ repos:
         exclude: .pre-commit-config.yaml
 
   - repo: https://github.com/abravalheri/validate-pyproject
-    rev: v0.20.2
+    rev: v0.22
     hooks:
       - id: validate-pyproject
 
   - repo: https://github.com/python-jsonschema/check-jsonschema
-    rev: 0.29.3
+    rev: 0.29.4
     hooks:
       - id: check-dependabot
       - id: check-github-workflows
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 155c38b..982d62a 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -58,6 +58,11 @@
 # sphinx-autodoc
 # Include __init__() docstring in class docstring
 autoclass_content = "both"
-autodoc_typehints = "both"
+autodoc_typehints = "description"
 autodoc_typehints_description_target = "documented_params"
 autodoc_typehints_format = "short"
+
+autodoc_type_aliases = {
+    "ArrayLike": "ArrayLike",
+    "NDArray": "NDArray",
+}
diff --git a/pyproject.toml b/pyproject.toml
index 6d9e9d0..cdf7a12 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -31,6 +31,7 @@ classifiers = [
 requires-python = ">=3.9"
 dependencies = [
     "numpy",
+    "numba",
     "pint != 0.24",
     "pyg4ometry",
     "pylegendmeta>=v0.9.0a2",
diff --git a/src/legendhpges/__init__.py b/src/legendhpges/__init__.py
index a82479e..141d0c8 100644
--- a/src/legendhpges/__init__.py
+++ b/src/legendhpges/__init__.py
@@ -18,6 +18,7 @@
     "PPC",
     "SemiCoax",
     "make_hpge",
+    "utils",
     "V07646A",
     "P00664B",
     "V02160A",
diff --git a/src/legendhpges/base.py b/src/legendhpges/base.py
index bb7d332..1112730 100644
--- a/src/legendhpges/base.py
+++ b/src/legendhpges/base.py
@@ -1,14 +1,18 @@
 from __future__ import annotations
 
 import json
+import logging
 import math
 from abc import ABC, abstractmethod
 from pathlib import Path
 
+import numpy as np
 from legendmeta import AttrsDict
+from numpy.typing import ArrayLike, NDArray
 from pint import Quantity
 from pyg4ometry import geant4
 
+from . import utils
 from .materials import make_natural_germanium
 from .registry import default_g4_registry
 from .registry import default_units_registry as u
@@ -63,6 +67,8 @@ def __init__(
 
         self.registry = registry
 
+        self.surfaces = []
+
         # build logical volume, default [mm]
         super().__init__(self._g4_solid(), material, self.name, self.registry)
 
@@ -70,16 +76,18 @@ def __repr__(self) -> str:
         return f"{self.__class__.__name__}({self.metadata})"
 
     def _g4_solid(self) -> geant4.solid.SolidBase:
-        """Build (by default) a :class:`pyg4ometry.solid.GenericPolycone` instance from the (r, z) information.
+        """Build a :class:`pyg4ometry.solid.GenericPolycone` instance from the ``(r, z)`` information.
 
         Returns
         -------
         g4_solid
-            A derived class of :class:`pyg4ometry.solid.SolidBase` to be used to construct the logical volume.
+            A derived class of :class:`pyg4ometry.solid.SolidBase` to be used
+            to construct the logical volume.
 
         Note
         ----
-            Detectors with a special geometry can have this method overridden in their class definition.
+            Detectors with a special geometry can have this method overridden
+            in their class definition.
         """
         # return ordered r,z lists, default unit [mm]
         r, z = self._decode_polycone_coord()
@@ -98,13 +106,128 @@ def _decode_polycone_coord(self) -> tuple[list[float], list[float]]:
         Returns
         -------
         (r, z)
-            two lists of r and z coordinates, respectively.
+            two lists of `r` and `z` coordinates, respectively.
 
         Note
         ----
         Must be overloaded by derived classes.
         """
 
+    def get_profile(self) -> tuple[list[float], list[float]]:
+        """Get the profile of the HPGe detector.
+
+        Returns
+        -------
+        (r.z)
+            two lists of `r` and `z` coordinates, respectively.
+
+        Note
+        -----
+            For V02160A and P00664A the detector profile is that of the solid without cut.
+        """
+        if not isinstance(self.solid, geant4.solid.GenericPolycone) and not isinstance(
+            self.solid.obj1, geant4.solid.GenericPolycone
+        ):
+            msg = "solid is not a polycone and neither is its primary consistient (obj1), thus no profile can be computed."
+            raise ValueError(msg)
+        if not isinstance(self.solid, geant4.solid.GenericPolycone):
+            r = self.solid.obj1.pR
+            z = self.solid.obj1.pZ
+        else:
+            r = self.solid.pR
+            z = self.solid.pZ
+
+        return r, z
+
+    def is_inside(self, coords: ArrayLike, tol: float = 1e-11) -> NDArray[bool]:
+        """Compute whether each point is inside the volume.
+
+        Parameters
+        ----------
+        coords
+            2D array of shape `(n,3)` of `(x,y,z)` coordinates for each of `n`
+            points, second index corresponds to `(x,y,z)`.
+        tol
+            distance outside the surface which is considered inside. Should be
+            on the order of numerical precision of the floating point representation.
+        """
+
+        if not isinstance(self.solid, geant4.solid.GenericPolycone):
+            msg = f"distance_to_surface is not implemented for {type(self.solid)} yet"
+            raise NotImplementedError(msg)
+
+        if not isinstance(coords, np.ndarray):
+            coords = np.array(coords)
+
+        if np.shape(coords)[1] != 3:
+            msg = "coords must be provided as a 2D array with x,y,z coordinates for each point."
+            raise ValueError(msg)
+
+        coords_rz = utils.convert_coords(coords)
+
+        # get the profile
+        r, z = self.get_profile()
+        s1, s2 = utils.get_line_segments(r, z, surface_indices=None)
+
+        # convert coords
+        coords_rz = utils.convert_coords(coords)
+
+        # compute shortest distances
+        dists = utils.shortest_distance(s1, s2, coords_rz, tol=tol)
+
+        # fnd the sign of the id with the lowest distance
+        ids = np.argmin(abs(dists), axis=1)
+        return np.where(dists[np.arange(dists.shape[0]), ids] > 0, True, False)
+
+    def distance_to_surface(
+        self,
+        coords: ArrayLike,
+        surface_indices: ArrayLike | None = None,
+        tol: float = 1e-11,
+    ) -> NDArray:
+        """Compute the distance of a set of points to the nearest detector surface.
+
+        Parameters
+        ----------
+        coords
+            2D array of shape `(n,3)` of `(x,y,z)` coordinates for each of `n`
+            points, second index corresponds to `(x,y,z)`.
+        surface_indices
+            list of indices of surfaces to consider. If ``None`` (the default)
+            all surfaces used.
+        tol
+            distance outside the surface which is considered inside. Should be
+            on the order of numerical precision of the floating point representation.
+
+        Note
+        ----
+        - Only implemented for solids based on :class:`geant4.solid.GenericPolycone`
+        - Coordinates should be relative to the origin of the polycone.
+        """
+        # check type of the solid
+        if not isinstance(self.solid, geant4.solid.GenericPolycone):
+            msg = f"distance_to_surface is not implemented for {type(self.solid)} yet"
+            raise NotImplementedError(msg)
+
+        if not isinstance(coords, np.ndarray):
+            coords = np.array(coords)
+
+        if np.shape(coords)[1] != 3:
+            msg = "coords must be provided as a 2D array with x,y,z coordinates for each point."
+            raise ValueError(msg)
+
+        # get the coordinates
+        r, z = self.get_profile()
+        s1, s2 = utils.get_line_segments(r, z, surface_indices=surface_indices)
+
+        # convert coords
+        coords_rz = utils.convert_coords(coords)
+
+        # compute shortest distances to every surface
+        dists = utils.shortest_distance(s1, s2, coords_rz, tol=tol, signed=False)
+
+        return np.min(abs(dists), axis=1)
+
     @property
     def volume(self) -> Quantity:
         """Volume of the HPGe."""
@@ -124,3 +247,47 @@ def volume(self) -> Quantity:
     def mass(self) -> Quantity:
         """Mass of the HPGe."""
         return (self.volume * (self.material.density * u.g / u.cm**3)).to(u.g)
+
+    def surface_area(
+        self, surface_indices: ArrayLike | None = None
+    ) -> NDArray[Quantity]:
+        """Surface area of the HPGe.
+
+        If a list of surface_indices is provided the area is computed only
+        considering these surfaces, from :math:`r_i` to :math:`r_{i+1}` and
+        :math:`z_i` to :math:`z_{i+1}`. Else the full area is computed.
+
+        Parameters
+        ----------
+        surface_indices
+            list of indices or ``None``.
+
+        Note
+        ----
+        Calculation is based on a polycone geometry so is incorrect for
+        asymmetric detectors.
+        """
+        if not isinstance(self.solid, geant4.solid.GenericPolycone):
+            logging.warning("The area is that of the solid without cut")
+
+        r, z = self.get_profile()
+
+        r = np.array(r)
+        z = np.array(z)
+
+        dr = np.array([r2 - r1 for r1, r2 in zip(r[:-1], r[1:])])
+        dz = np.array([z2 - z1 for z1, z2 in zip(z[:-1], z[1:])])
+        dl = np.sqrt(np.power(dr, 2) + np.power(dz, 2))
+        r0 = r[:-1]
+
+        if surface_indices is not None:
+            dr = dr[surface_indices]
+            dz = dz[surface_indices]
+            dl = dl[surface_indices]
+            r0 = r0[surface_indices]
+
+        return np.where(
+            dr == 0,
+            abs(dz) * r0 * 2 * np.pi * u.mm**2,
+            abs(dr) * dl * np.pi * u.mm**2,
+        )
diff --git a/src/legendhpges/bege.py b/src/legendhpges/bege.py
index 145111a..c951542 100644
--- a/src/legendhpges/bege.py
+++ b/src/legendhpges/bege.py
@@ -3,6 +3,7 @@
 import math
 
 from .base import HPGe
+from .build_utils import make_pplus
 
 
 class BEGe(HPGe):
@@ -16,22 +17,12 @@ def _tan(a):
 
         r = []
         z = []
+        surfaces = []
 
-        if c.pp_contact.depth_in_mm > 0:
-            r += [0, c.pp_contact.radius_in_mm, c.pp_contact.radius_in_mm]
-            z += [c.pp_contact.depth_in_mm, c.pp_contact.depth_in_mm, 0]
-        else:
-            r += [0]
-            z += [0]
-
-        r += [
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.outer,
-            c.groove.radius_in_mm.outer,
-        ]
-
-        z += [0, c.groove.depth_in_mm, c.groove.depth_in_mm, 0]
+        r_p, z_p, surface_p = make_pplus(c)
+        r += r_p
+        z += z_p
+        surfaces += surface_p
 
         if c.taper.bottom.height_in_mm > 0:
             r += [
@@ -40,9 +31,11 @@ def _tan(a):
                 c.radius_in_mm,
             ]
             z += [0, c.taper.bottom.height_in_mm]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [0]
+            surfaces += ["nplus"]
 
         if c.taper.top.height_in_mm > 0:
             r += [
@@ -51,11 +44,16 @@ def _tan(a):
                 - c.taper.top.height_in_mm * _tan(c.taper.top.angle_in_deg),
             ]
             z += [c.height_in_mm - c.taper.top.height_in_mm, c.height_in_mm]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         r += [0]
         z += [c.height_in_mm]
+        surfaces += ["nplus"]
+
+        self.surfaces = surfaces
 
         return r, z
diff --git a/src/legendhpges/build_utils.py b/src/legendhpges/build_utils.py
new file mode 100644
index 0000000..dd404c2
--- /dev/null
+++ b/src/legendhpges/build_utils.py
@@ -0,0 +1,66 @@
+from __future__ import annotations
+
+from collections.abc import Mapping
+
+
+def make_pplus(geometry: Mapping) -> tuple[list, list, list]:
+    """Make the pplus contact for BeGe and some ICPC.
+
+    Methods to avoid duplicating code.
+
+    Parameters
+    ----------
+    geometry
+        dictionary with the geometry information.
+
+    Returns
+    -------
+    (r,z,surfaces)
+        Tuple of lists of r,z coordinates and surface names.
+    """
+    r = []
+    z = []
+    surfaces = []
+
+    if geometry.pp_contact.depth_in_mm > 0:
+        r += [
+            0,
+            geometry.pp_contact.radius_in_mm,
+            geometry.pp_contact.radius_in_mm,
+            geometry.groove.radius_in_mm.inner,
+        ]
+        z += [
+            geometry.pp_contact.depth_in_mm,
+            geometry.pp_contact.depth_in_mm,
+            0,
+            0,
+        ]
+        surfaces += ["pplus", "passive", "passive"]
+
+    elif geometry.pp_contact.radius_in_mm < geometry.groove.radius_in_mm.inner:
+        r += [
+            0,
+            geometry.pp_contact.radius_in_mm,
+            geometry.groove.radius_in_mm.inner,
+        ]
+        z += [0, 0, 0]
+        surfaces += ["pplus", "passive"]
+    else:
+        r += [0, geometry.pp_contact.radius_in_mm]
+        z += [0, 0]
+        surfaces += ["pplus"]
+
+    r += [
+        geometry.groove.radius_in_mm.inner,
+        geometry.groove.radius_in_mm.outer,
+        geometry.groove.radius_in_mm.outer,
+    ]
+
+    z += [
+        geometry.groove.depth_in_mm,
+        geometry.groove.depth_in_mm,
+        0,
+    ]
+    surfaces += ["passive", "passive", "passive"]
+
+    return (r, z, surfaces)
diff --git a/src/legendhpges/draw.py b/src/legendhpges/draw.py
index 5cec622..50223ef 100644
--- a/src/legendhpges/draw.py
+++ b/src/legendhpges/draw.py
@@ -1,17 +1,14 @@
 from __future__ import annotations
 
-import logging
-
 import matplotlib.pyplot as plt
+import numpy as np
 from pyg4ometry.visualisation import VtkViewer
 
 from .base import HPGe
-from .p00664b import P00664B
-from .v02160a import V02160A
 
 
 def plot_profile(
-    hpge: HPGe, axes: plt.Axes = None, **kwargs
+    hpge: HPGe, axes: plt.Axes = None, split_by_type: bool = False, **kwargs
 ) -> tuple[plt.Figure, plt.Axes]:
     """Plot the HPGe profile with :mod:`matplotlib`.
 
@@ -21,22 +18,16 @@ def plot_profile(
         detector.
     axes
         pre-existing axes where the profile will be plotted.
+    split_by_type
+        boolean to separate surfaces of different types.
     **kwargs
         any keyword argument supported by :func:`matplotlib.pyplot.plot`.
 
     """
     # data
-    if isinstance(hpge, (V02160A, P00664B)):
-        r = hpge.solid.obj1.pR
-        z = hpge.solid.obj1.pZ
-        logging.warning("The detector profile is that of the solid without cut")
-    else:
-        r = hpge.solid.pR
-        z = hpge.solid.pZ
-
-    x = r + [-x for x in reversed(r)]
-    y = z + list(reversed(z))
+    r, z = hpge.get_profile()
 
+    # set options
     colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
 
     fig = None
@@ -50,15 +41,41 @@ def plot_profile(
 
     default_kwargs = {
         "marker": "o",
-        "markersize": 3,
+        "markersize": 2,
         "markeredgecolor": colors[1],
         "markerfacecolor": colors[1],
         "linewidth": 2,
     }
     default_kwargs |= kwargs
 
-    axes.plot(x, y, **default_kwargs)
+    if not split_by_type:
+        x = r + [-x for x in reversed(r)]
+        y = z + list(reversed(z))
 
+        axes.plot(x, y, **default_kwargs)
+    else:
+        surfaces = np.array(hpge.surfaces)
+        unique_surfaces = np.unique(surfaces)
+
+        dr = np.array([np.array([r1, r2]) for r1, r2 in zip(r[:-1], r[1:])])
+        dz = np.array([np.array([z1, z2]) for z1, z2 in zip(z[:-1], z[1:])])
+
+        for idx, u in enumerate(unique_surfaces):
+            drs_tmp = dr[surfaces == u]
+            dzs_tmp = dz[surfaces == u]
+
+            first = True
+            for r_tmp, z_tmp in zip(drs_tmp, dzs_tmp):
+                if first:
+                    axes.plot(
+                        r_tmp, z_tmp, color=colors[idx + 2], label=u, **default_kwargs
+                    )
+                    first = False
+                    axes.plot(-r_tmp, z_tmp, color=colors[idx + 2], **default_kwargs)
+                else:
+                    axes.plot(r_tmp, z_tmp, color=colors[idx + 2], **default_kwargs)
+                    axes.plot(-r_tmp, z_tmp, color=colors[idx + 2], **default_kwargs)
+        axes.legend(loc="upper right")
     return fig, axes
 
 
diff --git a/src/legendhpges/invcoax.py b/src/legendhpges/invcoax.py
index 7156bc6..66ffd4c 100644
--- a/src/legendhpges/invcoax.py
+++ b/src/legendhpges/invcoax.py
@@ -3,6 +3,7 @@
 import math
 
 from .base import HPGe
+from .build_utils import make_pplus
 
 
 class InvertedCoax(HPGe):
@@ -16,35 +17,12 @@ def _tan(a):
 
         r = []
         z = []
+        surfaces = []
 
-        if c.pp_contact.depth_in_mm > 0:
-            r += [
-                0,
-                c.pp_contact.radius_in_mm,
-                c.pp_contact.radius_in_mm,
-            ]
-            z += [
-                c.pp_contact.depth_in_mm,
-                c.pp_contact.depth_in_mm,
-                0,
-            ]
-        else:
-            r += [0]
-            z += [0]
-
-        r += [
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.outer,
-            c.groove.radius_in_mm.outer,
-        ]
-
-        z += [
-            0,
-            c.groove.depth_in_mm,
-            c.groove.depth_in_mm,
-            0,
-        ]
+        r_p, z_p, surface_p = make_pplus(c)
+        r += r_p
+        z += z_p
+        surfaces += surface_p
 
         if c.taper.bottom.height_in_mm > 0:
             r += [
@@ -57,9 +35,12 @@ def _tan(a):
                 0,
                 c.taper.bottom.height_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
+
         else:
             r += [c.radius_in_mm]
             z += [0]
+            surfaces += ["nplus"]
 
         if c.taper.top.height_in_mm > 0:
             r += [
@@ -72,9 +53,12 @@ def _tan(a):
                 c.height_in_mm - c.taper.top.height_in_mm,
                 c.height_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
+
         else:
             r += [c.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         if c.taper.borehole.height_in_mm > 0:
             r += [
@@ -87,9 +71,12 @@ def _tan(a):
                 c.height_in_mm,
                 c.height_in_mm - c.taper.borehole.height_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
+
         else:
             r += [c.borehole.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         if c.taper.borehole.height_in_mm != c.borehole.depth_in_mm:
             r += [
@@ -101,9 +88,14 @@ def _tan(a):
                 c.height_in_mm - c.borehole.depth_in_mm,
                 c.height_in_mm - c.borehole.depth_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
+
         else:
             r += [0]
 
             z += [c.height_in_mm - c.borehole.depth_in_mm]
+            surfaces += ["nplus"]
+
+        self.surfaces = surfaces
 
         return r, z
diff --git a/src/legendhpges/p00664b.py b/src/legendhpges/p00664b.py
index f96c370..534e863 100644
--- a/src/legendhpges/p00664b.py
+++ b/src/legendhpges/p00664b.py
@@ -56,13 +56,17 @@ def _tan(a):
 
         r = []
         z = []
+        surfaces = []
 
         if c.pp_contact.depth_in_mm > 0:
             r += [0, c.pp_contact.radius_in_mm, c.pp_contact.radius_in_mm]
             z += [c.pp_contact.depth_in_mm, c.pp_contact.depth_in_mm, 0]
+            surfaces += ["pplus", "passive"]
+
         else:
-            r += [0]
-            z += [0]
+            r += [0, c.pp_contact.radius_in_mm]
+            z += [0, 0]
+            surfaces += ["pplus"]
 
         if c.taper.bottom.height_in_mm > 0:
             r += [
@@ -71,9 +75,12 @@ def _tan(a):
                 c.radius_in_mm,
             ]
             z += [0, c.taper.bottom.height_in_mm]
+            surfaces += ["passive", "nplus"]
+
         else:
             r += [c.radius_in_mm]
             z += [0]
+            surfaces += ["passive"]
 
         if c.taper.top.height_in_mm > 0:
             r += [
@@ -82,13 +89,17 @@ def _tan(a):
                 - c.taper.top.height_in_mm * _tan(c.taper.top.angle_in_deg),
             ]
             z += [c.height_in_mm - c.taper.top.height_in_mm, c.height_in_mm]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         r += [0]
         z += [c.height_in_mm]
+        surfaces += ["nplus"]
 
+        self.surfaces = surfaces
         return r, z
 
     @property
diff --git a/src/legendhpges/ppc.py b/src/legendhpges/ppc.py
index a77deba..617b478 100644
--- a/src/legendhpges/ppc.py
+++ b/src/legendhpges/ppc.py
@@ -16,13 +16,17 @@ def _tan(a):
 
         r = []
         z = []
+        surfaces = []
 
         if c.pp_contact.depth_in_mm > 0:
             r += [0, c.pp_contact.radius_in_mm, c.pp_contact.radius_in_mm]
             z += [c.pp_contact.depth_in_mm, c.pp_contact.depth_in_mm, 0]
+            surfaces += ["pplus", "passive"]
+
         else:
-            r += [0]
-            z += [0]
+            r += [0, c.pp_contact.radius_in_mm]
+            z += [0, 0]
+            surfaces += ["pplus"]
 
         if c.taper.bottom.height_in_mm > 0:
             r += [
@@ -31,9 +35,11 @@ def _tan(a):
                 c.radius_in_mm,
             ]
             z += [0, c.taper.bottom.height_in_mm]
+            surfaces += ["passive", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [0]
+            surfaces += ["passive"]
 
         if c.taper.top.height_in_mm > 0:
             r += [
@@ -42,11 +48,16 @@ def _tan(a):
                 - c.taper.top.height_in_mm * _tan(c.taper.top.angle_in_deg),
             ]
             z += [c.height_in_mm - c.taper.top.height_in_mm, c.height_in_mm]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         r += [0]
         z += [c.height_in_mm]
+        surfaces += ["nplus"]
+
+        self.surfaces = surfaces
 
         return r, z
diff --git a/src/legendhpges/semicoax.py b/src/legendhpges/semicoax.py
index 2586d61..f73fe6d 100644
--- a/src/legendhpges/semicoax.py
+++ b/src/legendhpges/semicoax.py
@@ -16,10 +16,11 @@ def _tan(a):
 
         r = []
         z = []
+        surfaces = []
 
         r += [0, c.borehole.radius_in_mm]
-
         z += [c.borehole.depth_in_mm, c.borehole.depth_in_mm]
+        surfaces += ["pplus"]
 
         if c.taper.borehole.height_in_mm > 0:
             r += [
@@ -28,9 +29,11 @@ def _tan(a):
                 + c.taper.borehole.height_in_mm * _tan(c.taper.borehole.angle_in_deg),
             ]
             z += [c.taper.borehole.height_in_mm, 0]
+            surfaces += ["pplus", "pplus"]
         else:
             r += [c.borehole.radius_in_mm]
             z += [0]
+            surfaces += ["pplus"]
 
         r += [
             c.groove.radius_in_mm.inner,
@@ -40,6 +43,7 @@ def _tan(a):
         ]
 
         z += [0, c.groove.depth_in_mm, c.groove.depth_in_mm, 0]
+        surfaces += ["pplus", "passive", "passive", "passive"]
 
         if c.taper.bottom.height_in_mm > 0:
             r += [
@@ -48,9 +52,11 @@ def _tan(a):
                 c.radius_in_mm,
             ]
             z += [0, c.taper.bottom.height_in_mm]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [0]
+            surfaces += ["nplus"]
 
         if c.taper.top.height_in_mm > 0:
             r += [
@@ -59,11 +65,15 @@ def _tan(a):
                 - c.taper.top.height_in_mm * _tan(c.taper.top.angle_in_deg),
             ]
             z += [c.height_in_mm - c.taper.top.height_in_mm, c.height_in_mm]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         r += [0]
         z += [c.height_in_mm]
+        surfaces += ["nplus"]
 
+        self.surfaces = surfaces
         return r, z
diff --git a/src/legendhpges/utils.py b/src/legendhpges/utils.py
new file mode 100644
index 0000000..2712cf8
--- /dev/null
+++ b/src/legendhpges/utils.py
@@ -0,0 +1,221 @@
+from __future__ import annotations
+
+import numba
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+
+@numba.njit(cache=True)
+def convert_coords(coords: ArrayLike) -> NDArray:
+    """Converts (x,y,z) coordinates into (r,z)
+
+    Parameters
+    ----------
+    coords
+        numpy array of coordinates where the second index corresponds to (x,y,z) respectively
+
+    Returns
+    -------
+        numpy array of (r,z) coordinates for each point
+
+    """
+    r = np.sqrt(coords[:, 0] ** 2 + coords[:, 1] ** 2)
+    return np.column_stack((r, coords[:, 2]))
+
+
+def shortest_distance_to_plane(
+    a_vec: NDArray,
+    d: float,
+    points: NDArray,
+    rmax: float | None = None,
+    zrange: tuple[float, float] | None = None,
+) -> NDArray:
+    """Get the shortest distance from a plane (constrained in r and z) to each point.
+
+    The equation of the plane is given by :math:`a_1x+a_2y+a_3z=d`. Where
+    :math:`\\vec{a}=(a_1,a_2,a_3)`.
+    The closest point on the plane to the point (:math:`y`) is then given by:
+
+    .. math::
+        x =y-(y*a-d)*a/||a||^2
+
+    The distance is then given by the length of the vector :math:`x-y`. This
+    function also checks if the intersection point is above `rmax` or inside the
+    `zrange`, if not, ``numpy.nan`` is returned for that point.
+
+    Parameters
+    ----------
+    a_vec
+        3 coordinate array defining the plane.
+    d
+        scalar in the plane definition.
+    points
+        set of points to compute distance.
+    rmax
+        maximum radius for the plane.
+    zrange
+        range in z for the plane.
+    """
+
+    def _dot(a, b):
+        return np.sum(a * b, axis=1)
+
+    def _norm(a):
+        ax = 1 if a.ndim == 2 else 0
+        return np.sqrt(np.sum(a**2, axis=ax))
+
+    a_norm = np.sum(a_vec * a_vec)
+
+    proj_points = points - ((_dot(points, a_vec) - d)[:, np.newaxis]) * a_vec / a_norm
+
+    dist_vec = _norm((_dot(points, a_vec) - d)[:, np.newaxis] * a_vec / a_norm)
+
+    # check on r and z
+
+    proj_points_rz = convert_coords(proj_points)
+
+    if rmax is not None:
+        condition_r = proj_points_rz[:, 0] <= rmax
+    else:
+        condition_r = np.full(len(points), True)
+
+    if zrange is not None:
+        condition_z = (proj_points_rz[:, 1] > zrange[0]) & (
+            proj_points_rz[:, 1] < zrange[1]
+        )
+    else:
+        condition_z = np.full(len(points), True)
+
+    condition = condition_r * condition_z
+    return np.where(condition, dist_vec, np.full(len(points), np.nan))
+
+
+def get_line_segments(
+    r: ArrayLike, z: ArrayLike, surface_indices: ArrayLike = None
+) -> tuple[NDArray, NDArray]:
+    """Extracts the line segments from a shape.
+
+    Parameters
+    ---------
+    r
+        array or list of radial positions defining the polycone.
+    z
+        array or list of vertical positions defining the polycone.
+    surface_indices
+        list of integer indices of surfaces to consider. If ``None`` (the
+        default) all surfaces used.
+
+    Returns
+    -------
+        tuple of (s1,s2) arrays describing the line segments, both `s1` and
+        `s2` have shape `(n_segments,2)` where the first axis represents thhe
+        segment and the second `(r,z)`.
+    """
+    # build lists of pairs of coordinates
+    s1 = np.array([np.array([r1, z1]) for r1, z1 in zip(r[:-1], z[:-1])])
+    s2 = np.array([np.array([r2, z2]) for r2, z2 in zip(r[1:], z[1:])])
+
+    if surface_indices is not None:
+        s1 = s1[surface_indices]
+        s2 = s2[surface_indices]
+
+    return s1, s2
+
+
+@numba.njit(cache=True)
+def shortest_distance(
+    s1_list: NDArray,
+    s2_list: NDArray,
+    points: NDArray,
+    tol: float = 1e-11,
+    signed: bool = True,
+) -> tuple[NDArray, NDArray]:
+    """Get the shortest distance between each point and a line segment.
+
+    Based on vector algebra where the distance vector is given by:
+
+    .. math::
+        d = s_1 - p - ( (n · (s_1- p)) * n )
+
+    where:
+
+    - :math:`s_1` is a vector from which the distance is measured,
+    - `p` is a point vector,
+    - `n` is a unit direction vector from :math:`s_1` to :math:`s_2`,
+    - `a` is another point vector.
+
+    If the projection point lies inside the segment s1-s2. Else the closest
+    point is either :math:`s_1` or :math:`s_2`.  The distance is the modulus of
+    this vector and this calculation is performed for each point.  A sign is
+    attached based on the cross product of the line vector and the distance
+    vector.  To avoid numerical issues any point within the tolerance is
+    considered inside.
+
+    Parameters
+    ----------
+    s1_list
+        `(n_segments,2)` np.array of the first points in the line segment, for
+        the second axis indices `0,1` correspond to `r,z`.
+    s2_list
+        second points, same format as `s1_list`.
+    points
+        `(n_points,2)` array of points to compare, first axis corresponds to
+        the point index and the second to `(r,z)`.
+    tol
+        tolerance when computing sign, points within this distance to the
+        surface are pushed inside.
+    signed
+        boolean flag to attach a sign to the distance (positive if inside).
+
+    Returns
+    -------
+        ``(n_points,n_segments)`` numpy array of the shortest distances for each segment.
+    """
+
+    # helper functions
+    def _dot(a, b):
+        return np.sum(a * b, axis=1)
+
+    def _norm(a):
+        ax = 1 if a.ndim == 2 else 0
+        return np.sqrt(np.sum(a**2, axis=ax))
+
+    n_segments = len(s1_list)
+    dists = np.full((len(points), len(s1_list)), np.nan)
+
+    for segment in range(n_segments):
+        s1 = s1_list[segment]
+        s2 = s2_list[segment]
+
+        n = (s2 - s1) / _norm(s2 - s1)
+
+        proj_dist = -_dot(n, (n * _dot(s1 - points, n)[:, np.newaxis]))
+
+        dist_vec = np.empty_like(s1 - points)
+
+        condition1 = proj_dist < 0
+        condition2 = proj_dist > _norm(s2 - s1)
+        condition3 = (~condition1) & (~condition2)
+
+        diff_s1 = s1 - points
+        dist_vec[condition1] = diff_s1[condition1]
+        dist_vec[condition2] = s2 - points[condition2]
+        dist_vec[condition3] = (
+            diff_s1[condition3] - n * _dot(diff_s1, n)[condition3, np.newaxis]
+        )
+
+        # make this signed so inside is positive and outside negative
+        if signed:
+            sign_vec = n[0] * dist_vec[:, 1] - n[1] * dist_vec[:, 0]
+
+            # push points on surface inside
+            sign_vec = np.where(np.abs(sign_vec) < tol, -tol, sign_vec)
+            sign_vec_norm = -sign_vec / np.abs(sign_vec)
+
+        else:
+            sign_vec_norm = np.ones(len(dist_vec))
+
+        dists[:, segment] = np.where(
+            np.abs(_norm(dist_vec)) < tol, tol, np.abs(_norm(dist_vec)) * sign_vec_norm
+        )
+    return dists
diff --git a/src/legendhpges/v02160a.py b/src/legendhpges/v02160a.py
index 57bcea4..8f969cd 100644
--- a/src/legendhpges/v02160a.py
+++ b/src/legendhpges/v02160a.py
@@ -5,6 +5,7 @@
 from pyg4ometry import geant4
 
 from .base import HPGe
+from .build_utils import make_pplus
 from .registry import default_units_registry as u
 
 
@@ -59,35 +60,12 @@ def _tan(a):
 
         r = []
         z = []
+        surfaces = []
 
-        if c.pp_contact.depth_in_mm > 0:
-            r += [
-                0,
-                c.pp_contact.radius_in_mm,
-                c.pp_contact.radius_in_mm,
-            ]
-            z += [
-                c.pp_contact.depth_in_mm,
-                c.pp_contact.depth_in_mm,
-                0,
-            ]
-        else:
-            r += [0]
-            z += [0]
-
-        r += [
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.outer,
-            c.groove.radius_in_mm.outer,
-        ]
-
-        z += [
-            0,
-            c.groove.depth_in_mm,
-            c.groove.depth_in_mm,
-            0,
-        ]
+        r_p, z_p, surface_p = make_pplus(c)
+        r += r_p
+        z += z_p
+        surfaces += surface_p
 
         if c.taper.bottom.height_in_mm > 0:
             r += [
@@ -100,9 +78,12 @@ def _tan(a):
                 0,
                 c.taper.bottom.height_in_mm,
             ]
+
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [0]
+            surfaces += ["nplus"]
 
         if c.taper.top.height_in_mm > 0:
             r += [
@@ -115,9 +96,12 @@ def _tan(a):
                 c.height_in_mm - c.taper.top.height_in_mm,
                 c.height_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
+
         else:
             r += [c.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         if c.taper.borehole.height_in_mm > 0:
             r += [
@@ -130,9 +114,11 @@ def _tan(a):
                 c.height_in_mm,
                 c.height_in_mm - c.taper.borehole.height_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.borehole.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         if c.taper.borehole.height_in_mm != c.borehole.depth_in_mm:
             r += [
@@ -144,11 +130,16 @@ def _tan(a):
                 c.height_in_mm - c.borehole.depth_in_mm,
                 c.height_in_mm - c.borehole.depth_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [0]
 
             z += [c.height_in_mm - c.borehole.depth_in_mm]
 
+            surfaces += ["nplus"]
+
+        self.surfaces = surfaces
+
         return r, z
 
     @property
diff --git a/src/legendhpges/v02162b.py b/src/legendhpges/v02162b.py
index 3d85899..a006f45 100644
--- a/src/legendhpges/v02162b.py
+++ b/src/legendhpges/v02162b.py
@@ -3,6 +3,7 @@
 import math
 
 from .base import HPGe
+from .build_utils import make_pplus
 
 
 class V02162B(HPGe):
@@ -16,35 +17,12 @@ def _tan(a):
 
         r = []
         z = []
+        surfaces = []
 
-        if c.pp_contact.depth_in_mm > 0:
-            r += [
-                0,
-                c.pp_contact.radius_in_mm,
-                c.pp_contact.radius_in_mm,
-            ]
-            z += [
-                c.pp_contact.depth_in_mm,
-                c.pp_contact.depth_in_mm,
-                0,
-            ]
-        else:
-            r += [0]
-            z += [0]
-
-        r += [
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.outer,
-            c.groove.radius_in_mm.outer,
-        ]
-
-        z += [
-            0,
-            c.groove.depth_in_mm,
-            c.groove.depth_in_mm,
-            0,
-        ]
+        r_p, z_p, surface_p = make_pplus(c)
+        r += r_p
+        z += z_p
+        surfaces += surface_p
 
         if c.taper.bottom.height_in_mm > 0:
             r += [
@@ -57,9 +35,12 @@ def _tan(a):
                 0,
                 c.taper.bottom.height_in_mm,
             ]
+
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [0]
+            surfaces += ["nplus"]
 
         if c.taper.top.height_in_mm > 0:
             r += [
@@ -72,15 +53,18 @@ def _tan(a):
                 c.height_in_mm - c.taper.top.height_in_mm,
                 c.height_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [c.height_in_mm]
-
+            surfaces += ["nplus"]
         # top groove
         r += [c.extra.topgroove.radius_in_mm, c.extra.topgroove.radius_in_mm]
 
         z += [c.height_in_mm, c.height_in_mm - c.extra.topgroove.depth_in_mm]
 
+        surfaces += ["nplus", "nplus"]
+
         # borehole
         r += [c.borehole.radius_in_mm, c.borehole.radius_in_mm, 0]
 
@@ -89,5 +73,8 @@ def _tan(a):
             c.height_in_mm - c.borehole.depth_in_mm,
             c.height_in_mm - c.borehole.depth_in_mm,
         ]
+        surfaces += ["nplus", "nplus", "nplus"]
+
+        self.surfaces = surfaces
 
         return r, z
diff --git a/src/legendhpges/v07646a.py b/src/legendhpges/v07646a.py
index f4ab64c..afc6a40 100644
--- a/src/legendhpges/v07646a.py
+++ b/src/legendhpges/v07646a.py
@@ -3,6 +3,7 @@
 import math
 
 from .base import HPGe
+from .build_utils import make_pplus
 
 
 class V07646A(HPGe):
@@ -16,35 +17,12 @@ def _tan(a):
 
         r = []
         z = []
+        surfaces = []
 
-        if c.pp_contact.depth_in_mm > 0:
-            r += [
-                0,
-                c.pp_contact.radius_in_mm,
-                c.pp_contact.radius_in_mm,
-            ]
-            z += [
-                c.pp_contact.depth_in_mm,
-                c.pp_contact.depth_in_mm,
-                0,
-            ]
-        else:
-            r += [0]
-            z += [0]
-
-        r += [
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.inner,
-            c.groove.radius_in_mm.outer,
-            c.groove.radius_in_mm.outer,
-        ]
-
-        z += [
-            0,
-            c.groove.depth_in_mm,
-            c.groove.depth_in_mm,
-            0,
-        ]
+        r_p, z_p, surface_p = make_pplus(c)
+        r += r_p
+        z += z_p
+        surfaces += surface_p
 
         bottom_cylinder = c.extra.bottom_cylinder
 
@@ -59,15 +37,18 @@ def _tan(a):
                 0,
                 c.taper.bottom.height_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [bottom_cylinder.radius_in_mm]
             z += [0]
+            surfaces += ["nplus"]
 
         r += [bottom_cylinder.radius_in_mm, c.radius_in_mm]
         z += [
             bottom_cylinder.height_in_mm,
             bottom_cylinder.height_in_mm + bottom_cylinder.transition_in_mm,
         ]
+        surfaces += ["nplus", "nplus"]
 
         if c.taper.top.height_in_mm > 0:
             r += [
@@ -80,9 +61,11 @@ def _tan(a):
                 c.height_in_mm - c.taper.top.height_in_mm,
                 c.height_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
         else:
             r += [c.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         if c.taper.borehole.height_in_mm > 0:
             r += [
@@ -92,9 +75,12 @@ def _tan(a):
             ]
 
             z += [c.height_in_mm, c.height_in_mm - c.taper.borehole.height_in_mm]
+            surfaces += ["nplus", "nplus"]
+
         else:
             r += [c.borehole.radius_in_mm]
             z += [c.height_in_mm]
+            surfaces += ["nplus"]
 
         if c.taper.borehole.height_in_mm != c.borehole.depth_in_mm:
             r += [
@@ -106,9 +92,15 @@ def _tan(a):
                 c.height_in_mm - c.borehole.depth_in_mm,
                 c.height_in_mm - c.borehole.depth_in_mm,
             ]
+            surfaces += ["nplus", "nplus"]
+
         else:
             r += [0]
 
             z += [c.height_in_mm - c.borehole.depth_in_mm]
 
+            surfaces += ["nplus"]
+
+        self.surfaces = surfaces
+
         return r, z
diff --git a/tests/test_det_profile.py b/tests/test_det_profile.py
index 12ab996..0243d92 100644
--- a/tests/test_det_profile.py
+++ b/tests/test_det_profile.py
@@ -72,44 +72,60 @@ def test_make_icpc(test_data_configs):
     gedet = make_hpge(test_data_configs + "/V99000A.json")
     assert isinstance(gedet, InvertedCoax)
 
+    assert len(gedet._decode_polycone_coord()[0]) == len(gedet.surfaces) + 1
+
 
 def test_make_bege(test_data_configs):
     gedet = make_hpge(test_data_configs + "/B99000A.json")
     assert isinstance(gedet, BEGe)
 
+    assert len(gedet._decode_polycone_coord()[0]) == len(gedet.surfaces) + 1
+
 
 def test_make_ppc(test_data_configs):
     gedet = make_hpge(test_data_configs + "/P99000A.json")
     assert isinstance(gedet, PPC)
 
+    assert len(gedet._decode_polycone_coord()[0]) == len(gedet.surfaces) + 1
+
 
 def test_make_semicoax(test_data_configs):
     gedet = make_hpge(test_data_configs + "/C99000A.json")
     assert isinstance(gedet, SemiCoax)
 
+    assert len(gedet._decode_polycone_coord()[0]) == len(gedet.surfaces) + 1
+
 
 def make_v07646a():
     gedet = make_hpge(configs.V07646A)
     assert isinstance(gedet, V07646A)
 
+    assert len(gedet._decode_polycone_coord()[0]) == len(gedet.surfaces) + 1
+
 
 def test_make_p00664b():
     gedet = make_hpge(configs.P00664B)
     assert gedet.mass
     assert isinstance(gedet, P00664B)
 
+    assert len(gedet._decode_polycone_coord()[0]) == len(gedet.surfaces) + 1
+
 
 def test_make_v02162b():
     gedet = make_hpge(configs.V02162B)
     assert gedet.mass
     assert isinstance(gedet, V02162B)
 
+    assert len(gedet._decode_polycone_coord()[0]) == len(gedet.surfaces) + 1
+
 
 def test_make_v02160a():
     gedet = make_hpge(configs.V02160A)
     assert gedet.mass
     assert isinstance(gedet, V02160A)
 
+    assert len(gedet._decode_polycone_coord()[0]) == len(gedet.surfaces) + 1
+
 
 def test_null_enrichment():
     metadata = configs.V07646A
diff --git a/tests/test_det_properties.py b/tests/test_det_properties.py
new file mode 100644
index 0000000..e2daa77
--- /dev/null
+++ b/tests/test_det_properties.py
@@ -0,0 +1,33 @@
+from __future__ import annotations
+
+import pathlib
+
+import numpy as np
+import pytest
+from legendmeta import TextDB
+from legendtestdata import LegendTestData
+from pyg4ometry import geant4
+
+from legendhpges import (
+    make_hpge,
+)
+from legendhpges.materials import make_natural_germanium
+
+reg = geant4.Registry()
+natural_germanium = make_natural_germanium(reg)
+configs = TextDB(pathlib.Path(__file__).parent.resolve() / "configs")
+
+
+@pytest.fixture(scope="session")
+def test_data_configs():
+    ldata = LegendTestData()
+    ldata.checkout("5f9b368")
+    return ldata.get_path("legend/metadata/hardware/detectors/germanium/diodes")
+
+
+def test_surface_area(test_data_configs):
+    reg = geant4.Registry()
+    gedet = make_hpge(test_data_configs + "/C99000A.json", registry=reg)
+
+    assert len(gedet.surface_area(surface_indices=[])) == 0
+    assert np.sum(gedet.surface_area(surface_indices=None)) > 0
diff --git a/tests/test_distance.py b/tests/test_distance.py
new file mode 100644
index 0000000..0f790c8
--- /dev/null
+++ b/tests/test_distance.py
@@ -0,0 +1,111 @@
+from __future__ import annotations
+
+import pathlib
+
+import numpy as np
+import pytest
+from legendmeta import TextDB
+from legendtestdata import LegendTestData
+from pyg4ometry import geant4
+
+from legendhpges import (
+    make_hpge,
+)
+from legendhpges.materials import make_natural_germanium
+
+reg = geant4.Registry()
+natural_germanium = make_natural_germanium(reg)
+configs = TextDB(pathlib.Path(__file__).parent.resolve() / "configs")
+
+
+@pytest.fixture(scope="session")
+def test_data_configs():
+    ldata = LegendTestData()
+    ldata.checkout("2553a28")
+    return ldata.get_path("legend/metadata/hardware/detectors/germanium/diodes")
+
+
+def test_not_implemented():
+    reg = geant4.Registry()
+    ppc = make_hpge(configs.P00664B, registry=reg)
+
+    with pytest.raises(NotImplementedError):
+        ppc.distance_to_surface([[1, 0, 0]])
+
+
+def test_bad_dimensions(test_data_configs):
+    reg = geant4.Registry()
+    gedet = make_hpge(test_data_configs + "/C99000A.json", registry=reg)
+
+    with pytest.raises(ValueError):
+        gedet.distance_to_surface([[1, 0, 0, 0]])
+
+
+def test_output(test_data_configs):
+    reg = geant4.Registry()
+    gedet = make_hpge(test_data_configs + "/C99000A.json", registry=reg)
+    dist = gedet.distance_to_surface([[0, 0, 0], [1, 3, 3], [0, 0, 0]])
+
+    assert np.shape(dist == (3,))
+    assert np.all(dist >= 0)
+
+    dist_indices = gedet.distance_to_surface(
+        [[0, 0, 0], [1, 3, 3], [0, 0, 0]], surface_indices=[0, 3]
+    )
+    assert np.all(dist_indices >= dist)
+
+
+def test_inside_not_implemented():
+    reg = geant4.Registry()
+    ppc = make_hpge(configs.P00664B, registry=reg)
+
+    with pytest.raises(NotImplementedError):
+        ppc.is_inside([[1, 0, 0]])
+
+
+def test_inside_bad_dimensions(test_data_configs):
+    reg = geant4.Registry()
+    gedet = make_hpge(test_data_configs + "/C99000A.json", registry=reg)
+
+    with pytest.raises(ValueError):
+        gedet.is_inside([[1, 0, 0, 0]])
+
+
+def test_inside_output(test_data_configs):
+    reg = geant4.Registry()
+
+    gedet = make_hpge(test_data_configs + "/B99000A.json", registry=reg)
+    r, z = gedet._decode_polycone_coord()
+
+    # detetor is a simple bege
+    # p+ at 0-7.5 mm
+    # groove at 10-12 mm 2mm deep
+    # radius and height 40mm
+
+    theta = 33
+    cos = np.cos(np.deg2rad(theta))
+    sin = np.sin(np.deg2rad(theta))
+
+    # 1) on axis inside
+    # 2) below outside
+    # 3) inside groove
+    # 4) above groove
+    # 5) outside side
+    # 6) exactly on the side
+    # 7) far above top
+
+    coords = np.array(
+        [
+            [0, 0, 10],
+            [5 * cos, 5 * sin, -0.1],
+            [11 * cos, 11 * sin, 1],
+            [11 * cos, 11 * sin, 4],
+            [50 * cos, 50 * sin, 10],
+            [40 * cos, 40 * sin, 10],
+            [20, 20, 1000],
+        ]
+    )
+
+    is_in = gedet.is_inside(coords, tol=1e-11)
+
+    assert np.all(is_in == np.array([True, False, False, True, False, True, False]))
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 0000000..442c031
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,117 @@
+from __future__ import annotations
+
+import numpy as np
+
+from legendhpges import utils
+
+
+def test_distances():
+    # test conversion of coordinates
+
+    coords = np.array([[1, 1, 1], [0, 0, 0]])
+    coords_rz = utils.convert_coords(coords)
+
+    assert coords_rz.ndim == 2
+    assert np.shape(coords_rz) == (2, 2)
+    assert np.allclose(coords_rz[1], np.array([0, 0]))
+
+    # test shortest distance
+    # in all these
+    s1 = np.array([[0, 0]])
+    s2 = np.array([[1, 0]])
+
+    # first point on segment (distance is 0)
+    # second directly above start (distance is 5)
+    # third above the segment (distance is 7)
+    # fourth outside the segment by 3 units in x and 4 in y (distance is 5)
+    # last the same but for the first point
+
+    points = np.array([[0.5, 0], [0, 5], [0.3, 7], [4, 4], [-3, 4]])
+    res = np.array([0.0, 5.0, 7.0, 5.0, 5.0])
+    assert np.allclose(utils.shortest_distance(s1, s2, points)[:, 0], res)
+
+    # all distances shouldn't be affected by a global offset and rotation
+    offset = np.array([107, -203])
+    rot = np.array(
+        [
+            [np.cos(np.deg2rad(37)), -np.sin(np.deg2rad(37))],
+            [np.sin(np.deg2rad(37)), np.cos(np.deg2rad(37))],
+        ]
+    )
+
+    points_new = np.array([rot @ (p_tmp + offset) for p_tmp in points])
+    s1_new = np.array([rot @ (s1_t + offset) for s1_t in s1])
+    s2_new = np.array([rot @ (s2_t + offset) for s2_t in s2])
+    assert np.allclose(
+        utils.shortest_distance(s1_new, s2_new, points_new)[:, 0],
+        res,
+    )
+
+
+def test_plane_distance_unconstrained():
+    # start with a plane on the x,y plane
+    a = np.array([0, 0, 1])
+    d = 0
+    points = np.array([[0, 0, 1], [1, 1, 7], [-5, 2, 9], [1, 0, 0]])
+    assert np.allclose(
+        utils.shortest_distance_to_plane(a, d, points), np.array([1, 7, 9, 0])
+    )
+
+    # a is the normal vector while d=p*a where p is a point on the plane i.e 0
+
+    offset = np.array([107, -203, 197])
+    d_new = np.sum(a * offset)
+
+    assert np.allclose(
+        utils.shortest_distance_to_plane(a, d_new, points + offset),
+        np.array([1, 7, 9, 0]),
+    )
+
+    # now arbitrary rotation
+
+    alpha = 35
+    beta = 18
+    gamma = 48
+
+    ca = np.cos(np.deg2rad(alpha))
+    cb = np.cos(np.deg2rad(beta))
+    cg = np.cos(np.deg2rad(gamma))
+
+    sa = np.sin(np.deg2rad(alpha))
+    sb = np.sin(np.deg2rad(beta))
+    sg = np.sin(np.deg2rad(gamma))
+
+    rot = np.array(
+        [
+            [cb * cg, sa * sb * cg - ca * sg, ca * sb * cg + sa * sg],
+            [cb * sg, sa * sb * sg + ca * cg, ca * sb * sg - sa * cg],
+            [-sb, sa * cb, ca * cb],
+        ]
+    )
+    p_new = rot @ offset
+    points_new = np.array([rot @ (p_tmp + offset) for p_tmp in points])
+    a_new = rot @ a
+    d_new = np.sum(p_new * a_new)
+
+    assert np.allclose(
+        utils.shortest_distance_to_plane(a_new, d_new, points_new),
+        np.array([1, 7, 9, 0]),
+    )
+
+
+def test_plane_distance_constrained():
+    # vertical plane at x=3
+    a = np.array([1, 0, 0])
+    d = 3
+    points = np.array([[0, 0, 10], [2, 0, -7], [3, 5, -17], [2.9, 4, -5]])
+
+    assert np.allclose(
+        utils.shortest_distance_to_plane(a, d, points, rmax=5),
+        np.array([3, 1, np.nan, 0.1]),
+        equal_nan=True,
+    )
+    assert np.allclose(
+        utils.shortest_distance_to_plane(a, d, points, rmax=5, zrange=[-100, 1]),
+        np.array([np.nan, 1, np.nan, 0.1]),
+        equal_nan=True,
+    )