Compare commits

...

6 Commits

Author SHA1 Message Date
591d81b8a3 rename test to example landscape. Fixes #51 2026-03-30 08:52:08 +02:00
78d877c9bc PEP8 2026-03-30 08:49:18 +02:00
60edbd1976 fix typos / PEP8 2026-03-30 08:48:37 +02:00
db12d573b2 improve out of bounds handling of source placement 2026-03-30 08:48:15 +02:00
c635c7f594 Add background functionality 2026-03-30 08:24:07 +02:00
0a60bb09e9 improve plotting functionality 2026-03-30 08:23:02 +02:00
9 changed files with 71 additions and 34 deletions

View File

@ -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()

View File

@ -0,0 +1,3 @@
Energy,Data,cps
0,0,0
10000,0,0
1 Energy Data cps
2 0 0 0
3 10000 0 0

View File

@ -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 "

View File

@ -29,9 +29,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",
@ -48,9 +48,9 @@ def main():
setup_logger(args.loglevel)
logger = logging.getLogger(__name__)
if args.test:
test_yaml = """
name: Test landscape
if args.example:
example_yaml = """
name: Example landscape
speed: 8.33
acquisition_time: 1
@ -66,14 +66,14 @@ 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
"""
cp = ConfigParser(test_yaml).parse()
cp = ConfigParser(example_yaml).parse()
landscape = LandscapeDirector.build_from_config(cp)
output = SimulationEngine(
landscape=landscape,
@ -111,7 +111,8 @@ def main():
OutOfBoundsError,
DimensionError,
InvalidIsotopeError,
InvalidConfigValueError
InvalidConfigValueError,
NotImplementedError
) as e:
logger.critical(e)
logger.critical(

View File

@ -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)

View File

@ -58,12 +58,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]")

View File

@ -37,7 +37,7 @@ 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)
@ -77,15 +77,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.acquisition_points[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 +104,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 +130,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 +207,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")

View File

@ -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]:

View File

@ -9,6 +9,7 @@ class CountRateOutput:
sub_points: List[float]
cps: List[float]
integrated_counts: List[float]
mean_bkg_cps: List[float]
@dataclass