diff --git a/src/pg_rad/background/background.py b/src/pg_rad/background/background.py new file mode 100644 index 0000000..1eaa550 --- /dev/null +++ b/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 diff --git a/src/pg_rad/configs/defaults.py b/src/pg_rad/configs/defaults.py index 9730f70..9ad866e 100644 --- a/src/pg_rad/configs/defaults.py +++ b/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) +} diff --git a/src/pg_rad/data/backgrounds/LU_NaI_3inch.csv b/src/pg_rad/data/backgrounds/LU_NaI_3inch.csv new file mode 100644 index 0000000..b75af7d --- /dev/null +++ b/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 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/landscape/builder.py b/src/pg_rad/landscape/builder.py index 56ca378..40eefbf 100644 --- a/src/pg_rad/landscape/builder.py +++ b/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 " diff --git a/src/pg_rad/main.py b/src/pg_rad/main.py index 4a81386..0380f83 100644 --- a/src/pg_rad/main.py +++ b/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__": diff --git a/src/pg_rad/physics/fluence.py b/src/pg_rad/physics/fluence.py index ff0a803..be1ff2f 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/plotting/landscape_plotter.py b/src/pg_rad/plotting/landscape_plotter.py index aebfc36..3c42d8c 100644 --- a/src/pg_rad/plotting/landscape_plotter.py +++ b/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]") diff --git a/src/pg_rad/plotting/result_plotter.py b/src/pg_rad/plotting/result_plotter.py index c7cd323..ad595f2 100644 --- a/src/pg_rad/plotting/result_plotter.py +++ b/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") diff --git a/src/pg_rad/simulator/engine.py b/src/pg_rad/simulator/engine.py index 8e1a8d7..7b8c950 100644 --- a/src/pg_rad/simulator/engine.py +++ b/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]: diff --git a/src/pg_rad/simulator/outputs.py b/src/pg_rad/simulator/outputs.py index 264553f..9aca77e 100644 --- a/src/pg_rad/simulator/outputs.py +++ b/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 diff --git a/src/pg_rad/utils/export.py b/src/pg_rad/utils/export.py new file mode 100644 index 0000000..483e29f --- /dev/null +++ b/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 diff --git a/tests/test_roi.py b/tests/test_roi.py new file mode 100644 index 0000000..f4a0122 --- /dev/null +++ b/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