diff --git a/docs/requirements.txt b/docs/requirements.txt index 5d23b8347..c47a35f6d 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -18,7 +18,7 @@ matplotlib GitPython colorlog numpy-stl -radioactivedecay +radioactivedecay>=0.6.1 jsonpickle pandas requests diff --git a/docs/source/user_guide/user_guide_reference_sources_generic_source.rst b/docs/source/user_guide/user_guide_reference_sources_generic_source.rst index 21a5c2a1a..769274d61 100644 --- a/docs/source/user_guide/user_guide_reference_sources_generic_source.rst +++ b/docs/source/user_guide/user_guide_reference_sources_generic_source.rst @@ -343,22 +343,28 @@ normal distribution with: Spectra """"""" -**Discrete for gamma spectrum** +**Discrete energy spectrum** One can configure a generic source to produce particles with energies depending on weights. To do so, one must provide two lists of the same size: one for energies, one for weights. Each energy is associated to the corresponding weight. Probabilities are derived from weights simply by normalizing the weights list. -Several spectra are provided through the `get_rad_gamma_spectrum` function: +1252 isotopes spectra are provided through the `get_spectrum` function: .. code:: python - spectrum = gate.sources.base.get_rad_gamma_spectrum("Lu177") + spectrum = gate.sources.utility.get_spectrum("Lu177", spectrum_type, database="icrp107") +where ``spectrum_type`` is one of "gamma", "beta-", "beta+", "alpha", "X", "neutron", +"auger", "IE", "alpha recoil", "anihilation", "fission", "betaD", "b-spectra". From this list, +only b-spectra is histogram based (see next section), the rest are discrete. ``database`` can be "icrp107" or "radar". -The source can be configured like this: +ICRP107 data comes from `[ICRP, 2008. Nuclear Decay Data for Dosimetric Calculations. ICRP Publication 107. Ann. ICRP 38] ` +with the data from the `[Supplemental material] `. +`[Direct link to the zipped data] ` +The source can be configured like this: .. code:: python @@ -388,12 +394,12 @@ The produced particles will follow this pattern: One can configure a generic source to produce particles with energies according to a given histogram. Histograms are defined in the same way as `numpy`, using bin edges and histogram values. -Several spectra are provided through the `get_rad_beta_spectrum` function. -This data comes from `[doseinfo-radar] `_ (`[direct link to the excel file] `_). +Several spectra are provided through the `get_spectrum` function. You can use icrp107 data, or radar data. +Radar data comes from `[doseinfo-radar] `_ (`[direct link to the excel file] `_). .. code:: python - spectrum = gate.sources.base.get_rad_beta_spectrum("Lu177") + spectrum = gate.sources.utility.get_spectrum("Lu177", "beta-", database="radar") The source can be configured like this: @@ -405,7 +411,7 @@ The source can be configured like this: source.energy.spectrum_energy_bin_edges = spectrum.energy_bin_edges source.energy.spectrum_weights = spectrum.weights -For example, using this (which is what you get from `get_rad_beta_spectrum("Lu177")`): +For example, using this (which is what you get from `get_spectrum("Lu177", "beta-", database="radar")`): .. code:: python diff --git a/opengate/actors/base.py b/opengate/actors/base.py index 06ca85430..491d88c80 100644 --- a/opengate/actors/base.py +++ b/opengate/actors/base.py @@ -305,7 +305,7 @@ def _get_docstring_for_interface(cls, output_name, interface_name): docstring += "\n" return docstring - def __init__(self, *args, **kwargs): + def __init__(self, *args, **kwargs) -> None: GateObject.__init__(self, *args, **kwargs) self.actor_engine = ( None # this is set by the actor engine during initialization diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index d97ddc717..ae957e0c9 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -1083,7 +1083,7 @@ def StartSimulationAction(self): DigitizerBase.StartSimulationAction(self) DigitizerAdderActor.set_group_by_depth(self) if self.user_info.discretize_volume is None: - fatal(f'Please, set the option "discretize_volume"') + fatal('Please, set the option "discretize_volume"') depth = self.simulation.volume_manager.get_volume( self.discretize_volume ).volume_depth_in_tree diff --git a/opengate/base.py b/opengate/base.py index 75484bab3..4795be2ec 100644 --- a/opengate/base.py +++ b/opengate/base.py @@ -450,7 +450,7 @@ def __new__(cls, *args, **kwargs): new_instance = super(GateObject, cls).__new__(cls) return new_instance - def __init__(self, *args, simulation=None, **kwargs): + def __init__(self, *args, simulation=None, **kwargs) -> None: self._simulation = simulation # keep internal number of raised warnings (for debug) self.number_of_warnings = 0 diff --git a/opengate/contrib/linacs/dicomrtplan.py b/opengate/contrib/linacs/dicomrtplan.py index fa0842b65..2c20708bd 100644 --- a/opengate/contrib/linacs/dicomrtplan.py +++ b/opengate/contrib/linacs/dicomrtplan.py @@ -17,7 +17,7 @@ def list_of_beam_sequence_ID(file): patient_info.value[0], [[("0x300C", "0x0004"), i], [("0x300A", "0x0086"), None]], ) - if number_of_mu != None: + if number_of_mu is not None: id_of_beam_sequence.append(i) return id_of_beam_sequence @@ -41,7 +41,7 @@ def read(file, cp_id="all_cp", arc_id=None): dose_weight_array = [] collimation_angle_array = [] leaf_array = [] - if arc_id == None: + if arc_id is None: id_of_beam_sequence = list_of_beam_sequence_ID(file) else: arc_id = int(arc_id) @@ -81,7 +81,7 @@ def read(file, cp_id="all_cp", arc_id=None): l_cp_id = np.arange(0, nb_cp_id, 1) for id in l_cp_id: for i, key in enumerate(key_list): - if data_set[i][key][id] != None: + if data_set[i][key][id] is not None: if key == "isocenter": l_parameters[i].append(np.array(data_set[i][key][id]) * mm) elif ( @@ -94,8 +94,8 @@ def read(file, cp_id="all_cp", arc_id=None): l_parameters[i].append(data_set[i][key][id]) else: if id != 0: - if data_set[i][key][id] == None: - if l_parameters[i][id - 1] != None: + if data_set[i][key][id] is None: + if l_parameters[i][id - 1] is not None: l_parameters[i].append(l_parameters[i][id - 1]) else: l_parameters[i].append(data_set[i][key][id]) @@ -120,13 +120,13 @@ def read(file, cp_id="all_cp", arc_id=None): leaf_blocks = [leaf_block_1, leaf_block_2] for i in range(len(data_set[leaves][0])): for side_ID in range(2): - if data_set[leaves][0][i]["leaves"][id] != None: + if data_set[leaves][0][i]["leaves"][id] is not None: leaf_blocks[side_ID].append( data_set[leaves][side_ID][i]["leaves"][id] ) else: if id != 0: - if data_set[leaves][0][i]["leaves"][id - 1] != None: + if data_set[leaves][0][i]["leaves"][id - 1] is not None: leaf_blocks[side_ID].append( data_set[leaves][side_ID][i]["leaves"][id - 1] ) @@ -217,7 +217,7 @@ def read(file, cp_id="all_cp", arc_id=None): def retrieve_value_in_DICOM_RTPlan(cp, list): if list[0][0] in cp: - if list[0][1] != None: + if list[0][1] is not None: value = cp[list[0][0]].value[list[0][1]] else: value = cp[list[0][0]].value @@ -254,7 +254,7 @@ def extract_dataset(file, beam_sequence_ID): data_set[dir_angle] = {"dir_angle": []} data_set[isocenter] = {"isocenter": []} data_set[dose_weights] = {"cumulative weight": []} - data_set[leaves] = [[{"leaves": []} for i in range(nb_leaf)] for j in range(2)] + data_set[leaves] = [[{"leaves": []} for _ in range(nb_leaf)] for _ in range(2)] data_set[MU_number] = {"MU number": 0} data_set[limiting_device_angle] = {"collimation angle": []} count = 0 diff --git a/opengate/contrib/phantoms/nemaiec.py b/opengate/contrib/phantoms/nemaiec.py index 43f29a4d8..358909d72 100644 --- a/opengate/contrib/phantoms/nemaiec.py +++ b/opengate/contrib/phantoms/nemaiec.py @@ -3,6 +3,7 @@ import numpy as np import math +from opengate import utility import opengate.geometry.volumes from opengate.utility import fatal, g4_units from opengate.geometry.volumes import unite_volumes @@ -52,7 +53,7 @@ def add_iec_phantom( # Inside space for the water, same as the shell, with 3 mm less thickness = 3 * mm thickness_z = 10 * mm - interior, top_interior, c = add_iec_body( + interior, top_interior, _ = add_iec_body( simulation, f"{name}_interior", thickness, thickness_z ) interior.mother = iec.name @@ -634,7 +635,7 @@ def add_iec_phantom_vox_FIXME_TO_REMOVE(sim, name, image_filename, labels_filena iec = sim.add_volume("Image", name) iec.image = image_filename iec.material = "IEC_PLASTIC" - labels = json.loads(open(labels_filename).read()) + labels = utility.read_json_file(labels_filename) iec.voxel_materials = [] create_material(sim) material_list = {} @@ -656,3 +657,32 @@ def add_iec_phantom_vox_FIXME_TO_REMOVE(sim, name, image_filename, labels_filena m = [labels[l]["label"], labels[l]["label"] + 1, mat] iec.voxel_materials.append(m) return iec, material_list + + +def create_iec_phantom_source_vox( + image_filename, labels_filename, source_filename, activities=None +): + if activities is None: + activities = { + "iec_sphere_10mm": 1.0, + "iec_sphere_13mm": 1.0, + "iec_sphere_17mm": 1.0, + "iec_sphere_22mm": 1.0, + "iec_sphere_28mm": 1.0, + "iec_sphere_37mm": 1.0, + } + + img = itk.imread(image_filename) + labels = utility.read_json_file(labels_filename) + img_arr = itk.GetArrayViewFromImage(img) + + for label in labels: + l = labels[label]["label"] + if "sphere" in label and "shell" not in label: + img_arr[img_arr == l] = activities[label] + else: + img_arr[img_arr == l] = 0 + + # The coordinate system is different from IEC analytical volume + # 35mm should be added in Y + itk.imwrite(img, source_filename) diff --git a/opengate/sources/base.py b/opengate/sources/base.py index 11ebbe6b8..d1f1cb774 100644 --- a/opengate/sources/base.py +++ b/opengate/sources/base.py @@ -1,141 +1,8 @@ -from box import Box -import os -import pathlib -import numpy as np -import json - from ..actors.base import _setter_hook_attached_to from ..base import GateObject, process_cls from ..utility import g4_units -from ..exception import fatal from ..definitions import __world_name__ -gate_source_path = pathlib.Path(__file__).parent.resolve() - -# http://www.lnhb.fr/nuclear-data/module-lara/ -all_beta_plus_radionuclides = [ - "F18", - "Ga68", - "Zr89", - "Na22", - "C11", - "N13", - "O15", - "Rb82", -] - - -def read_beta_plus_spectra(rad_name): - """ - read the file downloaded from LNHB - there are 15 lines-long header to skip - first column is E(keV) - second column is dNtot/dE b+ - WARNING : bins width is not uniform (need to scale for density) - """ - filename = ( - f"{gate_source_path}/beta_plus_spectra/{rad_name}/beta+_{rad_name}_tot.bs" - ) - data = np.genfromtxt(filename, usecols=(0, 1), skip_header=15, dtype=float) - return data - - -def compute_bins_density(bins): - """ - Given a list of (energy) bins center, compute the width of each bin. - """ - lower = np.roll(bins, 1) - lower[0] = 0 - upper = bins - dx = upper - lower - return dx - - -def get_rad_yield(rad_name): - if not rad_name in all_beta_plus_radionuclides: - return 1.0 - data = read_beta_plus_spectra(rad_name) - ene = data[:, 0] / 1000 # convert from KeV to MeV - proba = data[:, 1] - cdf, total = compute_cdf_and_total_yield(proba, ene) - total = total * 1000 # (because was in MeV) - return total - - -def compute_cdf_and_total_yield(data, bins): - """ - Compute the CDF (Cumulative Density Function) of a list of non-uniform energy bins - with associated probability. - Also return the total probability. - """ - dx = compute_bins_density(bins) - p = data * dx - total = p.sum() - cdf = np.cumsum(p) / total - return cdf, total - - -def get_rad_gamma_spectrum(rad): - path = ( - pathlib.Path(os.path.dirname(__file__)) - / ".." - / "data" - / "rad_gamma_spectrum.json" - ) - with open(path, "r") as f: - data = json.load(f) - - # consider lower case - data = {key.lower(): value for key, value in data.items()} - rad = rad.lower() - - if rad not in data: - fatal(f"get_rad_gamma_spectrum: {path} does not contain data for ion {rad}") - - # select data for specific ion - data = Box(data[rad]) - - data.energies = np.array(data.energies) * g4_units.MeV - data.weights = np.array(data.weights) - - return data - - -def get_rad_beta_spectrum(rad: str): - path = ( - pathlib.Path(os.path.dirname(__file__)) - / ".." - / "data" - / "rad_beta_spectrum.json" - ) - with open(path, "r") as f: - data = json.load(f) - - if rad not in data: - fatal(f"get_rad_beta_spectrum: {path} does not contain data for ion {rad}") - - # select data for specific ion - data = Box(data[rad]) - - bin_edges = data.energy_bin_edges - n = len(bin_edges) - 1 - energies = [(bin_edges[i] + bin_edges[i + 1]) / 2 for i in range(n)] - - data.energy_bin_edges = np.array(data.energy_bin_edges) * g4_units.MeV - data.weights = np.array(data.weights) - data.energies = np.array(energies) - - return data - - -def set_source_rad_energy_spectrum(source, rad): - rad_spectrum = get_rad_gamma_spectrum(rad) - - source.particle = "gamma" - source.energy.type = "spectrum_discrete" - source.energy.spectrum_weights = rad_spectrum.weights - source.energy.spectrum_energies = rad_spectrum.energies - class SourceBase(GateObject): """ diff --git a/opengate/sources/generic.py b/opengate/sources/generic.py index b9b3e5f74..c8209e643 100644 --- a/opengate/sources/generic.py +++ b/opengate/sources/generic.py @@ -2,10 +2,10 @@ from scipy.spatial.transform import Rotation import opengate_core as g4 -from .base import ( - SourceBase, +from .base import SourceBase +from .utility import ( all_beta_plus_radionuclides, - read_beta_plus_spectra, + get_spectrum, compute_cdf_and_total_yield, ) from ..base import process_cls @@ -228,7 +228,7 @@ def initialize(self, run_timing_intervals): "range", ] l.extend(all_beta_plus_radionuclides) - if not self.energy.type in l: + if self.energy.type not in l: fatal( f"Cannot find the energy type {self.energy.type} for the source {self.name}.\n" f"Available types are {l}" @@ -251,10 +251,10 @@ def initialize(self, run_timing_intervals): # FIXME put this elsewhere if self.particle == "e+": if self.energy.type in all_beta_plus_radionuclides: - data = read_beta_plus_spectra(self.user_info.energy.type) + data = get_spectrum(self.user_info.energy.type, "e+", "radar") ene = data[:, 0] / 1000 # convert from KeV to MeV proba = data[:, 1] - cdf, total = compute_cdf_and_total_yield(proba, ene) + cdf, _ = compute_cdf_and_total_yield(proba, ene) # total = total * 1000 # (because was in MeV) # self.user_info.activity *= total self.energy.is_cdf = True @@ -274,7 +274,7 @@ def initialize(self, run_timing_intervals): # check direction type l = ["iso", "histogram", "momentum", "focused", "beam2d"] - if not self.direction.type in l: + if self.direction.type not in l: fatal( f"Cannot find the direction type {self.direction.type} for the source {self.name}.\n" f"Available types are {l}" diff --git a/opengate/sources/utility.py b/opengate/sources/utility.py new file mode 100644 index 000000000..133e2585f --- /dev/null +++ b/opengate/sources/utility.py @@ -0,0 +1,322 @@ +from box import Box +import pathlib +import numpy as np +import json +import re +import icrp107_database + +from ..utility import g4_units +from ..exception import fatal + +gate_source_path = pathlib.Path(__file__).parent.resolve() + +# http://www.lnhb.fr/nuclear-data/module-lara/ +all_beta_plus_radionuclides = [ + "F18", + "Ga68", + "Zr89", + "Na22", + "C11", + "N13", + "O15", + "Rb82", +] + +icrp107_emissions = [ + "alpha", + "beta-", + "beta+", + "gamma", + "X", + "neutron", + "auger", + "IE", + "alpha recoil", + "anihilation", + "fission", + "betaD", + "b-spectra", # beta spectras, both beta+ and beta- +] + +DEFAULT_DATABASE = "icrp107" +DEFAULT_SPECTRUM_TYPE = "gamma" + + +def get_spectrum( + rad_name: str, spectrum_type=DEFAULT_SPECTRUM_TYPE, database=DEFAULT_DATABASE +) -> Box: + """ + Retrieve the spectrum of a given radionuclide from the specified database. + + Parameters + ---------- + rad_name : str + The name of the radionuclide in Gate format (e.g. "Tc99m", "Lu177"). + + spectrum_type : str, optional + The type of spectrum to retrieve. Default is "gamma". Must be one of + "gamma", "beta+", "beta-", or "e+". icrp107 allows also one of + "alpha", "X", "neutron", "auger", "IE", "alpha recoil", "anihilation", "fission", "betaD", "b-spectra". + In the case of beta spectras, make use of the + :py:func:`set_source_energy_spectrum` function instead. + + database : str, optional + The database to retrieve the spectrum from. Default is "icrp107". + If not "icrp107", we use radar data for gammas and beta-, and LNHB data for beta+. + + Returns + ------- + Box + A Box object containing the spectrum data. + + Raises + ------ + fatal + If the specified database does not contain the requested spectrum type. + """ + if database == "icrp107": + return __get_icrp107_spectrum(rad_name, spectrum_type) + else: + if spectrum_type == "beta+" or spectrum_type == "e+": + return __read_beta_plus_spectra(rad_name) + elif spectrum_type == "beta-" or spectrum_type == "e-": + return __get_rad_beta_spectrum(rad_name) + elif spectrum_type == "gamma": + return __get_rad_gamma_spectrum(rad_name) + else: # FIXME use icrp107 for missing spectrum types + fatal(f"databse {database} doesn't contain spectrum type {spectrum_type}") + + +def set_source_energy_spectrum(source, rad: str, database=DEFAULT_DATABASE) -> None: + """ + Set the energy spectrum of a source. + + Parameters + ---------- + source : a source + The source to set the energy spectrum for + rad : str + The name of the radionuclide to use + database : str, optional + The database to use. Default is "icrp107". + + Notes + ----- + The source particle must be set before calling this function. + If the source particle is "beta-/e-" or "beta+/e+", the function will use the + "b-spectra" spectrum from the database data. + Otherwise, the function will use the discrete spectrum for the given particle. + + """ + if ( + source.particle == "beta-" + or source.particle == "e-" + or source.particle == "beta+" + or source.particle == "e+" + ): + rad_spectrum = get_spectrum(rad, "b-spectra", database) + source.energy.type = "spectrum_histogram" + source.energy.spectrum_weights = rad_spectrum.weights[:-1] + source.energy.spectrum_energy_bin_edges = rad_spectrum.energies + else: + rad_spectrum = get_spectrum(rad, source.particle, database) + source.energy.type = "spectrum_discrete" + source.energy.spectrum_weights = rad_spectrum.weights + source.energy.spectrum_energies = rad_spectrum.energies + + +def __gate_radname_to_icrp107(rad_name: str) -> str: + """ + Convert a radionuclide name from GATE to ICRP 107 format. + + GATE format is for example "F18" and ICRP 107 format is "F-18" + + Parameters + ---------- + rad_name : str + radionuclide name in GATE format + + Returns + ------- + str + radionuclide name in ICRP 107 format + """ + excited = rad_name[-1] == "m" # Handle Tc-99m + at_num = re.findall(r"\d+", rad_name)[0] # Extract atomic number + name = rad_name[:-1] if excited else rad_name # Remove final m + elem = name.replace(at_num, "").replace("-", "") # Find element code + elem = elem[0].upper() + elem[1:] # Fix + return f'{elem}-{at_num}{"m" if excited else ""}' + + +icrp107_time_units = { + "ms": g4_units.millisecond, + "s": g4_units.second, + "m": g4_units.minute, + "h": g4_units.hour, + "d": g4_units.day, + "y": g4_units.year, +} + + +def __convert_icrp107_time_unit(icrp_time_unit: str) -> float: + """ + Convert an ICRP 107 time unit to its corresponding value in GATE units. + + Parameters + ---------- + icrp_time_unit : str + The time unit from ICRP 107. Valid options are "ms" (milliseconds), + "s" (seconds), "m" (minutes), "h" (hours), "d" (days), and "y" (years). + + Returns + ------- + float + The equivalent time value in GATE units. + + Raises + ------ + Exception + If the provided time unit is not recognized. + """ + if icrp_time_unit in icrp107_time_units: + return icrp107_time_units[icrp_time_unit] + else: + fatal(f"unit {icrp_time_unit} not recognized") + + +def __read_beta_plus_spectra(rad_name): + """ + read the file downloaded from LNHB + there are 15 lines-long header to skip + first column is E(keV) + second column is dNtot/dE b+ + WARNING : bins width is not uniform (need to scale for density) + """ + if rad_name not in all_beta_plus_radionuclides: + # FIXME use icrp107 for missing isotopes + fatal(f"rad_name {rad_name} not in {all_beta_plus_radionuclides}") + + filename = ( + f"{gate_source_path}/beta_plus_spectra/{rad_name}/beta+_{rad_name}_tot.bs" + ) + data = np.genfromtxt(filename, usecols=(0, 1), skip_header=15, dtype=float) + # FIXME convert to MeV before return + return data + + +def compute_bins_density(bins): + """ + Given a list of (energy) bins center, compute the width of each bin. + """ + lower = np.roll(bins, 1) + lower[0] = 0 + upper = bins + dx = upper - lower + return dx + + +def get_rad_yield(rad_name): + if rad_name not in all_beta_plus_radionuclides: + return 1.0 + data = __read_beta_plus_spectra(rad_name) + ene = data[:, 0] / 1000 # convert from KeV to MeV + proba = data[:, 1] + _, total = compute_cdf_and_total_yield(proba, ene) + total = total * 1000 # (because was in MeV) + return total + + +def compute_cdf_and_total_yield(data, bins): + """ + Compute the CDF (Cumulative Density Function) of a list of non-uniform energy bins + with associated probability. + Also return the total probability. + """ + dx = compute_bins_density(bins) + p = data * dx + total = p.sum() + cdf = np.cumsum(p) / total + return cdf, total + + +def __get_icrp107_spectrum(rad_name: str, spectrum_type=DEFAULT_SPECTRUM_TYPE) -> Box: + """ + Get the spectrum of a given radionuclide according to ICRP107 recommendations. + + Parameters + ---------- + rad_name : str + The name of the radionuclide in Gate format, e.g. "Tc99m", "Lu177" + + spectrum_type : str + The type of spectrum to retrieve. Must be one of "gamma", "beta-", "beta+", "alpha", "X", "neutron", "auger", "IE", "alpha recoil", "anihilation", "fission", "betaD", "b-spectra" + + Returns + ------- + Box + A Box with two keys: "energies" and "weights". The first contains the energy of each emission, the second contains the weight of each emission (summing to 1). + + Raises + ------ + fatal + If the radionuclide or spectrum type is not valid. + """ + rad = __gate_radname_to_icrp107(rad_name) + icrp107_data = icrp107_database.get_icrp107_spectrum(rad, spectrum_type) + + # convert to Box with unit + gate_data = {} + gate_data["energies"] = np.array( + [v * g4_units.MeV for v in icrp107_data["energies"]] + ) + gate_data["weights"] = np.array([v for v in icrp107_data["weights"]]) + gate_data["half_life"] = icrp107_data["half_life"] * __convert_icrp107_time_unit( + icrp107_data["time_unit"] + ) + return Box(gate_data) + + +def __get_rad_gamma_spectrum(rad): + path = gate_source_path.parent / "data" / "rad_gamma_spectrum.json" + with open(path, "r") as f: + data = json.load(f) + + # consider lower case + data = {key.lower(): value for key, value in data.items()} + rad = rad.lower() + + if rad not in data: + # FIXME use icrp107 for missing isotopes + fatal(f"get_rad_gamma_spectrum: {path} does not contain data for ion {rad}") + + # select data for specific ion + data = Box(data[rad]) + + data.energies = np.array(data.energies) * g4_units.MeV + data.weights = np.array(data.weights) + + return data + + +def __get_rad_beta_spectrum(rad: str): + path = gate_source_path.parent / "data" / "rad_beta_spectrum.json" + with open(path, "r") as f: + data = json.load(f) + + if rad not in data: + # FIXME use icrp107 for missing isotopes + fatal(f"get_rad_beta_spectrum: {path} does not contain data for ion {rad}") + + # select data for specific ion + data = Box(data[rad]) + + bin_edges = data.energy_bin_edges + n = len(bin_edges) - 1 + energies = [(bin_edges[i] + bin_edges[i + 1]) / 2 for i in range(n)] + + data.energy_bin_edges = np.array(data.energy_bin_edges) * g4_units.MeV + data.weights = np.array(data.weights) + data.energies = np.array(energies) + + return data diff --git a/opengate/tests/src/test010_generic_source_energy_spectrum_discrete.py b/opengate/tests/src/test010_generic_source_energy_spectrum_discrete.py index 87231d829..319ab97c4 100755 --- a/opengate/tests/src/test010_generic_source_energy_spectrum_discrete.py +++ b/opengate/tests/src/test010_generic_source_energy_spectrum_discrete.py @@ -3,60 +3,13 @@ import opengate as gate from opengate.tests import utility -from opengate.sources.base import get_rad_gamma_spectrum +from opengate.sources.utility import get_spectrum import numpy as np import gatetools import matplotlib.pyplot as plt -def plot(output_file, ekin, data_x, data_y, relerrs): - bins = len(data_x) - hist_y, hist_x = np.histogram(ekin, bins=bins) - - fig, ax = plt.subplots(figsize=(8.5, 6)) - ax.set_xlabel("Energy (MeV)") - ax.set_ylabel("Number of particles") - ax.scatter(data_x, data_y, label="input energy spectrum", marker="o") - ax.scatter( - data_x, hist_y, label="simulated energy spectrum", marker=".", linewidth=0.9 - ) - ax.legend() - - ax2 = ax.twinx() - ax2.set_ylabel("Relative uncertainty") - ax2.set_ylim([0, 0.05]) - ax2.plot(data_x, relerrs, label="relative uncertainty", linewidth=0.4) - - fig.savefig(output_file, dpi=300) - plt.close(fig) - - -def root_load_ekin(root_file: str): - data_ref, keys_ref, m_ref = gatetools.phsp.load(root_file) - - index_ekin = keys_ref.index("KineticEnergy") - ekin = [data_ref_i[index_ekin] for data_ref_i in data_ref] - - return ekin - - -def add_source_energy_spectrum_discrete(sim, phsp): - spectrum = get_rad_gamma_spectrum("Lu177") - - source = sim.add_source("GenericSource", "beam") - source.attached_to = phsp.name - source.particle = "gamma" - source.n = 5e5 / sim.number_of_threads - source.position.type = "point" - source.direction.type = "iso" - source.energy.type = "spectrum_discrete" - source.energy.spectrum_energies = spectrum.energies - source.energy.spectrum_weights = spectrum.weights - - return source - - def run_simulation(paths): # create the simulation sim = gate.Simulation() @@ -107,40 +60,84 @@ def run_simulation(paths): ekin = root_load_ekin(str(phsp_actor.get_output_path())) # test - data_x = source.energy.spectrum_energies + energies = source.energy.spectrum_energies + weights = source.energy.spectrum_weights - data_y = ( - np.array(source.energy.spectrum_weights) - * source.n - / np.sum(source.energy.spectrum_weights) - ) + src_data = np.array(weights) + src_data = src_data / np.sum(src_data) - energy_counts = {} + energy_counts = {energy: 0 for energy in energies} for energy in ekin: if energy not in energy_counts: energy_counts[energy] = 0 energy_counts[energy] += 1 - output_y = [energy_counts[x] for x in data_x] + sim_data = [energy_counts[x] for x in energies] + sim_data = sim_data / np.sum(sim_data) - relerrs = [abs(output_y[i] - data_y[i]) / data_y[i] for i in range(len(data_y))] + relerrs = (sim_data - src_data) / src_data oks = [ - utility.check_diff_abs(0, relerr, tolerance=0.05, txt="relative error") + utility.check_diff_abs(0, abs(relerr), tolerance=0.05, txt="relative error") for relerr in relerrs ] is_ok = all(oks) plot( paths.output / "test010_plot_energy_spectrum_discrete.png", - ekin, - data_x, - data_y, + energies, + src_data, + sim_data, relerrs, ) utility.test_ok(is_ok) +def plot(output_file, energies, src_data, sim_data, relerrs): + fig, ax = plt.subplots(figsize=(8.5, 6)) + ax.set_xlabel("Energy (MeV)") + ax.set_ylabel("Number of particles") + ax.scatter(energies, src_data, label="input energy spectrum", marker="o") + ax.scatter( + energies, sim_data, label="simulated energy spectrum", marker=".", linewidth=0.9 + ) + + ax2 = ax.twinx() + ax2.set_ylabel("Relative error") + ax2.set_ylim([-0.05, +0.05]) + ax2.plot(energies, relerrs, label="relative error", linewidth=0.4, color="C2") + + ax.legend() + + fig.savefig(output_file, dpi=300) + plt.close(fig) + + +def root_load_ekin(root_file: str): + data_ref, keys_ref, _ = gatetools.phsp.load(root_file) + + index_ekin = keys_ref.index("KineticEnergy") + ekin = [data_ref_i[index_ekin] for data_ref_i in data_ref] + + return ekin + + +def add_source_energy_spectrum_discrete(sim, phsp): + spectrum = get_spectrum("Lu177", "gamma") + + source = sim.add_source("GenericSource", "beam") + source.attached_to = phsp.name + source.particle = "gamma" + source.n = 5e5 / sim.number_of_threads + source.position.type = "point" + source.direction.type = "iso" + source.energy.type = "spectrum_discrete" + source.energy.spectrum_energies = spectrum.energies + source.energy.spectrum_weights = spectrum.weights + + return source + + if __name__ == "__main__": paths = utility.get_default_test_paths( __file__, "test010_generic_source_energy_spectrum", output_folder="test010" diff --git a/opengate/tests/src/test010_generic_source_energy_spectrum_histogram.py b/opengate/tests/src/test010_generic_source_energy_spectrum_histogram.py index dfde34954..0c3cc0ecd 100755 --- a/opengate/tests/src/test010_generic_source_energy_spectrum_histogram.py +++ b/opengate/tests/src/test010_generic_source_energy_spectrum_histogram.py @@ -3,61 +3,13 @@ import opengate as gate from opengate.tests import utility -from opengate.sources.base import get_rad_beta_spectrum +from opengate.sources.utility import set_source_energy_spectrum import numpy as np import gatetools import matplotlib.pyplot as plt -def plot(output_file, ekin, data_x, data_y, relerrs): - bins = len(data_x) - hist_y, hist_x = np.histogram(ekin, bins=bins) - - fig, ax = plt.subplots(figsize=(8.5, 6)) - ax.set_xlabel("Energy (MeV)") - ax.set_ylabel("Number of particles") - ax.plot(data_x, data_y, label="input energy spectrum", marker="o") - ax.plot( - data_x, hist_y, label="simulated energy spectrum", marker=".", linewidth=0.9 - ) - ax.legend() - - ax2 = ax.twinx() - ax2.set_ylabel("Relative uncertainty") - ax2.set_ylim([0, 0.05]) - ax2.plot(data_x, relerrs, label="relative uncertainty", linewidth=0.4) - - fig.savefig(output_file, dpi=300) - plt.close(fig) - - -def root_load_ekin(root_file: str): - data_ref, keys_ref, m_ref = gatetools.phsp.load(root_file) - - index_ekin = keys_ref.index("KineticEnergy") - ekin = [data_ref_i[index_ekin] for data_ref_i in data_ref] - - return ekin - - -def add_source_energy_spectrum_histogram(sim, phsp, interpolation: str = None): - spectrum = get_rad_beta_spectrum("Lu177") - - source = sim.add_source("GenericSource", "beam") - source.attached_to = phsp.name - source.particle = "gamma" - source.n = 5e5 / sim.number_of_threads - source.position.type = "point" - source.direction.type = "iso" - source.energy.type = "spectrum_histogram" - source.energy.spectrum_energy_bin_edges = spectrum.energy_bin_edges - source.energy.spectrum_weights = spectrum.weights - source.energy.spectrum_histogram_interpolation = interpolation - - return source - - def run_simulation(paths, interpolation: str = None): # create the simulation sim = gate.Simulation() @@ -71,7 +23,7 @@ def run_simulation(paths, interpolation: str = None): sim.output_dir = paths.output # units - m = gate.g4_units.m + mm = gate.g4_units.mm g_cm3 = gate.g4_units.g_cm3 # materials @@ -80,24 +32,19 @@ def run_simulation(paths, interpolation: str = None): ) world = sim.world - world.size = [2 * m, 2 * m, 2 * m] + world.size = [10 * mm, 10 * mm, 10 * mm] world.material = "Vacuum" - phsp = sim.add_volume("Sphere", "phsp") - phsp.rmin = 0 * m - phsp.rmax = 1 * m - phsp.material = world.material - - source = add_source_energy_spectrum_histogram(sim, phsp, interpolation) + source = add_source_energy_spectrum_histogram(sim, interpolation) # actors - stats = sim.add_actor("SimulationStatisticsActor", "Stats") + # stats = sim.add_actor("SimulationStatisticsActor", "Stats") phsp_actor = sim.add_actor("PhaseSpaceActor", "phsp_actor") phsp_actor.output_filename = ( f"test010_energy_spectrum_histogram_{interpolation}.root" ) - # phsp_actor.attached_to = phsp + phsp_actor.steps_to_store = "first" phsp_actor.attributes = [ "KineticEnergy", ] @@ -113,35 +60,81 @@ def run_simulation(paths, interpolation: str = None): # test bin_edges = source.energy.spectrum_energy_bin_edges - data_x = [(bin_edges[i] + bin_edges[i + 1]) / 2 for i in range(len(bin_edges) - 1)] + weights = source.energy.spectrum_weights - data_y = ( - np.array(source.energy.spectrum_weights) - * source.n - / np.sum(source.energy.spectrum_weights) - ) + src_data = np.array(weights) + src_data = src_data / np.sum(src_data) - bins = len(data_x) - hist_y, hist_x = np.histogram(ekin, bins=bins) + sim_data, _ = np.histogram(ekin, bins=bin_edges, density=True) + sim_data = sim_data / np.sum(sim_data) - relerrs = [abs(hist_y[i] - data_y[i]) / data_y[i] for i in range(len(data_y))] + relerrs = (sim_data - src_data) / src_data oks = [ - utility.check_diff_abs(0, relerr, tolerance=0.05, txt="relative error") - for relerr in relerrs + utility.check_diff_abs(0, abs(relerrs[i]), tolerance=0.20, txt="relative error") + for i in range(len(relerrs)) + if src_data[i] > 0.001 ] is_ok = all(oks) plot( paths.output / f"test010_plot_energy_spectrum_histogram_{interpolation}.png", - ekin, - data_x, - data_y, + bin_edges, + src_data, + sim_data, relerrs, ) utility.test_ok(is_ok) +def plot(output_file, bin_edges, src_data, sim_data, relerrs): + fig, ax = plt.subplots(figsize=(8.5, 6)) + ax.set_xlabel("Energy (MeV)") + ax.set_ylabel("Probability") + ax.stairs(src_data, bin_edges, fill=True, label="input data") + ax.stairs(sim_data, bin_edges, fill=True, alpha=0.6, label="simulation") + + ax.set_xscale("log") + ax.set_yscale("log") + + ax2 = ax.twinx() + ax2.set_ylabel("Relative error") + ax2.set_ylim([-0.20, +0.20]) + ax2.plot( + (bin_edges[:-1] + bin_edges[1:]) / 2, + relerrs, + label="relative error", + color="C2", + ) + + ax.legend() + + fig.savefig(output_file, dpi=300) + plt.close(fig) + + +def root_load_ekin(root_file: str): + data_ref, keys_ref, _ = gatetools.phsp.load(root_file) + + index_ekin = keys_ref.index("KineticEnergy") + ekin = [data_ref_i[index_ekin] for data_ref_i in data_ref] + + return ekin + + +def add_source_energy_spectrum_histogram(sim, interpolation: str = None): + source = sim.add_source("GenericSource", "beam") + source.particle = "e-" + source.n = 2e6 / sim.number_of_threads + source.position.type = "point" + source.direction.type = "iso" + + set_source_energy_spectrum(source, "Lu177") # After defining the particle!! + source.energy.spectrum_histogram_interpolation = interpolation + + return source + + if __name__ == "__main__": paths = utility.get_default_test_paths( __file__, "test010_generic_source_energy_spectrum", output_folder="test010" diff --git a/opengate/tests/src/test013_half_life.py b/opengate/tests/src/test013_half_life.py index 6c677a65a..d30f1984a 100755 --- a/opengate/tests/src/test013_half_life.py +++ b/opengate/tests/src/test013_half_life.py @@ -3,7 +3,7 @@ import opengate as gate from opengate.tests import utility -from opengate.sources.base import get_rad_yield +from opengate.sources.utility import get_rad_yield if __name__ == "__main__": paths = utility.get_default_test_paths(__file__, output_folder="test013_hl") diff --git a/opengate/tests/src/test031_beta_plus_spectrum.py b/opengate/tests/src/test031_beta_plus_spectrum.py index 097b1bd34..5f8b0fcbb 100755 --- a/opengate/tests/src/test031_beta_plus_spectrum.py +++ b/opengate/tests/src/test031_beta_plus_spectrum.py @@ -7,7 +7,7 @@ import opengate as gate from opengate.tests import utility -from opengate.sources.base import compute_bins_density, get_rad_yield +from opengate.sources.utility import compute_bins_density, get_rad_yield, get_spectrum if __name__ == "__main__": paths = utility.get_default_test_paths(__file__, "", "test031") @@ -40,7 +40,7 @@ i += 1 for rad in l: - data = gate.sources.generic.read_beta_plus_spectra(rad) + data = get_spectrum(rad, "e+", "radar") x = data[:, 0] # energy E(keV) y = data[:, 1] # proba dNtot/dE b+ # normalize taking into account the bins density diff --git a/opengate/tests/src/test032_voxel_vs_volume.py b/opengate/tests/src/test032_voxel_vs_volume.py index 730cc81fe..f271a585e 100755 --- a/opengate/tests/src/test032_voxel_vs_volume.py +++ b/opengate/tests/src/test032_voxel_vs_volume.py @@ -1,8 +1,6 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -import json -import itk from scipy.spatial.transform import Rotation import opengate as gate import opengate.contrib.phantoms.nemaiec as gate_iec @@ -52,7 +50,7 @@ iec2.material = "G4_AIR" iec2.translation = [-40 * cm, 0 * cm, 0 * cm] iec2.dump_label_image = paths.output / "test032_iec_label.mhd" - labels = json.loads(open(paths.output_ref / "test032_labels.json").read()) + labels = utility.read_json_file(paths.output_ref / "test032_labels.json") iec2.voxel_materials = [] for l in labels: mat = "IEC_PLASTIC" @@ -76,7 +74,7 @@ iec2.voxel_materials.append(m) pMin, pMax = sim.volume_manager.volumes["iec1"].bounding_limits - print(f"pMin and pMax of iec1", pMin, pMax) + print("pMin and pMax of iec1", pMin, pMax) # the origin of iec1 is different from the origin of iec2 # we create fake images to be able to convert from diff --git a/opengate/tests/src/test037_pet_hits_singles_helpers.py b/opengate/tests/src/test037_pet_hits_singles_helpers.py index c29820c5d..436254e9f 100644 --- a/opengate/tests/src/test037_pet_hits_singles_helpers.py +++ b/opengate/tests/src/test037_pet_hits_singles_helpers.py @@ -6,7 +6,7 @@ import opengate.contrib.phantoms.necr as phantom_necr from opengate.tests import utility from opengate.userhooks import check_production_cuts -from opengate.sources.base import get_rad_yield +from opengate.sources.utility import get_rad_yield def create_pet_simulation(sim, paths, debug=False, create_mat=False): diff --git a/opengate/tests/src/test043_garf_helpers.py b/opengate/tests/src/test043_garf_helpers.py index f859d245a..779eb34ac 100644 --- a/opengate/tests/src/test043_garf_helpers.py +++ b/opengate/tests/src/test043_garf_helpers.py @@ -4,7 +4,7 @@ import opengate as gate import opengate.contrib.spect.ge_discovery_nm670 as gate_spect from opengate.tests import utility -from opengate.sources.base import get_rad_gamma_spectrum +from opengate.sources.utility import get_spectrum paths = utility.get_default_test_paths(__file__, "gate_test043_garf", "test043") @@ -47,7 +47,7 @@ def sim_phys(sim): def sim_source_test(sim, activity): - tc99m = get_rad_gamma_spectrum("Tc99m") + tc99m = get_spectrum("Tc99m", "gamma") e = tc99m.energies w = tc99m.weights diff --git a/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py b/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py index 65ab0a5b5..a37abe983 100644 --- a/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py +++ b/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py @@ -5,7 +5,6 @@ import opengate.contrib.phantoms.nemaiec as gate_iec import gatetools as gt import itk -import json import opengate.contrib.pet.philipsvereos as pet_vereos from opengate.tests import utility from box import Box @@ -105,7 +104,7 @@ def add_voxelized_phantom(sim, param): gate_iec.create_material(sim) iec.image = param.iec_vox_mhd iec.material = "G4_AIR" - labels = json.loads(open(param.iec_vox_json).read()) + labels = utility.read_json_file(param.iec_vox_json) # labels are not material, we assign the material belows # all spheres are water ; central hole is lung ; shell are plastic iec.voxel_materials = [] diff --git a/opengate/tests/src/test046_rad.py b/opengate/tests/src/test046_rad.py index 2d9a19e19..600497bc7 100755 --- a/opengate/tests/src/test046_rad.py +++ b/opengate/tests/src/test046_rad.py @@ -2,11 +2,10 @@ # -*- coding: utf-8 -*- from box import Box -import json import opengate.contrib.spect.ge_discovery_nm670 as gate_spect import opengate as gate from opengate.tests import utility -from opengate.sources.base import get_rad_gamma_spectrum +from opengate.sources.utility import get_spectrum if __name__ == "__main__": paths = utility.get_default_test_paths(__file__, "", output_folder="test046") @@ -53,8 +52,8 @@ json.dump(digit_ns, outfile, indent=4)""" # check - ref_digit = json.loads(open(paths.output_ref / "t046_digitizer.json").read()) - ref_digit_ns = json.loads(open(paths.output_ref / "t046_digitizer_ns.json").read()) + ref_digit = utility.read_json_file(paths.output_ref / "t046_digitizer.json") + ref_digit_ns = utility.read_json_file(paths.output_ref / "t046_digitizer_ns.json") ok = digit == ref_digit utility.print_test(ok, f"Test channels (with scatter): {ok}") is_ok = is_ok and ok @@ -64,13 +63,18 @@ # Test 3 print() - yields = [0.885, 0.172168, 1.847315, 1.0024600000000004] + yields = [0.8907654364665489, 0.18033, 1.8474, 1.0077] + tolerance = 0.01 + st = f"(tol = {tolerance * 100:.2f} %)" i = 0 for rad in radionuclides: - rad_spectrum = get_rad_gamma_spectrum(rad) + rad_spectrum = get_spectrum(rad, "gamma") tw = rad_spectrum["weights"].sum() - ok = tw == yields[i] - utility.print_test(ok, f"Test yield {rad}: {tw} {yields[i]} {ok}") + tw_d = tw / yields[i] * 100 - 100 + ok = abs(tw_d) <= tolerance * 100 + utility.print_test( + ok, f"Test yield {rad}: {tw} {yields[i]}: {tw_d:+.2f} % {st} {ok}" + ) is_ok = is_ok and ok i += 1 diff --git a/opengate/tests/src/test049_pet_digit_blurring_helpers.py b/opengate/tests/src/test049_pet_digit_blurring_helpers.py index 27bed3471..01c4bff13 100644 --- a/opengate/tests/src/test049_pet_digit_blurring_helpers.py +++ b/opengate/tests/src/test049_pet_digit_blurring_helpers.py @@ -12,7 +12,7 @@ ) from opengate.userhooks import check_production_cuts from opengate.tests import utility -from opengate.sources.base import get_rad_yield +from opengate.sources.utility import get_rad_yield paths = utility.get_default_test_paths(__file__, "gate_test049_pet_blur", "test049") diff --git a/opengate/tests/src/test073_helpers.py b/opengate/tests/src/test073_helpers.py index f6de36dc6..315271dbf 100644 --- a/opengate/tests/src/test073_helpers.py +++ b/opengate/tests/src/test073_helpers.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- import opengate as gate -from opengate.sources.base import set_source_rad_energy_spectrum +from opengate.sources.utility import set_source_energy_spectrum from opengate.exception import warning from opengate.tests import utility import opengate.contrib.spect.siemens_intevo as intevo @@ -33,7 +33,7 @@ def create_sim_tests(sim, threads=1, digitizer=1, debug=False): world.material = "G4_AIR" # spect head - head, colli, crystal = intevo.add_spect_head( + head, _, crystal = intevo.add_spect_head( sim, "spect", collimator_type="melp", debug=debug ) head.translation = [0, 0, -280 * mm] @@ -58,7 +58,7 @@ def create_sim_tests(sim, threads=1, digitizer=1, debug=False): s1.position.type = "sphere" s1.position.radius = 30 * mm s1.position.translation = [0, 0, 0] - set_source_rad_energy_spectrum(s1, "Lu177") + set_source_energy_spectrum(s1, "Lu177", "radar") s1.direction.type = "iso" s1.activity = activity @@ -67,7 +67,7 @@ def create_sim_tests(sim, threads=1, digitizer=1, debug=False): s2.position.type = "sphere" s2.position.radius = 60 * mm s2.position.translation = [0, 200 * mm, 0] - set_source_rad_energy_spectrum(s2, "Lu177") + set_source_energy_spectrum(s2, "Lu177", "radar") s2.direction.type = "iso" s2.activity = activity @@ -76,7 +76,7 @@ def create_sim_tests(sim, threads=1, digitizer=1, debug=False): s3.position.type = "sphere" s3.position.radius = 25 * mm s3.position.translation = [100, 0, 0 * mm] - set_source_rad_energy_spectrum(s3, "Lu177") + set_source_energy_spectrum(s3, "Lu177", "radar") s3.direction.type = "iso" s3.activity = activity @@ -204,13 +204,13 @@ def test073_setup_sim(sim, spect_type, collimator_type): # spect ? if spect_type == "intevo": - head, colli, crystal = intevo.add_spect_head( + head, _, _ = intevo.add_spect_head( sim, "spect", collimator_type=collimator_type, debug=sim.visu ) head.translation = [0, 0, -280 * mm] head.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() if spect_type == "discovery": - head, colli, crystal = discovery.add_spect_head( + head, _, _ = discovery.add_spect_head( sim, "spect", collimator_type=collimator_type, debug=sim.visu ) head.translation = [0, 0, -280 * mm] diff --git a/opengate/tests/src/test073_test3_intevo_tc99m_mt.py b/opengate/tests/src/test073_test3_intevo_tc99m_mt.py index d87fbb35c..31a632163 100755 --- a/opengate/tests/src/test073_test3_intevo_tc99m_mt.py +++ b/opengate/tests/src/test073_test3_intevo_tc99m_mt.py @@ -32,7 +32,7 @@ # source Bq = gate.g4_units.Bq - set_source_rad_energy_spectrum(source, "Tc99m") + set_source_energy_spectrum(source, "Tc99m") source.activity = 2e7 * Bq / sim.number_of_threads # start simulation diff --git a/opengate/tests/src/test073_test4_intevo_lu177_mt.py b/opengate/tests/src/test073_test4_intevo_lu177_mt.py index 89c606155..839acbc26 100755 --- a/opengate/tests/src/test073_test4_intevo_lu177_mt.py +++ b/opengate/tests/src/test073_test4_intevo_lu177_mt.py @@ -3,7 +3,7 @@ from test073_helpers import * from opengate.tests import utility -from opengate.sources.base import set_source_rad_energy_spectrum +from opengate.sources.utility import set_source_energy_spectrum if __name__ == "__main__": paths = utility.get_default_test_paths( @@ -33,7 +33,7 @@ # source Bq = gate.g4_units.Bq - set_source_rad_energy_spectrum(source, "Lu177") + set_source_energy_spectrum(source, "Lu177", "radar") source.activity = 2e8 * Bq / sim.number_of_threads # start simulation diff --git a/opengate/tests/src/test073_test5_discovery_lu177_mt.py b/opengate/tests/src/test073_test5_discovery_lu177_mt.py index 1ad43f320..cb4345927 100755 --- a/opengate/tests/src/test073_test5_discovery_lu177_mt.py +++ b/opengate/tests/src/test073_test5_discovery_lu177_mt.py @@ -3,7 +3,7 @@ from test073_helpers import * from opengate.tests import utility -from opengate.sources.base import set_source_rad_energy_spectrum +from opengate.sources.utility import set_source_energy_spectrum if __name__ == "__main__": paths = utility.get_default_test_paths( @@ -29,7 +29,7 @@ # source Bq = gate.g4_units.Bq - set_source_rad_energy_spectrum(source, "Lu177") + set_source_energy_spectrum(source, "Lu177", "radar") source.activity = 3e8 * Bq / sim.number_of_threads # start simulation diff --git a/opengate/tests/src/test073_test5_discovery_tc99m_mt.py b/opengate/tests/src/test073_test5_discovery_tc99m_mt.py index 7af48805c..90f5dbe80 100755 --- a/opengate/tests/src/test073_test5_discovery_tc99m_mt.py +++ b/opengate/tests/src/test073_test5_discovery_tc99m_mt.py @@ -28,7 +28,7 @@ # source Bq = gate.g4_units.Bq - set_source_rad_energy_spectrum(source, "Tc99m") + set_source_energy_spectrum(source, "Tc99m") source.activity = 4e7 * Bq / sim.number_of_threads # start simulation diff --git a/opengate/tests/src/test084_attenuation_map2_hu.py b/opengate/tests/src/test084_attenuation_map2_hu.py index 0e9297d41..470c6bd08 100755 --- a/opengate/tests/src/test084_attenuation_map2_hu.py +++ b/opengate/tests/src/test084_attenuation_map2_hu.py @@ -3,7 +3,7 @@ import opengate as gate from opengate.tests import utility -from opengate.sources.base import get_rad_gamma_spectrum +from opengate.sources.utility import get_spectrum if __name__ == "__main__": paths = utility.get_default_test_paths(__file__, "", output_folder="test084") @@ -47,7 +47,7 @@ mumap = sim.add_actor("AttenuationImageActor", "mumap") mumap.image_volume = patient # FIXME volume for the moment, not the name mumap.output_filename = "mumap2.mhd" - energies = get_rad_gamma_spectrum("Lu177").energies + energies = get_spectrum("Lu177", "gamma", "radar").energies mumap.energy = energies[3] # 0.208 keV mumap.database = "EPDL" # NIST or EPDL print(f"Energy is {mumap.energy/keV} keV") diff --git a/opengate/tests/src/test085_free_flight_helpers.py b/opengate/tests/src/test085_free_flight_helpers.py index 385c995eb..8dbc53765 100755 --- a/opengate/tests/src/test085_free_flight_helpers.py +++ b/opengate/tests/src/test085_free_flight_helpers.py @@ -5,7 +5,7 @@ import opengate.contrib.spect.ge_discovery_nm670 as nm670 import opengate.contrib.phantoms.nemaiec as nemaiec from opengate.image import get_translation_to_isocenter -from opengate.sources.base import set_source_rad_energy_spectrum +from opengate.sources.utility import set_source_energy_spectrum from pathlib import Path import numpy as np import opengate_core as g4 @@ -111,8 +111,8 @@ def create_simulation_test085(sim, paths, ac=1e5, angle_tolerance=None): source = sim.add_source("VoxelSource", "src") source.image = iec_source_filename source.position.translation = [0, 35 * mm, 0] - set_source_rad_energy_spectrum(source, "tc99m") source.particle = "gamma" + set_source_energy_spectrum(source, "tc99m", "radar") # After particle definition source.direction.acceptance_angle.volumes = [h.name for h in det_planes] source.direction.acceptance_angle.skip_policy = "SkipEvents" source.direction.acceptance_angle.intersection_flag = True @@ -125,7 +125,7 @@ def create_simulation_test085(sim, paths, ac=1e5, angle_tolerance=None): # add stat actor stats = sim.add_actor("SimulationStatisticsActor", "stats") - stats.output_filename = f"stats.txt" + stats.output_filename = "stats.txt" # set the gantry orientation nm670.rotate_gantry(det_plane1, radius, 0, 0, 1) diff --git a/opengate/tests/src/test085_free_flight_mt.py b/opengate/tests/src/test085_free_flight_mt.py index e01db9dce..731426ad7 100755 --- a/opengate/tests/src/test085_free_flight_mt.py +++ b/opengate/tests/src/test085_free_flight_mt.py @@ -1,8 +1,9 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- +import opengate as gate from opengate.tests import utility -from test085_free_flight_helpers import * +from test085_free_flight_helpers import create_simulation_test085 if __name__ == "__main__": @@ -16,8 +17,8 @@ arf1 = sim.get_actor("detector_arf_1") arf2 = sim.get_actor("detector_arf_2") - arf1.output_filename = f"projection_ff_1.mhd" - arf2.output_filename = f"projection_ff_2.mhd" + arf1.output_filename = "projection_ff_1.mhd" + arf2.output_filename = "projection_ff_2.mhd" stats = sim.get_actor("stats") stats.output_filename = "stats_ff.txt" diff --git a/opengate/tests/src/test085_free_flight_ref_mt.py b/opengate/tests/src/test085_free_flight_ref_mt.py index 056e4d3c3..ea333210e 100755 --- a/opengate/tests/src/test085_free_flight_ref_mt.py +++ b/opengate/tests/src/test085_free_flight_ref_mt.py @@ -1,7 +1,8 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -from test085_free_flight_helpers import * +import opengate as gate +from test085_free_flight_helpers import create_simulation_test085 from opengate.tests import utility diff --git a/opengate/tests/src/test085_free_flight_rotation.py b/opengate/tests/src/test085_free_flight_rotation.py index 1dc392b9f..9c478894f 100755 --- a/opengate/tests/src/test085_free_flight_rotation.py +++ b/opengate/tests/src/test085_free_flight_rotation.py @@ -1,10 +1,12 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- +import opengate as gate import opengate.logger +import opengate.contrib.spect.ge_discovery_nm670 as nm670 from opengate import g4_units from opengate.tests import utility -from test085_free_flight_helpers import * +from test085_free_flight_helpers import create_simulation_test085 if __name__ == "__main__": @@ -19,8 +21,8 @@ arf1 = sim.get_actor("detector_arf_1") arf2 = sim.get_actor("detector_arf_2") - arf1.output_filename = f"projection_ff_1.mhd" - arf2.output_filename = f"projection_ff_2.mhd" + arf1.output_filename = "projection_ff_1.mhd" + arf2.output_filename = "projection_ff_2.mhd" stats = sim.get_actor("stats") stats.output_filename = "stats_ff.txt" diff --git a/opengate/tests/utility.py b/opengate/tests/utility.py index dd32459a0..2510791c7 100644 --- a/opengate/tests/utility.py +++ b/opengate/tests/utility.py @@ -49,6 +49,22 @@ def test_ok(is_ok=False, exceptions=None): sys.exit(-1) +def read_json_file(filename: Path) -> dict: + """ + Read a JSON file into a Python dictionary. + + :param filename: Path object + The filename of the JSON file to read. + :return: dict + The data from the JSON file. + """ + if not filename.is_file(): + fatal(f"File {filename} does not exist.") + + with open(filename, "rb") as f: + return json.load(f) + + def read_stat_file(filename, encoder=None): if encoder == "json": return read_stat_file_json(filename) diff --git a/opengate/utility.py b/opengate/utility.py index ed12f8c68..5f4510104 100644 --- a/opengate/utility.py +++ b/opengate/utility.py @@ -513,7 +513,7 @@ def get_library_path(): files = os.listdir(path) lib_ext = "dll" if os.name == "nt" else "so" - libs = list(filter(lambda file: file.endswith(f".{lib_ext}"), files)) + libs = [file for file in files if file.endswith(f".{lib_ext}")] if len(libs) == 0: return "unknown" elif len(libs) > 1: @@ -604,3 +604,19 @@ def standard_error_c4_correction(n): return ( np.sqrt(2 / (n - 1)) * sc.special.gamma(n / 2) / sc.special.gamma((n - 1) / 2) ) + + +def read_json_file(filename: Path) -> dict: + """ + Read a JSON file into a Python dictionary. + + :param filename: Path object + The filename of the JSON file to read. + :return: dict + The data from the JSON file. + """ + if not filename.is_file(): + fatal(f"File {filename} does not exist.") + + with open(filename, "rb") as f: + return json.load(f) diff --git a/setup.py b/setup.py index 59b4fdf5b..dd9af365a 100644 --- a/setup.py +++ b/setup.py @@ -32,6 +32,7 @@ "PyYAML", "SimpleITK", "spekpy", + "icrp107-database", ] + install_requires_windows, )