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/exceptions/exceptions.py

@@ -10,6 +10,10 @@ class InvalidCSVError(DataLoadError):
     """Raised when a file is not a valid CSV."""
 
 
+class InvalidYAMLError(DataLoadError):
+    """Raised when a file is not a valid YAML."""
+
+
 class OutOfBoundsError(Exception):
     """Raised when an object is attempted to be placed out of bounds."""
 
@@ -18,7 +22,11 @@ class MissingConfigKeyError(KeyError):
     """Raised when a (nested) config key is missing in the config."""
     def __init__(self, key, subkey=None):
         if subkey:
-            self.message = f"Missing key in {key}: {', '.join(list(subkey))}"
+            if isinstance(subkey, str):
+                pass
+            elif isinstance(subkey, set):
+                subkey = ', '.join(list(subkey))
+            self.message = f"Missing key in {key}: {subkey}"
         else:
             self.message = f"Missing key: {key}"
 

src/pg_rad/inputparser/parser.py

@@ -7,9 +7,11 @@ import yaml
 from pg_rad.exceptions.exceptions import (
     MissingConfigKeyError,
     DimensionError,
-    InvalidConfigValueError
+    InvalidConfigValueError,
+    InvalidYAMLError
 )
 from pg_rad.configs import defaults
+from pg_rad.isotopes.isotope import get_isotope
 
 from .specs import (
     MetadataSpec,
@@ -45,7 +47,7 @@ class ConfigParser:
 
     def parse(self) -> SimulationSpec:
         self._warn_unknown_keys(
-            section="global",
+            section="root",
             provided=set(self.config.keys()),
             allowed=self._ALLOWED_ROOT_KEYS,
         )
@@ -74,7 +76,7 @@ class ConfigParser:
             data = yaml.safe_load(config_source)
 
         if not isinstance(data, dict):
-            raise ValueError(
+            raise InvalidYAMLError(
                 "Provided path or string is not a valid YAML representation."
             )
 
@@ -86,13 +88,13 @@ class ConfigParser:
                 name=self.config["name"]
             )
         except KeyError as e:
-            raise MissingConfigKeyError("global", {"name"}) from e
+            raise MissingConfigKeyError("root", {"name"}) from e
 
     def _parse_runtime(self) -> RuntimeSpec:
         required = {"speed", "acquisition_time"}
         missing = required - self.config.keys()
         if missing:
-            raise MissingConfigKeyError("global", missing)
+            raise MissingConfigKeyError("root", missing)
 
         return RuntimeSpec(
             speed=float(self.config["speed"]),
@@ -116,7 +118,7 @@ class ConfigParser:
         seed = options.get("seed")
 
         if not isinstance(air_density, float) or air_density <= 0:
-            raise ValueError(
+            raise InvalidConfigValueError(
                 "options.air_density_kg_per_m3 must be a positive float "
                 "in kg/m^3."
             )
@@ -124,7 +126,9 @@ class ConfigParser:
             seed is not None or
             (isinstance(seed, int) and seed <= 0)
         ):
-            raise ValueError("Seed must be a positive integer value.")
+            raise InvalidConfigValueError(
+                "Seed must be a positive integer value."
+            )
 
         return SimulationOptionsSpec(
             air_density=air_density,
@@ -133,14 +137,26 @@ class ConfigParser:
 
     def _parse_path(self) -> PathSpec:
         allowed_csv = {"file", "east_col_name", "north_col_name", "z"}
-        allowed_proc = {"segments", "length", "z", "alpha"}
+        allowed_proc = {"segments", "length", "z", "alpha", "direction"}
 
         path = self.config.get("path")
+
+        direction = path.get("direction", 'positive')
+
+        if direction == 'positive':
+            opposite_direction = False
+        elif direction == 'negative':
+            opposite_direction = True
+        else:
+            raise InvalidConfigValueError(
+                "Direction must be positive or negative."
+            )
+
         if path is None:
-            raise MissingConfigKeyError("global", {"path"})
+            raise MissingConfigKeyError("path")
 
         if not isinstance(path, dict):
-            raise ValueError("Path must be a dictionary.")
+            raise InvalidConfigValueError("Path must be a dictionary.")
 
         if "file" in path:
             self._warn_unknown_keys(
@@ -154,26 +170,31 @@ class ConfigParser:
                 east_col_name=path["east_col_name"],
                 north_col_name=path["north_col_name"],
                 z=path.get("z", 0),
+                opposite_direction=opposite_direction
             )
 
-        if "segments" in path:
+        elif "segments" in path:
             self._warn_unknown_keys(
                 section="path (procedural)",
                 provided=set(path.keys()),
                 allowed=allowed_proc,
             )
-            return self._parse_procedural_path(path)
+            return self._parse_procedural_path(path, opposite_direction)
 
-        raise ValueError("Invalid path configuration.")
+        else:
+            raise InvalidConfigValueError("Invalid path configuration.")
 
     def _parse_procedural_path(
         self,
-        path: Dict[str, Any]
+        path: Dict[str, Any],
+        opposite_direction: bool
     ) -> ProceduralPathSpec:
         raw_segments = path.get("segments")
 
         if not isinstance(raw_segments, List):
-            raise ValueError("path.segments must be a list of segments.")
+            raise InvalidConfigValueError(
+                "path.segments must be a list of segments."
+            )
 
         raw_length = path.get("length")
 
@@ -191,7 +212,8 @@ class ConfigParser:
             angles=angles,
             lengths=lengths,
             z=path.get("z", defaults.DEFAULT_PATH_HEIGHT),
-            alpha=path.get("alpha", defaults.DEFAULT_ALPHA)
+            alpha=path.get("alpha", defaults.DEFAULT_ALPHA),
+            opposite_direction=opposite_direction
         )
 
     def _process_segment_angles(
@@ -209,13 +231,17 @@ class ConfigParser:
 
             elif isinstance(segment, dict):
                 if len(segment) != 1:
-                    raise ValueError("Invalid segment definition.")
+                    raise InvalidConfigValueError(
+                        "Invalid segment definition."
+                    )
                 seg_type, angle = list(segment.items())[0]
                 segments.append(seg_type)
                 angles.append(angle)
 
             else:
-                raise ValueError("Invalid segment entry format.")
+                raise InvalidConfigValueError(
+                    "Invalid segment entry format."
+                )
 
         return segments, angles
 
@@ -239,18 +265,36 @@ class ConfigParser:
         return segment_type in {"turn_left", "turn_right"}
 
     def _parse_point_sources(self) -> List[PointSourceSpec]:
-        source_dict = self.config.get("sources", {})
+        source_dict = self.config.get("sources")
+        if source_dict is None:
+            raise MissingConfigKeyError("sources")
+
+        if not isinstance(source_dict, dict):
+            raise InvalidConfigValueError(
+                "sources must have subkeys representing point source names."
+            )
+
         specs: List[PointSourceSpec] = []
 
         for name, params in source_dict.items():
+            required = {
+                "activity_MBq", "isotope", "position", "gamma_energy_keV"
+            }
+            if not isinstance(params, dict):
+                raise InvalidConfigValueError(
+                    f"sources.{name} is not defined correctly."
+                    f" Must have subkeys {required}"
+                )
 
-            required = {"activity_MBq", "isotope", "position"}
             missing = required - params.keys()
             if missing:
                 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(
@@ -309,36 +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:
+            raise MissingConfigKeyError("detector")
 
-        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 ValueError(
-                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

@@ -22,7 +22,8 @@ class SimulationOptionsSpec:
 
 @dataclass
 class PathSpec(ABC):
-    pass
+    z: int | float
+    opposite_direction: bool
 
 
 @dataclass
@@ -30,7 +31,6 @@ class ProceduralPathSpec(PathSpec):
     segments: list[str]
     angles: list[float]
     lengths: list[int | None]
-    z: int | float
     alpha: float
 
 
@@ -39,7 +39,6 @@ class CSVPathSpec(PathSpec):
     file: str
     east_col_name: str
     north_col_name: str
-    z: int | float
 
 
 @dataclass
@@ -68,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

@@ -12,9 +12,10 @@ from pg_rad.inputparser.specs import (
     SimulationSpec,
     CSVPathSpec,
     AbsolutePointSourceSpec,
-    RelativePointSourceSpec
+    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
@@ -30,6 +31,7 @@ class LandscapeBuilder:
         self._point_sources = []
         self._size = None
         self._air_density = 1.243
+        self._detector = None
 
         logger.debug(f"LandscapeBuilder initialized: {self.name}")
 
@@ -58,10 +60,12 @@ class LandscapeBuilder:
         self,
         sim_spec: SimulationSpec
     ):
-        segments = sim_spec.path.segments
+        path = sim_spec.path
-        lengths = sim_spec.path.lengths
+        segments = path.segments
-        angles = sim_spec.path.angles
+        lengths = path.lengths
-        alpha = sim_spec.path.alpha
+        angles = path.angles
+        alpha = path.alpha
+        z = path.z
 
         sg = SegmentedRoadGenerator(
             ds=sim_spec.runtime.speed * sim_spec.runtime.acquisition_time,
@@ -76,7 +80,11 @@ class LandscapeBuilder:
             alpha=alpha
         )
 
-        self._path = Path(list(zip(x, y)))
+        self._path = Path(
+            list(zip(x, y)),
+            z=z,
+            opposite_direction=path.opposite_direction
+        )
         self._fit_landscape_to_path()
         return self
 
@@ -89,7 +97,8 @@ class LandscapeBuilder:
             df=df,
             east_col=spec.east_col_name,
             north_col=spec.north_col_name,
-            z=spec.z
+            z=spec.z,
+            opposite_direction=spec.opposite_direction
         )
 
         self._fit_landscape_to_path()
@@ -151,6 +160,10 @@ class LandscapeBuilder:
                 name=s.name
             ))
 
+    def set_detector(self, spec: DetectorSpec) -> Self:
+        self._detector = load_detector(spec.name)
+        return self
+
     def _fit_landscape_to_path(self) -> None:
         """The size of the landscape will be updated if
         1) _size is not set, or
@@ -221,6 +234,7 @@ class LandscapeBuilder:
             name=self.name,
             path=self._path,
             point_sources=self._point_sources,
+            detector=self._detector,
             size=self._size,
             air_density=self._air_density
         )

src/pg_rad/landscape/director.py

@@ -45,5 +45,6 @@ class LandscapeDirector:
                 sim_spec=config,
             )
         lb.set_point_sources(*config.point_sources)
+        lb.set_detector(config.detector)
         landscape = lb.build()
         return landscape

src/pg_rad/landscape/landscape.py

@@ -1,6 +1,7 @@
 import logging
 
 from pg_rad.path.path import Path
+from pg_rad.detector.detector import Detector
 from pg_rad.objects.sources import PointSource
 
 
@@ -17,6 +18,7 @@ class Landscape:
             path: Path,
             air_density: float,
             point_sources: list[PointSource],
+            detector: Detector,
             size: tuple[int, int, int]
             ):
         """Initialize a landscape.
@@ -26,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.
         """
 
@@ -34,5 +37,6 @@ class Landscape:
         self.point_sources = point_sources
         self.size = size
         self.air_density = air_density
+        self.detector = detector
 
         logger.debug(f"Landscape created: {self.name}")

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

@@ -3,15 +3,14 @@ import logging
 import sys
 
 from pandas.errors import ParserError
-from yaml import YAMLError
 
-from pg_rad.detector.builder import DetectorBuilder
 from pg_rad.exceptions.exceptions import (
     MissingConfigKeyError,
     OutOfBoundsError,
     DimensionError,
     InvalidConfigValueError,
-    InvalidIsotopeError
+    InvalidIsotopeError,
+    InvalidYAMLError
 )
 from pg_rad.logger.logger import setup_logger
 from pg_rad.inputparser.parser import ConfigParser
@@ -56,29 +55,30 @@ def main():
         acquisition_time: 1
 
         path:
-            length: 1000
+            length:
+                - 500
+                - 500
             segments:
                 - straight
+                - turn_left: 45
+            direction: negative
 
         sources:
             test_source:
-                activity_MBq: 1000
+                activity_MBq: 100
-                position: [500, 100, 0]
+                position: [250, 100, 0]
-                isotope: CS137
+                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)
@@ -88,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()
@@ -100,8 +98,7 @@ def main():
             plotter.plot()
         except (
             MissingConfigKeyError,
-            KeyError,
+            KeyError
-            YAMLError,
         ) as e:
             logger.critical(e)
             logger.critical(
@@ -125,8 +122,10 @@ def main():
 
         except (
             FileNotFoundError,
-            ParserError
+            ParserError,
-        ):
+            InvalidYAMLError
+        ) as e:
+            logger.critical(e)
             sys.exit(1)
 
 

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/path/path.py

@@ -42,6 +42,7 @@ class Path:
     def __init__(
             self,
             coord_list: Sequence[tuple[float, float]],
+            opposite_direction: bool,
             z: float = 0.,
             z_box: float = 50.
                  ):
@@ -74,6 +75,8 @@ class Path:
             ]
 
         self.z = z
+        self.opposite_direction = opposite_direction
+
         self.size = (
             np.ceil(max(self.x_list)),
             np.ceil(max(self.y_list)),

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,58 +31,62 @@ 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
 
 
-def calculate_fluence_at(
+def calculate_count_rate_per_second(
     landscape: "Landscape",
     pos: np.ndarray,
-    detector: IsotropicDetector | AngularDetector,
+    detector: "Detector",
-    tangent_vectors: np.ndarray,
     scaling=1E6
 ):
-    """Compute fluence at an arbitrary position in the landscape.
+    """Compute count rate in s^-1 m^-2 at a position in the landscape.
 
     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:
-        total_phi (np.ndarray): (N,) array of fluences.
+        total_phi (np.ndarray): (N,) array of count rates per second.
     """
     pos = np.atleast_2d(pos)
     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,
@@ -100,44 +102,56 @@ def calculate_fluence_at(
     return total_phi
 
 
-def calculate_fluence_along_path(
+def calculate_counts_along_path(
     landscape: "Landscape",
-    detector: "IsotropicDetector | AngularDetector",
+    detector: "Detector",
-    points_per_segment: int = 10
+    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 (Detector): _description_
+        points_per_segment (int, optional): _description_. Defaults to 100.
+
+    Returns:
+        Tuple[np.ndarray, np.ndarray]: Array of acquisition points and
+            integrated count rates.
+    """
+
     path = landscape.path
     num_points = len(path.x_list)
+    num_segments = len(path.segments)
 
-    dx = np.diff(path.x_list)
+    segment_lengths = np.array([seg.length for seg in path.segments])
-    dy = np.diff(path.y_list)
-    segment_lengths = np.sqrt(dx**2 + dy**2)
-
     original_distances = np.zeros(num_points)
     original_distances[1:] = np.cumsum(segment_lengths)
 
     # arc lengths at which to evaluate the path
-    s = np.linspace(
+    total_subpoints = num_segments * points_per_segment
-        0,
+    s = np.linspace(0, original_distances[-1], total_subpoints)
-        original_distances[-1],
-        num=num_points * points_per_segment)
 
     # Interpolate x and y as functions of arc length
     xnew = np.interp(s, original_distances, path.x_list)
     ynew = np.interp(s, original_distances, path.y_list)
-    z = np.full(xnew.shape, path.z)
+    z = np.full_like(xnew, path.z)
     full_positions = np.c_[xnew, ynew, z]
 
-    # to compute the angle between sources and the direction of travel, we
+    if path.opposite_direction:
-    # compute tangent vectors along the path.
+        full_positions = np.flip(full_positions, axis=0)
-    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)
 
-    phi_result = calculate_fluence_at(
+    # [counts/s]
-        landscape,
+    cps = calculate_count_rate_per_second(
-        full_positions,
+        landscape, full_positions, detector
-        detector,
+    )
-        tangent_vectors)
 
-    return s, phi_result
+    # reshape so each segment is on a row
+    cps_per_seg = cps.reshape(num_segments, points_per_segment)
+
+    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
+
+    return original_distances, s, cps, int_counts_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/landscape_plotter.py

@@ -69,7 +69,7 @@ class LandscapeSlicePlotter:
         ax.set_ylabel("Y [m]")
         ax.set_title(f"Landscape (top-down, z = {self.z})")
 
-    def _draw_path(self, ax, landscape):
+    def _draw_path(self, ax, landscape: Landscape):
         if landscape.path.z <= self.z:
             ax.plot(
                 landscape.path.x_list,
@@ -79,6 +79,9 @@ class LandscapeSlicePlotter:
                 markersize=3,
                 linewidth=1
             )
+            if len(landscape.path.x_list) >= 2:
+                ax = self._draw_path_direction_arrow(ax, landscape.path)
+
         else:
             logger.warning(
                 "Path is above the slice height z."
@@ -113,3 +116,35 @@ class LandscapeSlicePlotter:
                     f"Source {s.name} is above slice height z."
                     "It will not show on the plot."
                     )
+
+    def _draw_path_direction_arrow(self, ax, path) -> Axes:
+        inset_ax = ax.inset_axes([0.8, 0.1, 0.15, 0.15])
+
+        x_start, y_start = path.x_list[0], path.y_list[0]
+        x_end, y_end = path.x_list[1], path.y_list[1]
+
+        dx = x_end - x_start
+        dy = y_end - y_start
+
+        if path.opposite_direction:
+            dx = -dx
+            dy = -dy
+
+        length = 10
+        dx_norm = dx / (dx**2 + dy**2)**0.5 * length
+        dy_norm = dy / (dx**2 + dy**2)**0.5 * length
+
+        inset_ax.arrow(
+            0, 0,
+            dx_norm, dy_norm,
+            head_width=5, head_length=5,
+            fc='red', ec='red',
+            zorder=4, linewidth=1
+        )
+
+        inset_ax.set_xlim(-2*length, 2*length)
+        inset_ax.set_ylim(-2*length, 2*length)
+        inset_ax.set_title("Direction", fontsize=8)
+        inset_ax.set_xticks([])
+        inset_ax.set_yticks([])
+        return ax

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,43 +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()
-            3,
-            2,
-            width_ratios=[0.5, 0.5],
-            height_ratios=[0.7, 0.15, 0.15],
-            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_source_table(ax4)
-
-        plt.tight_layout()
         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()
@@ -57,7 +97,42 @@ class ResultPlotter:
         cols = ('Parameter', 'Value')
         data = [
             ["Air density (kg/m^3)", round(self.landscape.air_density, 3)],
-            ["Total path length (m)", round(self.landscape.path.length, 3)]
+            ["Total path length (m)", round(self.landscape.path.length, 3)],
+            ["Readout points", len(self.count_rate_res.integrated_counts)],
+        ]
+
+        ax.table(
+            cellText=data,
+            colLabels=cols,
+            loc='center'
+        )
+
+        return ax
+
+    def _draw_detector_table(self, ax):
+        det = self.landscape.detector
+
+        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 = [
+            ["Detector", f"{det.name} ({det_type})"],
+            ["Field efficiency", formatted_effs],
         ]
 
         ax.table(
@@ -84,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)
@@ -98,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,8 +6,8 @@ from pg_rad.simulator.outputs import (
     SimulationOutput,
     SourceOutput
 )
-from pg_rad.detector.detectors import IsotropicDetector, AngularDetector
+
-from pg_rad.physics.fluence import calculate_fluence_along_path
+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,11 +39,13 @@ class SimulationEngine:
         )
 
     def _calculate_count_rate_along_path(self) -> CountRateOutput:
-        arc_length, phi = calculate_fluence_along_path(
+        acq_points, sub_points, cps, int_counts = calculate_counts_along_path(
             self.landscape,
-            self.detector
+            self.detector,
+            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]:
 
@@ -61,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", [
@@ -8,7 +8,7 @@ from pg_rad.physics import get_mass_attenuation_coeff
     (1.00000E-02, 5.120E+00),
     (1.00000E-01, 1.541E-01),
     (1.00000E+00, 6.358E-02),
-    (1.00000E+01, 2.045E-02) 
+    (1.00000E+01, 2.045E-02)
 ])
 def test_exact_attenuation_retrieval(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