Merge pull request #55 from pim-n/feature-detectors-roadgen

Feature detectors roadgen

src/pg_rad/configs/__init__.py

@@ -1 +0,0 @@
-__all__ = []

src/pg_rad/configs/filepaths.py

@@ -1,3 +1,4 @@
 ATTENUATION_TABLE = 'attenuation_table.csv'
+ISOTOPE_TABLE = 'isotopes.csv'
 TEST_EXP_DATA = 'test_path_coords.csv'
 LOGGING_CONFIG = 'logging.yml'

src/pg_rad/data/__init__.py

@@ -1 +0,0 @@
-__all__ = []

src/pg_rad/data/angular_efficiencies/LU_HPGe_90.csv

@@ -0,0 +1,37 @@
+angle,662,1173,1332
+0,0.015,0.030,0.033
+10,0.011,0.021,0.024
+20,0.086,0.127,0.146
+30,0.294,0.356,0.397
+40,0.661,0.700,0.734
+50,1.054,1.057,1.057
+60,1.154,1.140,1.137
+70,1.186,1.152,1.138
+80,1.151,1.114,1.097
+90,1.000,1.000,1.000
+100,1.020,1.040,1.047
+110,1.074,1.093,1.103
+120,1.113,1.092,1.102
+130,1.139,1.122,1.113
+140,1.146,1.152,1.140
+150,1.113,1.118,1.104
+160,1.113,1.096,1.099
+170,1.091,1.076,1.083
+180,1.076,1.066,1.078
+-170,1.102,1.091,1.093
+-160,1.122,1.100,1.102
+-150,1.128,1.105,1.093
+-140,1.144,1.112,1.123
+-130,1.140,1.117,1.095
+-120,1.146,1.127,1.098
+-110,1.068,1.068,1.045
+-100,1.013,1.025,1.016
+-90,1.004,1.018,1.021
+-80,1.150,1.137,1.132
+-70,1.184,1.167,1.164
+-60,1.158,1.140,1.138
+-50,1.090,1.068,1.064
+-40,0.595,0.620,0.631
+-30,0.332,0.430,0.430
+-20,0.055,0.081,0.096
+-10,0.009,0.018,0.019

src/pg_rad/data/detectors.csv

@@ -0,0 +1,4 @@
+name,type,is_isotropic
+dummy,NaI,true
+LU_NaI_3inch,NaI,true
+LU_HPGe_90,HPGe,false

src/pg_rad/data/field_efficiencies/LU_HPGe_90.csv

@@ -0,0 +1,17 @@
+energy_keV,field_efficiency_m2,uncertainty
+59.5,0.00140,0.00005
+81.0,0.00310,0.00010
+122.1,0.00420,0.00013
+136.5,0.00428,0.00017
+160.6,0.00426,0.00030
+223.2,0.00418,0.00024
+276.4,0.00383,0.00012
+302.9,0.00370,0.00012
+356.0,0.00338,0.00010
+383.8,0.00323,0.00010
+511.0,0.00276,0.00008
+661.7,0.00241,0.00007
+834.8,0.00214,0.00007
+1173.2,0.00179,0.00005
+1274.5,0.00168,0.00005
+1332.5,0.00166,0.00005

src/pg_rad/data/field_efficiencies/LU_NaI_3inch.csv

@@ -0,0 +1,64 @@
+energy_keV,field_efficiency_m2
+10,5.50129E-09
+11.22018454,2.88553E-07
+12.58925412,4.81878E-06
+14.12537545,3.55112E-05
+15.84893192,0.000146367
+17.7827941,0.000397029
+19.95262315,0.000803336
+22.38721139,0.00131657
+25.11886432,0.001862377
+28.18382931,0.002373449
+31.6227766,0.002811046
+33.164,0.00269554
+33.164,0.002698792
+35.48133892,0.002509993
+39.81071706,0.002801304
+44.66835922,0.003015877
+50.11872336,0.003227431
+56.23413252,0.00341077
+63.09573445,0.003562051
+70.79457844,0.00368852
+79.43282347,0.003788875
+89,0.003867423
+100,0.003925025
+112.2018454,0.003967222
+125.8925412,0.003991551
+141.2537545,0.004000729
+158.4893192,0.003993145
+177.827941,0.003969163
+199.5262315,0.003925289
+223.8721139,0.003856247
+251.1886432,0.00375596
+281.8382931,0.003619634
+316.227766,0.003446087
+354.8133892,0.003242691
+398.1071706,0.003021761
+446.6835922,0.002791816
+501.1872336,0.002568349
+562.3413252,0.002350052
+630.9573445,0.002147662
+707.9457844,0.001957893
+794.3282347,0.001785694
+891,0.001626634
+1000,0.001482571
+1122.018454,0.00135047
+1258.925412,0.001231358
+1412.537545,0.001116695
+1584.893192,0.001011833
+1778.27941,0.000917017
+1995.262315,0.000828435
+2238.721139,0.000746854
+2511.886432,0.000672573
+2818.382931,0.00060493
+3162.27766,0.000544458
+3548.133892,0.000488446
+3981.071706,0.000438438
+4466.835922,0.000392416
+5011.872336,0.00035092
+5623.413252,0.000313959
+6309.573445,0.000279409
+7079.457844,0.000247794
+7943.282347,0.000218768
+8913,0.000190209
+10000,0.000164309

src/pg_rad/data/field_efficiencies/dummy.csv

@@ -0,0 +1,3 @@
+energy_keV,field_efficiency_m2
+0,1.0
+10000,1.0

src/pg_rad/data/isotopes.csv

@@ -0,0 +1,6 @@
+isotope,gamma_energy_keV,branching_ratio
+Cs134,604.7,0.976
+Cs134,795.9,0.855
+Cs137,661.657,0.851
+Co60,1173.228,1.0
+Co60,1332.5,1.0

src/pg_rad/dataloader/__init__.py

@@ -1,8 +0,0 @@
-# do not expose internal logger when running mkinit
-__ignore__ = ["logger"]
-
-from pg_rad.dataloader import dataloader
-
-from pg_rad.dataloader.dataloader import (load_data,)
-
-__all__ = ['dataloader', 'load_data']

src/pg_rad/detector/builder.py

@@ -1,20 +0,0 @@
-from pg_rad.inputparser.specs import DetectorSpec
-
-from .detectors import IsotropicDetector, AngularDetector
-
-
-class DetectorBuilder:
-    def __init__(
-        self,
-        detector_spec: DetectorSpec,
-    ):
-        self.detector_spec = detector_spec
-
-    def build(self) -> IsotropicDetector | AngularDetector:
-        if self.detector_spec.is_isotropic:
-            return IsotropicDetector(
-                self.detector_spec.name,
-                self.detector_spec.efficiency
-            )
-        else:
-            raise NotImplementedError("Angular detector not supported yet.")

src/pg_rad/detector/detector.py

@@ -0,0 +1,43 @@
+from importlib.resources import files
+
+from pandas import read_csv
+
+from pg_rad.utils.interpolators import (
+    get_field_efficiency, get_angular_efficiency
+)
+
+
+class Detector:
+    def __init__(
+        self,
+        name: str,
+        type: str,
+        is_isotropic: bool
+    ):
+        self.name = name
+        self.type = type
+        self.is_isotropic = is_isotropic
+
+    def get_efficiency(self, energy_keV, angle=None):
+        f_eff = get_field_efficiency(self.name, energy_keV)
+
+        if self.is_isotropic or angle is None:
+            return f_eff
+        else:
+            f_eff = get_field_efficiency(self.name, energy_keV)
+            a_eff = get_angular_efficiency(self.name, energy_keV, *angle)
+            return f_eff * a_eff
+
+
+def load_detector(name) -> Detector:
+    csv = files('pg_rad.data').joinpath('detectors.csv')
+    data = read_csv(csv)
+    dets = data['name'].values
+    if name in dets:
+        row = data[data['name'] == name].iloc[0]
+        return Detector(row['name'], row['type'], row['is_isotropic'])
+    else:
+        raise NotImplementedError(
+            f"Detector with name '{name}' not in detector library. Available:"
+            f"{', '.join(dets)}"
+        )

src/pg_rad/detector/detectors.py

@@ -1,38 +0,0 @@
-from abc import ABC
-
-
-class BaseDetector(ABC):
-    def __init__(
-        self,
-        name: str,
-        efficiency: float
-    ):
-        self.name = name
-        self.efficiency = efficiency
-
-    def get_efficiency(self):
-        pass
-
-
-class IsotropicDetector(BaseDetector):
-    def __init__(
-        self,
-        name: str,
-        efficiency: float,
-    ):
-        super().__init__(name, efficiency)
-
-    def get_efficiency(self, energy):
-        return self.efficiency
-
-
-class AngularDetector(BaseDetector):
-    def __init__(
-        self,
-        name: str,
-        efficiency: float
-    ):
-        super().__init__(name, efficiency)
-
-    def get_efficiency(self, angle, energy):
-        pass

src/pg_rad/exceptions/__init__.py

@@ -1,10 +0,0 @@
-# do not expose internal logger when running mkinit
-__ignore__ = ["logger"]
-
-from pg_rad.exceptions import exceptions
-
-from pg_rad.exceptions.exceptions import (ConvergenceError, DataLoadError,
-                                          InvalidCSVError, OutOfBoundsError,)
-
-__all__ = ['ConvergenceError', 'DataLoadError', 'InvalidCSVError',
-           'OutOfBoundsError', 'exceptions']

src/pg_rad/inputparser/parser.py

@@ -11,6 +11,7 @@ from pg_rad.exceptions.exceptions import (
     InvalidYAMLError
 )
 from pg_rad.configs import defaults
+from pg_rad.isotopes.isotope import get_isotope
 
 from .specs import (
     MetadataSpec,
@@ -276,7 +277,9 @@ class ConfigParser:
         specs: List[PointSourceSpec] = []
 
         for name, params in source_dict.items():
-            required = {"activity_MBq", "isotope", "position"}
+            required = {
+                "activity_MBq", "isotope", "position", "gamma_energy_keV"
+            }
             if not isinstance(params, dict):
                 raise InvalidConfigValueError(
                     f"sources.{name} is not defined correctly."
@@ -288,7 +291,10 @@ class ConfigParser:
                 raise MissingConfigKeyError(name, missing)
 
             activity = params.get("activity_MBq")
-            isotope = params.get("isotope")
+            isotope_name = params.get("isotope")
+            gamma_energy_keV = params.get("gamma_energy_keV")
+
+            isotope = get_isotope(isotope_name, gamma_energy_keV)
 
             if not isinstance(activity, int | float) or activity <= 0:
                 raise InvalidConfigValueError(
@@ -347,41 +353,11 @@ class ConfigParser:
         return specs
 
     def _parse_detector(self) -> DetectorSpec:
-        det_dict = self.config.get("detector")
+        det_name = self.config.get("detector")
-        required = {"name", "is_isotropic"}
+        if not det_name:
-        if not isinstance(det_dict, dict):
+            raise MissingConfigKeyError("detector")
-            raise InvalidConfigValueError(
-                "detector is not specified correctly. Must contain at least"
-                f"the subkeys {required}"
-                )
 
-        missing = required - det_dict.keys()
+        return DetectorSpec(name=det_name)
-        if missing:
-            raise MissingConfigKeyError("detector", missing)
-
-        name = det_dict.get("name")
-        is_isotropic = det_dict.get("is_isotropic")
-        eff = det_dict.get("efficiency")
-
-        default_detectors = defaults.DETECTOR_EFFICIENCIES
-
-        if name in default_detectors.keys() and not eff:
-            eff = default_detectors[name]
-        elif eff:
-            pass
-        else:
-            raise InvalidConfigValueError(
-                f"The detector {name} not found in library. Either "
-                f"specify {name}.efficiency or "
-                "choose a detector from the following list: "
-                f"{default_detectors.keys()}."
-            )
-
-        return DetectorSpec(
-            name=name,
-            efficiency=eff,
-            is_isotropic=is_isotropic
-        )
 
     def _warn_unknown_keys(self, section: str, provided: set, allowed: set):
         unknown = provided - allowed

src/pg_rad/inputparser/specs.py

@@ -67,8 +67,6 @@ class RelativePointSourceSpec(PointSourceSpec):
 @dataclass
 class DetectorSpec:
     name: str
-    efficiency: float
-    is_isotropic: bool
 
 
 @dataclass

src/pg_rad/isotopes/__init__.py

@@ -1,9 +0,0 @@
-# do not expose internal logger when running mkinit
-__ignore__ = ["logger"]
-
-from pg_rad.isotopes import isotope
-
-from pg_rad.isotopes.isotope import (CS137, Isotope, get_isotope,
-                                     preset_isotopes,)
-
-__all__ = ['CS137', 'Isotope', 'get_isotope', 'isotope', 'preset_isotopes']

src/pg_rad/isotopes/isotope.py

@@ -1,7 +1,10 @@
-from typing import Dict, Type
+from importlib.resources import files
 
+from pandas import read_csv
+
+from pg_rad.configs.filepaths import ISOTOPE_TABLE
 from pg_rad.exceptions.exceptions import InvalidIsotopeError
-from pg_rad.physics.attenuation import get_mass_attenuation_coeff
+from pg_rad.utils.interpolators import get_mass_attenuation_coeff
 
 
 class Isotope:
@@ -30,22 +33,35 @@ class Isotope:
         self.mu_mass_air = get_mass_attenuation_coeff(E / 1000)
 
 
-class CS137(Isotope):
+def get_isotope(isotope: str, energy_gamma_keV: float) -> Isotope:
-    def __init__(self):
+    """Lazy factory function to create isotope objects."""
-        super().__init__(
+    path = files('pg_rad.data').joinpath(ISOTOPE_TABLE)
-            name="Cs-137",
+    df = read_csv(path)
-            E=661.66,
+
-            b=0.851
+    isotope_df = df[df['isotope'] == isotope]
+
+    if isotope_df.empty:
+        raise InvalidIsotopeError(f"No data available for isotope {isotope}.")
+
+    tol = 0.01 * energy_gamma_keV
+    closest_energy = isotope_df[
+        (isotope_df['gamma_energy_keV'] >= energy_gamma_keV - tol) &
+        (isotope_df['gamma_energy_keV'] <= energy_gamma_keV + tol)
+    ]
+
+    if closest_energy.empty:
+        available_gammas = ', '.join(
+            str(x)+' keV' for x in isotope_df['gamma_energy_keV'].to_list()
+        )
+        raise InvalidIsotopeError(
+            f"No gamma of {energy_gamma_keV}±{tol} keV "
+            f"found for isotope {isotope}. "
+            f"Available gammas are: {available_gammas}"
         )
 
-
+    matched_row = closest_energy.iloc[0]
-preset_isotopes: Dict[str, Type[Isotope]] = {
+    return Isotope(
-    "Cs137": CS137
+        name=isotope,
-}
+        E=matched_row['gamma_energy_keV'],
-
+        b=matched_row['branching_ratio']
-
+    )
-def get_isotope(isotope_str: str) -> Isotope:
-    """Lazy factory function to create isotope objects."""
-    if isotope_str not in preset_isotopes:
-        raise InvalidIsotopeError(f"Unknown isotope: {isotope_str}")
-    return preset_isotopes[isotope_str]()

src/pg_rad/landscape/builder.py

@@ -15,7 +15,7 @@ from pg_rad.inputparser.specs import (
     RelativePointSourceSpec,
     DetectorSpec
 )
-
+from pg_rad.detector.detector import load_detector
 from pg_rad.path.path import Path, path_from_RT90
 
 from road_gen.generators.segmented_road_generator import SegmentedRoadGenerator
@@ -161,7 +161,7 @@ class LandscapeBuilder:
             ))
 
     def set_detector(self, spec: DetectorSpec) -> Self:
-        self._detector = spec
+        self._detector = load_detector(spec.name)
         return self
 
     def _fit_landscape_to_path(self) -> None:

src/pg_rad/landscape/landscape.py

@@ -1,7 +1,7 @@
 import logging
 
 from pg_rad.path.path import Path
-from pg_rad.detector.detectors import AngularDetector, IsotropicDetector
+from pg_rad.detector.detector import Detector
 from pg_rad.objects.sources import PointSource
 
 
@@ -18,7 +18,7 @@ class Landscape:
             path: Path,
             air_density: float,
             point_sources: list[PointSource],
-            detector: IsotropicDetector | AngularDetector,
+            detector: Detector,
             size: tuple[int, int, int]
             ):
         """Initialize a landscape.
@@ -28,6 +28,7 @@ class Landscape:
             path (Path): The path of the detector.
             air_density (float): Air density in kg/m^3.
             point_sources (list[PointSource]): List of point sources.
+            detector (Detector): The detector object.
             size (tuple[int, int, int]): Size of the world.
         """
 

src/pg_rad/logger/__init__.py

@@ -1,5 +0,0 @@
-from pg_rad.logger import logger
-
-from pg_rad.logger.logger import (setup_logger,)
-
-__all__ = ['logger', 'setup_logger']

src/pg_rad/main.py

@@ -4,7 +4,6 @@ import sys
 
 from pandas.errors import ParserError
 
-from pg_rad.detector.builder import DetectorBuilder
 from pg_rad.exceptions.exceptions import (
     MissingConfigKeyError,
     OutOfBoundsError,
@@ -56,10 +55,12 @@ def main():
         acquisition_time: 1
 
         path:
-            length: 1000
+            length:
+                - 500
+                - 500
             segments:
                 - straight
-                - turn_left
+                - turn_left: 45
             direction: negative
 
         sources:
@@ -67,20 +68,17 @@ def main():
                 activity_MBq: 100
                 position: [250, 100, 0]
                 isotope: Cs137
+                gamma_energy_keV: 661
 
-        detector:
+        detector: LU_NaI_3inch
-            name: dummy
-            is_isotropic: True
         """
 
         cp = ConfigParser(test_yaml).parse()
         landscape = LandscapeDirector.build_from_config(cp)
-        detector = DetectorBuilder(cp.detector).build()
         output = SimulationEngine(
             landscape=landscape,
             runtime_spec=cp.runtime,
             sim_spec=cp.options,
-            detector=detector
         ).simulate()
 
         plotter = ResultPlotter(landscape, output)
@@ -90,10 +88,8 @@ def main():
         try:
             cp = ConfigParser(args.config).parse()
             landscape = LandscapeDirector.build_from_config(cp)
-            detector = DetectorBuilder(cp.detector).build()
             output = SimulationEngine(
                 landscape=landscape,
-                detector=detector,
                 runtime_spec=cp.runtime,
                 sim_spec=cp.options
             ).simulate()

src/pg_rad/objects/__init__.py

@@ -1,11 +0,0 @@
-# do not expose internal logger when running mkinit
-__ignore__ = ["logger"]
-
-from pg_rad.objects import objects
-from pg_rad.objects import sources
-
-from pg_rad.objects.objects import (BaseObject,)
-from pg_rad.objects.sources import (PointSource,)
-
-__all__ = ['BaseObject', 'PointSource', 'objects',
-           'sources']

src/pg_rad/objects/sources.py

@@ -1,7 +1,7 @@
 import logging
 
 from .objects import BaseObject
-from pg_rad.isotopes.isotope import Isotope, get_isotope
+from pg_rad.isotopes.isotope import Isotope
 
 logger = logging.getLogger(__name__)
 
@@ -12,7 +12,7 @@ class PointSource(BaseObject):
     def __init__(
             self,
             activity_MBq: int,
-            isotope: str,
+            isotope: Isotope,
             position: tuple[float, float, float] = (0, 0, 0),
             name: str | None = None,
             color: str = 'red'
@@ -40,7 +40,7 @@ class PointSource(BaseObject):
         super().__init__(position, name, color)
 
         self.activity = activity_MBq
-        self.isotope: Isotope = get_isotope(isotope)
+        self.isotope = isotope
 
         logger.debug(f"Source created: {self.name}")
 

src/pg_rad/path/__init__.py

@@ -1,8 +0,0 @@
-# do not expose internal logger when running mkinit
-__ignore__ = ["logger"]
-
-from pg_rad.path import path
-
-from pg_rad.path.path import (Path, PathSegment, path_from_RT90,)
-
-__all__ = ['Path', 'PathSegment', 'path', 'path_from_RT90']

src/pg_rad/physics/__init__.py

@@ -1,12 +0,0 @@
-# do not expose internal logger when running mkinit
-__ignore__ = ["logger"]
-from pg_rad.physics import attenuation
-from pg_rad.physics import fluence
-
-from pg_rad.physics.attenuation import (get_mass_attenuation_coeff,)
-from pg_rad.physics.fluence import (calculate_fluence_along_path,
-                                    calculate_fluence_at, phi,)
-
-__all__ = ['attenuation', 'calculate_fluence_along_path',
-           'calculate_fluence_at', 'fluence', 'get_mass_attenuation_coeff',
-           'phi']

src/pg_rad/physics/attenuation.py

@@ -1,17 +0,0 @@
-from importlib.resources import files
-
-from pandas import read_csv
-from scipy.interpolate import interp1d
-
-from pg_rad.configs.filepaths import ATTENUATION_TABLE
-
-
-def get_mass_attenuation_coeff(
-        *args
-        ) -> float:
-    csv = files('pg_rad.data').joinpath(ATTENUATION_TABLE)
-    data = read_csv(csv)
-    x = data["energy_mev"].to_numpy()
-    y = data["mu"].to_numpy()
-    f = interp1d(x, y)
-    return f(*args)

src/pg_rad/physics/fluence.py

@@ -2,11 +2,9 @@ from typing import Tuple, TYPE_CHECKING
 
 import numpy as np
 
-from pg_rad.detector.detectors import IsotropicDetector, AngularDetector
-
-
 if TYPE_CHECKING:
     from pg_rad.landscape.landscape import Landscape
+    from pg_rad.detector.detector import Detector
 
 
 def phi(
@@ -33,16 +31,14 @@ def phi(
     """
 
     # Linear photon attenuation coefficient in m^-1.
-    mu_mass_air *= 0.1
+    mu_air = 0.1 * mu_mass_air * air_density
-    mu_air = mu_mass_air * air_density
 
     phi_r = (
         activity
         * eff
         * branching_ratio
         * np.exp(-mu_air * r)
-        / (4 * np.pi * r**2)
+    ) / (4 * np.pi * r**2)
-    )
 
     return phi_r
 
@@ -50,8 +46,7 @@ def phi(
 def calculate_count_rate_per_second(
     landscape: "Landscape",
     pos: np.ndarray,
-    detector: IsotropicDetector | AngularDetector,
+    detector: "Detector",
-    tangent_vectors: np.ndarray,
     scaling=1E6
 ):
     """Compute count rate in s^-1 m^-2 at a position in the landscape.
@@ -59,7 +54,7 @@ def calculate_count_rate_per_second(
     Args:
         landscape (Landscape): The landscape to compute.
         pos (np.ndarray): (N, 3) array of positions.
-        detector (IsotropicDetector | AngularDetector):
+        detector (Detector):
             Detector object, needed to compute correct efficiency.
 
     Returns:
@@ -69,22 +64,29 @@ def calculate_count_rate_per_second(
     total_phi = np.zeros(pos.shape[0])
 
     for source in landscape.point_sources:
+        # See Bukartas (2021) page 25 for incidence angle math
         source_to_detector = pos - np.array(source.pos)
         r = np.linalg.norm(source_to_detector, axis=1)
         r = np.maximum(r, 1E-3)  # enforce minimum distance of 1cm
 
-        if isinstance(detector, AngularDetector):
+        if not detector.is_isotropic:
-            cos_theta = (
+            v = np.zeros_like(pos)
-                np.sum(tangent_vectors * source_to_detector, axis=1) / (
+            v[1:] = pos[1:] - pos[:-1]
-                    np.linalg.norm(source_to_detector, axis=1) *
+            v[0] = v[1]  # handle first point
-                    np.linalg.norm(tangent_vectors, axis=1)
+            vx, vy = v[:, 0], v[:, 1]
-                )
+
-            )
+            r_vec = pos - np.array(source.pos)
-            cos_theta = np.clip(cos_theta, -1, 1)
+            rx, ry = r_vec[:, 0], r_vec[:, 1]
-            theta = np.arccos(cos_theta)
+
-            eff = detector.get_efficiency(theta, energy=source.isotope.E)
+            theta = np.arctan2(vy, vx) - np.arctan2(ry, rx)
+
+            # normalise to [-pi, pi] and convert to degrees
+            theta = (theta + np.pi) % (2 * np.pi) - np.pi
+            theta_deg = np.degrees(theta)
+
+            eff = detector.get_efficiency(source.isotope.E, theta_deg)
         else:
-            eff = detector.get_efficiency(energy=source.isotope.E)
+            eff = detector.get_efficiency(source.isotope.E)
 
         phi_source = phi(
             r=r,
@@ -102,15 +104,15 @@ def calculate_count_rate_per_second(
 
 def calculate_counts_along_path(
     landscape: "Landscape",
-    detector: "IsotropicDetector | AngularDetector",
+    detector: "Detector",
-    acquisition_time: int,
+    velocity: float,
     points_per_segment: int = 10,
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Compute the counts recorded in each acquisition period in the landscape.
 
     Args:
         landscape (Landscape): _description_
-        detector (IsotropicDetector | AngularDetector): _description_
+        detector (Detector): _description_
         points_per_segment (int, optional): _description_. Defaults to 100.
 
     Returns:
@@ -123,7 +125,6 @@ def calculate_counts_along_path(
     num_segments = len(path.segments)
 
     segment_lengths = np.array([seg.length for seg in path.segments])
-    ds = segment_lengths[0]
     original_distances = np.zeros(num_points)
     original_distances[1:] = np.cumsum(segment_lengths)
 
@@ -140,26 +141,17 @@ def calculate_counts_along_path(
     if path.opposite_direction:
         full_positions = np.flip(full_positions, axis=0)
 
-    # to compute the angle between sources and the direction of travel, we
+    # [counts/s]
-    # compute tangent vectors along the path.
+    cps = calculate_count_rate_per_second(
-    dx_ds = np.gradient(xnew, s)
+        landscape, full_positions, detector
-    dy_ds = np.gradient(ynew, s)
-    tangent_vectors = np.c_[dx_ds, dy_ds, np.zeros_like(dx_ds)]
-    tangent_vectors /= np.linalg.norm(tangent_vectors, axis=1, keepdims=True)
-
-    count_rate = calculate_count_rate_per_second(
-        landscape, full_positions, detector, tangent_vectors
     )
 
-    count_rate *= (acquisition_time / points_per_segment)
+    # reshape so each segment is on a row
+    cps_per_seg = cps.reshape(num_segments, points_per_segment)
 
-    count_rate_segs = count_rate.reshape(num_segments, points_per_segment)
+    du = s[1] - s[0]
-    integrated = np.trapezoid(
+    integrated_counts = np.trapezoid(cps_per_seg, dx=du, axis=1) / velocity
-        count_rate_segs,
+    int_counts_result = np.zeros(num_points)
-        dx=ds/points_per_segment,
+    int_counts_result[1:] = integrated_counts
-        axis=1
-    )
 
-    result = np.zeros(num_points)
+    return original_distances, s, cps, int_counts_result
-    result[1:] = integrated
-    return original_distances, result

src/pg_rad/plotting/__init__.py

@@ -1,7 +0,0 @@
-# do not expose internal logger when running mkinit
-__ignore__ = ["logger"]
-from pg_rad.plotting import landscape_plotter
-
-from pg_rad.plotting.landscape_plotter import (LandscapeSlicePlotter,)
-
-__all__ = ['LandscapeSlicePlotter', 'landscape_plotter']

src/pg_rad/plotting/result_plotter.py

@@ -1,3 +1,7 @@
+from importlib.resources import files
+
+import numpy as np
+import pandas as pd
 from matplotlib import pyplot as plt
 from matplotlib.gridspec import GridSpec
 
@@ -13,45 +17,79 @@ class ResultPlotter:
         self.source_res = output.sources
 
     def plot(self, landscape_z: float = 0):
-        fig = plt.figure(figsize=(12, 10), constrained_layout=True)
+        self._plot_main(landscape_z)
-        fig.suptitle(self.landscape.name)
+        self._plot_detector()
-        gs = GridSpec(
+        self._plot_metadata()
-            4,
-            2,
-            width_ratios=[0.5, 0.5],
-            height_ratios=[0.7, 0.1, 0.1, 0.1],
-            hspace=0.2)
-
-        ax1 = fig.add_subplot(gs[0, 0])
-        self._draw_count_rate(ax1)
-
-        ax2 = fig.add_subplot(gs[0, 1])
-        self._plot_landscape(ax2, landscape_z)
-
-        ax3 = fig.add_subplot(gs[1, :])
-        self._draw_table(ax3)
-
-        ax4 = fig.add_subplot(gs[2, :])
-        self._draw_detector_table(ax4)
-
-        ax5 = fig.add_subplot(gs[3, :])
-        self._draw_source_table(ax5)
 
         plt.show()
 
+    def _plot_main(self, landscape_z):
+        fig = plt.figure(figsize=(12, 8))
+        fig.suptitle(self.landscape.name)
+        fig.canvas.manager.set_window_title("Main results")
+
+        gs = GridSpec(
+            2, 2, figure=fig,
+            height_ratios=[1, 2], width_ratios=[1, 1],
+            hspace=0.3, wspace=0.3)
+
+        ax_cps = fig.add_subplot(gs[0, 0])
+        self._draw_cps(ax_cps)
+
+        ax_counts = fig.add_subplot(gs[0, 1])
+        self._draw_count_rate(ax_counts)
+
+        ax_landscape = fig.add_subplot(gs[1, :])
+        self._plot_landscape(ax_landscape, landscape_z)
+
+    def _plot_detector(self):
+        det = self.landscape.detector
+
+        fig = plt.figure(figsize=(10, 4))
+        fig.canvas.manager.set_window_title("Detector")
+
+        gs = GridSpec(1, 2, figure=fig, width_ratios=[0.5, 0.5])
+
+        ax_table = fig.add_subplot(gs[0, 0])
+        self._draw_detector_table(ax_table)
+
+        if not det.is_isotropic:
+            ax_polar = fig.add_subplot(gs[0, 1], projection='polar')
+
+            energies = [
+                source.primary_gamma for source in self.source_res
+            ]
+
+            self._draw_angular_efficiency_polar(ax_polar, det, energies[0])
+
+    def _plot_metadata(self):
+        fig, axs = plt.subplots(2, 1, figsize=(10, 6))
+        fig.canvas.manager.set_window_title("Simulation Metadata")
+
+        self._draw_table(axs[0])
+        self._draw_source_table(axs[1])
+
     def _plot_landscape(self, ax, z):
         lp = LandscapeSlicePlotter()
         ax = lp.plot(landscape=self.landscape, z=z, ax=ax, show=False)
         return ax
 
-    def _draw_count_rate(self, ax):
+    def _draw_cps(self, ax):
-        x = self.count_rate_res.arc_length
+        x = self.count_rate_res.sub_points
-        y = self.count_rate_res.count_rate
+        y = self.count_rate_res.cps
-        ax.plot(x, y, label='Count rate', color='r')
+        ax.plot(x, y, color='b')
-        ax.set_title('Count rate')
+        ax.set_title('Counts per second (CPS)')
         ax.set_xlabel('Arc length s [m]')
-        ax.set_ylabel('Counts')
+        ax.set_ylabel('CPS [s$^{-1}$]')
-        ax.legend()
+
+    def _draw_count_rate(self, ax):
+        x = self.count_rate_res.acquisition_points
+        y = self.count_rate_res.integrated_counts
+        ax.plot(x, y, color='r', linestyle='--', alpha=0.2)
+        ax.scatter(x, y, color='r', marker='x')
+        ax.set_title('Integrated counts')
+        ax.set_xlabel('Arc length s [m]')
+        ax.set_ylabel('N')
 
     def _draw_table(self, ax):
         ax.set_axis_off()
@@ -60,6 +98,7 @@ class ResultPlotter:
         data = [
             ["Air density (kg/m^3)", round(self.landscape.air_density, 3)],
             ["Total path length (m)", round(self.landscape.path.length, 3)],
+            ["Readout points", len(self.count_rate_res.integrated_counts)],
         ]
 
         ax.table(
@@ -72,13 +111,28 @@ class ResultPlotter:
 
     def _draw_detector_table(self, ax):
         det = self.landscape.detector
-        print(det)
+
+        if det.is_isotropic:
+            det_type = "Isotropic"
+        else:
+            det_type = "Angular"
+
+        source_energies = [
+            source.primary_gamma for source in self.source_res
+        ]
+
+        # list field efficiencies for each primary gamma in the landscape
+        effs = {e: det.get_efficiency(e) for e in source_energies}
+        formatted_effs = ", ".join(
+            f"{value:.3f} @ {key:.1f} keV"
+            for key, value in effs.items()
+        )
         ax.set_axis_off()
         ax.set_title('Detector')
         cols = ('Parameter', 'Value')
         data = [
-            ["Name", det.name],
+            ["Detector", f"{det.name} ({det_type})"],
-            ["Type", det.efficiency],
+            ["Field efficiency", formatted_effs],
         ]
 
         ax.table(
@@ -105,7 +159,7 @@ class ResultPlotter:
         data = [
             [
                 s.name,
-                s.isotope,
+                s.isotope+f" ({s.primary_gamma} keV)",
                 s.activity,
                 "("+", ".join(f"{val:.2f}" for val in s.position)+")",
                 round(s.dist_from_path, 2)
@@ -119,3 +173,34 @@ class ResultPlotter:
         )
 
         return ax
+
+    def _draw_angular_efficiency_polar(self, ax, detector, energy_keV):
+        # find the energies available for this detector.
+        csv = files('pg_rad.data.angular_efficiencies').joinpath(
+            detector.name + '.csv'
+        )
+        data = pd.read_csv(csv)
+        energy_cols = [col for col in data.columns if col != "angle"]
+        energies = np.array([float(col) for col in energy_cols])
+
+        # take energy column that is within 1% tolerance of energy_keV
+        rel_diff = np.abs(energies - energy_keV) / energies
+        match_idx = np.where(rel_diff <= 0.01)[0]
+        best_idx = match_idx[np.argmin(rel_diff[match_idx])]
+        col = energy_cols[best_idx]
+
+        theta_deg = data["angle"].to_numpy()
+        eff = data[col].to_numpy()
+
+        idx = np.argsort(theta_deg)
+        theta_deg = theta_deg[idx]
+        eff = eff[idx]
+
+        theta_deg = np.append(theta_deg, theta_deg[0])
+        eff = np.append(eff, eff[0])
+
+        theta_rad = np.radians(theta_deg)
+
+        print(theta_rad)
+        ax.plot(theta_rad, eff)
+        ax.set_title(f"Rel. angular efficiency @ {energy_keV:.1f} keV")

src/pg_rad/simulator/engine.py

@@ -6,7 +6,7 @@ from pg_rad.simulator.outputs import (
     SimulationOutput,
     SourceOutput
 )
-from pg_rad.detector.detectors import IsotropicDetector, AngularDetector
+
 from pg_rad.physics.fluence import calculate_counts_along_path
 from pg_rad.utils.projection import minimal_distance_to_path
 from pg_rad.inputparser.specs import RuntimeSpec, SimulationOptionsSpec
@@ -17,13 +17,12 @@ class SimulationEngine:
     def __init__(
         self,
         landscape: Landscape,
-        detector: IsotropicDetector | AngularDetector,
         runtime_spec: RuntimeSpec,
         sim_spec: SimulationOptionsSpec,
     ):
 
         self.landscape = landscape
-        self.detector = detector
+        self.detector = self.landscape.detector
         self.runtime_spec = runtime_spec
         self.sim_spec = sim_spec
 
@@ -40,12 +39,13 @@ class SimulationEngine:
         )
 
     def _calculate_count_rate_along_path(self) -> CountRateOutput:
-        arc_length, phi = calculate_counts_along_path(
+        acq_points, sub_points, cps, int_counts = calculate_counts_along_path(
             self.landscape,
             self.detector,
-            acquisition_time=self.runtime_spec.acquisition_time
+            velocity=self.runtime_spec.speed
         )
-        return CountRateOutput(arc_length, phi)
+
+        return CountRateOutput(acq_points, sub_points, cps, int_counts)
 
     def _calculate_point_source_distance_to_path(self) -> List[SourceOutput]:
 
@@ -62,6 +62,7 @@ class SimulationEngine:
                 SourceOutput(
                     s.name,
                     s.isotope.name,
+                    s.isotope.E,
                     s.activity,
                     s.pos,
                     dist_to_path)

src/pg_rad/simulator/outputs.py

@@ -5,14 +5,17 @@ from dataclasses import dataclass
 
 @dataclass
 class CountRateOutput:
-    arc_length: List[float]
+    acquisition_points: List[float]
-    count_rate: List[float]
+    sub_points: List[float]
+    cps: List[float]
+    integrated_counts: List[float]
 
 
 @dataclass
 class SourceOutput:
     name: str
     isotope: str
+    primary_gamma: float
     activity: float
     position: Tuple[float, float, float]
     dist_from_path: float

src/pg_rad/utils/interpolators.py

@@ -0,0 +1,56 @@
+from importlib.resources import files
+
+import numpy as np
+from pandas import read_csv
+from scipy.interpolate import interp1d, CubicSpline
+
+from pg_rad.configs.filepaths import ATTENUATION_TABLE
+
+
+def get_mass_attenuation_coeff(*args) -> float:
+    csv = files('pg_rad.data').joinpath(ATTENUATION_TABLE)
+    data = read_csv(csv)
+    x = data["energy_mev"].to_numpy()
+    y = data["mu"].to_numpy()
+    f = interp1d(x, y)
+    return f(*args)
+
+
+def get_field_efficiency(name: str, energy_keV: float) -> float:
+    csv = files('pg_rad.data.field_efficiencies').joinpath(name+'.csv')
+    data = read_csv(csv)
+    data = data.groupby("energy_keV", as_index=False).mean()
+    x = data["energy_keV"].to_numpy()
+    y = data["field_efficiency_m2"].to_numpy()
+    f = CubicSpline(x, y)
+    return f(energy_keV)
+
+
+def get_angular_efficiency(name: str, energy_keV: float, *angle: float):
+    csv = files('pg_rad.data.angular_efficiencies').joinpath(name+'.csv')
+    data = read_csv(csv)
+
+    # check all energies at which angular eff. is available for this detector.
+    # this is done within 1% tolerance
+    energy_cols = [col for col in data.columns if col != "angle"]
+    energies = np.array([float(col) for col in energy_cols])
+    rel_diff = np.abs(energies - energy_keV) / energies
+    match_idx = np.where(rel_diff <= 0.01)[0]
+
+    if len(match_idx) == 0:
+        raise NotImplementedError(
+            f"No angular efficiency defined for {energy_keV} keV "
+            f"in detector '{name}'. Available: {energies}"
+        )
+
+    best_idx = match_idx[np.argmin(rel_diff[match_idx])]
+    selected_energy_col = energy_cols[best_idx]
+
+    x = data["angle"].to_numpy()
+    y = data[selected_energy_col].to_numpy()
+    idx = np.argsort(x)
+    x = x[idx]
+    y = y[idx]
+    f = interp1d(x, y)
+
+    return f(angle)

tests/test_attenuation_functions.py

@@ -1,6 +1,6 @@
 import pytest
 
-from pg_rad.physics import get_mass_attenuation_coeff
+from pg_rad.utils.interpolators import get_mass_attenuation_coeff
 
 
 @pytest.mark.parametrize("energy,mu", [

tests/test_fluence_rate.py

@@ -3,20 +3,20 @@ import pytest
 
 from pg_rad.inputparser.parser import ConfigParser
 from pg_rad.landscape.director import LandscapeDirector
-from pg_rad.physics import calculate_fluence_at
+from pg_rad.physics.fluence import phi
 
 
 @pytest.fixture
 def isotropic_detector():
-    from pg_rad.detector.detectors import IsotropicDetector
+    from pg_rad.detector.detector import load_detector
-    return IsotropicDetector(name="test_detector", eff=1.0)
+    return load_detector('dummy')
 
 
 @pytest.fixture
 def phi_ref(test_landscape, isotropic_detector):
     source = test_landscape.point_sources[0]
 
-    r = np.linalg.norm(np.array([10, 10, 0]) - np.array(source.pos))
+    r = np.linalg.norm(np.array([0, 0, 0]) - np.array(source.pos))
 
     A = source.activity * 1E6
     b = source.isotope.b
@@ -43,12 +43,11 @@ def test_landscape():
     sources:
         test_source:
             activity_MBq: 100
-            position: [0, 0, 0]
+            position: [0, 100, 0]
-            isotope: CS137
+            isotope: Cs137
+            gamma_energy_keV: 661
 
-    detector:
+    detector: dummy
-        name: dummy
-        is_isotropic: True
     """
 
     cp = ConfigParser(test_yaml).parse()
@@ -57,10 +56,15 @@ def test_landscape():
 
 
 def test_single_source_fluence(phi_ref, test_landscape, isotropic_detector):
-    phi = calculate_fluence_at(
+    s = test_landscape.point_sources[0]
-        test_landscape,
+    r = np.linalg.norm(np.array([0, 0, 0]) - np.array(s.pos))
-        np.array([10, 10, 0]),
+    phi_calc = phi(
-        isotropic_detector,
+        r,
-        tangent_vectors=None,
+        s.activity*1E6,
+        s.isotope.b,
+        s.isotope.mu_mass_air,
+        test_landscape.air_density,
+        isotropic_detector.get_efficiency(s.isotope.E)
         )
-    assert pytest.approx(phi[0], rel=1E-3) == phi_ref
+
+    assert pytest.approx(phi_calc, rel=1E-6) == phi_ref

tests/test_sources.py

@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from pg_rad.objects import PointSource
+from pg_rad.objects.sources import PointSource
 
 
 @pytest.fixture