Merge pull request #74 from pim-n/dev

BUGFIX: incorrect shape comparison causing all saves to fail

CITATION.cff

@@ -0,0 +1,26 @@
+# This CITATION.cff file was generated with cffinit.
+# Visit https://bit.ly/cffinit to generate yours today!
+
+cff-version: 1.2.0
+title: PG-RAD - Primary Gamma RADiation Simulator
+message: >-
+  If you use this software, please cite it using the
+  metadata from this file.
+type: software
+authors:
+  - given-names: Pim
+    family-names: Nelissen
+    email: pi0274ne-s@student.lu.se
+    affiliation: Lund University
+repository-code: 'https://github.com/pim-n/pg-rad'
+abstract: >-
+  Primary Gamma RADiation (PG-RAD) is a deterministic
+  landscape simulator to simulate vehicle-based gamma source
+  localization scenarios.
+keywords:
+  - gamma spectrometry
+  - mobile gamma spectrometry
+  - emergency preparedness
+  - environmental radiation
+license: MIT
+

README.md

@@ -1,59 +1,32 @@
 [![Python 3.12](https://img.shields.io/badge/python-3.12-blue.svg)](https://www.python.org/downloads/release/python-312/)
 [![Tests](https://github.com/pim-n/pg-rad/actions/workflows/ci-tests.yml/badge.svg)](https://github.com/pim-n/pg-rad/actions/workflows/ci-tests.yml)
-# pg-rad
+# PG-RAD - Primary Gamma RADiation landscape simulator
-Primary Gamma RADiation landscape - Development
 
-## Clone
+PG-RAD is a command-line software package for simulating localisation of orphan sources using mobile gamma spectrometry. PG-RAD provides a framework for construction road geometries, arbitrary distributions of point sources, modelling the detector response to primary gammas. PG-RAD provides a configurable framework for constructing detector trajectories, defining radioactive source distributions, modelling detector response, and generating synthetic acquisition data. User input is specified through YAML configuration files, which makes simulations reproducible and easily shared.
-```
+
-git clone https://github.com/pim-n/pg-rad
+## Installation
-cd pg-rad
+
-git checkout dev
+PG-RAD is a Python 3 package and is only tested on x86_64 Linux systems. The following installation instructions were testing for a fresh virtual machine running Ubuntu 26.04 LTS.
+
+## Conda installation (Tested and recommended)
+
+1. Install git by `sudo apt update && sudo apt install git`
+2. Install miniforge by following [these](https://conda-forge.org/download/) instructions.
+3. Create a file called `environment.yml`, and paste the following in there:
+
+```yaml
+name: my-pgrad-env
+channels:
+  - conda-forge
+dependencies:
+  - python=3.12
+  - pip:
+    - git+ssh://git@github.com/pim-n/pg-rad.git@v0.1.0
 ```
+4. Run `conda env create -f environment.yml` to create the environment
+5. Run `conda activate pg-rad-analysis`. You should now be in the conda environment `my-pgrad-env`.
+6. To test if installation was succesful, run `pgrad --example --show`. If this runs without errors and produces visual output, PG-RAD is correctly installed.
 
-or
+## Manual installation
 
-```
+If you prefer another virtual environment, you still need git installed. Ensure the Python version of the environment is `>3.12.4` and `<3.13`. Then, with your virtual environment activated, run `pip install git+ssh://git@github.com/pim-n/pg-rad.git@main`. Run `pgrad --example --show` to check if the installation was successful.
-git@github.com:pim-n/pg-rad.git
-cd pg-rad
-git checkout dev
-```
-
-## Dependencies / venv
-
-With Python verion `>=3.12.4` and `<3.13`, create a virtual environment and install pg-rad.
-
-```
-python3 -m venv .venv
-source .venv/bin/activate
-```
-
-With the virtual environment activated, run:
-
-```
-pip install -e .[dev]
-```
-
-## Running example landscape
-
-The example landscape can be generated using the command-line interface. Still in the virtual environment, run
-
-```
-pgrad --test --loglevel DEBUG
-```
-
-## Tests
-
-Tests can be run with `pytest` from the root directory of the repository. With the virtual environment activated, run:
-
-```
-pytest
-```
-
-## Local viewing of documentation
-
-PG-RAD uses [Material for MkDocs](https://squidfunk.github.io/mkdocs-material/) for generating documentation. It can be locally viewed by (in the venv) running:
-```
-mkdocs serve
-```
-
-where you can add the `--livereload` flag to automatically update the documentation as you write to the Markdown files.

docs/config-spec.md

@@ -2,10 +2,7 @@
     To get started quickly, you may copy and modify one of the example configs found [here](quickstart.md#example-configs).
 
 
-The config file must be a [YAML](https://yaml.org/) file. YAML is a serialization language that works with key-value pairs, but in a syntax more readable than some other alternatives. In YAML, the indentation matters. I
+The config file must be a [YAML](https://yaml.org/) file. YAML is a serialization language that works with key-value pairs, but in a syntax more readable than some other alternatives. The remainder of this chapter will explain the different required and optionals keys, what they represent, and allowed values.
-
-
-The remainder of this chapter will explain the different required and optionals keys, what they represent, and allowed values. 
 
 ## Required keys
 
@@ -124,11 +121,11 @@ Like with the lengths, if a turn segment has no angle specified, a random one (w
     Letting PG-RAD randomly assign lengths and angles can cause (expected) issues. That is because of physics restrictions. If the combination of length, angle (radius) and velocity of the vehicle is such that the centrifugal force makes it impossible to take this turn, PG-RAD will raise an error. To fix it, you can 1) reduce the speed; 2) define a smaller angle for the turn; or 3) assign more length to the turn segment.
 
 !!! info
-    For more information about how procedural roads are generated, including the random sampling of lengths and angles, see X
+    For more information about how procedural roads are generated, including the random sampling of lengths and angles, see [this](explainers/prefab_roads.ipynb) explainer.
 
 ### Sources
 
-Currently, the only type of source supported is a point source. Point sources can be added under the `sources` key, where the **subkey is the name** of the source:
+Currently, the only type of source supported is an isotropic point source. However, an arbitrary number of point sources can be added to the landscape. Point sources can be added under the `sources` key, where the **subkey is the name** of the source:
 
 ```yaml
 sources:
@@ -190,21 +187,10 @@ Note that side is relative to the direction of travel. The path will by default
 
 ### Detector
 
-The final required key is the `detector`. Currently, only isotropic detectors are supported. Nonetheless, you must specify it with `name`, `is_isotropic` and `efficiency`:
+The final required key is the `detector`. Currently, custom detectors are not yet supported and you must choose from a list of existing detectors:
 
 ```yaml
-detector:
+detector: LU_HPGe_90
-  name: test
-  is_isotropic: True
-  efficiency: 0.02
-```
-
-Note there are some existing detectors available, where efficiency is not required and will be looked up by PG-RAD itself:
-
-```yaml
-detector:
-  name: NaIR
-  is_isotropic: True
 ```
 
 ## Optional keys

docs/explainers/count_rate_along_path.ipynb

@@ -0,0 +1,121 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "a8d303ad",
+   "metadata": {},
+   "source": [
+    "# Gamma detectors along a path"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "08dda386",
+   "metadata": {},
+   "source": [
+    "## Fluence rate at $\\vec{r}$\n",
+    "\n",
+    "Let $\\vec{r}_{p} = (x_{p},y_{p},z_{p})$ denote the location of a point source $p$. Let $\\vec{r}_{i} = (x_{i},y_{i},z_{i})$ denote an arbitrary point in space. The primary photon fluence rate at $\\vec{r}$ is then given by\n",
+    "\n",
+    "$$\n",
+    "\\dot{\\phi}(r) = \\frac{A n_\\gamma \\exp(-\\mu_{air} r)}{4\\pi r^2}\n",
+    "$$\n",
+    "\n",
+    "where $r = ||\\vec{r}_p - \\vec{r}_i ||$. The units are $\\dot{\\phi} \\sim \\frac{\\text{photons}}{s \\cdot m^2}$\n",
+    "\n",
+    "## Count rate\n",
+    "\n",
+    "Gamma detectors are not perfectly efficient and efficiency is dependent on both photon energy $E_\\gamma$ and incident angle $\\theta$ [1].\n",
+    "\n",
+    "- the field efficiency $\\varepsilon_D (E_\\gamma) \\in [0, 1]$, in units of area $\\text{m}^2$,\n",
+    "- the relative angular efficiency $\\varepsilon_\\theta (E_\\gamma, \\theta) \\in [0, 1]$, dimensionless.\n",
+    "\n",
+    "The total efficiency of the detector is then defined as\n",
+    "\n",
+    "$$\n",
+    "\\varepsilon(E_\\gamma, \\theta) = \\varepsilon_D (E_\\gamma) \\varepsilon_\\theta (E_\\gamma, \\theta) \\; .\n",
+    "$$\n",
+    "\n",
+    "Where $\\varepsilon(E_\\gamma, \\theta) \\sim \\text{m}^2$.\n",
+    "\n",
+    "If the detector $D$ is positioned at $\\vec{r}_i$, the **count rate** becomes\n",
+    "\n",
+    "$$\n",
+    "\\dot{N}(r, E_\\gamma, \\theta) = \\varepsilon(E_\\gamma, \\theta) \\phi(r)\n",
+    "$$\n",
+    "\n",
+    "where $\\dot{N} \\sim \\frac{\\text{counts}}{s}$.\n",
+    "\n",
+    "## Acquisiton time \n",
+    "\n",
+    "The acquisition time window $t_{w}$ is the time during which counts are accumulated in the detector until readout into the digital system. A typical $t_{w}$ in mobile gamma spectrometry is 1 to 10 seconds [2]. \n",
+    "\n",
+    "## Integration of counts\n",
+    "\n",
+    "Suppose an acquisition time of $t_{w}$ seconds and a fixed velocity $v$ in meters per seconds. Let $R(u)$ describe a road of $L$ meters long in the xy-plane, described as a function of arc length $u$ in meters (distance traveled along the road), where $u \\in [0, L]$. The euclidian norm between the point $R(u)$ and point source $\\vec{r}_p$ is then\n",
+    "\n",
+    "$$\n",
+    "r(u) = || \\vec{r}_p - R(u) ||\n",
+    "$$\n",
+    "\n",
+    "Assuming a fixed velocity $v$, the distance traveled during one acquisition window $t_{w}$ is $\\Delta_s \\equiv vt_{w}$ meters. The path is divided into $K = L/\\Delta s$ segments, where the $k$-th segment represents the interval\n",
+    "\n",
+    "$$\n",
+    "u \\in [(k-1) \\Delta_s, k\\Delta_s] \\; , \\; k = 1, 2, \\dots, K\n",
+    "$$\n",
+    "\n",
+    "The total count rate acquired during segment $k$-th is then\n",
+    "\n",
+    "$$\n",
+    "N_{w}(k) = \\frac{1}{v} \\int_{(k-1)\\Delta_s}^{k\\Delta_s} \\underbrace{\\dot{N}(r(u), E_\\gamma, \\theta(u))}_{\\text{CPS}} du\n",
+    "$$\n",
+    "\n",
+    "## Numerical approximation\n",
+    "\n",
+    "Let us divide each segment into $N$ equally spaced points with step size $\\Delta u = \\Delta s / N$. Applying the trapezoidal rule then gives\n",
+    "\n",
+    "$$\n",
+    "N_w(k) \\approx \\frac{\\Delta u}{v}\n",
+    "\\left[\n",
+    "\\frac{\\dot{N}_0 + \\dot{N}_N}{2}\n",
+    "+ \\sum_{n=1}^{N-1} \\dot{N}_n\n",
+    "\\right],\n",
+    "$$\n",
+    "\n",
+    "where\n",
+    "\n",
+    "$$\n",
+    "\\dot{N}_n = \\dot{N}\\big(r(u_n), E_\\gamma, \\theta(u_n)\\big), \\quad\n",
+    "u_n = (k-1)\\Delta s + n \\Delta u.\n",
+    "$$\n",
+    "\n",
+    "## References\n",
+    "\n",
+    "[1] A. Bukartas, ‘Assessment of mobile radiometry data in radiological emergencies using Bayesian statistical methods’, thesis/doccomp, Lund University, 2021. Accessed: Jan. 19, 2026. [Online]. Available: http://lup.lub.lu.se/record/4c298e71-3278-42a7-818a-6f17a5121d56\n",
+    "\n",
+    "[2] R. Finck, A. Bukartas, M. Jönsson, and C. Rääf, ‘Maximum detection distances for gamma emitting point sources in mobile gamma spectrometry’, Applied Radiation and Isotopes, vol. 184, p. 110195, Jun. 2022, doi: 10.1016/j.apradiso.2022.110195.\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": ".venv",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/quickstart.md

@@ -22,13 +22,18 @@ Primary Gamma RADiation landscape tool
 
 If you get something like `pgrad: command not found`, please consult the [installation guide](installation.md).
 
-You can run a quick test scenario as follows:
+You can run a quick test by running the example landscape as follows:
 
 ```
-pgrad --test
+pgrad --example
 ```
 
-This should produce a plot of a scenario containing a single point source and a path.
+This should produce an output like
+
+```
+INFO: Landscape built successfully: Example landscape
+WARNING: No output produced. Use --save flag to save outputs and/or --showplots to display interactive plots.
+```
 
 ## Running PG-RAD
 
@@ -38,11 +43,11 @@ In order to use the CLI for your own simulations, you need to provide a *config
 pgrad --config path/to/my_config.yml
 ```
 
-where `path/to/my_config.yml` points to your config file.
+where `path/to/my_config.yml` points to your config file. To check the results live, add the `--showplots` flag. If you want to save the results directly, then add the `--save` flag (you can use them at the same time as well).
 
 ## Example configs
 
-The easiest way is to take one of these example configs, and adjust them as needed. Alternatively, there is a detailed guide on how to write your own config file [here](config-spec.md).
+The easiest way to get started is to take one of these example configs, and adjust them as needed. Alternatively, there is a detailed guide on how to write your own config file [here](config-spec.md).
 
 === "Example 1"
 
@@ -61,15 +66,14 @@ The easiest way is to take one of these example configs, and adjust them as need
     sources:
       source1:
       activity_MBq: 1000
-      isotope: CS137
+      isotope: Cs137
+      gamma_energy_keV: 662
       position:
         along_path: 100
         dist_from_path: 50
         side: left
     
-    detector:
+    detector: dummy
-      name: dummy
-      is_isotropic: True
     ```
 
 === "Example 2"
@@ -89,21 +93,21 @@ The easiest way is to take one of these example configs, and adjust them as need
     sources:
       source1:
         activity_MBq: 1000
-        isotope: CS137
+        isotope: Cs137
+        gamma_energy_keV: 662
         position: [104.3, 32.5, 0]
       source2:
         activity_MBq: 100
-        isotope: CS137
+        isotope: Cs137
+        gamma_energy_keV: 662
         position: [0, 0, 0]
         
-    detector:
+    detector: dummy
-      name: dummy
-      is_isotropic: True
     ```
 
 === "Example 3"
 
-    This is an example of a procedural path with random apportionment of total length and random angles being assigned to turns. The parameter `alpha` is optional, and is related to randomness. A higher value leads to more uniform apportionment of lengths and a lower value to more random apportionment. More information about `alpha` can be found [here](pg-rad-config-spec.md).
+    This is an example of a procedural path with random apportionment of total length and random angles being assigned to turns. The parameter `alpha` is optional, and is related to randomness. A higher value leads to more uniform apportionment of lengths and a lower value to more random apportionment. More information about `alpha` can be found [here](explainers/prefab_roads.ipynb).
 
     ``` yaml
     name: Example 3
@@ -121,12 +125,11 @@ The easiest way is to take one of these example configs, and adjust them as need
     sources:
       source1:
         activity_MBq: 1000
-        isotope: CS137
+        isotope: Cs137
+        gamma_energy_keV: 662
         position: [0, 0, 0]
 
-    detector:
+    detector: dummy
-      name: dummy
-      is_isotropic: True
     ```
 
 === "Example 4"
@@ -148,12 +151,11 @@ The easiest way is to take one of these example configs, and adjust them as need
     sources:
       source1:
         activity_MBq: 1000
-        isotope: CS137
+        isotope: Cs137
+        gamma_energy_keV: 662
         position: [0, 0, 0]
         
-    detector:
+    detector: dummy
-      name: dummy
-      is_isotropic: True
     ```
 
 === "Example 5"
@@ -178,10 +180,9 @@ The easiest way is to take one of these example configs, and adjust them as need
     sources:
       source1:
         activity_MBq: 1000
-        isotope: CS137
+        isotope: Cs137
+        gamma_energy_keV: 662
         position: [0, 0, 0]
         
-    detector:
+    detector: dummy
-      name: dummy
-      is_isotropic: True
     ```

pyproject.toml

@@ -6,16 +6,16 @@ build-backend = "setuptools.build_meta"
 where = ["src"]
 
 [tool.setuptools.package-data]
-"pg_rad.data" = ["*.csv"]
+"pg_rad.data" = ["**/*.csv"]
 "pg_rad.configs" = ["*.yml"]
 
 [project]
-name = "pg-rad"
+name = "pgrad"
-version = "0.2.1"
+version = "0.1.0"
 authors = [
   { name="Pim Nelissen", email="pi0274ne-s@student.lu.se" },
 ]
-description = "Primary Gamma RADiation Landscape"
+description = "PG-RAD - Primary Gamma RADiation Simulator"
 readme = "README.md"
 requires-python = ">=3.12.4,<3.13"
 dependencies = [
@@ -37,4 +37,4 @@ Issues = "https://github.com/pim-n/pg-rad/issues"
 
 [project.optional-dependencies]
 
-dev = ["pytest", "mkinit", "notebook", "mkdocs-material", "mkdocstrings-python", "mkdocs-jupyter", "flake8"]
+dev = ["pytest", "notebook", "mkdocs-material", "mkdocstrings-python", "mkdocs-jupyter", "flake8"]

requirements.txt

@@ -1,6 +0,0 @@
-matplotlib>=3.9.2
-notebook>=7.2.1
-numpy>=2
-pandas>=2.3.1
-piecewise_regression==1.5.0
-pyyaml>=6.0.2

src/pg_rad/background/background.py

@@ -0,0 +1,68 @@
+from importlib.resources import files
+from typing import Tuple
+import numpy as np
+from pandas import read_csv
+
+from pg_rad.configs.defaults import FWHM_PARAMS
+from pg_rad.detector.detector import Detector
+
+
+def generate_background(
+    cps_array: np.ndarray,
+    detector: Detector,
+    energy_keV: float,
+    lam_inp: int | None = None,
+    seed: int | None = None
+
+) -> np.ndarray:
+    """
+    Generate synthetic background cps for a given detector and energy.
+    """
+    if not lam_inp:
+        ROI_lo, ROI_hi = get_roi_from_fwhm(detector, energy_keV)
+        lam = get_cps_from_roi(detector, ROI_lo, ROI_hi)
+    else:
+        lam = lam_inp
+
+    rng = np.random.default_rng(seed=seed)
+    return rng.poisson(lam=lam, size=cps_array.shape)
+
+
+def fwhm(A: float, B: float, C: float, E: float) -> float:
+    return np.sqrt(A + B * E + C * E**2)
+
+
+def get_roi_from_fwhm(
+    detector: Detector,
+    energy_keV: float
+) -> Tuple[float, float]:
+    """
+    Get the region of interest for given primary gamma energy, using
+    the 3*FWHM rule.
+    """
+    A, B, C = FWHM_PARAMS.get(detector.type)
+    delta = 3*fwhm(A, B, C, energy_keV)
+    return (energy_keV-delta, energy_keV+delta)
+
+
+def get_cps_from_roi(
+    detector: Detector,
+    roi_lo: float,
+    roi_hi: float
+) -> float:
+
+    try:
+        csv = files('pg_rad.data.backgrounds').joinpath(detector.name+'.csv')
+        data = read_csv(csv)
+    except FileNotFoundError:
+        raise NotImplementedError(
+            f"Detector {detector.name} does not have backgrounds implemented."
+        )
+
+    # get indices of nearest bins
+    idx_min = (data["Energy"] - roi_lo).abs().idxmin()
+    idx_max = (data["Energy"] - roi_hi).abs().idxmin()
+    idx_start, idx_end = sorted([idx_min, idx_max])
+
+    cps_sum = data.loc[idx_start:idx_end, "cps"].sum()
+    return cps_sum

src/pg_rad/configs/defaults.py

@@ -23,3 +23,10 @@ DETECTOR_EFFICIENCIES = {
     "NaIR": 0.0216,
     "NaIF": 0.0254
 }
+
+# A, B, C parameters for different detector types
+# for formula FWHM(E) = sqrt(A + B*E + C*E^2)
+FWHM_PARAMS = {
+    "NaI": (-70.55, 2.678, 0.000602),
+    "HPGe": (1.778, 0.001843, 0.000001)
+}

src/pg_rad/data/angular_efficiencies/LU_HPGe_90.csv

@@ -0,0 +1,38 @@
+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
+-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/angular_efficiencies/LU_NaIR.csv

@@ -0,0 +1,38 @@
+angle,662
+0,0.027
+10,0.162
+20,0.346
+30,0.517
+40,0.662
+50,0.794
+60,0.882
+70,0.947
+80,0.995
+90,1.000
+100,0.970
+110,0.895
+120,0.778
+130,0.649
+140,0.546
+150,0.477
+160,0.387
+170,0.267
+180,0.205
+-180,0.205
+-170,0.266
+-160,0.385
+-150,0.527
+-140,0.671
+-130,0.764
+-120,0.838
+-110,0.763
+-100,0.838
+-90,0.903
+-80,0.904
+-70,0.898
+-60,0.862
+-50,0.717
+-40,0.337
+-30,0.154
+-20,0.253
+-10,0.125

src/pg_rad/data/backgrounds/LU_NaI_3inch.csv

@@ -0,0 +1,257 @@
+Energy,Data,cps
+0,0,0
+11.6866,0,0
+23.3732,1653,1.181558256
+35.0598,3827,2.735525375
+46.7464,5010,3.581129378
+58.433,7880,5.632594711
+70.1196,11068,7.911365261
+81.8062,13414,9.588277341
+93.4928,14536,10.39027877
+105.179,14580,10.42172981
+116.866,14138,10.10578985
+128.553,13559,9.691922802
+140.239,12725,9.095782702
+151.926,11597,8.289492495
+163.612,10696,7.645461044
+175.299,9847,7.038598999
+186.986,9063,6.478198713
+198.672,8192,5.855611151
+210.359,7627,5.451751251
+222.045,6941,4.961401001
+233.732,6450,4.610436026
+245.419,6046,4.321658327
+257.105,5560,3.974267334
+268.792,4845,3.463187991
+280.478,4267,3.05003574
+292.165,4066,2.906361687
+303.852,3711,2.652609006
+315.538,3413,2.439599714
+327.225,3080,2.201572552
+338.911,3035,2.169406719
+350.598,2789,1.993566833
+362.285,2766,1.977126519
+373.971,2432,1.73838456
+385.658,2182,1.55968549
+397.344,2010,1.436740529
+409.031,1911,1.365975697
+420.718,1871,1.337383846
+432.404,1678,1.199428163
+444.091,1601,1.144388849
+455.777,1538,1.099356683
+467.464,1473,1.052894925
+479.151,1390,0.993566833
+490.837,1369,0.978556112
+502.524,1358,0.970693352
+514.21,1356,0.96926376
+525.897,1323,0.945675482
+537.584,1256,0.897784132
+549.27,1171,0.837026447
+560.957,1072,0.766261615
+572.643,1094,0.781987134
+584.33,1157,0.827019299
+596.017,1073,0.766976412
+607.703,1083,0.774124375
+619.39,1042,0.744817727
+631.076,995,0.711222302
+642.763,938,0.670478914
+654.45,770,0.550393138
+666.136,787,0.562544675
+677.823,761,0.543959971
+689.509,703,0.502501787
+701.196,671,0.479628306
+712.883,611,0.436740529
+724.569,599,0.428162974
+736.256,573,0.40957827
+747.942,623,0.445318084
+759.629,636,0.454610436
+771.316,634,0.453180843
+783.002,602,0.430307362
+794.689,525,0.375268049
+806.375,544,0.388849178
+818.062,557,0.39814153
+829.749,527,0.376697641
+841.435,508,0.363116512
+853.122,509,0.363831308
+864.808,478,0.341672623
+876.495,498,0.355968549
+888.182,515,0.368120086
+899.868,521,0.372408863
+911.555,539,0.385275197
+923.241,546,0.390278771
+934.928,545,0.389563974
+946.615,529,0.378127234
+958.301,489,0.349535382
+969.988,490,0.350250179
+981.674,447,0.319513939
+993.361,437,0.312365976
+1005.05,406,0.290207291
+1016.73,380,0.271622588
+1028.42,357,0.255182273
+1040.11,346,0.247319514
+1051.79,348,0.248749107
+1063.48,337,0.240886347
+1075.17,376,0.268763402
+1086.85,376,0.268763402
+1098.54,363,0.259471051
+1110.23,388,0.277340958
+1121.91,330,0.235882773
+1133.6,364,0.260185847
+1145.29,323,0.230879199
+1156.97,359,0.256611866
+1168.66,330,0.235882773
+1180.35,284,0.203002144
+1192.03,256,0.182987848
+1203.72,252,0.180128663
+1215.41,243,0.173695497
+1227.09,264,0.188706219
+1238.78,226,0.16154396
+1250.47,214,0.152966405
+1262.15,184,0.131522516
+1273.84,179,0.127948535
+1285.53,168,0.120085776
+1297.21,185,0.132237312
+1308.9,177,0.126518942
+1320.59,175,0.12508935
+1332.27,171,0.122230164
+1343.96,187,0.133666905
+1355.65,206,0.147248034
+1367.33,255,0.182273052
+1379.02,331,0.23659757
+1390.71,395,0.282344532
+1402.39,540,0.385989993
+1414.08,608,0.43459614
+1425.77,612,0.437455325
+1437.45,575,0.411007863
+1449.14,516,0.368834882
+1460.82,467,0.333809864
+1472.51,369,0.263759828
+1484.2,254,0.181558256
+1495.88,220,0.157255182
+1507.57,125,0.089349535
+1519.26,136,0.097212294
+1530.94,92,0.065761258
+1542.63,74,0.052894925
+1554.32,67,0.047891351
+1566,68,0.048606147
+1577.69,70,0.05003574
+1589.38,58,0.041458184
+1601.06,66,0.047176555
+1612.75,59,0.042172981
+1624.44,47,0.033595425
+1636.12,60,0.042887777
+1647.81,74,0.052894925
+1659.5,74,0.052894925
+1671.18,84,0.060042888
+1682.87,88,0.062902073
+1694.56,88,0.062902073
+1706.24,65,0.046461758
+1717.93,78,0.05575411
+1729.62,84,0.060042888
+1741.3,62,0.04431737
+1752.99,71,0.050750536
+1764.68,54,0.038598999
+1776.36,53,0.037884203
+1788.05,53,0.037884203
+1799.74,34,0.024303074
+1811.42,36,0.025732666
+1823.11,46,0.032880629
+1834.8,45,0.032165833
+1846.48,41,0.029306648
+1858.17,33,0.023588277
+1869.86,33,0.023588277
+1881.54,29,0.020729092
+1893.23,39,0.027877055
+1904.92,37,0.026447462
+1916.6,28,0.020014296
+1928.29,41,0.029306648
+1939.98,39,0.027877055
+1951.66,46,0.032880629
+1963.35,46,0.032880629
+1975.04,36,0.025732666
+1986.72,39,0.027877055
+1998.41,47,0.033595425
+2010.1,59,0.042172981
+2021.78,56,0.040028592
+2033.47,44,0.031451036
+2045.15,57,0.040743388
+2056.84,50,0.035739814
+2068.53,68,0.048606147
+2080.21,56,0.040028592
+2091.9,53,0.037884203
+2103.59,51,0.03645461
+2115.27,30,0.021443888
+2126.96,40,0.028591851
+2138.65,47,0.033595425
+2150.33,42,0.030021444
+2162.02,49,0.035025018
+2173.71,31,0.022158685
+2185.39,29,0.020729092
+2197.08,49,0.035025018
+2208.77,26,0.018584703
+2220.45,32,0.022873481
+2232.14,23,0.016440315
+2243.83,19,0.013581129
+2255.51,31,0.022158685
+2267.2,16,0.011436741
+2278.89,21,0.015010722
+2290.57,19,0.013581129
+2302.26,20,0.014295926
+2313.95,14,0.010007148
+2325.63,18,0.012866333
+2337.32,18,0.012866333
+2349.01,22,0.015725518
+2360.69,27,0.0192995
+2372.38,21,0.015010722
+2384.07,21,0.015010722
+2395.75,35,0.02501787
+2407.44,49,0.035025018
+2419.13,55,0.039313796
+2430.81,47,0.033595425
+2442.5,68,0.048606147
+2454.19,56,0.040028592
+2465.87,57,0.040743388
+2477.56,63,0.045032166
+2489.25,52,0.037169407
+2500.93,61,0.043602573
+2512.62,37,0.026447462
+2524.31,34,0.024303074
+2535.99,34,0.024303074
+2547.68,26,0.018584703
+2559.37,15,0.010721944
+2571.05,16,0.011436741
+2582.74,9,0.006433167
+2594.43,12,0.008577555
+2606.11,8,0.00571837
+2617.8,11,0.007862759
+2629.48,2,0.001429593
+2641.17,7,0.005003574
+2652.86,5,0.003573981
+2664.54,2,0.001429593
+2676.23,2,0.001429593
+2687.92,0,0
+2699.6,6,0.004288778
+2711.29,2,0.001429593
+2722.98,6,0.004288778
+2734.66,5,0.003573981
+2746.35,5,0.003573981
+2758.04,2,0.001429593
+2769.72,5,0.003573981
+2781.41,6,0.004288778
+2793.1,7,0.005003574
+2804.78,2,0.001429593
+2816.47,2,0.001429593
+2828.16,6,0.004288778
+2839.84,5,0.003573981
+2851.53,4,0.002859185
+2863.22,1,0.000714796
+2874.9,4,0.002859185
+2886.59,4,0.002859185
+2898.28,5,0.003573981
+2909.96,2,0.001429593
+2921.65,3,0.002144389
+2933.34,4,0.002859185
+2945.02,5,0.003573981
+2956.71,3,0.002144389
+2968.4,0,0
+2980.08,0,0

src/pg_rad/data/backgrounds/dummy.csv

@@ -0,0 +1,3 @@
+Energy,Data,cps
+0,0,0
+10000,0,0

src/pg_rad/data/detectors.csv

@@ -0,0 +1,5 @@
+name,type,is_isotropic
+dummy,NaI,true
+LU_NaI_3inch,NaI,true
+LU_HPGe_90,HPGe,false
+LU_NaIR,NaI,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_NaIR.csv

@@ -0,0 +1,5 @@
+energy_keV,field_efficiency_m2
+0,0
+661.657,0.0261
+1173.228,0.0203
+1332.492,0.0166

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

@@ -104,7 +104,7 @@ class ConfigParser:
     def _parse_options(self) -> SimulationOptionsSpec:
         options = self.config.get("options", {})
 
-        allowed = {"air_density_kg_per_m3", "seed"}
+        allowed = {"air_density_kg_per_m3", "seed", "bkg_cps"}
         self._warn_unknown_keys(
             section="options",
             provided=set(options.keys()),
@@ -116,23 +116,33 @@ class ConfigParser:
             defaults.DEFAULT_AIR_DENSITY
         )
         seed = options.get("seed")
+        bkg_cps = options.get("bkg_cps")
 
         if not isinstance(air_density, float) or air_density <= 0:
             raise InvalidConfigValueError(
                 "options.air_density_kg_per_m3 must be a positive float "
                 "in kg/m^3."
             )
+
         if (
-            seed is not None or
+            seed is not None and
             (isinstance(seed, int) and seed <= 0)
         ):
             raise InvalidConfigValueError(
                 "Seed must be a positive integer value."
             )
 
+        if bkg_cps is not None and (
+            not isinstance(bkg_cps, int) or bkg_cps < 0
+        ):
+            raise InvalidConfigValueError(
+                "Background CPS must be an integer >= 0."
+            )
+
         return SimulationOptionsSpec(
             air_density=air_density,
             seed=seed,
+            bkg_cps=bkg_cps
         )
 
     def _parse_path(self) -> PathSpec:
@@ -353,41 +363,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

@@ -18,6 +18,7 @@ class RuntimeSpec:
 class SimulationOptionsSpec:
     air_density: float
     seed: int | None = None
+    bkg_cps: int | None = None
 
 
 @dataclass
@@ -67,8 +68,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
@@ -107,7 +107,8 @@ class LandscapeBuilder:
 
     def set_point_sources(
         self,
-        *sources: AbsolutePointSourceSpec | RelativePointSourceSpec
+        *sources: AbsolutePointSourceSpec | RelativePointSourceSpec,
+        bounds_check: bool = False
     ):
         """Add one or more point sources to the world.
 
@@ -145,8 +146,12 @@ class LandscapeBuilder:
                     along_path=along_path,
                     side=s.side,
                     dist_from_path=s.dist_from_path)
-            if any(
+
-                p < 0 or p >= s for p, s in zip(pos, self._size)
+            # we dont support -x values, but negative y values are possible as
+            # the path is centered in the y direction.
+            if bounds_check and not (
+                (0 <= pos[0] <= self._size[0]) and
+                (-0.5 * self._size[1] <= pos[1] <= 0.5 * self._size[1])
             ):
                 raise OutOfBoundsError(
                     "One or more sources attempted to "
@@ -161,7 +166,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

@@ -2,9 +2,9 @@ import argparse
 import logging
 import sys
 
+from numpy.random import SeedSequence
 from pandas.errors import ParserError
 
-from pg_rad.detector.builder import DetectorBuilder
 from pg_rad.exceptions.exceptions import (
     MissingConfigKeyError,
     OutOfBoundsError,
@@ -18,6 +18,7 @@ from pg_rad.inputparser.parser import ConfigParser
 from pg_rad.landscape.director import LandscapeDirector
 from pg_rad.plotting.result_plotter import ResultPlotter
 from pg_rad.simulator.engine import SimulationEngine
+from pg_rad.utils.export import generate_folder_name, save_results
 
 
 def main():
@@ -30,9 +31,9 @@ def main():
         help="Build from a config file."
     )
     parser.add_argument(
-        "--test",
+        "--example",
         action="store_true",
-        help="Load and run the test landscape."
+        help="Load and run an example landscape."
     )
     parser.add_argument(
         "--loglevel",
@@ -40,23 +41,30 @@ def main():
         choices=["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"],
     )
     parser.add_argument(
-        "--saveplot",
+        "--showplots",
         action="store_true",
-        help="Save the plot or not."
+        help="Show the plots immediately."
+    )
+    parser.add_argument(
+        "--save",
+        action="store_true",
+        help="Save the outputs"
     )
 
     args = parser.parse_args()
     setup_logger(args.loglevel)
     logger = logging.getLogger(__name__)
 
-    if args.test:
+    if args.example:
-        test_yaml = """
+        input_config = """
-        name: Test landscape
+        name: Example landscape
         speed: 8.33
         acquisition_time: 1
 
         path:
-            length: 1000
+            length:
+                - 500
+                - 500
             segments:
                 - straight
                 - turn_left: 45
@@ -65,42 +73,47 @@ def main():
         sources:
             test_source:
                 activity_MBq: 100
-                position: [250, 100, 0]
+                position: [250, 30, 0]
                 isotope: Cs137
                 gamma_energy_keV: 661
 
-        detector:
+        detector: LU_NaI_3inch
-            name: dummy
+
-            is_isotropic: True
+        options:
+          seed: 1234
         """
-
-        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)
-        plotter.plot()
-
     elif args.config:
+        input_config = args.config
+    else:
+        logger.warning(
+            "No input provided. Try --example or --config path/to/config.yml. "
+        )
+        sys.exit(1)
     try:
-            cp = ConfigParser(args.config).parse()
+        cp = ConfigParser(input_config).parse()
+        if cp.options.seed is None:
+            entr = SeedSequence().entropy
+            cp.options.seed = int(str(entr)[:6])
         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()
 
         plotter = ResultPlotter(landscape, output)
+        if args.save:
+            folder_name = generate_folder_name(output)
+            save_results(output, folder_name)
+            plotter.save(folder_name)
+        if args.showplots:
             plotter.plot()
+
+        if not (args.save or args.showplots):
+            logger.warning(
+                "No output produced. Use --save flag to save outputs and/or "
+                "--showplots to display interactive plots."
+            )
     except (
         MissingConfigKeyError,
         KeyError
@@ -116,7 +129,8 @@ def main():
         OutOfBoundsError,
         DimensionError,
         InvalidIsotopeError,
-            InvalidConfigValueError
+        InvalidConfigValueError,
+        NotImplementedError
     ) as e:
         logger.critical(e)
         logger.critical(

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,11 @@ from typing import Tuple, TYPE_CHECKING
 
 import numpy as np
 
-from pg_rad.detector.detectors import IsotropicDetector, AngularDetector
+from pg_rad.background.background import generate_background
-
 
 if TYPE_CHECKING:
     from pg_rad.landscape.landscape import Landscape
+    from pg_rad.detector.detector import Detector
 
 
 def phi(
@@ -33,16 +33,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 +48,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 +56,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 +66,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,16 +106,20 @@ def calculate_count_rate_per_second(
 
 def calculate_counts_along_path(
     landscape: "Landscape",
-    detector: "IsotropicDetector | AngularDetector",
+    detector: "Detector",
-    acquisition_time: int,
+    velocity: float,
+    t_acq: float,
     points_per_segment: int = 10,
+    bkg_cps_input: int | None = None,
+    seed: int | None = None
 ) -> 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.
+        bkg_cps_input (int | None, optional): Optional background CPS.
 
     Returns:
         Tuple[np.ndarray, np.ndarray]: Array of acquisition points and
@@ -123,7 +131,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 +147,46 @@ 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)
+    if bkg_cps_input is None:
-
+        bkg = generate_background(
-    count_rate_segs = count_rate.reshape(num_segments, points_per_segment)
+            cps, detector, landscape.point_sources[0].isotope.E,
-    integrated = np.trapezoid(
+            seed=seed
-        count_rate_segs,
+        )
-        dx=ds/points_per_segment,
+    elif bkg_cps_input == 0:
-        axis=1
+        bkg = bkg_cps_input
+    else:
+        bkg = generate_background(
+            cps, detector, landscape.point_sources[0].isotope.E,
+            lam_inp=bkg_cps_input,
+            seed=seed
         )
 
-    result = np.zeros(num_points)
+    cps_with_bg = cps + bkg
-    result[1:] = integrated
+
-    return original_distances, result
+    # Integrate along time dimension. see e.g. Bukartas (2021 doccomp.)
+    t = s / abs(velocity)
+
+    # acquisition bins
+    t_bins = np.arange(0, t[-1] + t_acq, t_acq)
+
+    int_counts = np.zeros(len(t_bins) - 1)
+    acq_points = np.zeros(len(t_bins) - 1)
+
+    for i in range(len(t_bins) - 1):
+        mask = (t >= t_bins[i]) & (t < t_bins[i + 1])
+
+        int_counts[i] = np.trapezoid(cps_with_bg[mask], t[mask])
+        acq_points[i] = np.mean(s[mask])
+
+    return (
+        acq_points[:-1],
+        s[:-1],
+        cps_with_bg[:-1],
+        int_counts[:-1],
+        np.mean(bkg)
+        )

src/pg_rad/plotting/landscape_plotter.py

@@ -4,11 +4,10 @@ from matplotlib import pyplot as plt
 from matplotlib.axes import Axes
 from matplotlib.patches import Circle
 
-from numpy import median
-
 from pg_rad.landscape.landscape import Landscape
 
 logger = logging.getLogger(__name__)
+plt.set_loglevel(level='warning')
 
 
 class LandscapeSlicePlotter:
@@ -58,12 +57,7 @@ class LandscapeSlicePlotter:
         ax.set_xlim(right=max(width, .5*height))
 
         # if the road is very flat, we center it vertically (looks better)
-        if median(landscape.path.y_list) == 0:
+        ax.set_ylim(bottom=-.5*width, top=.5*width)
-            h = max(height, .5*width)
-            ax.set_ylim(bottom=-h//2,
-                        top=h//2)
-        else:
-            ax.set_ylim(top=max(height, .5*width))
 
         ax.set_xlabel("X [m]")
         ax.set_ylabel("Y [m]")

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
 
@@ -6,6 +10,9 @@ from pg_rad.simulator.outputs import SimulationOutput
 from pg_rad.landscape.landscape import Landscape
 
 
+plt.set_loglevel(level='warning')
+
+
 class ResultPlotter:
     def __init__(self, landscape: Landscape, output: SimulationOutput):
         self.landscape = landscape
@@ -13,45 +20,109 @@ 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 save(self, path: str, landscape_z: float = 0) -> None:
+        fig_1 = self._plot_main(landscape_z)
+        fig_2 = self._plot_detector()
+        fig_3 = self._plot_metadata()
+
+        fig_1.savefig(path+"/main.jpg")
+        fig_2.savefig(path+"/detector.jpg")
+        fig_3.savefig(path+"/metadata.jpg")
+
+    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_counts(ax_counts)
+
+        ax_landscape = fig.add_subplot(gs[1, :])
+        self._plot_landscape(ax_landscape, landscape_z)
+        return fig
+
+    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])
+        return fig
+
+    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])
+        return fig
+
     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', label=f'max(CPS) = {y.max():.2f}')
-        ax.set_title('Count rate')
+        ax.legend(handlelength=0, handletextpad=0, fancybox=True)
+        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_counts(self, ax):
+        x = self.count_rate_res.distance[1:]
+        y = self.count_rate_res.integrated_counts[1:]
+
+        yerr = np.sqrt(y)
+        peak_idx = y.argmax()
+
+        ax.fill_between(
+            x,
+            y - yerr,
+            y + yerr,
+            color='red',
+            alpha=0.15,
+        )
+
+        ax.plot(
+            x, y,
+            color='r', linestyle='--', alpha=0.2,
+            label=rf'max(counts)={y[peak_idx]:.2f} $\pm$ {yerr[peak_idx]:.2f}'
+        )
+        ax.legend(handlelength=0, handletextpad=0, fancybox=True)
+        ax.scatter(x, y, color='r', marker='x')
+        ax.set_title('Counts')
+        ax.set_xlabel('Arc length s [m]')
+        ax.set_ylabel('N')
 
     def _draw_table(self, ax):
         ax.set_axis_off()
@@ -60,6 +131,9 @@ 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)],
+            ["Seed", self.count_rate_res.seed],
+            ["Mean background cps", round(self.count_rate_res.mean_bkg_cps, 3)]
         ]
 
         ax.table(
@@ -72,13 +146,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:.5f} @ {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(
@@ -119,3 +208,33 @@ 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)
+
+        ax.plot(theta_rad, eff)
+        ax.set_title(f"Rel. angular efficiency @ {energy_keV:.1f} keV")

src/pg_rad/simulator/engine.py

@@ -3,10 +3,11 @@ from typing import List
 from pg_rad.landscape.landscape import Landscape
 from pg_rad.simulator.outputs import (
     CountRateOutput,
+    DetectorOutput,
     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 +18,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
 
@@ -32,20 +32,38 @@ class SimulationEngine:
 
         count_rate_results = self._calculate_count_rate_along_path()
         source_results = self._calculate_point_source_distance_to_path()
+        detector_results = self._generate_detector_output()
 
         return SimulationOutput(
             name=self.landscape.name,
+            size=self.landscape.size,
             count_rate=count_rate_results,
+            detector=detector_results,
             sources=source_results
         )
 
     def _calculate_count_rate_along_path(self) -> CountRateOutput:
-        arc_length, phi = calculate_counts_along_path(
+        acq_points, sub_points, cps, int_counts, mean_bkg_counts = (
+            calculate_counts_along_path(
                 self.landscape,
                 self.detector,
-            acquisition_time=self.runtime_spec.acquisition_time
+                velocity=self.runtime_spec.speed,
+                bkg_cps_input=self.sim_spec.bkg_cps,
+                t_acq=self.runtime_spec.acquisition_time,
+                seed=self.sim_spec.seed
+            )
+        )
+
+        return CountRateOutput(
+            self.landscape.path.x_list,
+            self.landscape.path.y_list,
+            acq_points,
+            sub_points,
+            cps,
+            int_counts,
+            mean_bkg_counts,
+            self.sim_spec.seed
         )
-        return CountRateOutput(arc_length, phi)
 
     def _calculate_point_source_distance_to_path(self) -> List[SourceOutput]:
 
@@ -69,3 +87,13 @@ class SimulationEngine:
             )
 
         return source_output
+
+    def _generate_detector_output(self) -> DetectorOutput:
+        return DetectorOutput(
+            name=self.detector.name,
+            type=self.detector.type,
+            is_isotropic=self.detector.is_isotropic,
+            field_eff=self.detector.get_efficiency(
+                self.landscape.point_sources[0].isotope.E
+            )
+        )

src/pg_rad/simulator/outputs.py

@@ -5,8 +5,14 @@ from dataclasses import dataclass
 
 @dataclass
 class CountRateOutput:
-    arc_length: List[float]
+    x: List[float]
-    count_rate: List[float]
+    y: List[float]
+    distance: List[float]
+    sub_points: List[float]
+    cps: List[float]
+    integrated_counts: List[float]
+    mean_bkg_cps: List[float]
+    seed: int
 
 
 @dataclass
@@ -19,8 +25,18 @@ class SourceOutput:
     dist_from_path: float
 
 
+@dataclass
+class DetectorOutput:
+    name: str
+    type: str
+    is_isotropic: bool
+    field_eff: float
+
+
 @dataclass
 class SimulationOutput:
     name: str
+    size: tuple
+    detector: DetectorOutput
     count_rate: CountRateOutput
     sources: List[SourceOutput]

src/pg_rad/utils/export.py

@@ -0,0 +1,127 @@
+from dataclasses import asdict
+from datetime import datetime as dt
+import json
+import os
+import logging
+import re
+
+from numpy import array, full_like, ndarray, bool_
+from pandas import DataFrame
+
+from pg_rad.simulator.outputs import SimulationOutput
+
+
+logger = logging.getLogger(__name__)
+logging.getLogger("PIL").setLevel(logging.WARNING)
+
+
+class NumpyEncoder(json.JSONEncoder):
+    def default(self, obj):
+        if isinstance(obj, ndarray):
+            return obj.tolist()
+        elif isinstance(obj, bool_):
+            return bool(obj)
+        return super().default(obj)
+
+
+def generate_folder_name(sim: SimulationOutput) -> str:
+    formatted_sim_name = re.sub(r"\s+", '_',  sim.name.lower())
+    folder_name = (
+        formatted_sim_name +
+        '_result_' +
+        dt.today().strftime('%Y%m%d_%H%M')
+    )
+    return folder_name
+
+
+def save_results(sim: SimulationOutput, folder_name: str) -> None:
+    """Parse all simulation output and save to a folder."""
+    if not os.path.exists(folder_name):
+        os.makedirs(folder_name)
+    else:
+        logger.warning("Folder already exists. Overwrite?")
+        ans = input("[type 'n' to cancel overwrite] ")
+        if ans.lower() == 'n':
+            return
+
+    logger.debug(
+        f"Integrated counts: {list(sim.count_rate.integrated_counts)}"
+    )
+    logger.debug(
+        f"Distances: {list(sim.count_rate.distance)}"
+    )
+
+    df = generate_df(sim)
+    csv_name = generate_csv_name(sim)
+    df.to_csv(f"{folder_name}/{csv_name}.csv", index=False)
+    param_dict = generate_sim_param_dict(sim)
+    with open(f"{folder_name}/parameters.json", 'w') as f:
+        json.dump(param_dict, f, cls=NumpyEncoder)
+    logger.info(f"Simulation output saved to {folder_name}!")
+
+
+def generate_sim_param_dict(sim: SimulationOutput) -> dict:
+    """Parse simulation parameters and hyperparameters to dictionary."""
+    d = asdict(sim)
+    d.pop('count_rate')
+    return d
+
+
+def generate_df(sim: SimulationOutput) -> DataFrame:
+    """Parse simulation output to CSV format and the name of CSV."""
+
+    br_array = full_like(
+        sim.count_rate.integrated_counts,
+        sim.count_rate.mean_bkg_cps
+    )
+
+    east_coords = sim.count_rate.x[1:]
+    north_coords = sim.count_rate.y[1:]
+
+    if len(east_coords) != sim.count_rate.integrated_counts.shape[0]:
+        east_coords = None
+        north_coords = None
+        logger.warning(
+            "PG-RAD currently does not support this experimental path."
+            " Only ROI_P, ROI_BR and Dist will be saved."
+        )
+
+    result_df = DataFrame(
+        {
+            "East": east_coords,
+            "North": north_coords,
+            "ROI_P": sim.count_rate.integrated_counts,
+            "ROI_BR": br_array,
+            "Dist": sim.count_rate.distance
+        }
+    )
+
+    return result_df
+
+
+def generate_csv_name(sim: SimulationOutput) -> str:
+    """Generate CSV name according to Alex' specification"""
+    num_src = len(sim.sources)
+    src_ids = [str(i+1) for i in range(num_src)]
+    bkg_cps = round(sim.count_rate.mean_bkg_cps)
+    source_param_strings = [
+        [
+            str(round(s.activity))+"MBq",
+            str(round(s.dist_from_path))+"m",
+            str(round(s.position[0]))+'_'+str(round(s.position[1]))
+        ]
+        for s in sim.sources
+    ]
+
+    if num_src == 1:
+        src_str = "_".join(source_param_strings[0])
+    else:
+        src_str_array = array(
+            [list(item) for item in zip(*source_param_strings)]
+        )
+
+        src_str = "_".join(src_str_array.flat)
+
+    src_ids_str = "_".join(src_ids)
+    csv_name = f"{src_ids_str}_src_{bkg_cps}_cps_bkg_{src_str}"
+    return csv_name

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_roi.py

@@ -0,0 +1,42 @@
+import pytest
+
+from pg_rad.background.background import get_roi_from_fwhm, get_cps_from_roi
+from pg_rad.detector.detector import load_detector
+
+E_gamma = 661.7  # keV
+
+
+@pytest.fixture
+def detector():
+    return load_detector("LU_NaI_3inch")
+
+
+@pytest.fixture
+def roi_Cs137_NaI(detector):
+    return get_roi_from_fwhm(detector, E_gamma)
+
+
+@pytest.mark.parametrize("ref_roi_lo, ref_roi_hi", [
+    (537.58, 792.69)
+])
+def test_roi_acquisition_Cs137_NaI(ref_roi_lo, ref_roi_hi, roi_Cs137_NaI):
+    """Test against ROI used in GammaVision"""
+    est_roi_lo, est_roi_hi = roi_Cs137_NaI
+
+    assert pytest.approx(est_roi_lo, rel=0.05) == ref_roi_lo
+    assert pytest.approx(est_roi_hi, rel=0.05) == ref_roi_hi
+
+
+@pytest.mark.parametrize("ref_cps", [
+    13.61
+])
+def test_cps_acquisition_Cs137_NaI(
+    ref_cps,
+    roi_Cs137_NaI,
+    detector
+):
+    """Test against cps from GammaVision"""
+    est_roi_lo, est_roi_hi = roi_Cs137_NaI
+    est_cps = get_cps_from_roi(detector, est_roi_lo, est_roi_hi)
+
+    assert pytest.approx(est_cps, rel=0.1) == ref_cps

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