mirror of
https://github.com/pim-n/pg-rad
synced 2026-05-14 00:58:10 +02:00
Compare commits
6 Commits
8098ab9329
...
591d81b8a3
| Author | SHA1 | Date | |
|---|---|---|---|
| 591d81b8a3 | |||
| 78d877c9bc | |||
| 60edbd1976 | |||
| db12d573b2 | |||
| c635c7f594 | |||
| 0a60bb09e9 |
@ -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()
|
||||
|
||||
3
src/pg_rad/data/backgrounds/dummy.csv
Normal file
3
src/pg_rad/data/backgrounds/dummy.csv
Normal file
@ -0,0 +1,3 @@
|
||||
Energy,Data,cps
|
||||
0,0,0
|
||||
10000,0,0
|
||||
|
@ -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 "
|
||||
|
||||
@ -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(
|
||||
|
||||
@ -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)
|
||||
|
||||
@ -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]")
|
||||
|
||||
@ -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")
|
||||
|
||||
@ -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]:
|
||||
|
||||
|
||||
@ -9,6 +9,7 @@ class CountRateOutput:
|
||||
sub_points: List[float]
|
||||
cps: List[float]
|
||||
integrated_counts: List[float]
|
||||
mean_bkg_cps: List[float]
|
||||
|
||||
|
||||
@dataclass
|
||||
|
||||
Reference in New Issue
Block a user