Merge pull request #51 from tdixon97/main

implement distance to surface and add the surface types to HPGe object

docs/source/conf.py

@@ -61,3 +61,8 @@
 autodoc_typehints = "description"
 autodoc_typehints_description_target = "documented_params"
 autodoc_typehints_format = "short"
+
+autodoc_type_aliases = {
+    "ArrayLike": "ArrayLike",
+    "NDArray": "NDArray",
+}

pyproject.toml

@@ -31,6 +31,7 @@ classifiers = [
 requires-python = ">=3.9"
 dependencies = [
     "numpy",
+    "numba",
     "pint != 0.24",
     "pyg4ometry",
     "pylegendmeta>=v0.9.0a2",

src/legendhpges/__init__.py

@@ -18,6 +18,7 @@
     "PPC",
     "SemiCoax",
     "make_hpge",
+    "utils",
     "V07646A",
     "P00664B",
     "V02160A",

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,23 +67,27 @@ def __init__(
 
         self.registry = registry
 
+        self.surfaces = []
+
         # build logical volume, default [mm]
         super().__init__(self._g4_solid(), material, self.name, self.registry)
 
     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,
+        )

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

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)