mirror of
https://github.com/pim-n/pg-rad
synced 2026-03-23 21:58:12 +01:00
Compare commits
4 Commits
b882f20358
...
09e74051f0
| Author | SHA1 | Date | |
|---|---|---|---|
| 09e74051f0 | |||
| 7dec54fa1c | |||
| 938b3a7afc | |||
| b82196e431 |
@ -10,6 +10,10 @@ class InvalidCSVError(DataLoadError):
|
|||||||
"""Raised when a file is not a valid CSV."""
|
"""Raised when a file is not a valid CSV."""
|
||||||
|
|
||||||
|
|
||||||
|
class InvalidYAMLError(DataLoadError):
|
||||||
|
"""Raised when a file is not a valid YAML."""
|
||||||
|
|
||||||
|
|
||||||
class OutOfBoundsError(Exception):
|
class OutOfBoundsError(Exception):
|
||||||
"""Raised when an object is attempted to be placed out of bounds."""
|
"""Raised when an object is attempted to be placed out of bounds."""
|
||||||
|
|
||||||
@ -18,7 +22,11 @@ class MissingConfigKeyError(KeyError):
|
|||||||
"""Raised when a (nested) config key is missing in the config."""
|
"""Raised when a (nested) config key is missing in the config."""
|
||||||
def __init__(self, key, subkey=None):
|
def __init__(self, key, subkey=None):
|
||||||
if subkey:
|
if subkey:
|
||||||
self.message = f"Missing key in {key}: {', '.join(list(subkey))}"
|
if isinstance(subkey, str):
|
||||||
|
pass
|
||||||
|
elif isinstance(subkey, set):
|
||||||
|
subkey = ', '.join(list(subkey))
|
||||||
|
self.message = f"Missing key in {key}: {subkey}"
|
||||||
else:
|
else:
|
||||||
self.message = f"Missing key: {key}"
|
self.message = f"Missing key: {key}"
|
||||||
|
|
||||||
|
|||||||
@ -7,7 +7,8 @@ import yaml
|
|||||||
from pg_rad.exceptions.exceptions import (
|
from pg_rad.exceptions.exceptions import (
|
||||||
MissingConfigKeyError,
|
MissingConfigKeyError,
|
||||||
DimensionError,
|
DimensionError,
|
||||||
InvalidConfigValueError
|
InvalidConfigValueError,
|
||||||
|
InvalidYAMLError
|
||||||
)
|
)
|
||||||
from pg_rad.configs import defaults
|
from pg_rad.configs import defaults
|
||||||
|
|
||||||
@ -45,7 +46,7 @@ class ConfigParser:
|
|||||||
|
|
||||||
def parse(self) -> SimulationSpec:
|
def parse(self) -> SimulationSpec:
|
||||||
self._warn_unknown_keys(
|
self._warn_unknown_keys(
|
||||||
section="global",
|
section="root",
|
||||||
provided=set(self.config.keys()),
|
provided=set(self.config.keys()),
|
||||||
allowed=self._ALLOWED_ROOT_KEYS,
|
allowed=self._ALLOWED_ROOT_KEYS,
|
||||||
)
|
)
|
||||||
@ -74,7 +75,7 @@ class ConfigParser:
|
|||||||
data = yaml.safe_load(config_source)
|
data = yaml.safe_load(config_source)
|
||||||
|
|
||||||
if not isinstance(data, dict):
|
if not isinstance(data, dict):
|
||||||
raise ValueError(
|
raise InvalidYAMLError(
|
||||||
"Provided path or string is not a valid YAML representation."
|
"Provided path or string is not a valid YAML representation."
|
||||||
)
|
)
|
||||||
|
|
||||||
@ -86,13 +87,13 @@ class ConfigParser:
|
|||||||
name=self.config["name"]
|
name=self.config["name"]
|
||||||
)
|
)
|
||||||
except KeyError as e:
|
except KeyError as e:
|
||||||
raise MissingConfigKeyError("global", {"name"}) from e
|
raise MissingConfigKeyError("root", {"name"}) from e
|
||||||
|
|
||||||
def _parse_runtime(self) -> RuntimeSpec:
|
def _parse_runtime(self) -> RuntimeSpec:
|
||||||
required = {"speed", "acquisition_time"}
|
required = {"speed", "acquisition_time"}
|
||||||
missing = required - self.config.keys()
|
missing = required - self.config.keys()
|
||||||
if missing:
|
if missing:
|
||||||
raise MissingConfigKeyError("global", missing)
|
raise MissingConfigKeyError("root", missing)
|
||||||
|
|
||||||
return RuntimeSpec(
|
return RuntimeSpec(
|
||||||
speed=float(self.config["speed"]),
|
speed=float(self.config["speed"]),
|
||||||
@ -116,7 +117,7 @@ class ConfigParser:
|
|||||||
seed = options.get("seed")
|
seed = options.get("seed")
|
||||||
|
|
||||||
if not isinstance(air_density, float) or air_density <= 0:
|
if not isinstance(air_density, float) or air_density <= 0:
|
||||||
raise ValueError(
|
raise InvalidConfigValueError(
|
||||||
"options.air_density_kg_per_m3 must be a positive float "
|
"options.air_density_kg_per_m3 must be a positive float "
|
||||||
"in kg/m^3."
|
"in kg/m^3."
|
||||||
)
|
)
|
||||||
@ -124,7 +125,9 @@ class ConfigParser:
|
|||||||
seed is not None or
|
seed is not None or
|
||||||
(isinstance(seed, int) and seed <= 0)
|
(isinstance(seed, int) and seed <= 0)
|
||||||
):
|
):
|
||||||
raise ValueError("Seed must be a positive integer value.")
|
raise InvalidConfigValueError(
|
||||||
|
"Seed must be a positive integer value."
|
||||||
|
)
|
||||||
|
|
||||||
return SimulationOptionsSpec(
|
return SimulationOptionsSpec(
|
||||||
air_density=air_density,
|
air_density=air_density,
|
||||||
@ -133,14 +136,26 @@ class ConfigParser:
|
|||||||
|
|
||||||
def _parse_path(self) -> PathSpec:
|
def _parse_path(self) -> PathSpec:
|
||||||
allowed_csv = {"file", "east_col_name", "north_col_name", "z"}
|
allowed_csv = {"file", "east_col_name", "north_col_name", "z"}
|
||||||
allowed_proc = {"segments", "length", "z", "alpha"}
|
allowed_proc = {"segments", "length", "z", "alpha", "direction"}
|
||||||
|
|
||||||
path = self.config.get("path")
|
path = self.config.get("path")
|
||||||
|
|
||||||
|
direction = path.get("direction", 'positive')
|
||||||
|
|
||||||
|
if direction == 'positive':
|
||||||
|
opposite_direction = False
|
||||||
|
elif direction == 'negative':
|
||||||
|
opposite_direction = True
|
||||||
|
else:
|
||||||
|
raise InvalidConfigValueError(
|
||||||
|
"Direction must be positive or negative."
|
||||||
|
)
|
||||||
|
|
||||||
if path is None:
|
if path is None:
|
||||||
raise MissingConfigKeyError("global", {"path"})
|
raise MissingConfigKeyError("path")
|
||||||
|
|
||||||
if not isinstance(path, dict):
|
if not isinstance(path, dict):
|
||||||
raise ValueError("Path must be a dictionary.")
|
raise InvalidConfigValueError("Path must be a dictionary.")
|
||||||
|
|
||||||
if "file" in path:
|
if "file" in path:
|
||||||
self._warn_unknown_keys(
|
self._warn_unknown_keys(
|
||||||
@ -154,26 +169,31 @@ class ConfigParser:
|
|||||||
east_col_name=path["east_col_name"],
|
east_col_name=path["east_col_name"],
|
||||||
north_col_name=path["north_col_name"],
|
north_col_name=path["north_col_name"],
|
||||||
z=path.get("z", 0),
|
z=path.get("z", 0),
|
||||||
|
opposite_direction=opposite_direction
|
||||||
)
|
)
|
||||||
|
|
||||||
if "segments" in path:
|
elif "segments" in path:
|
||||||
self._warn_unknown_keys(
|
self._warn_unknown_keys(
|
||||||
section="path (procedural)",
|
section="path (procedural)",
|
||||||
provided=set(path.keys()),
|
provided=set(path.keys()),
|
||||||
allowed=allowed_proc,
|
allowed=allowed_proc,
|
||||||
)
|
)
|
||||||
return self._parse_procedural_path(path)
|
return self._parse_procedural_path(path, opposite_direction)
|
||||||
|
|
||||||
raise ValueError("Invalid path configuration.")
|
else:
|
||||||
|
raise InvalidConfigValueError("Invalid path configuration.")
|
||||||
|
|
||||||
def _parse_procedural_path(
|
def _parse_procedural_path(
|
||||||
self,
|
self,
|
||||||
path: Dict[str, Any]
|
path: Dict[str, Any],
|
||||||
|
opposite_direction: bool
|
||||||
) -> ProceduralPathSpec:
|
) -> ProceduralPathSpec:
|
||||||
raw_segments = path.get("segments")
|
raw_segments = path.get("segments")
|
||||||
|
|
||||||
if not isinstance(raw_segments, List):
|
if not isinstance(raw_segments, List):
|
||||||
raise ValueError("path.segments must be a list of segments.")
|
raise InvalidConfigValueError(
|
||||||
|
"path.segments must be a list of segments."
|
||||||
|
)
|
||||||
|
|
||||||
raw_length = path.get("length")
|
raw_length = path.get("length")
|
||||||
|
|
||||||
@ -191,7 +211,8 @@ class ConfigParser:
|
|||||||
angles=angles,
|
angles=angles,
|
||||||
lengths=lengths,
|
lengths=lengths,
|
||||||
z=path.get("z", defaults.DEFAULT_PATH_HEIGHT),
|
z=path.get("z", defaults.DEFAULT_PATH_HEIGHT),
|
||||||
alpha=path.get("alpha", defaults.DEFAULT_ALPHA)
|
alpha=path.get("alpha", defaults.DEFAULT_ALPHA),
|
||||||
|
opposite_direction=opposite_direction
|
||||||
)
|
)
|
||||||
|
|
||||||
def _process_segment_angles(
|
def _process_segment_angles(
|
||||||
@ -209,13 +230,17 @@ class ConfigParser:
|
|||||||
|
|
||||||
elif isinstance(segment, dict):
|
elif isinstance(segment, dict):
|
||||||
if len(segment) != 1:
|
if len(segment) != 1:
|
||||||
raise ValueError("Invalid segment definition.")
|
raise InvalidConfigValueError(
|
||||||
|
"Invalid segment definition."
|
||||||
|
)
|
||||||
seg_type, angle = list(segment.items())[0]
|
seg_type, angle = list(segment.items())[0]
|
||||||
segments.append(seg_type)
|
segments.append(seg_type)
|
||||||
angles.append(angle)
|
angles.append(angle)
|
||||||
|
|
||||||
else:
|
else:
|
||||||
raise ValueError("Invalid segment entry format.")
|
raise InvalidConfigValueError(
|
||||||
|
"Invalid segment entry format."
|
||||||
|
)
|
||||||
|
|
||||||
return segments, angles
|
return segments, angles
|
||||||
|
|
||||||
@ -239,12 +264,25 @@ class ConfigParser:
|
|||||||
return segment_type in {"turn_left", "turn_right"}
|
return segment_type in {"turn_left", "turn_right"}
|
||||||
|
|
||||||
def _parse_point_sources(self) -> List[PointSourceSpec]:
|
def _parse_point_sources(self) -> List[PointSourceSpec]:
|
||||||
source_dict = self.config.get("sources", {})
|
source_dict = self.config.get("sources")
|
||||||
|
if source_dict is None:
|
||||||
|
raise MissingConfigKeyError("sources")
|
||||||
|
|
||||||
|
if not isinstance(source_dict, dict):
|
||||||
|
raise InvalidConfigValueError(
|
||||||
|
"sources must have subkeys representing point source names."
|
||||||
|
)
|
||||||
|
|
||||||
specs: List[PointSourceSpec] = []
|
specs: List[PointSourceSpec] = []
|
||||||
|
|
||||||
for name, params in source_dict.items():
|
for name, params in source_dict.items():
|
||||||
|
|
||||||
required = {"activity_MBq", "isotope", "position"}
|
required = {"activity_MBq", "isotope", "position"}
|
||||||
|
if not isinstance(params, dict):
|
||||||
|
raise InvalidConfigValueError(
|
||||||
|
f"sources.{name} is not defined correctly."
|
||||||
|
f" Must have subkeys {required}"
|
||||||
|
)
|
||||||
|
|
||||||
missing = required - params.keys()
|
missing = required - params.keys()
|
||||||
if missing:
|
if missing:
|
||||||
raise MissingConfigKeyError(name, missing)
|
raise MissingConfigKeyError(name, missing)
|
||||||
@ -309,8 +347,13 @@ class ConfigParser:
|
|||||||
return specs
|
return specs
|
||||||
|
|
||||||
def _parse_detector(self) -> DetectorSpec:
|
def _parse_detector(self) -> DetectorSpec:
|
||||||
det_dict = self.config.get("detector", {})
|
det_dict = self.config.get("detector")
|
||||||
required = {"name", "is_isotropic"}
|
required = {"name", "is_isotropic"}
|
||||||
|
if not isinstance(det_dict, dict):
|
||||||
|
raise InvalidConfigValueError(
|
||||||
|
"detector is not specified correctly. Must contain at least"
|
||||||
|
f"the subkeys {required}"
|
||||||
|
)
|
||||||
|
|
||||||
missing = required - det_dict.keys()
|
missing = required - det_dict.keys()
|
||||||
if missing:
|
if missing:
|
||||||
@ -327,7 +370,7 @@ class ConfigParser:
|
|||||||
elif eff:
|
elif eff:
|
||||||
pass
|
pass
|
||||||
else:
|
else:
|
||||||
raise ValueError(
|
raise InvalidConfigValueError(
|
||||||
f"The detector {name} not found in library. Either "
|
f"The detector {name} not found in library. Either "
|
||||||
f"specify {name}.efficiency or "
|
f"specify {name}.efficiency or "
|
||||||
"choose a detector from the following list: "
|
"choose a detector from the following list: "
|
||||||
|
|||||||
@ -22,7 +22,8 @@ class SimulationOptionsSpec:
|
|||||||
|
|
||||||
@dataclass
|
@dataclass
|
||||||
class PathSpec(ABC):
|
class PathSpec(ABC):
|
||||||
pass
|
z: int | float
|
||||||
|
opposite_direction: bool
|
||||||
|
|
||||||
|
|
||||||
@dataclass
|
@dataclass
|
||||||
@ -30,7 +31,6 @@ class ProceduralPathSpec(PathSpec):
|
|||||||
segments: list[str]
|
segments: list[str]
|
||||||
angles: list[float]
|
angles: list[float]
|
||||||
lengths: list[int | None]
|
lengths: list[int | None]
|
||||||
z: int | float
|
|
||||||
alpha: float
|
alpha: float
|
||||||
|
|
||||||
|
|
||||||
@ -39,7 +39,6 @@ class CSVPathSpec(PathSpec):
|
|||||||
file: str
|
file: str
|
||||||
east_col_name: str
|
east_col_name: str
|
||||||
north_col_name: str
|
north_col_name: str
|
||||||
z: int | float
|
|
||||||
|
|
||||||
|
|
||||||
@dataclass
|
@dataclass
|
||||||
|
|||||||
@ -12,7 +12,8 @@ from pg_rad.inputparser.specs import (
|
|||||||
SimulationSpec,
|
SimulationSpec,
|
||||||
CSVPathSpec,
|
CSVPathSpec,
|
||||||
AbsolutePointSourceSpec,
|
AbsolutePointSourceSpec,
|
||||||
RelativePointSourceSpec
|
RelativePointSourceSpec,
|
||||||
|
DetectorSpec
|
||||||
)
|
)
|
||||||
|
|
||||||
from pg_rad.path.path import Path, path_from_RT90
|
from pg_rad.path.path import Path, path_from_RT90
|
||||||
@ -30,6 +31,7 @@ class LandscapeBuilder:
|
|||||||
self._point_sources = []
|
self._point_sources = []
|
||||||
self._size = None
|
self._size = None
|
||||||
self._air_density = 1.243
|
self._air_density = 1.243
|
||||||
|
self._detector = None
|
||||||
|
|
||||||
logger.debug(f"LandscapeBuilder initialized: {self.name}")
|
logger.debug(f"LandscapeBuilder initialized: {self.name}")
|
||||||
|
|
||||||
@ -58,10 +60,12 @@ class LandscapeBuilder:
|
|||||||
self,
|
self,
|
||||||
sim_spec: SimulationSpec
|
sim_spec: SimulationSpec
|
||||||
):
|
):
|
||||||
segments = sim_spec.path.segments
|
path = sim_spec.path
|
||||||
lengths = sim_spec.path.lengths
|
segments = path.segments
|
||||||
angles = sim_spec.path.angles
|
lengths = path.lengths
|
||||||
alpha = sim_spec.path.alpha
|
angles = path.angles
|
||||||
|
alpha = path.alpha
|
||||||
|
z = path.z
|
||||||
|
|
||||||
sg = SegmentedRoadGenerator(
|
sg = SegmentedRoadGenerator(
|
||||||
ds=sim_spec.runtime.speed * sim_spec.runtime.acquisition_time,
|
ds=sim_spec.runtime.speed * sim_spec.runtime.acquisition_time,
|
||||||
@ -76,7 +80,11 @@ class LandscapeBuilder:
|
|||||||
alpha=alpha
|
alpha=alpha
|
||||||
)
|
)
|
||||||
|
|
||||||
self._path = Path(list(zip(x, y)))
|
self._path = Path(
|
||||||
|
list(zip(x, y)),
|
||||||
|
z=z,
|
||||||
|
opposite_direction=path.opposite_direction
|
||||||
|
)
|
||||||
self._fit_landscape_to_path()
|
self._fit_landscape_to_path()
|
||||||
return self
|
return self
|
||||||
|
|
||||||
@ -89,7 +97,8 @@ class LandscapeBuilder:
|
|||||||
df=df,
|
df=df,
|
||||||
east_col=spec.east_col_name,
|
east_col=spec.east_col_name,
|
||||||
north_col=spec.north_col_name,
|
north_col=spec.north_col_name,
|
||||||
z=spec.z
|
z=spec.z,
|
||||||
|
opposite_direction=spec.opposite_direction
|
||||||
)
|
)
|
||||||
|
|
||||||
self._fit_landscape_to_path()
|
self._fit_landscape_to_path()
|
||||||
@ -151,6 +160,10 @@ class LandscapeBuilder:
|
|||||||
name=s.name
|
name=s.name
|
||||||
))
|
))
|
||||||
|
|
||||||
|
def set_detector(self, spec: DetectorSpec) -> Self:
|
||||||
|
self._detector = spec
|
||||||
|
return self
|
||||||
|
|
||||||
def _fit_landscape_to_path(self) -> None:
|
def _fit_landscape_to_path(self) -> None:
|
||||||
"""The size of the landscape will be updated if
|
"""The size of the landscape will be updated if
|
||||||
1) _size is not set, or
|
1) _size is not set, or
|
||||||
@ -221,6 +234,7 @@ class LandscapeBuilder:
|
|||||||
name=self.name,
|
name=self.name,
|
||||||
path=self._path,
|
path=self._path,
|
||||||
point_sources=self._point_sources,
|
point_sources=self._point_sources,
|
||||||
|
detector=self._detector,
|
||||||
size=self._size,
|
size=self._size,
|
||||||
air_density=self._air_density
|
air_density=self._air_density
|
||||||
)
|
)
|
||||||
|
|||||||
@ -45,5 +45,6 @@ class LandscapeDirector:
|
|||||||
sim_spec=config,
|
sim_spec=config,
|
||||||
)
|
)
|
||||||
lb.set_point_sources(*config.point_sources)
|
lb.set_point_sources(*config.point_sources)
|
||||||
|
lb.set_detector(config.detector)
|
||||||
landscape = lb.build()
|
landscape = lb.build()
|
||||||
return landscape
|
return landscape
|
||||||
|
|||||||
@ -1,6 +1,7 @@
|
|||||||
import logging
|
import logging
|
||||||
|
|
||||||
from pg_rad.path.path import Path
|
from pg_rad.path.path import Path
|
||||||
|
from pg_rad.detector.detectors import AngularDetector, IsotropicDetector
|
||||||
from pg_rad.objects.sources import PointSource
|
from pg_rad.objects.sources import PointSource
|
||||||
|
|
||||||
|
|
||||||
@ -17,6 +18,7 @@ class Landscape:
|
|||||||
path: Path,
|
path: Path,
|
||||||
air_density: float,
|
air_density: float,
|
||||||
point_sources: list[PointSource],
|
point_sources: list[PointSource],
|
||||||
|
detector: IsotropicDetector | AngularDetector,
|
||||||
size: tuple[int, int, int]
|
size: tuple[int, int, int]
|
||||||
):
|
):
|
||||||
"""Initialize a landscape.
|
"""Initialize a landscape.
|
||||||
@ -34,5 +36,6 @@ class Landscape:
|
|||||||
self.point_sources = point_sources
|
self.point_sources = point_sources
|
||||||
self.size = size
|
self.size = size
|
||||||
self.air_density = air_density
|
self.air_density = air_density
|
||||||
|
self.detector = detector
|
||||||
|
|
||||||
logger.debug(f"Landscape created: {self.name}")
|
logger.debug(f"Landscape created: {self.name}")
|
||||||
|
|||||||
@ -3,7 +3,6 @@ import logging
|
|||||||
import sys
|
import sys
|
||||||
|
|
||||||
from pandas.errors import ParserError
|
from pandas.errors import ParserError
|
||||||
from yaml import YAMLError
|
|
||||||
|
|
||||||
from pg_rad.detector.builder import DetectorBuilder
|
from pg_rad.detector.builder import DetectorBuilder
|
||||||
from pg_rad.exceptions.exceptions import (
|
from pg_rad.exceptions.exceptions import (
|
||||||
@ -11,7 +10,8 @@ from pg_rad.exceptions.exceptions import (
|
|||||||
OutOfBoundsError,
|
OutOfBoundsError,
|
||||||
DimensionError,
|
DimensionError,
|
||||||
InvalidConfigValueError,
|
InvalidConfigValueError,
|
||||||
InvalidIsotopeError
|
InvalidIsotopeError,
|
||||||
|
InvalidYAMLError
|
||||||
)
|
)
|
||||||
from pg_rad.logger.logger import setup_logger
|
from pg_rad.logger.logger import setup_logger
|
||||||
from pg_rad.inputparser.parser import ConfigParser
|
from pg_rad.inputparser.parser import ConfigParser
|
||||||
@ -59,12 +59,14 @@ def main():
|
|||||||
length: 1000
|
length: 1000
|
||||||
segments:
|
segments:
|
||||||
- straight
|
- straight
|
||||||
|
- turn_left
|
||||||
|
direction: negative
|
||||||
|
|
||||||
sources:
|
sources:
|
||||||
test_source:
|
test_source:
|
||||||
activity_MBq: 1000
|
activity_MBq: 100
|
||||||
position: [500, 100, 0]
|
position: [250, 100, 0]
|
||||||
isotope: CS137
|
isotope: Cs137
|
||||||
|
|
||||||
detector:
|
detector:
|
||||||
name: dummy
|
name: dummy
|
||||||
@ -100,8 +102,7 @@ def main():
|
|||||||
plotter.plot()
|
plotter.plot()
|
||||||
except (
|
except (
|
||||||
MissingConfigKeyError,
|
MissingConfigKeyError,
|
||||||
KeyError,
|
KeyError
|
||||||
YAMLError,
|
|
||||||
) as e:
|
) as e:
|
||||||
logger.critical(e)
|
logger.critical(e)
|
||||||
logger.critical(
|
logger.critical(
|
||||||
@ -125,8 +126,10 @@ def main():
|
|||||||
|
|
||||||
except (
|
except (
|
||||||
FileNotFoundError,
|
FileNotFoundError,
|
||||||
ParserError
|
ParserError,
|
||||||
):
|
InvalidYAMLError
|
||||||
|
) as e:
|
||||||
|
logger.critical(e)
|
||||||
sys.exit(1)
|
sys.exit(1)
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
@ -42,6 +42,7 @@ class Path:
|
|||||||
def __init__(
|
def __init__(
|
||||||
self,
|
self,
|
||||||
coord_list: Sequence[tuple[float, float]],
|
coord_list: Sequence[tuple[float, float]],
|
||||||
|
opposite_direction: bool,
|
||||||
z: float = 0.,
|
z: float = 0.,
|
||||||
z_box: float = 50.
|
z_box: float = 50.
|
||||||
):
|
):
|
||||||
@ -74,6 +75,8 @@ class Path:
|
|||||||
]
|
]
|
||||||
|
|
||||||
self.z = z
|
self.z = z
|
||||||
|
self.opposite_direction = opposite_direction
|
||||||
|
|
||||||
self.size = (
|
self.size = (
|
||||||
np.ceil(max(self.x_list)),
|
np.ceil(max(self.x_list)),
|
||||||
np.ceil(max(self.y_list)),
|
np.ceil(max(self.y_list)),
|
||||||
|
|||||||
@ -47,14 +47,14 @@ def phi(
|
|||||||
return phi_r
|
return phi_r
|
||||||
|
|
||||||
|
|
||||||
def calculate_fluence_at(
|
def calculate_count_rate_per_second(
|
||||||
landscape: "Landscape",
|
landscape: "Landscape",
|
||||||
pos: np.ndarray,
|
pos: np.ndarray,
|
||||||
detector: IsotropicDetector | AngularDetector,
|
detector: IsotropicDetector | AngularDetector,
|
||||||
tangent_vectors: np.ndarray,
|
tangent_vectors: np.ndarray,
|
||||||
scaling=1E6
|
scaling=1E6
|
||||||
):
|
):
|
||||||
"""Compute fluence at an arbitrary position in the landscape.
|
"""Compute count rate in s^-1 m^-2 at a position in the landscape.
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
landscape (Landscape): The landscape to compute.
|
landscape (Landscape): The landscape to compute.
|
||||||
@ -63,7 +63,7 @@ def calculate_fluence_at(
|
|||||||
Detector object, needed to compute correct efficiency.
|
Detector object, needed to compute correct efficiency.
|
||||||
|
|
||||||
Returns:
|
Returns:
|
||||||
total_phi (np.ndarray): (N,) array of fluences.
|
total_phi (np.ndarray): (N,) array of count rates per second.
|
||||||
"""
|
"""
|
||||||
pos = np.atleast_2d(pos)
|
pos = np.atleast_2d(pos)
|
||||||
total_phi = np.zeros(pos.shape[0])
|
total_phi = np.zeros(pos.shape[0])
|
||||||
@ -100,33 +100,46 @@ def calculate_fluence_at(
|
|||||||
return total_phi
|
return total_phi
|
||||||
|
|
||||||
|
|
||||||
def calculate_fluence_along_path(
|
def calculate_counts_along_path(
|
||||||
landscape: "Landscape",
|
landscape: "Landscape",
|
||||||
detector: "IsotropicDetector | AngularDetector",
|
detector: "IsotropicDetector | AngularDetector",
|
||||||
points_per_segment: int = 10
|
acquisition_time: int,
|
||||||
|
points_per_segment: int = 10,
|
||||||
) -> Tuple[np.ndarray, np.ndarray]:
|
) -> Tuple[np.ndarray, np.ndarray]:
|
||||||
|
"""Compute the counts recorded in each acquisition period in the landscape.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
landscape (Landscape): _description_
|
||||||
|
detector (IsotropicDetector | AngularDetector): _description_
|
||||||
|
points_per_segment (int, optional): _description_. Defaults to 100.
|
||||||
|
|
||||||
|
Returns:
|
||||||
|
Tuple[np.ndarray, np.ndarray]: Array of acquisition points and
|
||||||
|
integrated count rates.
|
||||||
|
"""
|
||||||
|
|
||||||
path = landscape.path
|
path = landscape.path
|
||||||
num_points = len(path.x_list)
|
num_points = len(path.x_list)
|
||||||
|
num_segments = len(path.segments)
|
||||||
|
|
||||||
dx = np.diff(path.x_list)
|
segment_lengths = np.array([seg.length for seg in path.segments])
|
||||||
dy = np.diff(path.y_list)
|
ds = segment_lengths[0]
|
||||||
segment_lengths = np.sqrt(dx**2 + dy**2)
|
|
||||||
|
|
||||||
original_distances = np.zeros(num_points)
|
original_distances = np.zeros(num_points)
|
||||||
original_distances[1:] = np.cumsum(segment_lengths)
|
original_distances[1:] = np.cumsum(segment_lengths)
|
||||||
|
|
||||||
# arc lengths at which to evaluate the path
|
# arc lengths at which to evaluate the path
|
||||||
s = np.linspace(
|
total_subpoints = num_segments * points_per_segment
|
||||||
0,
|
s = np.linspace(0, original_distances[-1], total_subpoints)
|
||||||
original_distances[-1],
|
|
||||||
num=num_points * points_per_segment)
|
|
||||||
|
|
||||||
# Interpolate x and y as functions of arc length
|
# Interpolate x and y as functions of arc length
|
||||||
xnew = np.interp(s, original_distances, path.x_list)
|
xnew = np.interp(s, original_distances, path.x_list)
|
||||||
ynew = np.interp(s, original_distances, path.y_list)
|
ynew = np.interp(s, original_distances, path.y_list)
|
||||||
z = np.full(xnew.shape, path.z)
|
z = np.full_like(xnew, path.z)
|
||||||
full_positions = np.c_[xnew, ynew, z]
|
full_positions = np.c_[xnew, ynew, z]
|
||||||
|
|
||||||
|
if path.opposite_direction:
|
||||||
|
full_positions = np.flip(full_positions, axis=0)
|
||||||
|
|
||||||
# to compute the angle between sources and the direction of travel, we
|
# to compute the angle between sources and the direction of travel, we
|
||||||
# compute tangent vectors along the path.
|
# compute tangent vectors along the path.
|
||||||
dx_ds = np.gradient(xnew, s)
|
dx_ds = np.gradient(xnew, s)
|
||||||
@ -134,10 +147,19 @@ def calculate_fluence_along_path(
|
|||||||
tangent_vectors = np.c_[dx_ds, dy_ds, np.zeros_like(dx_ds)]
|
tangent_vectors = np.c_[dx_ds, dy_ds, np.zeros_like(dx_ds)]
|
||||||
tangent_vectors /= np.linalg.norm(tangent_vectors, axis=1, keepdims=True)
|
tangent_vectors /= np.linalg.norm(tangent_vectors, axis=1, keepdims=True)
|
||||||
|
|
||||||
phi_result = calculate_fluence_at(
|
count_rate = calculate_count_rate_per_second(
|
||||||
landscape,
|
landscape, full_positions, detector, tangent_vectors
|
||||||
full_positions,
|
)
|
||||||
detector,
|
|
||||||
tangent_vectors)
|
|
||||||
|
|
||||||
return s, phi_result
|
count_rate *= (acquisition_time / points_per_segment)
|
||||||
|
|
||||||
|
count_rate_segs = count_rate.reshape(num_segments, points_per_segment)
|
||||||
|
integrated = np.trapezoid(
|
||||||
|
count_rate_segs,
|
||||||
|
dx=ds/points_per_segment,
|
||||||
|
axis=1
|
||||||
|
)
|
||||||
|
|
||||||
|
result = np.zeros(num_points)
|
||||||
|
result[1:] = integrated
|
||||||
|
return original_distances, result
|
||||||
|
|||||||
@ -69,7 +69,7 @@ class LandscapeSlicePlotter:
|
|||||||
ax.set_ylabel("Y [m]")
|
ax.set_ylabel("Y [m]")
|
||||||
ax.set_title(f"Landscape (top-down, z = {self.z})")
|
ax.set_title(f"Landscape (top-down, z = {self.z})")
|
||||||
|
|
||||||
def _draw_path(self, ax, landscape):
|
def _draw_path(self, ax, landscape: Landscape):
|
||||||
if landscape.path.z <= self.z:
|
if landscape.path.z <= self.z:
|
||||||
ax.plot(
|
ax.plot(
|
||||||
landscape.path.x_list,
|
landscape.path.x_list,
|
||||||
@ -79,6 +79,9 @@ class LandscapeSlicePlotter:
|
|||||||
markersize=3,
|
markersize=3,
|
||||||
linewidth=1
|
linewidth=1
|
||||||
)
|
)
|
||||||
|
if len(landscape.path.x_list) >= 2:
|
||||||
|
ax = self._draw_path_direction_arrow(ax, landscape.path)
|
||||||
|
|
||||||
else:
|
else:
|
||||||
logger.warning(
|
logger.warning(
|
||||||
"Path is above the slice height z."
|
"Path is above the slice height z."
|
||||||
@ -113,3 +116,35 @@ class LandscapeSlicePlotter:
|
|||||||
f"Source {s.name} is above slice height z."
|
f"Source {s.name} is above slice height z."
|
||||||
"It will not show on the plot."
|
"It will not show on the plot."
|
||||||
)
|
)
|
||||||
|
|
||||||
|
def _draw_path_direction_arrow(self, ax, path) -> Axes:
|
||||||
|
inset_ax = ax.inset_axes([0.8, 0.1, 0.15, 0.15])
|
||||||
|
|
||||||
|
x_start, y_start = path.x_list[0], path.y_list[0]
|
||||||
|
x_end, y_end = path.x_list[1], path.y_list[1]
|
||||||
|
|
||||||
|
dx = x_end - x_start
|
||||||
|
dy = y_end - y_start
|
||||||
|
|
||||||
|
if path.opposite_direction:
|
||||||
|
dx = -dx
|
||||||
|
dy = -dy
|
||||||
|
|
||||||
|
length = 10
|
||||||
|
dx_norm = dx / (dx**2 + dy**2)**0.5 * length
|
||||||
|
dy_norm = dy / (dx**2 + dy**2)**0.5 * length
|
||||||
|
|
||||||
|
inset_ax.arrow(
|
||||||
|
0, 0,
|
||||||
|
dx_norm, dy_norm,
|
||||||
|
head_width=5, head_length=5,
|
||||||
|
fc='red', ec='red',
|
||||||
|
zorder=4, linewidth=1
|
||||||
|
)
|
||||||
|
|
||||||
|
inset_ax.set_xlim(-2*length, 2*length)
|
||||||
|
inset_ax.set_ylim(-2*length, 2*length)
|
||||||
|
inset_ax.set_title("Direction", fontsize=8)
|
||||||
|
inset_ax.set_xticks([])
|
||||||
|
inset_ax.set_yticks([])
|
||||||
|
return ax
|
||||||
|
|||||||
@ -16,10 +16,10 @@ class ResultPlotter:
|
|||||||
fig = plt.figure(figsize=(12, 10), constrained_layout=True)
|
fig = plt.figure(figsize=(12, 10), constrained_layout=True)
|
||||||
fig.suptitle(self.landscape.name)
|
fig.suptitle(self.landscape.name)
|
||||||
gs = GridSpec(
|
gs = GridSpec(
|
||||||
3,
|
4,
|
||||||
2,
|
2,
|
||||||
width_ratios=[0.5, 0.5],
|
width_ratios=[0.5, 0.5],
|
||||||
height_ratios=[0.7, 0.15, 0.15],
|
height_ratios=[0.7, 0.1, 0.1, 0.1],
|
||||||
hspace=0.2)
|
hspace=0.2)
|
||||||
|
|
||||||
ax1 = fig.add_subplot(gs[0, 0])
|
ax1 = fig.add_subplot(gs[0, 0])
|
||||||
@ -32,9 +32,11 @@ class ResultPlotter:
|
|||||||
self._draw_table(ax3)
|
self._draw_table(ax3)
|
||||||
|
|
||||||
ax4 = fig.add_subplot(gs[2, :])
|
ax4 = fig.add_subplot(gs[2, :])
|
||||||
self._draw_source_table(ax4)
|
self._draw_detector_table(ax4)
|
||||||
|
|
||||||
|
ax5 = fig.add_subplot(gs[3, :])
|
||||||
|
self._draw_source_table(ax5)
|
||||||
|
|
||||||
plt.tight_layout()
|
|
||||||
plt.show()
|
plt.show()
|
||||||
|
|
||||||
def _plot_landscape(self, ax, z):
|
def _plot_landscape(self, ax, z):
|
||||||
@ -57,7 +59,26 @@ class ResultPlotter:
|
|||||||
cols = ('Parameter', 'Value')
|
cols = ('Parameter', 'Value')
|
||||||
data = [
|
data = [
|
||||||
["Air density (kg/m^3)", round(self.landscape.air_density, 3)],
|
["Air density (kg/m^3)", round(self.landscape.air_density, 3)],
|
||||||
["Total path length (m)", round(self.landscape.path.length, 3)]
|
["Total path length (m)", round(self.landscape.path.length, 3)],
|
||||||
|
]
|
||||||
|
|
||||||
|
ax.table(
|
||||||
|
cellText=data,
|
||||||
|
colLabels=cols,
|
||||||
|
loc='center'
|
||||||
|
)
|
||||||
|
|
||||||
|
return ax
|
||||||
|
|
||||||
|
def _draw_detector_table(self, ax):
|
||||||
|
det = self.landscape.detector
|
||||||
|
print(det)
|
||||||
|
ax.set_axis_off()
|
||||||
|
ax.set_title('Detector')
|
||||||
|
cols = ('Parameter', 'Value')
|
||||||
|
data = [
|
||||||
|
["Name", det.name],
|
||||||
|
["Type", det.efficiency],
|
||||||
]
|
]
|
||||||
|
|
||||||
ax.table(
|
ax.table(
|
||||||
|
|||||||
@ -7,7 +7,7 @@ from pg_rad.simulator.outputs import (
|
|||||||
SourceOutput
|
SourceOutput
|
||||||
)
|
)
|
||||||
from pg_rad.detector.detectors import IsotropicDetector, AngularDetector
|
from pg_rad.detector.detectors import IsotropicDetector, AngularDetector
|
||||||
from pg_rad.physics.fluence import calculate_fluence_along_path
|
from pg_rad.physics.fluence import calculate_counts_along_path
|
||||||
from pg_rad.utils.projection import minimal_distance_to_path
|
from pg_rad.utils.projection import minimal_distance_to_path
|
||||||
from pg_rad.inputparser.specs import RuntimeSpec, SimulationOptionsSpec
|
from pg_rad.inputparser.specs import RuntimeSpec, SimulationOptionsSpec
|
||||||
|
|
||||||
@ -40,9 +40,10 @@ class SimulationEngine:
|
|||||||
)
|
)
|
||||||
|
|
||||||
def _calculate_count_rate_along_path(self) -> CountRateOutput:
|
def _calculate_count_rate_along_path(self) -> CountRateOutput:
|
||||||
arc_length, phi = calculate_fluence_along_path(
|
arc_length, phi = calculate_counts_along_path(
|
||||||
self.landscape,
|
self.landscape,
|
||||||
self.detector
|
self.detector,
|
||||||
|
acquisition_time=self.runtime_spec.acquisition_time
|
||||||
)
|
)
|
||||||
return CountRateOutput(arc_length, phi)
|
return CountRateOutput(arc_length, phi)
|
||||||
|
|
||||||
|
|||||||
Reference in New Issue
Block a user