From c635c7f5947753eaf9af8b524bde664f4ef67020 Mon Sep 17 00:00:00 2001
From: Pim Nelissen <pim.nelissen@protonmail.com>
Date: Mon, 30 Mar 2026 08:24:07 +0200
Subject: [PATCH] Add background functionality

---
 src/pg_rad/background/background.py   | 18 +++++++++++++++---
 src/pg_rad/data/backgrounds/dummy.csv |  3 +++
 src/pg_rad/main.py                    |  3 ++-
 src/pg_rad/physics/fluence.py         | 11 +++++++++--
 src/pg_rad/simulator/engine.py        | 18 +++++++++++++-----
 src/pg_rad/simulator/outputs.py       |  1 +
 6 files changed, 43 insertions(+), 11 deletions(-)
 create mode 100644 src/pg_rad/data/backgrounds/dummy.csv

diff --git a/src/pg_rad/background/background.py b/src/pg_rad/background/background.py
index e37e875..0075052 100644
--- a/src/pg_rad/background/background.py
+++ b/src/pg_rad/background/background.py
@@ -9,11 +9,18 @@ from pg_rad.detector.detector import Detector
 
 def generate_background(
     cps_array: np.ndarray,
+    detector: Detector,
+    energy_keV: float,
+    
 ) -> np.ndarray:
     """
     Generate synthetic background cps for a given detector and energy.
     """
-    pass
+    ROI_lo, ROI_hi = get_roi_from_fwhm(detector, energy_keV)
+    lam = get_cps_from_roi(detector, ROI_lo, ROI_hi)
+
+    rng = np.random.default_rng()
+    return rng.poisson(lam=lam, size=cps_array.shape)
 
 
 def fwhm(A: float, B: float, C: float, E: float) -> float:
@@ -39,8 +46,13 @@ def get_cps_from_roi(
     roi_hi: float
 ) -> float:
 
-    csv = files('pg_rad.data.backgrounds').joinpath(detector.name+'.csv')
-    data = read_csv(csv)
+    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()
diff --git a/src/pg_rad/data/backgrounds/dummy.csv b/src/pg_rad/data/backgrounds/dummy.csv
new file mode 100644
index 0000000..d01376b
--- /dev/null
+++ b/src/pg_rad/data/backgrounds/dummy.csv
@@ -0,0 +1,3 @@
+Energy,Data,cps
+0,0,0
+10000,0,0
\ No newline at end of file
diff --git a/src/pg_rad/main.py b/src/pg_rad/main.py
index 4a81386..9d2759c 100644
--- a/src/pg_rad/main.py
+++ b/src/pg_rad/main.py
@@ -111,7 +111,8 @@ def main():
             OutOfBoundsError,
             DimensionError,
             InvalidIsotopeError,
-            InvalidConfigValueError
+            InvalidConfigValueError,
+            NotImplementedError
         ) as e:
             logger.critical(e)
             logger.critical(
diff --git a/src/pg_rad/physics/fluence.py b/src/pg_rad/physics/fluence.py
index ff0a803..07b4279 100644
--- a/src/pg_rad/physics/fluence.py
+++ b/src/pg_rad/physics/fluence.py
@@ -2,6 +2,8 @@ from typing import Tuple, TYPE_CHECKING
 
 import numpy as np
 
+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
@@ -146,12 +148,17 @@ def calculate_counts_along_path(
         landscape, full_positions, detector
     )
 
+    bkg = generate_background(
+        cps, detector, landscape.point_sources[0].isotope.E
+    )
+
+    cps_with_bg = cps#  + bkg
     # reshape so each segment is on a row
-    cps_per_seg = cps.reshape(num_segments, points_per_segment)
+    cps_per_seg = cps_with_bg.reshape(num_segments, points_per_segment)
 
     du = s[1] - s[0]
     integrated_counts = np.trapezoid(cps_per_seg, dx=du, axis=1) / velocity
     int_counts_result = np.zeros(num_points)
     int_counts_result[1:] = integrated_counts
 
-    return original_distances, s, cps, int_counts_result
+    return original_distances, s, cps_with_bg, int_counts_result, np.mean(bkg)
diff --git a/src/pg_rad/simulator/engine.py b/src/pg_rad/simulator/engine.py
index 8e1a8d7..7562206 100644
--- a/src/pg_rad/simulator/engine.py
+++ b/src/pg_rad/simulator/engine.py
@@ -39,13 +39,21 @@ class SimulationEngine:
         )
 
     def _calculate_count_rate_along_path(self) -> CountRateOutput:
-        acq_points, sub_points, cps, int_counts = calculate_counts_along_path(
-            self.landscape,
-            self.detector,
-            velocity=self.runtime_spec.speed
+        acq_points, sub_points, cps, int_counts, mean_bkg_counts = (
+            calculate_counts_along_path(
+                self.landscape,
+                self.detector,
+                velocity=self.runtime_spec.speed
+            )
         )
 
-        return CountRateOutput(acq_points, sub_points, cps, int_counts)
+        return CountRateOutput(
+            acq_points,
+            sub_points,
+            cps,
+            int_counts,
+            mean_bkg_counts
+        )
 
     def _calculate_point_source_distance_to_path(self) -> List[SourceOutput]:
 
diff --git a/src/pg_rad/simulator/outputs.py b/src/pg_rad/simulator/outputs.py
index 264553f..89c597c 100644
--- a/src/pg_rad/simulator/outputs.py
+++ b/src/pg_rad/simulator/outputs.py
@@ -9,6 +9,7 @@ class CountRateOutput:
     sub_points: List[float]
     cps: List[float]
     integrated_counts: List[float]
+    mean_bkg_cps: List[float]
 
 
 @dataclass