Merge pull request #60 from pim-n/feature-background-activity

background CPS for 3x3 inch NaI detector based on FWHM lookup method
export functionality (to CSV)

src/pg_rad/background/background.py

@@ -0,0 +1,63 @@
+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,
+
+) -> np.ndarray:
+    """
+    Generate synthetic background cps for a given detector and energy.
+    """
+    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:
+    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/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/landscape/builder.py

@@ -145,8 +145,13 @@ 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.
+            print(pos)
+            if 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 "

src/pg_rad/main.py

@@ -17,6 +17,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():
@@ -29,9 +30,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",
@@ -39,18 +40,23 @@ 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:
-        test_yaml = """
-        name: Test landscape
+    if args.example:
+        input_config = """
+        name: Example landscape
         speed: 8.33
         acquisition_time: 1
 
@@ -66,67 +72,64 @@ def main():
         sources:
             test_source:
                 activity_MBq: 100
-                position: [250, 100, 0]
+                position: [250, 30, 0]
                 isotope: Cs137
                 gamma_energy_keV: 661
 
         detector: LU_NaI_3inch
         """
+    elif args.config:
+        input_config = args.config
 
-        cp = ConfigParser(test_yaml).parse()
+    try:
+        cp = ConfigParser(input_config).parse()
         landscape = LandscapeDirector.build_from_config(cp)
         output = SimulationEngine(
             landscape=landscape,
             runtime_spec=cp.runtime,
-            sim_spec=cp.options,
+            sim_spec=cp.options
         ).simulate()
 
         plotter = ResultPlotter(landscape, output)
-        plotter.plot()
-
-    elif args.config:
-        try:
-            cp = ConfigParser(args.config).parse()
-            landscape = LandscapeDirector.build_from_config(cp)
-            output = SimulationEngine(
-                landscape=landscape,
-                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()
-        except (
-            MissingConfigKeyError,
-            KeyError
-        ) as e:
-            logger.critical(e)
-            logger.critical(
-                "The config file is missing required keys or may be an "
-                "invalid YAML file. Check the log above. Consult the "
-                "documentation for examples of how to write a config file."
-                )
-            sys.exit(1)
-        except (
-            OutOfBoundsError,
-            DimensionError,
-            InvalidIsotopeError,
-            InvalidConfigValueError
-        ) as e:
-            logger.critical(e)
-            logger.critical(
-                "One or more items in config are not specified correctly. "
-                "Please consult this log and fix the problem."
-                )
-            sys.exit(1)
 
-        except (
-            FileNotFoundError,
-            ParserError,
-            InvalidYAMLError
-        ) as e:
-            logger.critical(e)
-            sys.exit(1)
+    except (
+        MissingConfigKeyError,
+        KeyError
+    ) as e:
+        logger.critical(e)
+        logger.critical(
+            "The config file is missing required keys or may be an "
+            "invalid YAML file. Check the log above. Consult the "
+            "documentation for examples of how to write a config file."
+            )
+        sys.exit(1)
+    except (
+        OutOfBoundsError,
+        DimensionError,
+        InvalidIsotopeError,
+        InvalidConfigValueError,
+        NotImplementedError
+    ) as e:
+        logger.critical(e)
+        logger.critical(
+            "One or more items in config are not specified correctly. "
+            "Please consult this log and fix the problem."
+            )
+        sys.exit(1)
+
+    except (
+        FileNotFoundError,
+        ParserError,
+        InvalidYAMLError
+    ) as e:
+        logger.critical(e)
+        sys.exit(1)
 
 
 if __name__ == "__main__":

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)

src/pg_rad/plotting/landscape_plotter.py

@@ -4,8 +4,6 @@ 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__)
@@ -58,12 +56,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:
-            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_ylim(bottom=-.5*width, top=.5*width)
 
         ax.set_xlabel("X [m]")
         ax.set_ylabel("Y [m]")

src/pg_rad/plotting/result_plotter.py

@@ -23,6 +23,15 @@ class ResultPlotter:
 
         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)
@@ -37,10 +46,11 @@ class ResultPlotter:
         self._draw_cps(ax_cps)
 
         ax_counts = fig.add_subplot(gs[0, 1])
-        self._draw_count_rate(ax_counts)
+        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
@@ -61,6 +71,7 @@ class ResultPlotter:
             ]
 
             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))
@@ -68,6 +79,7 @@ class ResultPlotter:
 
         self._draw_table(axs[0])
         self._draw_source_table(axs[1])
+        return fig
 
     def _plot_landscape(self, ax, z):
         lp = LandscapeSlicePlotter()
@@ -77,15 +89,20 @@ class ResultPlotter:
     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.plot(x, y, color='b', label=f'max(CPS) = {y.max():.2f}')
+        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('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)
+    def _draw_counts(self, ax):
+        x = self.count_rate_res.distance[1:]
+        y = self.count_rate_res.integrated_counts[1:]
+        ax.plot(
+            x, y, color='r', linestyle='--',
+            alpha=0.2, label=f'max(counts) = {y.max():.2f}'
+        )
+        ax.legend(handlelength=0, handletextpad=0, fancybox=True)
         ax.scatter(x, y, color='r', marker='x')
         ax.set_title('Integrated counts')
         ax.set_xlabel('Arc length s [m]')
@@ -99,6 +116,7 @@ class ResultPlotter:
             ["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)],
+            ["Mean background cps", round(self.count_rate_res.mean_bkg_cps, 3)]
         ]
 
         ax.table(
@@ -124,7 +142,7 @@ class ResultPlotter:
         # 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"
+            f"{value:.5f} @ {key:.1f} keV"
             for key, value in effs.items()
         )
         ax.set_axis_off()
@@ -201,6 +219,5 @@ class ResultPlotter:
 
         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

@@ -39,13 +39,23 @@ 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(
+            self.landscape.path.x_list,
+            self.landscape.path.y_list,
+            acq_points,
+            sub_points,
+            cps,
+            int_counts,
+            mean_bkg_counts
+        )
 
     def _calculate_point_source_distance_to_path(self) -> List[SourceOutput]:
 

src/pg_rad/simulator/outputs.py

@@ -5,10 +5,13 @@ from dataclasses import dataclass
 
 @dataclass
 class CountRateOutput:
-    acquisition_points: List[float]
+    x: List[float]
+    y: List[float]
+    distance: List[float]
     sub_points: List[float]
     cps: List[float]
     integrated_counts: List[float]
+    mean_bkg_cps: List[float]
 
 
 @dataclass

src/pg_rad/utils/export.py

@@ -0,0 +1,86 @@
+from datetime import datetime as dt
+import os
+import logging
+import re
+
+from numpy import array, full_like
+from pandas import DataFrame
+
+from pg_rad.simulator.outputs import SimulationOutput
+
+
+logger = logging.getLogger(__name__)
+
+
+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
+
+    df = generate_df(sim)
+    csv_name = generate_csv_name(sim)
+    df.to_csv(f"{folder_name}/{csv_name}.csv", index=False)
+    logger.info(f"Simulation output saved to {folder_name}!")
+
+
+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
+    )
+
+    result_df = DataFrame(
+        {
+            "East": sim.count_rate.x,
+            "North": sim.count_rate.y,
+            "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)
+    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)
+
+    csv_name = f"{num_src}_src_{bkg_cps}_cps_bkg_{src_str}"
+    return csv_name

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