Add background counts for 3x3 inch NaI. Add auto ROI lookup function from the FWHM. Add ROI test.

Feature detectors roadgen

src/pg_rad/background/background.py

@@ -0,0 +1,51 @@
+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,
+) -> np.ndarray:
+    """
+    Generate synthetic background cps for a given detector and energy.
+    """
+    pass
+
+
+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:
+
+    csv = files('pg_rad.data.backgrounds').joinpath(detector.name+'.csv')
+    data = read_csv(csv)
+
+    # 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,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/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/detectors.csv

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

src/pg_rad/data/field_efficiencies/LU_HPGe_90.csv

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

src/pg_rad/data/field_efficiencies/LU_NaI_3inch.csv

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

src/pg_rad/data/field_efficiencies/dummy.csv

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

src/pg_rad/detector/builder.py

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

src/pg_rad/detector/detector.py

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

src/pg_rad/detector/detectors.py

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

src/pg_rad/inputparser/parser.py

@@ -353,41 +353,11 @@ class ConfigParser:
         return specs
 
     def _parse_detector(self) -> DetectorSpec:
-        det_dict = self.config.get("detector")
-        required = {"name", "is_isotropic"}
-        if not isinstance(det_dict, dict):
-            raise InvalidConfigValueError(
-                "detector is not specified correctly. Must contain at least"
-                f"the subkeys {required}"
-                )
+        det_name = self.config.get("detector")
+        if not det_name:
+            raise MissingConfigKeyError("detector")
 
-        missing = required - det_dict.keys()
-        if missing:
-            raise MissingConfigKeyError("detector", missing)
-
-        name = det_dict.get("name")
-        is_isotropic = det_dict.get("is_isotropic")
-        eff = det_dict.get("efficiency")
-
-        default_detectors = defaults.DETECTOR_EFFICIENCIES
-
-        if name in default_detectors.keys() and not eff:
-            eff = default_detectors[name]
-        elif eff:
-            pass
-        else:
-            raise InvalidConfigValueError(
-                f"The detector {name} not found in library. Either "
-                f"specify {name}.efficiency or "
-                "choose a detector from the following list: "
-                f"{default_detectors.keys()}."
-            )
-
-        return DetectorSpec(
-            name=name,
-            efficiency=eff,
-            is_isotropic=is_isotropic
-        )
+        return DetectorSpec(name=det_name)
 
     def _warn_unknown_keys(self, section: str, provided: set, allowed: set):
         unknown = provided - allowed

src/pg_rad/inputparser/specs.py

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

src/pg_rad/isotopes/isotope.py

@@ -4,7 +4,7 @@ from pandas import read_csv
 
 from pg_rad.configs.filepaths import ISOTOPE_TABLE
 from pg_rad.exceptions.exceptions import InvalidIsotopeError
-from pg_rad.physics.attenuation import get_mass_attenuation_coeff
+from pg_rad.utils.interpolators import get_mass_attenuation_coeff
 
 
 class Isotope:

src/pg_rad/landscape/builder.py

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

src/pg_rad/landscape/landscape.py

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

src/pg_rad/main.py

@@ -4,7 +4,6 @@ import sys
 
 from pandas.errors import ParserError
 
-from pg_rad.detector.builder import DetectorBuilder
 from pg_rad.exceptions.exceptions import (
     MissingConfigKeyError,
     OutOfBoundsError,
@@ -56,7 +55,9 @@ def main():
         acquisition_time: 1
 
         path:
-            length: 1000
+            length:
+                - 500
+                - 500
             segments:
                 - straight
                 - turn_left: 45
@@ -69,19 +70,15 @@ def main():
                 isotope: Cs137
                 gamma_energy_keV: 661
 
-        detector:
-            name: dummy
-            is_isotropic: True
+        detector: LU_NaI_3inch
         """
 
         cp = ConfigParser(test_yaml).parse()
         landscape = LandscapeDirector.build_from_config(cp)
-        detector = DetectorBuilder(cp.detector).build()
         output = SimulationEngine(
             landscape=landscape,
             runtime_spec=cp.runtime,
             sim_spec=cp.options,
-            detector=detector
         ).simulate()
 
         plotter = ResultPlotter(landscape, output)
@@ -91,10 +88,8 @@ def main():
         try:
             cp = ConfigParser(args.config).parse()
             landscape = LandscapeDirector.build_from_config(cp)
-            detector = DetectorBuilder(cp.detector).build()
             output = SimulationEngine(
                 landscape=landscape,
-                detector=detector,
                 runtime_spec=cp.runtime,
                 sim_spec=cp.options
             ).simulate()

src/pg_rad/physics/attenuation.py

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

src/pg_rad/physics/fluence.py

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

src/pg_rad/plotting/result_plotter.py

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

src/pg_rad/simulator/engine.py

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

src/pg_rad/simulator/outputs.py

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

src/pg_rad/utils/interpolators.py

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

tests/test_attenuation_functions.py

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

tests/test_fluence_rate.py

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

tests/test_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