diff --git a/src/pg_rad/plotting/result_plotter.py b/src/pg_rad/plotting/result_plotter.py index 5163913..c7cd323 100644 --- a/src/pg_rad/plotting/result_plotter.py +++ b/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") diff --git a/tests/test_attenuation_functions.py b/tests/test_attenuation_functions.py index 385ace5..75de49c 100644 --- a/tests/test_attenuation_functions.py +++ b/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", [ @@ -8,7 +8,7 @@ from pg_rad.physics import get_mass_attenuation_coeff (1.00000E-02, 5.120E+00), (1.00000E-01, 1.541E-01), (1.00000E+00, 6.358E-02), - (1.00000E+01, 2.045E-02) + (1.00000E+01, 2.045E-02) ]) def test_exact_attenuation_retrieval(energy, mu): """ diff --git a/tests/test_fluence_rate.py b/tests/test_fluence_rate.py index 7220612..bc1e6cb 100644 --- a/tests/test_fluence_rate.py +++ b/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 diff --git a/tests/test_sources.py b/tests/test_sources.py index 5395269..b878da1 100644 --- a/tests/test_sources.py +++ b/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