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

Feature detectors roadgen

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/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/inputparser/parser.py

@@ -353,41 +353,11 @@ class ConfigParser:
         return specs
 
     def _parse_detector(self) -> DetectorSpec:
-        det_dict = self.config.get("detector")
-        required = {"name", "is_isotropic"}
-        if not isinstance(det_dict, dict):
-            raise InvalidConfigValueError(
-                "detector is not specified correctly. Must contain at least"
-                f"the subkeys {required}"
-                )
+        det_name = self.config.get("detector")
+        if not det_name:
+            raise MissingConfigKeyError("detector")
 
-        missing = required - det_dict.keys()
-        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
-        )
+        return DetectorSpec(name=det_name)
 
     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/isotope.py

@@ -4,7 +4,7 @@ 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:

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/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,7 +55,9 @@ def main():
         acquisition_time: 1
 
         path:
-            length: 1000
+            length:
+                - 500
+                - 500
             segments:
                 - straight
                 - turn_left: 45
@@ -69,19 +70,15 @@ def main():
                 isotope: Cs137
                 gamma_energy_keV: 661
 
-        detector:
-            name: dummy
-            is_isotropic: True
+        detector: LU_NaI_3inch
         """
 
         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)
@@ -91,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/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 = mu_mass_air * air_density
+    mu_air = 0.1 * 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,
-    tangent_vectors: np.ndarray,
+    detector: "Detector",
     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):
-            cos_theta = (
-                np.sum(tangent_vectors * source_to_detector, axis=1) / (
-                    np.linalg.norm(source_to_detector, axis=1) *
-                    np.linalg.norm(tangent_vectors, axis=1)
-                )
-            )
-            cos_theta = np.clip(cos_theta, -1, 1)
-            theta = np.arccos(cos_theta)
-            eff = detector.get_efficiency(theta, energy=source.isotope.E)
+        if not detector.is_isotropic:
+            v = np.zeros_like(pos)
+            v[1:] = pos[1:] - pos[:-1]
+            v[0] = v[1]  # handle first point
+            vx, vy = v[:, 0], v[:, 1]
+
+            r_vec = pos - np.array(source.pos)
+            rx, ry = r_vec[:, 0], r_vec[:, 1]
+
+            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",
-    acquisition_time: int,
+    detector: "Detector",
+    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
-    # compute tangent vectors along the path.
-    dx_ds = np.gradient(xnew, s)
-    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
+    # [counts/s]
+    cps = calculate_count_rate_per_second(
+        landscape, full_positions, detector
     )
 
-    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)
-    integrated = np.trapezoid(
-        count_rate_segs,
-        dx=ds/points_per_segment,
-        axis=1
-    )
+    du = s[1] - s[0]
+    integrated_counts = np.trapezoid(cps_per_seg, dx=du, axis=1) / velocity
+    int_counts_result = np.zeros(num_points)
+    int_counts_result[1:] = integrated_counts
 
-    result = np.zeros(num_points)
-    result[1:] = integrated
-    return original_distances, result
+    return original_distances, s, cps, int_counts_result

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)
-        fig.suptitle(self.landscape.name)
-        gs = GridSpec(
-            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)
+        self._plot_main(landscape_z)
+        self._plot_detector()
+        self._plot_metadata()
 
         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):
-        x = self.count_rate_res.arc_length
-        y = self.count_rate_res.count_rate
-        ax.plot(x, y, label='Count rate', color='r')
-        ax.set_title('Count rate')
+    def _draw_cps(self, ax):
+        x = self.count_rate_res.sub_points
+        y = self.count_rate_res.cps
+        ax.plot(x, y, color='b')
+        ax.set_title('Counts per second (CPS)')
         ax.set_xlabel('Arc length s [m]')
-        ax.set_ylabel('Counts')
-        ax.legend()
+        ax.set_ylabel('CPS [s$^{-1}$]')
+
+    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],
-            ["Type", det.efficiency],
+            ["Detector", f"{det.name} ({det_type})"],
+            ["Field efficiency", formatted_effs],
         ]
 
         ax.table(
@@ -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]:
 

src/pg_rad/simulator/outputs.py

@@ -5,8 +5,10 @@ from dataclasses import dataclass
 
 @dataclass
 class CountRateOutput:
-    arc_length: List[float]
-    count_rate: List[float]
+    acquisition_points: List[float]
+    sub_points: List[float]
+    cps: List[float]
+    integrated_counts: List[float]
 
 
 @dataclass

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
-    return IsotropicDetector(name="test_detector", eff=1.0)
+    from pg_rad.detector.detector import load_detector
+    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]
-            isotope: CS137
+            position: [0, 100, 0]
+            isotope: Cs137
+            gamma_energy_keV: 661
 
-    detector:
-        name: dummy
-        is_isotropic: True
+    detector: dummy
     """
 
     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(
-        test_landscape,
-        np.array([10, 10, 0]),
-        isotropic_detector,
-        tangent_vectors=None,
+    s = test_landscape.point_sources[0]
+    r = np.linalg.norm(np.array([0, 0, 0]) - np.array(s.pos))
+    phi_calc = phi(
+        r,
+        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