From 359f7a6df6966ec6180ca5f19f69f621760b565a Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Thu, 20 Jul 2023 13:29:09 +0200 Subject: [PATCH 01/46] fixed comment error --- opengate/__init__.py | 1 + opengate/source/TreatmentPlanPhsSource.py | 280 ++++++++++++++++++++++ 2 files changed, 281 insertions(+) create mode 100644 opengate/source/TreatmentPlanPhsSource.py diff --git a/opengate/__init__.py b/opengate/__init__.py index 78efc3abc..1ab2dc82f 100644 --- a/opengate/__init__.py +++ b/opengate/__init__.py @@ -20,6 +20,7 @@ from .actor.ActorBase import * from .geometry.VolumeBase import * from .source.TreatmentPlanSource import * +from .source.TreatmentPlanPhsSource import * # main object from .Simulation import * diff --git a/opengate/source/TreatmentPlanPhsSource.py b/opengate/source/TreatmentPlanPhsSource.py new file mode 100644 index 000000000..7705f5e53 --- /dev/null +++ b/opengate/source/TreatmentPlanPhsSource.py @@ -0,0 +1,280 @@ +import numpy as np +from scipy.spatial.transform import Rotation +import opengate as gate +from .TreatmentPlanSource import * +import os +import pathlib +import uproot +from pathlib import Path + + +class TreatmentPlanPhsSource(TreatmentPlanSource): + def __init__(self, name, sim): + self.name = name + # self.mother = None + self.rotation = Rotation.identity() + self.translation = [0, 0, 0] + self.spots = None + # self.path_to_phaseSpaceFiles = "" + self.phaseSpaceList_file_name = "" + self.phaseSpaceList = {} + self.n_sim = 0 + self.sim = sim # simulation obj to which we want to add the tpPhS source + self.distance_source_to_isocenter = None + # SMX to Isocenter distance + self.distance_stearmag_to_isocenter_x = None + # SMY to Isocenter distance + self.distance_stearmag_to_isocenter_y = None + self.batch_size = None + + def __del__(self): + pass + + def set_particles_to_simulate(self, n_sim): + self.n_sim = n_sim + + def set_distance_source_to_isocenter(self, distance): + self.distance_source_to_isocenter = distance + + def set_distance_stearmag_to_isocenter(self, distance_x, distance_y): + self.distance_stearmag_to_isocenter_x = distance_x + self.distance_stearmag_to_isocenter_y = distance_y + + def set_spots(self, spots): + self.spots = spots + + def set_spots_from_rtplan(self, rt_plan_path): + beamset = gate.beamset_info(rt_plan_path) + gantry_angle = beamset.beam_angles[0] + spots = gate.get_spots_from_beamset(beamset) + self.spots = spots + self.rotation = Rotation.from_euler("z", gantry_angle, degrees=True) + + def set_phaseSpaceList_file_name(self, file_name): + self.phaseSpaceList_file_name = file_name + + def initialize_tpsource(self): + if self.batch_size is None: + self.batch_size = 30000 + # verify that all necessary parameters are set + self.verify_necessary_parameters_are_set() + + # read in the phase space list + self.phaseSpaceList = self.read_list_of_Phs(self.phaseSpaceList_file_name) + # verify the phase space list + self.verify_phs_files_exist(self.phaseSpaceList) + + spots_array = self.spots + sim = self.sim + nSim = self.n_sim + + # mapping factors between iso center plane and nozzle plane (due to steering magnets) + cal_proportion_factor = ( + lambda d_magnet_iso: 1 + if (d_magnet_iso == float("inf")) + else (d_magnet_iso - self.distance_source_to_isocenter) / d_magnet_iso + ) + self.proportion_factor_x = cal_proportion_factor( + self.distance_stearmag_to_isocenter_x + ) + self.proportion_factor_y = cal_proportion_factor( + self.distance_stearmag_to_isocenter_y + ) + tot_sim_particles = 0 + # initialize a pencil beam for each spot + for i, spot in enumerate(spots_array): + # simulate a fraction of the beam particles for this spot + nspot = np.round(spot.beamFraction * nSim) + if nspot == 0: + continue + tot_sim_particles += nspot + # create a source + source = sim.add_source("PhaseSpaceSource", f"{self.name}_spot_{i}") + + # set energy + # find corresponding phase space file + if self.phaseSpaceList.get(spot.energy) is not None: + source.phsp_file = self.phaseSpaceList.get(spot.energy) + + else: + print( + "ERROR in TreatmentPlanPhsSource: Energy requested from plan file does not exist. Aborting." + ) + print("Requested energy was: ", spot.energy) + exit(-1) + + # use the local positions in phase space file + source.position_key = "PrePositionLocal" + source.direction_key = "PreDirectionLocal" + if self.batch_size is not None: + source.batch_size = self.batch_size + else: + source.batch_size = 30000 + + source.particle = spot.particle_name + + # # set mother + # if self.mother is not None: + # source.mother = self.mother + + # POSITION: + source.override_position = True + source.position.translation = self._get_pbs_position(spot) + + # ROTATION: + source.override_direction = True + source.position.rotation = self._get_pbs_rotation(spot) + + # add weight + # source.weight = -1 + source.n = nspot + + self.actual_sim_particles = tot_sim_particles + + def _get_pbs_position(self, spot): + # (x,y) referr to isocenter plane. + # Need to be corrected to referr to nozzle plane + pos = [ + (spot.xiec) * self.proportion_factor_x, + (spot.yiec) * self.proportion_factor_y, + self.distance_source_to_isocenter, + ] + # Gantry angle = 0 -> source comes from +y and is positioned along negative side of y-axis + # https://opengate.readthedocs.io/en/latest/source_and_particle_management.html + + position = (self.rotation * Rotation.from_euler("x", np.pi / 2)).apply( + pos + ) + self.translation + + return position + + def _get_pbs_rotation(self, spot): + # by default the source points in direction z+. + # Need to account for SM direction deviation and rotation thoward isocenter (270 deg around x) + # then rotate of gantry angle + rotation = [0.0, 0.0, 0.0] + beta = np.arctan(spot.yiec / self.distance_stearmag_to_isocenter_y) + alpha = np.arctan(spot.xiec / self.distance_stearmag_to_isocenter_x) + rotation[0] = -np.pi / 2 + beta + rotation[2] = -alpha + + # apply gantry angle + spot_rotation = ( + self.rotation * Rotation.from_euler("xyz", rotation) + ).as_matrix() + + return spot_rotation + + def verify_phs_files_exist(self, phs_dict): + """Check if all the files in the dictionary exist. + Returns True if all the files exist, False otherwise + If one file does not exist, the program is aborted.""" + + for key, phs_file in phs_dict.items(): + # print( + # "key: ", + # key, + # " phs_file: ", + # phs_file, + # " type: ", + # type(phs_file), + # ) + + if not os.path.isfile(phs_file): + print( + "ERROR in ThreatmenPlanPhsSource: File {} does not exist".format( + phs_file + ) + ) + print("Error: File in Phase space dictionary does not exist. Aborting.") + exit(-1) + return True + + def read_list_of_Phs(self, file_name: str, path_to_phsp=""): + """Read a list of Phs from a file. + + Parameters + ---------- + file_name : str + The name of the file to read from. + File needs to have at least two columns + First column is the energy in MeV which is the energy label + Second column is the phase space file name + Returns + ------- + dictionary of Phs""" + + phs_dict = {} + try: + input_arr = np.genfromtxt( + file_name, + delimiter="\t", + comments="#", + usecols=(0, 1), + dtype=None, + encoding=None, + ) + + # convert to dictionary + if input_arr.shape == (0,): + print( + "Error in TreatmentPlanPhsSource: No data found in file: ", + file_name, + " Aborting.", + ) + exit(-1) + if input_arr.ndim == 0: + # only single line read, convert to array + input_arr = np.array([input_arr]) + + # normal case, multiple lines + if path_to_phsp != "": + phs_dict = { + float(i[0]): str(path_to_phsp + "/" + (i[1])) for i in input_arr + } + else: + phs_dict = {float(i[0]): str(i[1]) for i in input_arr} + print("phs_dict read: ", phs_dict) + except Exception as e: + print( + "Error in TreatmentPlanPhsSource: could not read the phase space file list. Aborting." + ) + print("The error was: ", e) + exit(-1) + if len(phs_dict) == 0: + print( + "Error in TreatmentPlanPhsSource: the phase space file list is empty. Aborting." + ) + exit(-1) + return phs_dict + + def verify_necessary_parameters_are_set(self): + """This function checks whether all necessary parameters were set. + E.g. not None + It does not check sensibility""" + if self.phaseSpaceList_file_name is None: + print( + "Error in TreatmentPlanPhsSource: phaseSpaceList_file_name is None. Aborting." + ) + exit(-1) + if self.distance_source_to_isocenter is None: + print( + "Error in TreatmentPlanPhsSource: distance_source_to_isocenter is None. Aborting." + ) + exit(-1) + if self.distance_stearmag_to_isocenter_x is None: + print( + "Error in TreatmentPlanPhsSource: distance_stearmag_to_isocenter_x is None. Aborting." + ) + exit(-1) + if self.distance_stearmag_to_isocenter_y is None: + print( + "Error in TreatmentPlanPhsSource: distance_stearmag_to_isocenter_y is None. Aborting." + ) + exit(-1) + if self.batch_size is None: + print("Error in TreatmentPlanPhsSource: batch_size is None. Aborting.") + exit(-1) + if self.spots is None: + print("Error in TreatmentPlanPhsSource: No spots have been set. Aborting.") + exit(-1) From 0f262353a1eedc53a98c28d668e96a922f1d5c20 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Thu, 20 Jul 2023 16:06:49 +0200 Subject: [PATCH 02/46] Test for TPPhsSource without treatment plan --- opengate/tests/src/test061_TPPhsSource.py | 75 +++++ .../tests/src/test061_TPPhsSource_helpers.py | 311 ++++++++++++++++++ 2 files changed, 386 insertions(+) create mode 100644 opengate/tests/src/test061_TPPhsSource.py create mode 100644 opengate/tests/src/test061_TPPhsSource_helpers.py diff --git a/opengate/tests/src/test061_TPPhsSource.py b/opengate/tests/src/test061_TPPhsSource.py new file mode 100644 index 000000000..75911886d --- /dev/null +++ b/opengate/tests/src/test061_TPPhsSource.py @@ -0,0 +1,75 @@ +import test061_TPPhsSource_helpers as t +import opengate as gate +from scipy.spatial.transform import Rotation + + +paths = gate.get_default_test_paths( + __file__, "test061_TPPhsSource", output_folder="test061" +) + +ref_path = paths.data / "test061_ref" + +# units +m = gate.g4_units("m") +mm = gate.g4_units("mm") +cm = gate.g4_units("cm") +nm = gate.g4_units("nm") +Bq = gate.g4_units("Bq") +MeV = gate.g4_units("MeV") +deg: float = gate.g4_units("deg") + + +def main(): + # print("create reference PhS file") + # t.create_test_Phs( + # particle="proton", + # phs_name=paths.output / "test_proton", + # number_of_particles=1, + # translation=[0 * cm, 0 * cm, 0 * mm], + # ) + print("Testing TPPhS source rotations") + + t.test_source_rotation_A( + plan_file_name=ref_path / "PlanSpot.txt", + phs_list_file_name=ref_path / "PhsList.txt", + phs_file_name_out=paths.output / "output.root", + ) + + a = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="KineticEnergy", + ref_value=150 * MeV, + ) + b = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PrePositionLocal_X", + ref_value=-60, + ) + c = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PrePositionLocal_Y", + ref_value=50, + ) + d = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PrePositionLocal_Z", + ref_value=0, + ) + + e = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PreDirectionLocal_X", + ref_value=-0.012, + ) + f = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PreDirectionLocal_Y", + ref_value=0.01, + ) + + # this is the end, my friend + gate.test_ok(a and b and c and d and e and f) + + +if __name__ == "__main__": + main() diff --git a/opengate/tests/src/test061_TPPhsSource_helpers.py b/opengate/tests/src/test061_TPPhsSource_helpers.py new file mode 100644 index 000000000..d74933a69 --- /dev/null +++ b/opengate/tests/src/test061_TPPhsSource_helpers.py @@ -0,0 +1,311 @@ +import opengate as gate + +# import numpy as np +# import matplotlib.pyplot as plot +from scipy.spatial.transform import Rotation +import gatetools.phsp as phsp +import os + + +# paths = gate.get_default_test_paths( +# __file__, "gate_test019_linac_phsp", output_folder="test019" +# ) +# paths = {output: "output"} + +# units +m = gate.g4_units("m") +mm = gate.g4_units("mm") +cm = gate.g4_units("cm") +nm = gate.g4_units("nm") +Bq = gate.g4_units("Bq") +MeV = gate.g4_units("MeV") +deg: float = gate.g4_units("deg") + + +def create_test_Phs( + particle="proton", + phs_name="output/test_proton.root", + number_of_particles=1, + translation=[0 * mm, 0 * mm, 0 * mm], +): + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + # ui.visu = True + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + # ui.running_verbose_level = gate.EVENT + ui.number_of_threads = 1 + ui.random_seed = "auto" + + # units + m = gate.g4_units("m") + mm = gate.g4_units("mm") + cm = gate.g4_units("cm") + nm = gate.g4_units("nm") + Bq = gate.g4_units("Bq") + MeV = gate.g4_units("MeV") + + ########################################################################################## + # geometry + ########################################################################################## + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + world.material = "G4_Galactic" + + # virtual plane for phase space + plane = sim.add_volume("Tubs", "phase_space_plane") + # plane.material = "G4_AIR" + plane.material = "G4_Galactic" + plane.rmin = 0 + plane.rmax = 30 * cm + plane.dz = 1 * nm # half height + # plane.rotation = Rotation.from_euler("xy", [180, 30], degrees=True).as_matrix() + # plane.translation = [0 * mm, 0 * mm, 0 * mm] + plane.translation = translation + + plane.color = [1, 0, 0, 1] # red + + ########################################################################################## + # Actors + ########################################################################################## + # Split the joined path into the directory path and filename + directory_path, filename = os.path.split(phs_name) + # Extract the base filename and extension + base_filename, extension = os.path.splitext(filename) + new_extension = ".root" + + # PhaseSpace Actor + ta1 = sim.add_actor("PhaseSpaceActor", "PhaseSpace1") + ta1.mother = plane.name + ta1.attributes = [ + "KineticEnergy", + "Weight", + "PrePosition", + "PrePositionLocal", + "ParticleName", + "PreDirection", + "PreDirectionLocal", + "PDGCode", + ] + new_joined_path = os.path.join(directory_path, base_filename + new_extension) + ta1.output = new_joined_path + ta1.debug = False + f = sim.add_filter("ParticleFilter", "f") + f.particle = particle + ta1.filters.append(f) + + phys = sim.get_physics_user_info() + # ~ phys.physics_list_name = "FTFP_BERT" + phys.physics_list_name = "QGSP_BIC_EMZ" + + ########################################################################################## + # Source + ########################################################################################## + source = sim.add_source("GenericSource", "particle_source") + source.mother = "world" + # source.particle = "ion 6 12" # Carbon ions + source.particle = "proton" + source.position.type = "point" + source.direction.type = "momentum" + source.direction.momentum = [0, 0, 1] + source.position.translation = [0 * cm, 0 * cm, -0.1 * cm] + source.energy.type = "mono" + source.energy.mono = 150 * MeV + source.n = number_of_particles + + # output = sim.run() + output = sim.start(start_new_process=True) + + +def create_PhS_withoutSource( + phs_name="output/test_proton.root", +): + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + # ui.visu = True + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + # ui.running_verbose_level = gate.EVENT + ui.number_of_threads = 1 + ui.random_seed = "auto" + + # units + m = gate.g4_units("m") + mm = gate.g4_units("mm") + cm = gate.g4_units("cm") + nm = gate.g4_units("nm") + Bq = gate.g4_units("Bq") + MeV = gate.g4_units("MeV") + deg: float = gate.g4_units("deg") + + ########################################################################################## + # geometry + ########################################################################################## + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + world.material = "G4_Galactic" + print(world.name) + + # virtual plane for phase space + plane = sim.add_volume("Tubs", "phase_space_plane") + # plane.material = "G4_AIR" + plane.material = "G4_Galactic" + plane.rmin = 0 + plane.rmax = 30 * cm + plane.dz = 1 * nm # half height + # plane.rotation = Rotation.from_euler("xy", [180, 30], degrees=True).as_matrix() + plane.translation = [0 * mm, 0 * mm, 0 * mm] + # plane.translation = translation + + plane.color = [1, 0, 0, 1] # red + + ########################################################################################## + # Actors + ########################################################################################## + + # PhaseSpace Actor + ta1 = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + ta1.mother = plane.name + ta1.attributes = [ + "KineticEnergy", + "Weight", + "PrePosition", + "PrePositionLocal", + "ParticleName", + "PreDirection", + "PreDirectionLocal", + "PDGCode", + ] + ta1.output = phs_name + ta1.debug = False + + phys = sim.get_physics_user_info() + # ~ phys.physics_list_name = "FTFP_BERT" + phys.physics_list_name = "QGSP_BIC_EMZ" + + # ########################################################################################## + # # Source + # ########################################################################################## + # # phsp source + # source = sim.add_source("PhaseSpaceSource", "phsp_source_global") + # source.mother = world.name + # source.phsp_file = source_name + # source.position_key = "PrePosition" + # source.direction_key = "PreDirection" + # source.global_flag = True + # source.particle = particle + # source.batch_size = 3000 + # source.n = number_of_particles / ui.number_of_threads + # # source.position.translation = [0 * cm, 0 * cm, -35 * cm] + + # output = sim.run() + # output = sim.start() + # output = sim.start(start_new_process=True) + return sim, plane + + +def test_source_rotation_A( + plan_file_name="output/test_proton_offset.root", + phs_list_file_name="output/PhsList.txt", + phs_file_name_out="output/output/test_source_electron.root", +) -> None: + sim, plane = create_PhS_withoutSource( + phs_name=phs_file_name_out, + ) + number_of_particles = 1 + ########################################################################################## + # Source + ########################################################################################## + # TreatmentPlanPhsSource source + tpPhSs = gate.TreatmentPlanPhsSource("RT_plan", sim) + tpPhSs.set_phaseSpaceList_file_name(phs_list_file_name) + spots, ntot, energies, G = gate.spots_info_from_txt(plan_file_name, "") + # print(spots, ntot, energies, G) + tpPhSs.set_spots(spots) + tpPhSs.set_particles_to_simulate(number_of_particles) + tpPhSs.set_distance_source_to_isocenter(100 * cm) + tpPhSs.set_distance_stearmag_to_isocenter(5 * m, 5 * m) + tpPhSs.rotation = Rotation.from_euler("z", G, degrees=True) + tpPhSs.initialize_tpsource() + + # depending on the rotation of the gantry, the rotation of the phase space to catch the particles is different + plane.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() + + output = sim.start() + + +def test_source_rotation_B( + plan_file_name="output/test_proton_offset.root", + phs_list_file_name="output/PhsList.txt", + phs_file_name_out="output/output/test_source_electron.root", +) -> None: + sim, plane = create_PhS_withoutSource( + phs_name=phs_file_name_out, + ) + number_of_particles = 1 + ########################################################################################## + # Source + ########################################################################################## + # TreatmentPlanPhsSource source + tpPhSs = gate.TreatmentPlanPhsSource("RT_plan", sim) + tpPhSs.set_phaseSpaceList_file_name(phs_list_file_name) + spots, ntot, energies, G = gate.spots_info_from_txt(plan_file_name, "") + # print(spots, ntot, energies, G) + tpPhSs.set_spots(spots) + tpPhSs.set_particles_to_simulate(number_of_particles) + tpPhSs.set_distance_source_to_isocenter(100 * cm) + tpPhSs.set_distance_stearmag_to_isocenter(5 * m, 5 * m) + # tpPhSs.rotation = Rotation.from_euler("z", G, degrees=True) + tpPhSs.initialize_tpsource() + + # depending on the rotation of the gantry, the rotation of the phase space to catch the particles is different + # plane.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() + + output = sim.start() + + +def get_first_entry_of_key( + file_name_root="output/test_source_electron.root", key="ParticleName" +) -> None: + # read root file + data_ref, keys_ref, m_ref = phsp.load(file_name_root) + # print(data_ref) + # print(keys_ref) + index = keys_ref.index(key) + # print(index) + # print(data_ref[index][0]) + return data_ref[0][index] + + +def check_value_from_root_file( + file_name_root="output/test_source_electron.root", + key="ParticleName", + ref_value="e-", +): + # read root file + value = get_first_entry_of_key(file_name_root=file_name_root, key=key) + if (type(ref_value) != str) and (type(value) != str): + is_ok = gate.check_diff_abs( + float(value), float(ref_value), tolerance=1e-3, txt=key + ) + # gate.check_diff_abs(float(value), float(ref_value), tolerance=1e-6, txt=key) + else: + if value == ref_value: + # print("Is correct") + is_ok = True + else: + # print("Is not correct") + is_ok = False + # print("ref_value: ", ref_value) + # print("value: ", value) + return is_ok From 3a22c23b8fcaf50920ab2f340c1ba36bc2ee84ce Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Fri, 21 Jul 2023 09:37:57 +0200 Subject: [PATCH 03/46] removal of duplicated code from inherited function TreatmenPlanSource.py --- opengate/source/TreatmentPlanPhsSource.py | 74 +++---------------- .../tests/src/test061_TPPhsSource_helpers.py | 4 +- 2 files changed, 14 insertions(+), 64 deletions(-) diff --git a/opengate/source/TreatmentPlanPhsSource.py b/opengate/source/TreatmentPlanPhsSource.py index 7705f5e53..3d821327a 100644 --- a/opengate/source/TreatmentPlanPhsSource.py +++ b/opengate/source/TreatmentPlanPhsSource.py @@ -3,19 +3,18 @@ import opengate as gate from .TreatmentPlanSource import * import os -import pathlib -import uproot -from pathlib import Path + +# import pathlib +# import uproot +# from pathlib import Path class TreatmentPlanPhsSource(TreatmentPlanSource): def __init__(self, name, sim): self.name = name - # self.mother = None self.rotation = Rotation.identity() self.translation = [0, 0, 0] self.spots = None - # self.path_to_phaseSpaceFiles = "" self.phaseSpaceList_file_name = "" self.phaseSpaceList = {} self.n_sim = 0 @@ -30,40 +29,28 @@ def __init__(self, name, sim): def __del__(self): pass - def set_particles_to_simulate(self, n_sim): - self.n_sim = n_sim - def set_distance_source_to_isocenter(self, distance): self.distance_source_to_isocenter = distance + self.d_nozzle_to_iso = distance def set_distance_stearmag_to_isocenter(self, distance_x, distance_y): + self.d_stearMag_to_iso_x = distance_x + self.d_stearMag_to_iso_y = distance_y self.distance_stearmag_to_isocenter_x = distance_x self.distance_stearmag_to_isocenter_y = distance_y - def set_spots(self, spots): - self.spots = spots - - def set_spots_from_rtplan(self, rt_plan_path): - beamset = gate.beamset_info(rt_plan_path) - gantry_angle = beamset.beam_angles[0] - spots = gate.get_spots_from_beamset(beamset) - self.spots = spots - self.rotation = Rotation.from_euler("z", gantry_angle, degrees=True) - def set_phaseSpaceList_file_name(self, file_name): self.phaseSpaceList_file_name = file_name - def initialize_tpsource(self): + def initialize_tpPhssource(self): if self.batch_size is None: self.batch_size = 30000 # verify that all necessary parameters are set self.verify_necessary_parameters_are_set() - # read in the phase space list self.phaseSpaceList = self.read_list_of_Phs(self.phaseSpaceList_file_name) # verify the phase space list self.verify_phs_files_exist(self.phaseSpaceList) - spots_array = self.spots sim = self.sim nSim = self.n_sim @@ -72,14 +59,11 @@ def initialize_tpsource(self): cal_proportion_factor = ( lambda d_magnet_iso: 1 if (d_magnet_iso == float("inf")) - else (d_magnet_iso - self.distance_source_to_isocenter) / d_magnet_iso - ) - self.proportion_factor_x = cal_proportion_factor( - self.distance_stearmag_to_isocenter_x - ) - self.proportion_factor_y = cal_proportion_factor( - self.distance_stearmag_to_isocenter_y + else (d_magnet_iso - self.d_nozzle_to_iso) / d_magnet_iso ) + self.proportion_factor_x = cal_proportion_factor(self.d_stearMag_to_iso_x) + self.proportion_factor_y = cal_proportion_factor(self.d_stearMag_to_iso_y) + tot_sim_particles = 0 # initialize a pencil beam for each spot for i, spot in enumerate(spots_array): @@ -131,40 +115,6 @@ def initialize_tpsource(self): self.actual_sim_particles = tot_sim_particles - def _get_pbs_position(self, spot): - # (x,y) referr to isocenter plane. - # Need to be corrected to referr to nozzle plane - pos = [ - (spot.xiec) * self.proportion_factor_x, - (spot.yiec) * self.proportion_factor_y, - self.distance_source_to_isocenter, - ] - # Gantry angle = 0 -> source comes from +y and is positioned along negative side of y-axis - # https://opengate.readthedocs.io/en/latest/source_and_particle_management.html - - position = (self.rotation * Rotation.from_euler("x", np.pi / 2)).apply( - pos - ) + self.translation - - return position - - def _get_pbs_rotation(self, spot): - # by default the source points in direction z+. - # Need to account for SM direction deviation and rotation thoward isocenter (270 deg around x) - # then rotate of gantry angle - rotation = [0.0, 0.0, 0.0] - beta = np.arctan(spot.yiec / self.distance_stearmag_to_isocenter_y) - alpha = np.arctan(spot.xiec / self.distance_stearmag_to_isocenter_x) - rotation[0] = -np.pi / 2 + beta - rotation[2] = -alpha - - # apply gantry angle - spot_rotation = ( - self.rotation * Rotation.from_euler("xyz", rotation) - ).as_matrix() - - return spot_rotation - def verify_phs_files_exist(self, phs_dict): """Check if all the files in the dictionary exist. Returns True if all the files exist, False otherwise diff --git a/opengate/tests/src/test061_TPPhsSource_helpers.py b/opengate/tests/src/test061_TPPhsSource_helpers.py index d74933a69..992688092 100644 --- a/opengate/tests/src/test061_TPPhsSource_helpers.py +++ b/opengate/tests/src/test061_TPPhsSource_helpers.py @@ -236,7 +236,7 @@ def test_source_rotation_A( tpPhSs.set_distance_source_to_isocenter(100 * cm) tpPhSs.set_distance_stearmag_to_isocenter(5 * m, 5 * m) tpPhSs.rotation = Rotation.from_euler("z", G, degrees=True) - tpPhSs.initialize_tpsource() + tpPhSs.initialize_tpPhssource() # depending on the rotation of the gantry, the rotation of the phase space to catch the particles is different plane.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() @@ -266,7 +266,7 @@ def test_source_rotation_B( tpPhSs.set_distance_source_to_isocenter(100 * cm) tpPhSs.set_distance_stearmag_to_isocenter(5 * m, 5 * m) # tpPhSs.rotation = Rotation.from_euler("z", G, degrees=True) - tpPhSs.initialize_tpsource() + tpPhSs.initialize_tpPhssource() # depending on the rotation of the gantry, the rotation of the phase space to catch the particles is different # plane.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() From 91d3d52d0f7c785c1626de6b4a22ba04758efd30 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Fri, 21 Jul 2023 09:39:05 +0200 Subject: [PATCH 04/46] small code changes --- opengate/source/TreatmentPlanPhsSource.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/opengate/source/TreatmentPlanPhsSource.py b/opengate/source/TreatmentPlanPhsSource.py index 3d821327a..78e048389 100644 --- a/opengate/source/TreatmentPlanPhsSource.py +++ b/opengate/source/TreatmentPlanPhsSource.py @@ -4,10 +4,6 @@ from .TreatmentPlanSource import * import os -# import pathlib -# import uproot -# from pathlib import Path - class TreatmentPlanPhsSource(TreatmentPlanSource): def __init__(self, name, sim): From 6cdd52c4c8cc3ac77e596b8cb56cd83789073816 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Fri, 21 Jul 2023 10:15:03 +0200 Subject: [PATCH 05/46] Add data for test061 --- opengate/tests/data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opengate/tests/data b/opengate/tests/data index d09897ab4..bdac91904 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit d09897ab4445299141011a772deb0d08e7014958 +Subproject commit bdac91904851981cdc12826619df2d6de09b0a26 From 812ee1284e1fdb2e5e581fe606c581f61ce0a3cf Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Fri, 21 Jul 2023 11:21:50 +0200 Subject: [PATCH 06/46] Change output ref folder path --- opengate/tests/src/test061_TPPhsSource.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opengate/tests/src/test061_TPPhsSource.py b/opengate/tests/src/test061_TPPhsSource.py index 75911886d..1ccf1161a 100644 --- a/opengate/tests/src/test061_TPPhsSource.py +++ b/opengate/tests/src/test061_TPPhsSource.py @@ -6,8 +6,9 @@ paths = gate.get_default_test_paths( __file__, "test061_TPPhsSource", output_folder="test061" ) +paths.output_ref = paths.output_ref / "test061_ref" -ref_path = paths.data / "test061_ref" +ref_path = paths.output_ref # units m = gate.g4_units("m") From 3985374a5b542f5b5184b4b1b6dfaa8255966f89 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Fri, 21 Jul 2023 13:13:39 +0200 Subject: [PATCH 07/46] Correct path for test061 --- opengate/tests/data | 2 +- opengate/tests/src/test061_TPPhsSource.py | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/opengate/tests/data b/opengate/tests/data index bdac91904..6d32057de 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit bdac91904851981cdc12826619df2d6de09b0a26 +Subproject commit 6d32057dea5640cb4a55113ecc8d175c57eefbc1 diff --git a/opengate/tests/src/test061_TPPhsSource.py b/opengate/tests/src/test061_TPPhsSource.py index 1ccf1161a..548f9dfe1 100644 --- a/opengate/tests/src/test061_TPPhsSource.py +++ b/opengate/tests/src/test061_TPPhsSource.py @@ -3,9 +3,7 @@ from scipy.spatial.transform import Rotation -paths = gate.get_default_test_paths( - __file__, "test061_TPPhsSource", output_folder="test061" -) +paths = gate.get_default_test_paths(__file__, "test061_TPPhsSource") paths.output_ref = paths.output_ref / "test061_ref" ref_path = paths.output_ref From debc8a07aa9dcab3b766296002c5a770df01b48a Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Tue, 25 Jul 2023 16:02:26 +0200 Subject: [PATCH 08/46] added PhS path folder to TPPhs source --- opengate/source/TreatmentPlanPhsSource.py | 13 +++++-- opengate/tests/src/test061_TPPhsSource.py | 3 +- .../tests/src/test061_TPPhsSource_helpers.py | 34 ++----------------- 3 files changed, 16 insertions(+), 34 deletions(-) diff --git a/opengate/source/TreatmentPlanPhsSource.py b/opengate/source/TreatmentPlanPhsSource.py index 78e048389..d7e73a479 100644 --- a/opengate/source/TreatmentPlanPhsSource.py +++ b/opengate/source/TreatmentPlanPhsSource.py @@ -11,6 +11,7 @@ def __init__(self, name, sim): self.rotation = Rotation.identity() self.translation = [0, 0, 0] self.spots = None + self.phaseSpaceFolder = "" self.phaseSpaceList_file_name = "" self.phaseSpaceList = {} self.n_sim = 0 @@ -36,7 +37,10 @@ def set_distance_stearmag_to_isocenter(self, distance_x, distance_y): self.distance_stearmag_to_isocenter_y = distance_y def set_phaseSpaceList_file_name(self, file_name): - self.phaseSpaceList_file_name = file_name + self.phaseSpaceList_file_name = str(file_name) + + def set_phaseSpaceFolder(self, folder_name): + self.phaseSpaceFolder = str(folder_name) def initialize_tpPhssource(self): if self.batch_size is None: @@ -44,7 +48,9 @@ def initialize_tpPhssource(self): # verify that all necessary parameters are set self.verify_necessary_parameters_are_set() # read in the phase space list - self.phaseSpaceList = self.read_list_of_Phs(self.phaseSpaceList_file_name) + self.phaseSpaceList = self.read_list_of_Phs( + self.phaseSpaceList_file_name, self.phaseSpaceFolder + ) # verify the phase space list self.verify_phs_files_exist(self.phaseSpaceList) spots_array = self.spots @@ -150,6 +156,9 @@ def read_list_of_Phs(self, file_name: str, path_to_phsp=""): ------- dictionary of Phs""" + # prepend the path to the phase space files + file_name = path_to_phsp + "/" + file_name + phs_dict = {} try: input_arr = np.genfromtxt( diff --git a/opengate/tests/src/test061_TPPhsSource.py b/opengate/tests/src/test061_TPPhsSource.py index 75911886d..44b4b993d 100644 --- a/opengate/tests/src/test061_TPPhsSource.py +++ b/opengate/tests/src/test061_TPPhsSource.py @@ -31,7 +31,8 @@ def main(): t.test_source_rotation_A( plan_file_name=ref_path / "PlanSpot.txt", - phs_list_file_name=ref_path / "PhsList.txt", + phs_list_file_name="PhsList.txt", + phs_folder_name=ref_path, phs_file_name_out=paths.output / "output.root", ) diff --git a/opengate/tests/src/test061_TPPhsSource_helpers.py b/opengate/tests/src/test061_TPPhsSource_helpers.py index 992688092..c6866d0ce 100644 --- a/opengate/tests/src/test061_TPPhsSource_helpers.py +++ b/opengate/tests/src/test061_TPPhsSource_helpers.py @@ -216,7 +216,8 @@ def create_PhS_withoutSource( def test_source_rotation_A( plan_file_name="output/test_proton_offset.root", - phs_list_file_name="output/PhsList.txt", + phs_list_file_name="PhsList.txt", + phs_folder_name="", phs_file_name_out="output/output/test_source_electron.root", ) -> None: sim, plane = create_PhS_withoutSource( @@ -229,6 +230,7 @@ def test_source_rotation_A( # TreatmentPlanPhsSource source tpPhSs = gate.TreatmentPlanPhsSource("RT_plan", sim) tpPhSs.set_phaseSpaceList_file_name(phs_list_file_name) + tpPhSs.set_phaseSpaceFolder(phs_folder_name) spots, ntot, energies, G = gate.spots_info_from_txt(plan_file_name, "") # print(spots, ntot, energies, G) tpPhSs.set_spots(spots) @@ -244,36 +246,6 @@ def test_source_rotation_A( output = sim.start() -def test_source_rotation_B( - plan_file_name="output/test_proton_offset.root", - phs_list_file_name="output/PhsList.txt", - phs_file_name_out="output/output/test_source_electron.root", -) -> None: - sim, plane = create_PhS_withoutSource( - phs_name=phs_file_name_out, - ) - number_of_particles = 1 - ########################################################################################## - # Source - ########################################################################################## - # TreatmentPlanPhsSource source - tpPhSs = gate.TreatmentPlanPhsSource("RT_plan", sim) - tpPhSs.set_phaseSpaceList_file_name(phs_list_file_name) - spots, ntot, energies, G = gate.spots_info_from_txt(plan_file_name, "") - # print(spots, ntot, energies, G) - tpPhSs.set_spots(spots) - tpPhSs.set_particles_to_simulate(number_of_particles) - tpPhSs.set_distance_source_to_isocenter(100 * cm) - tpPhSs.set_distance_stearmag_to_isocenter(5 * m, 5 * m) - # tpPhSs.rotation = Rotation.from_euler("z", G, degrees=True) - tpPhSs.initialize_tpPhssource() - - # depending on the rotation of the gantry, the rotation of the phase space to catch the particles is different - # plane.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() - - output = sim.start() - - def get_first_entry_of_key( file_name_root="output/test_source_electron.root", key="ParticleName" ) -> None: From 5398250620ca9fbc1c6cc43719ff2e5c77912159 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Tue, 1 Aug 2023 12:12:51 +0200 Subject: [PATCH 09/46] added function to allow the selection of spots of one beam of a dicom beamset --- opengate/helpers_rt_plan.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/opengate/helpers_rt_plan.py b/opengate/helpers_rt_plan.py index 0d8ba9492..a90346b95 100644 --- a/opengate/helpers_rt_plan.py +++ b/opengate/helpers_rt_plan.py @@ -695,4 +695,27 @@ def get_spots_from_beamset(beamset): return spots_array +def get_spots_from_beam(beamset, beam_name): + """Get the spots from a beamset with a given name + the beamFraction is weighted by the total number of particles in the beam""" + rad_type = beamset.bs_info["Radiation Type Opengate"] + spots_array = [] + # mswtot = beamset.mswtot + # find the beam name and get its index + beam_index = beamset.beam_names.index(beam_name) + print(f"Beam index: {beam_index}") + beam = beamset.beams[beam_index] + mswtot = beam.FinalCumulativeMetersetWeight + # mswtot = beam.mswtot + for energy_layer in beam.layers: + for spot in energy_layer.spots: + nPlannedSpot = spot.w + spot.beamFraction = ( + nPlannedSpot / mswtot + ) # nr particles planned for the spot/tot particles planned for the beam + spot.particle_name = rad_type + spots_array.append(spot) + return spots_array + + # vim: set et softtabstop=4 sw=4 smartindent: From b810dee64423f3cac6d24a601342287706b62ddc Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Wed, 2 Aug 2023 08:18:13 +0200 Subject: [PATCH 10/46] made TpsPhsSource compatible with other naming conventions in phase space files, e.g. old root files from Gate v9 --- opengate/source/PhaseSpaceSource.py | 26 +++++++++--------- opengate/source/TreatmentPlanPhsSource.py | 33 ++++++++++++++--------- 2 files changed, 35 insertions(+), 24 deletions(-) diff --git a/opengate/source/PhaseSpaceSource.py b/opengate/source/PhaseSpaceSource.py index e877d9cbf..ca99115bb 100644 --- a/opengate/source/PhaseSpaceSource.py +++ b/opengate/source/PhaseSpaceSource.py @@ -41,13 +41,13 @@ def set_default_user_info(user_info): user_info.batch_size = 10000 # branch name in the phsp file user_info.position_key = "PrePositionLocal" - user_info.position_key_x = None - user_info.position_key_y = None - user_info.position_key_z = None + user_info.position_key_x = "" + user_info.position_key_y = "" + user_info.position_key_z = "" user_info.direction_key = "PreDirectionLocal" - user_info.direction_key_x = None - user_info.direction_key_y = None - user_info.direction_key_z = None + user_info.direction_key_x = "" + user_info.direction_key_y = "" + user_info.direction_key_z = "" user_info.energy_key = "KineticEnergy" user_info.weight_key = "Weight" user_info.particle_name_key = "ParticleName" @@ -83,18 +83,20 @@ def initialize(self, run_timing_intervals): ) # check user info + # only if position_key_x etc is not set, then set then based in position_key ui = self.user_info - if ui.position_key_x is None: + if ui.position_key_x == "": ui.position_key_x = f"{ui.position_key}_X" - if ui.position_key_y is None: + if ui.position_key_y == "": ui.position_key_y = f"{ui.position_key}_Y" - if ui.position_key_z is None: + if ui.position_key_z == "": ui.position_key_z = f"{ui.position_key}_Z" - if ui.direction_key_x is None: + # only if direction_key_x etc is not set, then set then based in direction_key + if ui.direction_key_x == "": ui.direction_key_x = f"{ui.direction_key}_X" - if ui.direction_key_y is None: + if ui.direction_key_y == "": ui.direction_key_y = f"{ui.direction_key}_Y" - if ui.direction_key_z is None: + if ui.direction_key_z == "": ui.direction_key_z = f"{ui.direction_key}_Z" # initialize the generator (read the phsp file) diff --git a/opengate/source/TreatmentPlanPhsSource.py b/opengate/source/TreatmentPlanPhsSource.py index d7e73a479..8a88ba019 100644 --- a/opengate/source/TreatmentPlanPhsSource.py +++ b/opengate/source/TreatmentPlanPhsSource.py @@ -14,6 +14,15 @@ def __init__(self, name, sim): self.phaseSpaceFolder = "" self.phaseSpaceList_file_name = "" self.phaseSpaceList = {} + self.position_key_x = "PrePositionLocal_X" + self.position_key_y = "PrePositionLocal_Y" + self.position_key_z = "PrePositionLocal_Z" + self.direction_key_x = "PreDirectionLocal_X" + self.direction_key_y = "PreDirectionLocal_Y" + self.direction_key_z = "PreDirectionLocal_Z" + self.energy_key = "KineticEnergy" + self.weight_key = "Weight" + self.PDGCode_key = "PDGCode" self.n_sim = 0 self.sim = sim # simulation obj to which we want to add the tpPhS source self.distance_source_to_isocenter = None @@ -81,7 +90,6 @@ def initialize_tpPhssource(self): # find corresponding phase space file if self.phaseSpaceList.get(spot.energy) is not None: source.phsp_file = self.phaseSpaceList.get(spot.energy) - else: print( "ERROR in TreatmentPlanPhsSource: Energy requested from plan file does not exist. Aborting." @@ -89,20 +97,22 @@ def initialize_tpPhssource(self): print("Requested energy was: ", spot.energy) exit(-1) - # use the local positions in phase space file - source.position_key = "PrePositionLocal" - source.direction_key = "PreDirectionLocal" + # set keys of phase space file to use + source.position_key_x = self.position_key_x + source.position_key_y = self.position_key_y + source.position_key_z = self.position_key_z + source.direction_key_x = self.direction_key_x + source.direction_key_y = self.direction_key_y + source.direction_key_z = self.direction_key_z + source.energy_key = self.energy_key + source.weight_key = self.weight_key + source.PDGCode_key = self.PDGCode_key + if self.batch_size is not None: source.batch_size = self.batch_size else: source.batch_size = 30000 - source.particle = spot.particle_name - - # # set mother - # if self.mother is not None: - # source.mother = self.mother - # POSITION: source.override_position = True source.position.translation = self._get_pbs_position(spot) @@ -112,7 +122,6 @@ def initialize_tpPhssource(self): source.position.rotation = self._get_pbs_rotation(spot) # add weight - # source.weight = -1 source.n = nspot self.actual_sim_particles = tot_sim_particles @@ -189,7 +198,7 @@ def read_list_of_Phs(self, file_name: str, path_to_phsp=""): } else: phs_dict = {float(i[0]): str(i[1]) for i in input_arr} - print("phs_dict read: ", phs_dict) + # print("phs_dict read: ", phs_dict) except Exception as e: print( "Error in TreatmentPlanPhsSource: could not read the phase space file list. Aborting." From 843866eafcac9e9959e567720dc63b4d602ffdca Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Fri, 4 Aug 2023 12:04:25 +0200 Subject: [PATCH 11/46] extended phase space source and TP phase space source to count particles by a given primary and energy threshold, instead of counting every particle --- .../opengate_lib/GatePhaseSpaceSource.cpp | 98 ++++++++++++++---- .../opengate_lib/GatePhaseSpaceSource.h | 7 ++ .../opengate_lib/pyGatePhaseSpaceSource.cpp | 1 + opengate/source/PhaseSpaceSource.py | 21 ++++ opengate/source/TreatmentPlanPhsSource.py | 10 ++ .../src/test060_PhsSource_UntilNexPrimary.py | 99 +++++++++++++++++++ .../tests/src/test060_PhsSource_helpers.py | 29 ++++++ 7 files changed, 245 insertions(+), 20 deletions(-) create mode 100644 opengate/tests/src/test060_PhsSource_UntilNexPrimary.py diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp index 32761eab6..223464588 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp @@ -49,6 +49,13 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) { fMass = fParticleDefinition->GetPDGMass(); } + fgenerate_until_next_primary = + DictGetBool(user_info, "generate_until_next_primary"); + fprimary_lower_energy_threshold = + DictGetDouble(user_info, "primary_lower_energy_threshold"); + fprimary_PDGCode = DictGetInt(user_info, "primary_PDGCode"); + fprimary_pname = DictGetStr(user_info, "primary_particle_name"); + // Init fNumberOfGeneratedEvents = 0; fCurrentIndex = -1; @@ -90,19 +97,51 @@ void GatePhaseSpaceSource::GenerateBatchOfParticles() { void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event, double current_simulation_time) { + // check if we should simulate until next primary + // in this case, generate until a second primary is in the list, excluding the + // second primary + if (fgenerate_until_next_primary) { + int num_primaries = 0; + while (num_primaries <= 2) { + // If batch is empty, we generate some particles + if (fCurrentIndex >= fCurrentBatchSize) + GenerateBatchOfParticles(); + + // check if next particle is primary + if (ParticleIsPrimary()) + num_primaries++; + + // don't generate the second primary + if (num_primaries < 2) { + // Go + GenerateOnePrimary(event, current_simulation_time); + + // update the index; + fCurrentIndex++; + + // // update the number of generated event + // fNumberOfGeneratedEvents++; + } else + break; + } + // update the number of generated event + fNumberOfGeneratedEvents++; + } - // If batch is empty, we generate some particles - if (fCurrentIndex >= fCurrentBatchSize) - GenerateBatchOfParticles(); + else { + // If batch is empty, we generate some particles + if (fCurrentIndex >= fCurrentBatchSize) + GenerateBatchOfParticles(); - // Go - GenerateOnePrimary(event, current_simulation_time); + // Go + GenerateOnePrimary(event, current_simulation_time); - // update the index; - fCurrentIndex++; + // update the index; + fCurrentIndex++; - // update the number of generated event - fNumberOfGeneratedEvents++; + // update the number of generated event + fNumberOfGeneratedEvents++; + } } void GatePhaseSpaceSource::GenerateOnePrimary(G4Event *event, @@ -139,17 +178,6 @@ void GatePhaseSpaceSource::AddOnePrimaryVertex(G4Event *event, const G4ThreeVector &direction, double energy, double time, double w) { - // create primary particle - // SetPDGcode(G4int c); - // G4PrimaryParticle(const G4ParticleDefinition *Gcode, G4double px, G4double - // py, G4double pz); else - // { - // fUseParticleTypeFromFile = false; - // auto *particle_table = G4ParticleTable::GetParticleTable(); - // fParticleDefinition = particle_table->FindParticle(pname); - // fCharge = fParticleDefinition->GetPDGCharge(); - // fMass = fParticleDefinition->GetPDGMass(); - // } auto *particle = new G4PrimaryParticle(); // if we use the particle type from the file, there are two options @@ -195,3 +223,33 @@ void GatePhaseSpaceSource::AddOnePrimaryVertex(G4Event *event, // weights event->GetPrimaryVertex(0)->SetWeight(w); } + +bool GatePhaseSpaceSource::ParticleIsPrimary() { + // check if particle is primary + bool is_primary = false; + + // if PDGCode exists in file + if ((fPDGCode[fCurrentIndex] != 0) && (fprimary_PDGCode != 0)) { + if ((fprimary_PDGCode == fPDGCode[fCurrentIndex]) && + (fprimary_lower_energy_threshold <= fEnergy[fCurrentIndex])) { + is_primary = true; + } + } + // if PDGCode does not exist in file, but particle name does and is not + // empty + else if ((fParticleName[fCurrentIndex].length() != 0) && + (fprimary_pname.length() != 0)) { + if ((fprimary_pname == fParticleName[fCurrentIndex]) && + (fprimary_lower_energy_threshold <= fEnergy[fCurrentIndex])) { + is_primary = true; + } + } else { + G4Exception("GatePhaseSpaceSource::ParticleIsPrimary", "Error", + FatalException, "Particle type not defined in file"); + std::cout << "ERROR: Particle name nor PDGCode defined in file. Aborting." + << std::endl; + exit(1); + } + + return is_primary; +} diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h index a3d0c6555..68c7d5081 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h @@ -47,6 +47,8 @@ class GatePhaseSpaceSource : public GateVSource { void SetGeneratorFunction(ParticleGeneratorType &f); + bool ParticleIsPrimary(); + // virtual void SetGeneratorInfo(py::dict &user_info); void GenerateBatchOfParticles(); @@ -77,6 +79,11 @@ class GatePhaseSpaceSource : public GateVSource { std::vector fWeight; // std::vector fTime; + bool fgenerate_until_next_primary; + string fprimary_pname; + int fprimary_PDGCode; + double fprimary_lower_energy_threshold; + ParticleGeneratorType fGenerator; size_t fCurrentIndex; }; diff --git a/core/opengate_core/opengate_lib/pyGatePhaseSpaceSource.cpp b/core/opengate_core/opengate_lib/pyGatePhaseSpaceSource.cpp index 8382a62c8..8a8e062b5 100644 --- a/core/opengate_core/opengate_lib/pyGatePhaseSpaceSource.cpp +++ b/core/opengate_core/opengate_lib/pyGatePhaseSpaceSource.cpp @@ -34,6 +34,7 @@ void init_GatePhaseSpaceSource(py::module &m) { .def_readwrite("fEnergy", &GatePhaseSpaceSource::fEnergy) .def_readwrite("fWeight", &GatePhaseSpaceSource::fWeight) + // .def_readwrite("fTime", &GatePhaseSpaceSource::fTime) ; } diff --git a/opengate/source/PhaseSpaceSource.py b/opengate/source/PhaseSpaceSource.py index ca99115bb..9f307d6a8 100644 --- a/opengate/source/PhaseSpaceSource.py +++ b/opengate/source/PhaseSpaceSource.py @@ -60,6 +60,10 @@ def set_default_user_info(user_info): user_info.position = Box() user_info.position.translation = [0, 0, 0] user_info.position.rotation = Rotation.identity().as_matrix() + user_info.generate_until_next_primary = False + user_info.primary_particle_name = "" + user_info.primary_lower_energy_threshold = 0 + user_info.primary_PDGCode = 0 # user_info.time_key = None # FIXME later def __del__(self): @@ -99,6 +103,23 @@ def initialize(self, run_timing_intervals): if ui.direction_key_z == "": ui.direction_key_z = f"{ui.direction_key}_Z" + # check if the source should generate particles until the second one + # which is identified as primary by name, PDGCode and above a threshold + if ui.generate_until_next_primary == True: + if len(ui.primary_particle_name) <= 0 and ui.primary_PDGCode == 0: + gate.fatal( + f"PhaseSpaceSource: generate_until_next_primary is True but no primary particle is defined" + ) + if ui.primary_lower_energy_threshold <= 0: + gate.fatal( + f"PhaseSpaceSource: generate_until_next_primary is True but no primary_lower_energy_threshold is defined" + ) + + print( + "PhaseSpaceSource primary_lower_energy_threshold: ", + ui.primary_lower_energy_threshold, + ) + # initialize the generator (read the phsp file) self.particle_generator.initialize(self.user_info) diff --git a/opengate/source/TreatmentPlanPhsSource.py b/opengate/source/TreatmentPlanPhsSource.py index 8a88ba019..7a796b5a1 100644 --- a/opengate/source/TreatmentPlanPhsSource.py +++ b/opengate/source/TreatmentPlanPhsSource.py @@ -23,6 +23,10 @@ def __init__(self, name, sim): self.energy_key = "KineticEnergy" self.weight_key = "Weight" self.PDGCode_key = "PDGCode" + self.generate_until_next_primary = False + self.primary_particle_name = "" + self.primary_lower_energy_threshold = 0 + self.primary_PDGCode = 0 self.n_sim = 0 self.sim = sim # simulation obj to which we want to add the tpPhS source self.distance_source_to_isocenter = None @@ -124,6 +128,12 @@ def initialize_tpPhssource(self): # add weight source.n = nspot + # allow the possibility to count primaries + source.generate_until_next_primary = self.generate_until_next_primary + source.primary_particle_name = self.primary_particle_name + source.primary_lower_energy_threshold = self.primary_lower_energy_threshold + source.primary_PDGCode = self.primary_PDGCode + self.actual_sim_particles = tot_sim_particles def verify_phs_files_exist(self, phs_dict): diff --git a/opengate/tests/src/test060_PhsSource_UntilNexPrimary.py b/opengate/tests/src/test060_PhsSource_UntilNexPrimary.py new file mode 100644 index 000000000..538eca807 --- /dev/null +++ b/opengate/tests/src/test060_PhsSource_UntilNexPrimary.py @@ -0,0 +1,99 @@ +import test060_PhsSource_helpers as t +import opengate as gate +import uproot +import numpy as np +import pandas as pd +import gatetools.phsp as phsp + + +paths = gate.get_default_test_paths( + __file__, "test060_PhsSource_ParticleName_direct", output_folder="test060" +) + +# units +m = gate.g4_units("m") +mm = gate.g4_units("mm") +cm = gate.g4_units("cm") +nm = gate.g4_units("nm") +Bq = gate.g4_units("Bq") +MeV = gate.g4_units("MeV") +deg: float = gate.g4_units("deg") + + +def createRootFile(file_name): + # # create pandas dataframe + df = pd.DataFrame( + { + "KineticEnergy": [100, 2.0, 3.0, 80, 2.0, 3.0, 2.0, 110.0, 3.0, 111.0, 2], + "PDGCode": [2212, 11, 11, 2212, 11, 11, 11, 2212, 11, 2212, 11], + "PreDirectionLocal_X": np.zeros(11), + "PreDirectionLocal_Y": np.zeros(11), + "PreDirectionLocal_Z": np.ones(11), + "PrePositionLocal_X": np.zeros(11), + "PrePositionLocal_Y": np.zeros(11), + "PrePositionLocal_Z": np.zeros(11), + "PreDirection_X": np.zeros(11), + "PreDirection_Y": np.zeros(11), + "PreDirection_Z": np.ones(11), + "PrePosition_X": np.zeros(11), + "PrePosition_Y": np.zeros(11), + "PrePosition_Z": np.zeros(11), + "Weight": np.ones(11), + } + ) + # PDGCode proton = 2212 + # PDGCode electron = 11 + # PDGCode positron = -11 + # PDGCode photon = 22 + # PDGCode neutron = 2112 + # # "ParticleName": ["gamma", "e-", "e+"], + # # df["ParticleName"] = df["ParticleName"].astype(str) + # df = df.astype({"ParticleName": "char"}) + + # print(df, df.dtypes) + + # generate root file + tfile = uproot.recreate(file_name) + + tfile["PhaseSpace1"] = df + tfile.close() + + +def main(): + print("create reference PhS file") + createRootFile(paths.output / "test_phs.root") + + print("testing until primary") + t.test_source_untilPrimary( + source_file_name=paths.output / "test_phs.root", + phs_file_name_out=paths.output / "test_source_untilPrimary.root", + ) + + is_ok = t.check_value_from_root_file( + file_name_root=paths.output / "test_source_translation.root", + key="PrePositionLocal_X", + ref_value=30 * mm, + ) + # load data from root file + data_ref, keys_ref, m_ref = phsp.load( + paths.output / "test_source_untilPrimary.root" + ) + # the root file should contain 9 particles + # print(m_ref) + + # there should be 3 protons in the root file + index = keys_ref.index("PDGCode") + particle_list = data_ref[:, index] + # print(particle_list) + count = np.count_nonzero(particle_list == 2212) + is_ok = False + if m_ref == 9 and count == 3: + is_ok = True + print("test ok") + + # this is the end, my friend + gate.test_ok(is_ok) + + +if __name__ == "__main__": + main() diff --git a/opengate/tests/src/test060_PhsSource_helpers.py b/opengate/tests/src/test060_PhsSource_helpers.py index 1506972fb..50b0839a8 100644 --- a/opengate/tests/src/test060_PhsSource_helpers.py +++ b/opengate/tests/src/test060_PhsSource_helpers.py @@ -383,6 +383,35 @@ def test_source_rotation( output = sim.start() +def test_source_untilPrimary( + source_file_name="output/test_proton_offset.root", + phs_file_name_out="output/output/test_source_electron.root", +) -> None: + sim = create_PhS_withoutSource( + phs_name=phs_file_name_out, + ) + number_of_particles = 2 + ########################################################################################## + # Source + ########################################################################################## + # phsp source + source = sim.add_source("PhaseSpaceSource", "phsp_source_global") + source.mother = "world" + source.phsp_file = source_file_name + source.position_key = "PrePosition" + source.direction_key = "PreDirection" + source.global_flag = True + source.particle = "" + source.batch_size = 3000 + source.n = number_of_particles + source.generate_until_next_primary = True + source.primary_lower_energy_threshold = 90.0 * MeV + source.primary_PDGCode = 2212 + print(source) + + output = sim.start() + + def get_first_entry_of_key( file_name_root="output/test_source_electron.root", key="ParticleName" ) -> None: From 0da87ff0ec5e159a2a47c953ab77b3ae2fe9ddcb Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Fri, 4 Aug 2023 12:33:59 +0200 Subject: [PATCH 12/46] removed unneccessary terminal output --- opengate/source/PhaseSpaceSource.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/opengate/source/PhaseSpaceSource.py b/opengate/source/PhaseSpaceSource.py index 9f307d6a8..3762a5b07 100644 --- a/opengate/source/PhaseSpaceSource.py +++ b/opengate/source/PhaseSpaceSource.py @@ -115,11 +115,6 @@ def initialize(self, run_timing_intervals): f"PhaseSpaceSource: generate_until_next_primary is True but no primary_lower_energy_threshold is defined" ) - print( - "PhaseSpaceSource primary_lower_energy_threshold: ", - ui.primary_lower_energy_threshold, - ) - # initialize the generator (read the phsp file) self.particle_generator.initialize(self.user_info) From 1207652d63f6566f08743db14d0697478b011cc6 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Mon, 7 Aug 2023 10:23:48 +0200 Subject: [PATCH 13/46] Added safety check for batch size --- opengate/source/PhaseSpaceSourceGenerator.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/opengate/source/PhaseSpaceSourceGenerator.py b/opengate/source/PhaseSpaceSourceGenerator.py index 6e086e632..aa0731009 100644 --- a/opengate/source/PhaseSpaceSourceGenerator.py +++ b/opengate/source/PhaseSpaceSourceGenerator.py @@ -33,6 +33,8 @@ def initialize(self, user_info): def read_phsp_and_keys(self): # convert str like 1e5 to int self.user_info.batch_size = int(float(self.user_info.batch_size)) + if self.user_info.batch_size < 1: + gate.fatal("PhaseSpaceSourceGenerator: Batch size should be > 0") # open root file and get the first branch # FIXME could have an option to select the branch From 7a1213a04915e67c670aec8a19c77bf7c66147d1 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Thu, 10 Aug 2023 16:16:06 +0200 Subject: [PATCH 14/46] added a check to verify that there are branches in the root file --- opengate/source/PhaseSpaceSourceGenerator.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/opengate/source/PhaseSpaceSourceGenerator.py b/opengate/source/PhaseSpaceSourceGenerator.py index aa0731009..74b169241 100644 --- a/opengate/source/PhaseSpaceSourceGenerator.py +++ b/opengate/source/PhaseSpaceSourceGenerator.py @@ -40,7 +40,13 @@ def read_phsp_and_keys(self): # FIXME could have an option to select the branch self.root_file = uproot.open(self.user_info.phsp_file) branches = self.root_file.keys() - self.root_file = self.root_file[branches[0]] + if len(branches) > 0: + self.root_file = self.root_file[branches[0]] + else: + gate.fatal( + f"PhaseSpaceSourceGenerator: No useable branches in the root file {self.user_info.phsp_file}. Aborting." + ) + exit() # initialize the iterator self.iter = self.root_file.iterate(step_size=self.user_info.batch_size) From bebe7f2d6ed28a8566f004b4a4056c57b65f69c0 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Wed, 27 Sep 2023 13:37:26 +0200 Subject: [PATCH 15/46] small changes --- opengate/source/TreatmentPlanPhsSource.py | 1 + 1 file changed, 1 insertion(+) diff --git a/opengate/source/TreatmentPlanPhsSource.py b/opengate/source/TreatmentPlanPhsSource.py index 7a796b5a1..910fb75e9 100644 --- a/opengate/source/TreatmentPlanPhsSource.py +++ b/opengate/source/TreatmentPlanPhsSource.py @@ -120,6 +120,7 @@ def initialize_tpPhssource(self): # POSITION: source.override_position = True source.position.translation = self._get_pbs_position(spot) + print("source.position.translation: ", source.position.translation) # ROTATION: source.override_direction = True From 6152de3c73668cfe6074e90585568ca3d7999062 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Fri, 24 Nov 2023 12:50:28 +0100 Subject: [PATCH 16/46] fixed phase space source files in Gate core --- .../opengate_lib/GatePhaseSpaceSource.cpp | 260 +++++++++--------- .../opengate_lib/GatePhaseSpaceSource.h | 2 +- 2 files changed, 128 insertions(+), 134 deletions(-) diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp index f19a7311a..789fe17a5 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp @@ -50,17 +50,17 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) { fMass = fParticleDefinition->GetPDGMass(); } - fgenerate_until_next_primary = - DictGetBool(user_info, "generate_until_next_primary"); - fprimary_lower_energy_threshold = - DictGetDouble(user_info, "primary_lower_energy_threshold"); - fprimary_PDGCode = DictGetInt(user_info, "primary_PDGCode"); - fprimary_pname = DictGetStr(user_info, "primary_particle_name"); - // Init l.fNumberOfGeneratedEvents = 0; l.fCurrentIndex = 0; l.fCurrentBatchSize = 0; + + l.fgenerate_until_next_primary = + DictGetBool(user_info, "generate_until_next_primary"); + l.fprimary_lower_energy_threshold = + DictGetDouble(user_info, "primary_lower_energy_threshold"); + l.fprimary_PDGCode = DictGetInt(user_info, "primary_PDGCode"); + l.fprimary_pname = DictGetStr(user_info, "primary_particle_name"); } void GatePhaseSpaceSource::PrepareNextRun() { @@ -103,11 +103,11 @@ void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event, // check if we should simulate until next primary // in this case, generate until a second primary is in the list, excluding the // second primary - if (fgenerate_until_next_primary) { + if (l.fgenerate_until_next_primary) { int num_primaries = 0; while (num_primaries <= 2) { // If batch is empty, we generate some particles - if (fCurrentIndex >= fCurrentBatchSize) + if (l.fCurrentIndex >= l.fCurrentBatchSize) GenerateBatchOfParticles(); // check if next particle is primary @@ -120,7 +120,7 @@ void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event, GenerateOnePrimary(event, current_simulation_time); // update the index; - fCurrentIndex++; + l.fCurrentIndex++; // // update the number of generated event // fNumberOfGeneratedEvents++; @@ -128,7 +128,7 @@ void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event, break; } // update the number of generated event - fNumberOfGeneratedEvents++; + l.fNumberOfGeneratedEvents++; } else { @@ -145,145 +145,139 @@ void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event, // update the number of generated event l.fNumberOfGeneratedEvents++; } +} - void GatePhaseSpaceSource::GenerateOnePrimary( - G4Event * event, double current_simulation_time) { - auto &l = fThreadLocalDataPhsp.Get(); - - auto position = G4ThreeVector(l.fPositionX[l.fCurrentIndex], - l.fPositionY[l.fCurrentIndex], - l.fPositionZ[l.fCurrentIndex]); - auto direction = G4ParticleMomentum(l.fDirectionX[l.fCurrentIndex], - l.fDirectionY[l.fCurrentIndex], - l.fDirectionZ[l.fCurrentIndex]); - auto energy = l.fEnergy[l.fCurrentIndex]; - auto weight = l.fWeight[l.fCurrentIndex]; - - // FIXME auto time = fTime[l.fCurrentIndex]; - - // transform according to mother - if (!fGlobalFag) { - auto &ls = fThreadLocalData.Get(); - position = ls.fGlobalRotation * position + ls.fGlobalTranslation; - direction = direction / direction.mag(); - direction = ls.fGlobalRotation * direction; - } +void GatePhaseSpaceSource::GenerateOnePrimary(G4Event *event, + double current_simulation_time) { + auto &l = fThreadLocalDataPhsp.Get(); - // Create the final vertex - AddOnePrimaryVertex(event, position, direction, energy, - current_simulation_time, weight); + auto position = G4ThreeVector(l.fPositionX[l.fCurrentIndex], + l.fPositionY[l.fCurrentIndex], + l.fPositionZ[l.fCurrentIndex]); + auto direction = G4ParticleMomentum(l.fDirectionX[l.fCurrentIndex], + l.fDirectionY[l.fCurrentIndex], + l.fDirectionZ[l.fCurrentIndex]); + auto energy = l.fEnergy[l.fCurrentIndex]; + auto weight = l.fWeight[l.fCurrentIndex]; + + // FIXME auto time = fTime[l.fCurrentIndex]; + + // transform according to mother + if (!fGlobalFag) { + auto &ls = fThreadLocalData.Get(); + position = ls.fGlobalRotation * position + ls.fGlobalTranslation; + direction = direction / direction.mag(); + direction = ls.fGlobalRotation * direction; } - void GatePhaseSpaceSource::AddOnePrimaryVertex( - G4Event * event, const G4ThreeVector &position, - const G4ThreeVector &direction, double energy, double time, double w) { - auto *particle = new G4PrimaryParticle(); - auto &l = fThreadLocalDataPhsp.Get(); - if (fUseParticleTypeFromFile) { - if (l.fPDGCode[l.fCurrentIndex] != 0) { - fParticleDefinition = - fParticleTable->FindParticle(l.fPDGCode[l.fCurrentIndex]); - particle->SetParticleDefinition(fParticleDefinition); - } else { - Fatal("GatePhaseSpaceSource: PDGCode not available. Aborting."); - } - } else { + // Create the final vertex + AddOnePrimaryVertex(event, position, direction, energy, + current_simulation_time, weight); +} + +void GatePhaseSpaceSource::AddOnePrimaryVertex(G4Event *event, + const G4ThreeVector &position, + const G4ThreeVector &direction, + double energy, double time, + double w) { + auto *particle = new G4PrimaryParticle(); + auto &l = fThreadLocalDataPhsp.Get(); + if (fUseParticleTypeFromFile) { + if (l.fPDGCode[l.fCurrentIndex] != 0) { + fParticleDefinition = + fParticleTable->FindParticle(l.fPDGCode[l.fCurrentIndex]); particle->SetParticleDefinition(fParticleDefinition); - particle->SetMass(fMass); - particle->SetCharge(fCharge); + } else { + Fatal("GatePhaseSpaceSource: PDGCode not available. Aborting."); } - particle->SetKineticEnergy(energy); - particle->SetMomentumDirection(direction); - - // set vertex - auto *vertex = new G4PrimaryVertex(position, time); - vertex->SetPrimary(particle); - event->AddPrimaryVertex(vertex); - - // weights - event->GetPrimaryVertex(0)->SetWeight(w); + } else { + particle->SetParticleDefinition(fParticleDefinition); + particle->SetMass(fMass); + particle->SetCharge(fCharge); } + particle->SetKineticEnergy(energy); + particle->SetMomentumDirection(direction); - void GatePhaseSpaceSource::SetPDGCodeBatch( - const py::array_t &fPDGCode) const { - auto &l = fThreadLocalDataPhsp.Get(); - l.fPDGCode = PyBindGetVector(fPDGCode); - } + // set vertex + auto *vertex = new G4PrimaryVertex(position, time); + vertex->SetPrimary(particle); + event->AddPrimaryVertex(vertex); - void GatePhaseSpaceSource::SetEnergyBatch(const py::array_t &fEnergy) - const { - auto &l = fThreadLocalDataPhsp.Get(); - l.fEnergy = PyBindGetVector(fEnergy); - } + // weights + event->GetPrimaryVertex(0)->SetWeight(w); +} - void GatePhaseSpaceSource::SetWeightBatch(const py::array_t &fWeight) - const { - auto &l = fThreadLocalDataPhsp.Get(); - l.fWeight = PyBindGetVector(fWeight); - } +void GatePhaseSpaceSource::SetPDGCodeBatch( + const py::array_t &fPDGCode) const { + auto &l = fThreadLocalDataPhsp.Get(); + l.fPDGCode = PyBindGetVector(fPDGCode); +} - void GatePhaseSpaceSource::SetPositionXBatch( - const py::array_t &fPositionX) const { - auto &l = fThreadLocalDataPhsp.Get(); - l.fPositionX = PyBindGetVector(fPositionX); - } +void GatePhaseSpaceSource::SetEnergyBatch( + const py::array_t &fEnergy) const { + auto &l = fThreadLocalDataPhsp.Get(); + l.fEnergy = PyBindGetVector(fEnergy); +} - void GatePhaseSpaceSource::SetPositionYBatch( - const py::array_t &fPositionY) const { - auto &l = fThreadLocalDataPhsp.Get(); - l.fPositionY = PyBindGetVector(fPositionY); - } +void GatePhaseSpaceSource::SetWeightBatch( + const py::array_t &fWeight) const { + auto &l = fThreadLocalDataPhsp.Get(); + l.fWeight = PyBindGetVector(fWeight); +} - void GatePhaseSpaceSource::SetPositionZBatch( - const py::array_t &fPositionZ) const { - auto &l = fThreadLocalDataPhsp.Get(); - l.fPositionZ = PyBindGetVector(fPositionZ); - } +void GatePhaseSpaceSource::SetPositionXBatch( + const py::array_t &fPositionX) const { + auto &l = fThreadLocalDataPhsp.Get(); + l.fPositionX = PyBindGetVector(fPositionX); +} - void GatePhaseSpaceSource::SetDirectionXBatch( - const py::array_t &fDirectionX) const { - auto &l = fThreadLocalDataPhsp.Get(); - l.fDirectionX = PyBindGetVector(fDirectionX); - } +void GatePhaseSpaceSource::SetPositionYBatch( + const py::array_t &fPositionY) const { + auto &l = fThreadLocalDataPhsp.Get(); + l.fPositionY = PyBindGetVector(fPositionY); +} - void GatePhaseSpaceSource::SetDirectionYBatch( - const py::array_t &fDirectionY) const { - auto &l = fThreadLocalDataPhsp.Get(); - l.fDirectionY = PyBindGetVector(fDirectionY); - } +void GatePhaseSpaceSource::SetPositionZBatch( + const py::array_t &fPositionZ) const { + auto &l = fThreadLocalDataPhsp.Get(); + l.fPositionZ = PyBindGetVector(fPositionZ); +} - void GatePhaseSpaceSource::SetDirectionZBatch( - const py::array_t &fDirectionZ) const { - auto &l = fThreadLocalDataPhsp.Get(); - l.fDirectionZ = PyBindGetVector(fDirectionZ); - } +void GatePhaseSpaceSource::SetDirectionXBatch( + const py::array_t &fDirectionX) const { + auto &l = fThreadLocalDataPhsp.Get(); + l.fDirectionX = PyBindGetVector(fDirectionX); +} - bool GatePhaseSpaceSource::ParticleIsPrimary() { - // check if particle is primary - bool is_primary = false; +void GatePhaseSpaceSource::SetDirectionYBatch( + const py::array_t &fDirectionY) const { + auto &l = fThreadLocalDataPhsp.Get(); + l.fDirectionY = PyBindGetVector(fDirectionY); +} - // if PDGCode exists in file - if ((fPDGCode[fCurrentIndex] != 0) && (fprimary_PDGCode != 0)) { - if ((fprimary_PDGCode == fPDGCode[fCurrentIndex]) && - (fprimary_lower_energy_threshold <= fEnergy[fCurrentIndex])) { - is_primary = true; - } - } - // if PDGCode does not exist in file, but particle name does and is not - // empty - else if ((fParticleName[fCurrentIndex].length() != 0) && - (fprimary_pname.length() != 0)) { - if ((fprimary_pname == fParticleName[fCurrentIndex]) && - (fprimary_lower_energy_threshold <= fEnergy[fCurrentIndex])) { - is_primary = true; - } - } else { - G4Exception("GatePhaseSpaceSource::ParticleIsPrimary", "Error", - FatalException, "Particle type not defined in file"); - std::cout << "ERROR: Particle name nor PDGCode defined in file. Aborting." - << std::endl; - exit(1); - } +void GatePhaseSpaceSource::SetDirectionZBatch( + const py::array_t &fDirectionZ) const { + auto &l = fThreadLocalDataPhsp.Get(); + l.fDirectionZ = PyBindGetVector(fDirectionZ); +} - return is_primary; +bool GatePhaseSpaceSource::ParticleIsPrimary() { + auto &l = fThreadLocalDataPhsp.Get(); + // check if particle is primary + bool is_primary = false; + + // if PDGCode exists in file + if ((l.fPDGCode[l.fCurrentIndex] != 0) && (l.fprimary_PDGCode != 0)) { + if ((l.fprimary_PDGCode == l.fPDGCode[l.fCurrentIndex]) && + (l.fprimary_lower_energy_threshold <= l.fEnergy[l.fCurrentIndex])) { + is_primary = true; + } + } else { + G4Exception("GatePhaseSpaceSource::ParticleIsPrimary", "Error", + FatalException, "Particle type not defined in file"); + std::cout << "ERROR: PDGCode not defined in file. Aborting." << std::endl; + exit(1); } + + return is_primary; +} diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h index fe25a48bc..4d4268412 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h @@ -43,7 +43,7 @@ class GatePhaseSpaceSource : public GateVSource { const G4ThreeVector &momentum_direction, double energy, double time, double w); - void SetGeneratorFunction(ParticleGeneratorType &f); + void SetGeneratorFunction(ParticleGeneratorType &f) const; bool ParticleIsPrimary(); From 4a2884d672da1dd1ba19230c032ea97fcb058600 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Fri, 24 Nov 2023 13:38:18 +0100 Subject: [PATCH 17/46] fixed test060 --- .../opengate_lib/GatePhaseSpaceSource.cpp | 1 - .../opengate_lib/GatePhaseSpaceSource.h | 1 - opengate/source/PhaseSpaceSource.py | 122 -------------- opengate/source/PhaseSpaceSourceGenerator.py | 157 ------------------ opengate/sources/TreatmentPlanPhsSource.py | 2 - opengate/sources/phspsources.py | 26 ++- .../tests/src/test060_phsp_source_helpers.py | 2 +- ...> test060_phsp_source_untilNextPrimary.py} | 21 +-- 8 files changed, 37 insertions(+), 295 deletions(-) delete mode 100644 opengate/source/PhaseSpaceSource.py delete mode 100644 opengate/source/PhaseSpaceSourceGenerator.py rename opengate/tests/src/{test060_PhsSource_UntilNexPrimary.py => test060_phsp_source_untilNextPrimary.py} (89%) diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp index 789fe17a5..b888601fa 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp @@ -60,7 +60,6 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) { l.fprimary_lower_energy_threshold = DictGetDouble(user_info, "primary_lower_energy_threshold"); l.fprimary_PDGCode = DictGetInt(user_info, "primary_PDGCode"); - l.fprimary_pname = DictGetStr(user_info, "primary_particle_name"); } void GatePhaseSpaceSource::PrepareNextRun() { diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h index 4d4268412..b5998b90e 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h @@ -82,7 +82,6 @@ class GatePhaseSpaceSource : public GateVSource { struct threadLocalTPhsp { bool fgenerate_until_next_primary; - string fprimary_pname; int fprimary_PDGCode; double fprimary_lower_energy_threshold; diff --git a/opengate/source/PhaseSpaceSource.py b/opengate/source/PhaseSpaceSource.py deleted file mode 100644 index 3762a5b07..000000000 --- a/opengate/source/PhaseSpaceSource.py +++ /dev/null @@ -1,122 +0,0 @@ -from .SourceBase import * -import opengate_core as g4 -from .PhaseSpaceSourceGenerator import * -from scipy.spatial.transform import Rotation -from box import Box - - -class PhaseSpaceSource(SourceBase): - """ - Source of particles from a (root) phase space. - Read position + direction + energy + weight from the root and use them as event. - - if global flag is True, the position/direction are global, not - in the coordinate system of the mother volume. - if the global flag is False, the position/direction are relative - to the mother volume - - - Timing is not used (yet) - - NOT ready for multithread (yet) - - type of particle not read in the phase space but set by user - - """ - - type_name = "PhaseSpaceSource" - - @staticmethod - def set_default_user_info(user_info): - gate.SourceBase.set_default_user_info(user_info) - # initial user info - user_info.phsp_file = None - user_info.n = 1 - user_info.particle = "" # FIXME later as key - # if a particle name is supplied, the particle type is set to it - # otherwise, information from the phase space is used - - # if global flag is True, the position/direction are global, not - # in the coordinate system of the mother volume. - # if the global flag is False, the position/direction are relative - # to the mother volume - user_info.global_flag = False - user_info.batch_size = 10000 - # branch name in the phsp file - user_info.position_key = "PrePositionLocal" - user_info.position_key_x = "" - user_info.position_key_y = "" - user_info.position_key_z = "" - user_info.direction_key = "PreDirectionLocal" - user_info.direction_key_x = "" - user_info.direction_key_y = "" - user_info.direction_key_z = "" - user_info.energy_key = "KineticEnergy" - user_info.weight_key = "Weight" - user_info.particle_name_key = "ParticleName" - user_info.PDGCode_key = "PDGCode" - # change position and direction of the source - # position is relative to the stored coordinates - # direction is a rotation of the stored direction - user_info.override_position = False - user_info.override_direction = False - user_info.position = Box() - user_info.position.translation = [0, 0, 0] - user_info.position.rotation = Rotation.identity().as_matrix() - user_info.generate_until_next_primary = False - user_info.primary_particle_name = "" - user_info.primary_lower_energy_threshold = 0 - user_info.primary_PDGCode = 0 - # user_info.time_key = None # FIXME later - - def __del__(self): - pass - - def create_g4_source(self): - return g4.GatePhaseSpaceSource() - - def __init__(self, user_info): - super().__init__(user_info) - self.particle_generator = gate.PhaseSpaceSourceGenerator() - - def initialize(self, run_timing_intervals): - # initialize the mother class generic source - - gate.SourceBase.initialize(self, run_timing_intervals) - if self.simulation.use_multithread: - gate.fatal( - f"Cannot use phsp source in MT mode for the moment" - f" (need to create a generator that read the root tree randomly" - ) - - # check user info - # only if position_key_x etc is not set, then set then based in position_key - ui = self.user_info - if ui.position_key_x == "": - ui.position_key_x = f"{ui.position_key}_X" - if ui.position_key_y == "": - ui.position_key_y = f"{ui.position_key}_Y" - if ui.position_key_z == "": - ui.position_key_z = f"{ui.position_key}_Z" - # only if direction_key_x etc is not set, then set then based in direction_key - if ui.direction_key_x == "": - ui.direction_key_x = f"{ui.direction_key}_X" - if ui.direction_key_y == "": - ui.direction_key_y = f"{ui.direction_key}_Y" - if ui.direction_key_z == "": - ui.direction_key_z = f"{ui.direction_key}_Z" - - # check if the source should generate particles until the second one - # which is identified as primary by name, PDGCode and above a threshold - if ui.generate_until_next_primary == True: - if len(ui.primary_particle_name) <= 0 and ui.primary_PDGCode == 0: - gate.fatal( - f"PhaseSpaceSource: generate_until_next_primary is True but no primary particle is defined" - ) - if ui.primary_lower_energy_threshold <= 0: - gate.fatal( - f"PhaseSpaceSource: generate_until_next_primary is True but no primary_lower_energy_threshold is defined" - ) - - # initialize the generator (read the phsp file) - self.particle_generator.initialize(self.user_info) - - # set the function pointer to the cpp side - self.g4_source.SetGeneratorFunction(self.particle_generator.generate) diff --git a/opengate/source/PhaseSpaceSourceGenerator.py b/opengate/source/PhaseSpaceSourceGenerator.py deleted file mode 100644 index 74b169241..000000000 --- a/opengate/source/PhaseSpaceSourceGenerator.py +++ /dev/null @@ -1,157 +0,0 @@ -import threading -import uproot -import opengate as gate -from scipy.spatial.transform import Rotation -import numpy as np - - -class PhaseSpaceSourceGenerator: - """ - Class that read phase space root file and extract position/direction/energy/weights of particles. - Particles information will be copied to the c++ side to be used as a source - """ - - def __init__(self): - self.lock = threading.Lock() - self.initialize_is_done = False - self.user_info = None - self.root_file = None - self.num_entries = 0 - self.cycle_count = 0 - - def __getstate__(self): - self.lock = None - return self.__dict__ - - def initialize(self, user_info): - with self.lock: - if not self.initialize_is_done: - self.user_info = user_info - self.read_phsp_and_keys() - self.initialize_is_done = True - - def read_phsp_and_keys(self): - # convert str like 1e5 to int - self.user_info.batch_size = int(float(self.user_info.batch_size)) - if self.user_info.batch_size < 1: - gate.fatal("PhaseSpaceSourceGenerator: Batch size should be > 0") - - # open root file and get the first branch - # FIXME could have an option to select the branch - self.root_file = uproot.open(self.user_info.phsp_file) - branches = self.root_file.keys() - if len(branches) > 0: - self.root_file = self.root_file[branches[0]] - else: - gate.fatal( - f"PhaseSpaceSourceGenerator: No useable branches in the root file {self.user_info.phsp_file}. Aborting." - ) - exit() - - # initialize the iterator - self.iter = self.root_file.iterate(step_size=self.user_info.batch_size) - - # initialize counters - self.num_entries = int(self.root_file.num_entries) - self.cycle_count = 0 - - def generate(self, source): - """ - Main function that will be called from the cpp side every time a batch - of particles should be created. - Once created here, the particles are copied to cpp. - (Yes maybe the copy could be avoided, but I did not manage to do it) - """ - - # read data from root tree - try: - batch = next(self.iter) - except: - self.cycle_count += 1 - gate.warning( - f"End of the phase-space {self.num_entries} elements, " - f"restart from beginning. Cycle count = {self.cycle_count}" - ) - self.iter = self.root_file.iterate(step_size=self.user_info.batch_size) - batch = next(self.iter) - # print("batch", batch) - # print("type of batch", type(batch)) - # print("type of batch[0]", type(batch[0])) - # print("batch fields", batch.fields) - - # copy to cpp - ui = self.user_info - - # # check if the keys for particle name or PDGCode are in the root file - - if ui.PDGCode_key in batch.fields: - source.fPDGCode = batch[ui.PDGCode_key] - # print("source.fPDGCode", source.fPDGCode) - # print("type of source.fPDGCode", type(source.fPDGCode)) - # print("source.fPDGCode in there") - else: - source.fPDGCode = np.zeros(len(batch), dtype=int) - # print("type of source.fPDGCode", type(source.fPDGCode)) - # print("source.fPDGCode", source.fPDGCode) - if ui.particle_name_key in batch.fields: - # print("source.fParticleName in there") - source.fParticleName = batch[ui.particle_name_key] - else: - source.fParticleName = [""] * len(batch) - - # if neither particle name, nor PDGCode are in the root file, - # nor the particle type was set, raise an error - if ( - ui.particle == "" - and ui.PDGCode_key not in batch.fields - and ui.particle_name_key not in batch.fields - ): - print( - "ERROR: PhaseSpaceSource: No particle name or PDGCode key in the root file, " - "and no particle type was set. " - "Please set the particle type or add the particle name or PDGCode key in the root file. Aborting." - ) - exit(-1) - - # change position and direction of the source - - # if override_position is set to True, the position - # supplied will be added to the root file position - if ui.override_position: - source.fPositionX = batch[ui.position_key_x] + ui.position.translation[0] - source.fPositionY = batch[ui.position_key_y] + ui.position.translation[1] - source.fPositionZ = batch[ui.position_key_z] + ui.position.translation[2] - else: - source.fPositionX = batch[ui.position_key_x] - source.fPositionY = batch[ui.position_key_y] - source.fPositionZ = batch[ui.position_key_z] - - # direction is a rotation of the stored direction - # if override_direction is set to True, the direction - # in the root file will be rotated based on the supplied rotation matrix - if ui.override_direction: - # create point vectors - points = np.column_stack( - ( - batch[ui.direction_key_x], - batch[ui.direction_key_y], - batch[ui.direction_key_z], - ) - ) - # create rotation matrix - r = Rotation.from_matrix(ui.position.rotation) - # rotate vector with rotation matrix - points = r.apply(points) - # assign rotated vector to direction - # source.fDirectionX = points[:, 0] - # source.fDirectionY = points[:, 1] - # source.fDirectionZ = points[:, 2] - source.fDirectionX, source.fDirectionY, source.fDirectionZ = points.T - else: - source.fDirectionX = batch[ui.direction_key_x] - source.fDirectionY = batch[ui.direction_key_y] - source.fDirectionZ = batch[ui.direction_key_z] - - # pass energy and weight - source.fEnergy = batch[ui.energy_key] - source.fWeight = batch[ui.weight_key] diff --git a/opengate/sources/TreatmentPlanPhsSource.py b/opengate/sources/TreatmentPlanPhsSource.py index 910fb75e9..b8410720d 100644 --- a/opengate/sources/TreatmentPlanPhsSource.py +++ b/opengate/sources/TreatmentPlanPhsSource.py @@ -24,7 +24,6 @@ def __init__(self, name, sim): self.weight_key = "Weight" self.PDGCode_key = "PDGCode" self.generate_until_next_primary = False - self.primary_particle_name = "" self.primary_lower_energy_threshold = 0 self.primary_PDGCode = 0 self.n_sim = 0 @@ -131,7 +130,6 @@ def initialize_tpPhssource(self): # allow the possibility to count primaries source.generate_until_next_primary = self.generate_until_next_primary - source.primary_particle_name = self.primary_particle_name source.primary_lower_energy_threshold = self.primary_lower_energy_threshold source.primary_PDGCode = self.primary_PDGCode diff --git a/opengate/sources/phspsources.py b/opengate/sources/phspsources.py index 084cebe4a..5a5353da6 100644 --- a/opengate/sources/phspsources.py +++ b/opengate/sources/phspsources.py @@ -38,6 +38,8 @@ def initialize(self, user_info): def read_phsp_and_keys(self): # convert str like 1e5 to int self.user_info.batch_size = int(float(self.user_info.batch_size)) + if self.user_info.batch_size < 1: + gate.fatal("PhaseSpaceSourceGenerator: Batch size should be > 0") if ( opengate_core.IsMultithreadedApplication() @@ -50,7 +52,14 @@ def read_phsp_and_keys(self): # FIXME could have an option to select the branch self.root_file = uproot.open(self.user_info.phsp_file) branches = self.root_file.keys() - self.root_file = self.root_file[branches[0]] + if len(branches) > 0: + self.root_file = self.root_file[branches[0]] + else: + gate.fatal( + f"PhaseSpaceSourceGenerator: No useable branches in the root file {self.user_info.phsp_file}. Aborting." + ) + exit() + self.num_entries = int(self.root_file.num_entries) # initialize the index to start @@ -256,6 +265,9 @@ def set_default_user_info(user_info): user_info.position = Box() user_info.position.translation = [0, 0, 0] user_info.position.rotation = Rotation.identity().as_matrix() + user_info.generate_until_next_primary = False + user_info.primary_lower_energy_threshold = 0 + user_info.primary_PDGCode = 0 # user_info.time_key = None # FIXME TODO later # for debug user_info.verbose_batch = False @@ -290,6 +302,18 @@ def initialize(self, run_timing_intervals): if ui.direction_key_z is None: ui.direction_key_z = f"{ui.direction_key}_Z" + # check if the source should generate particles until the second one + # which is identified as primary by name, PDGCode and above a threshold + if ui.generate_until_next_primary == True: + if ui.primary_PDGCode == 0: + gate.fatal( + f"PhaseSpaceSource: generate_until_next_primary is True but no primary particle is defined" + ) + if ui.primary_lower_energy_threshold <= 0: + gate.fatal( + f"PhaseSpaceSource: generate_until_next_primary is True but no primary_lower_energy_threshold is defined" + ) + # initialize the generator (read the phsp file) self.particle_generator.initialize(self.user_info) diff --git a/opengate/tests/src/test060_phsp_source_helpers.py b/opengate/tests/src/test060_phsp_source_helpers.py index 82b8545f7..4112e1a8c 100644 --- a/opengate/tests/src/test060_phsp_source_helpers.py +++ b/opengate/tests/src/test060_phsp_source_helpers.py @@ -377,7 +377,7 @@ def test_source_untilPrimary( source_file_name="output/test_proton_offset.root", phs_file_name_out="output/output/test_source_electron.root", ) -> None: - sim = create_PhS_withoutSource( + sim = create_phs_without_source( phs_name=phs_file_name_out, ) number_of_particles = 2 diff --git a/opengate/tests/src/test060_PhsSource_UntilNexPrimary.py b/opengate/tests/src/test060_phsp_source_untilNextPrimary.py similarity index 89% rename from opengate/tests/src/test060_PhsSource_UntilNexPrimary.py rename to opengate/tests/src/test060_phsp_source_untilNextPrimary.py index 538eca807..3bb1506b6 100644 --- a/opengate/tests/src/test060_PhsSource_UntilNexPrimary.py +++ b/opengate/tests/src/test060_phsp_source_untilNextPrimary.py @@ -1,23 +1,24 @@ -import test060_PhsSource_helpers as t +import test060_phsp_source_helpers as t import opengate as gate import uproot import numpy as np import pandas as pd import gatetools.phsp as phsp +from opengate.tests import utility -paths = gate.get_default_test_paths( +paths = utility.get_default_test_paths( __file__, "test060_PhsSource_ParticleName_direct", output_folder="test060" ) # units -m = gate.g4_units("m") -mm = gate.g4_units("mm") -cm = gate.g4_units("cm") -nm = gate.g4_units("nm") -Bq = gate.g4_units("Bq") -MeV = gate.g4_units("MeV") -deg: float = gate.g4_units("deg") +m = gate.g4_units.m +mm = gate.g4_units.mm +cm = gate.g4_units.cm +nm = gate.g4_units.nm +Bq = gate.g4_units.Bq +MeV = gate.g4_units.MeV +deg: float = gate.g4_units.deg def createRootFile(file_name): @@ -92,7 +93,7 @@ def main(): print("test ok") # this is the end, my friend - gate.test_ok(is_ok) + utility.test_ok(is_ok) if __name__ == "__main__": From 58072a490581d386e820fa9ceebe23298aff6eaa Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Fri, 24 Nov 2023 14:14:12 +0100 Subject: [PATCH 18/46] moved opengate/sources/treatmentPlanPhsSource.py to contrib/tps and fixed associated tests061 --- .../tps/treatmentPlanPhsSource.py} | 6 ++- opengate/tests/src/test061_TPPhsSource.py | 22 ++++++---- .../tests/src/test061_TPPhsSource_helpers.py | 42 +++++++------------ 3 files changed, 32 insertions(+), 38 deletions(-) rename opengate/{sources/TreatmentPlanPhsSource.py => contrib/tps/treatmentPlanPhsSource.py} (99%) diff --git a/opengate/sources/TreatmentPlanPhsSource.py b/opengate/contrib/tps/treatmentPlanPhsSource.py similarity index 99% rename from opengate/sources/TreatmentPlanPhsSource.py rename to opengate/contrib/tps/treatmentPlanPhsSource.py index b8410720d..d5190252a 100644 --- a/opengate/sources/TreatmentPlanPhsSource.py +++ b/opengate/contrib/tps/treatmentPlanPhsSource.py @@ -1,9 +1,13 @@ import numpy as np from scipy.spatial.transform import Rotation import opengate as gate -from .TreatmentPlanSource import * + +from opengate.contrib.tps.ionbeamtherapy import * + import os +# TreatmentPlanSource + class TreatmentPlanPhsSource(TreatmentPlanSource): def __init__(self, name, sim): diff --git a/opengate/tests/src/test061_TPPhsSource.py b/opengate/tests/src/test061_TPPhsSource.py index 4b5f5f7e6..dc2d3281a 100644 --- a/opengate/tests/src/test061_TPPhsSource.py +++ b/opengate/tests/src/test061_TPPhsSource.py @@ -1,21 +1,25 @@ import test061_TPPhsSource_helpers as t import opengate as gate from scipy.spatial.transform import Rotation +from opengate.tests import utility +from opengate.contrib.tps.treatmentPlanPhsSource import TreatmentPlanPhsSource +from opengate.contrib.tps.ionbeamtherapy import spots_info_from_txt, TreatmentPlanSource -paths = gate.get_default_test_paths(__file__, "test061_TPPhsSource") + +paths = utility.get_default_test_paths(__file__, "test061_TPPhsSource") paths.output_ref = paths.output_ref / "test061_ref" ref_path = paths.output_ref # units -m = gate.g4_units("m") -mm = gate.g4_units("mm") -cm = gate.g4_units("cm") -nm = gate.g4_units("nm") -Bq = gate.g4_units("Bq") -MeV = gate.g4_units("MeV") -deg: float = gate.g4_units("deg") +m = gate.g4_units.m +mm = gate.g4_units.mm +cm = gate.g4_units.cm +nm = gate.g4_units.nm +Bq = gate.g4_units.Bq +MeV = gate.g4_units.MeV +deg: float = gate.g4_units.deg def main(): @@ -68,7 +72,7 @@ def main(): ) # this is the end, my friend - gate.test_ok(a and b and c and d and e and f) + utility.test_ok(a and b and c and d and e and f) if __name__ == "__main__": diff --git a/opengate/tests/src/test061_TPPhsSource_helpers.py b/opengate/tests/src/test061_TPPhsSource_helpers.py index c6866d0ce..ae468d78c 100644 --- a/opengate/tests/src/test061_TPPhsSource_helpers.py +++ b/opengate/tests/src/test061_TPPhsSource_helpers.py @@ -5,6 +5,9 @@ from scipy.spatial.transform import Rotation import gatetools.phsp as phsp import os +from opengate.tests import utility +from opengate.contrib.tps.treatmentPlanPhsSource import TreatmentPlanPhsSource +from opengate.contrib.tps.ionbeamtherapy import spots_info_from_txt, TreatmentPlanSource # paths = gate.get_default_test_paths( @@ -13,13 +16,13 @@ # paths = {output: "output"} # units -m = gate.g4_units("m") -mm = gate.g4_units("mm") -cm = gate.g4_units("cm") -nm = gate.g4_units("nm") -Bq = gate.g4_units("Bq") -MeV = gate.g4_units("MeV") -deg: float = gate.g4_units("deg") +m = gate.g4_units.m +mm = gate.g4_units.mm +cm = gate.g4_units.cm +nm = gate.g4_units.nm +Bq = gate.g4_units.Bq +MeV = gate.g4_units.MeV +deg: float = gate.g4_units.deg def create_test_Phs( @@ -41,14 +44,6 @@ def create_test_Phs( ui.number_of_threads = 1 ui.random_seed = "auto" - # units - m = gate.g4_units("m") - mm = gate.g4_units("mm") - cm = gate.g4_units("cm") - nm = gate.g4_units("nm") - Bq = gate.g4_units("Bq") - MeV = gate.g4_units("MeV") - ########################################################################################## # geometry ########################################################################################## @@ -138,15 +133,6 @@ def create_PhS_withoutSource( ui.number_of_threads = 1 ui.random_seed = "auto" - # units - m = gate.g4_units("m") - mm = gate.g4_units("mm") - cm = gate.g4_units("cm") - nm = gate.g4_units("nm") - Bq = gate.g4_units("Bq") - MeV = gate.g4_units("MeV") - deg: float = gate.g4_units("deg") - ########################################################################################## # geometry ########################################################################################## @@ -228,10 +214,10 @@ def test_source_rotation_A( # Source ########################################################################################## # TreatmentPlanPhsSource source - tpPhSs = gate.TreatmentPlanPhsSource("RT_plan", sim) + tpPhSs = TreatmentPlanPhsSource("RT_plan", sim) tpPhSs.set_phaseSpaceList_file_name(phs_list_file_name) tpPhSs.set_phaseSpaceFolder(phs_folder_name) - spots, ntot, energies, G = gate.spots_info_from_txt(plan_file_name, "") + spots, ntot, energies, G = spots_info_from_txt(plan_file_name, "") # print(spots, ntot, energies, G) tpPhSs.set_spots(spots) tpPhSs.set_particles_to_simulate(number_of_particles) @@ -267,10 +253,10 @@ def check_value_from_root_file( # read root file value = get_first_entry_of_key(file_name_root=file_name_root, key=key) if (type(ref_value) != str) and (type(value) != str): - is_ok = gate.check_diff_abs( + is_ok = utility.check_diff_abs( float(value), float(ref_value), tolerance=1e-3, txt=key ) - # gate.check_diff_abs(float(value), float(ref_value), tolerance=1e-6, txt=key) + # utility.check_diff_abs(float(value), float(ref_value), tolerance=1e-6, txt=key) else: if value == ref_value: # print("Is correct") From 812bf62f1f116649b2ea34b018296169be3f24dc Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Mon, 27 Nov 2023 09:41:01 +0100 Subject: [PATCH 19/46] removed pandas requirement from test060_phsp_source_untilNextPrimary.py --- .../test060_phsp_source_untilNextPrimary.py | 66 +++++++++---------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/opengate/tests/src/test060_phsp_source_untilNextPrimary.py b/opengate/tests/src/test060_phsp_source_untilNextPrimary.py index 3bb1506b6..de675743c 100644 --- a/opengate/tests/src/test060_phsp_source_untilNextPrimary.py +++ b/opengate/tests/src/test060_phsp_source_untilNextPrimary.py @@ -2,7 +2,6 @@ import opengate as gate import uproot import numpy as np -import pandas as pd import gatetools.phsp as phsp from opengate.tests import utility @@ -22,41 +21,42 @@ def createRootFile(file_name): - # # create pandas dataframe - df = pd.DataFrame( - { - "KineticEnergy": [100, 2.0, 3.0, 80, 2.0, 3.0, 2.0, 110.0, 3.0, 111.0, 2], - "PDGCode": [2212, 11, 11, 2212, 11, 11, 11, 2212, 11, 2212, 11], - "PreDirectionLocal_X": np.zeros(11), - "PreDirectionLocal_Y": np.zeros(11), - "PreDirectionLocal_Z": np.ones(11), - "PrePositionLocal_X": np.zeros(11), - "PrePositionLocal_Y": np.zeros(11), - "PrePositionLocal_Z": np.zeros(11), - "PreDirection_X": np.zeros(11), - "PreDirection_Y": np.zeros(11), - "PreDirection_Z": np.ones(11), - "PrePosition_X": np.zeros(11), - "PrePosition_Y": np.zeros(11), - "PrePosition_Z": np.zeros(11), - "Weight": np.ones(11), - } - ) - # PDGCode proton = 2212 - # PDGCode electron = 11 - # PDGCode positron = -11 - # PDGCode photon = 22 - # PDGCode neutron = 2112 - # # "ParticleName": ["gamma", "e-", "e+"], - # # df["ParticleName"] = df["ParticleName"].astype(str) - # df = df.astype({"ParticleName": "char"}) - - # print(df, df.dtypes) + # create numpy arrays + kinetic_energy = np.array([100, 2.0, 3.0, 80, 2.0, 3.0, 2.0, 110.0, 3.0, 111.0, 2]) + pdg_code = np.array([2212, 11, 11, 2212, 11, 11, 11, 2212, 11, 2212, 11]) + pre_direction_local_x = np.zeros(11) + pre_direction_local_y = np.zeros(11) + pre_direction_local_z = np.ones(11) + pre_position_local_x = np.zeros(11) + pre_position_local_y = np.zeros(11) + pre_position_local_z = np.zeros(11) + pre_direction_x = np.zeros(11) + pre_direction_y = np.zeros(11) + pre_direction_z = np.ones(11) + pre_position_x = np.zeros(11) + pre_position_y = np.zeros(11) + pre_position_z = np.zeros(11) + weight = np.ones(11) # generate root file tfile = uproot.recreate(file_name) - - tfile["PhaseSpace1"] = df + tfile["PhaseSpace1"] = { + "KineticEnergy": kinetic_energy, + "PDGCode": pdg_code, + "PreDirectionLocal_X": pre_direction_local_x, + "PreDirectionLocal_Y": pre_direction_local_y, + "PreDirectionLocal_Z": pre_direction_local_z, + "PrePositionLocal_X": pre_position_local_x, + "PrePositionLocal_Y": pre_position_local_y, + "PrePositionLocal_Z": pre_position_local_z, + "PreDirection_X": pre_direction_x, + "PreDirection_Y": pre_direction_y, + "PreDirection_Z": pre_direction_z, + "PrePosition_X": pre_position_x, + "PrePosition_Y": pre_position_y, + "PrePosition_Z": pre_position_z, + "Weight": weight, + } tfile.close() From 927c17a433ed6ff5b8d983792b52ce7dbaf49292 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Tue, 28 Nov 2023 08:35:55 +0100 Subject: [PATCH 20/46] changed override_position to translate_position and override_direction to rotate_direction. --- opengate/contrib/tps/treatmentPlanPhsSource.py | 4 ++-- opengate/sources/phspsources.py | 12 ++++++------ opengate/tests/src/test060_phsp_source_helpers.py | 6 +++--- ...py => test060_phsp_source_particletype_direct.py} | 0 ..._source_particletype_fromphs_particlename_wip.py} | 0 ...t060_phsp_source_particletype_fromphs_pdgcode.py} | 0 6 files changed, 11 insertions(+), 11 deletions(-) rename opengate/tests/src/{test060_phsp_source_particlename_direct.py => test060_phsp_source_particletype_direct.py} (100%) rename opengate/tests/src/{test060_phsp_source_particlename_fromphs_particlename_wip.py => test060_phsp_source_particletype_fromphs_particlename_wip.py} (100%) rename opengate/tests/src/{test060_phsp_source_particlename_fromphs_pdgcode.py => test060_phsp_source_particletype_fromphs_pdgcode.py} (100%) diff --git a/opengate/contrib/tps/treatmentPlanPhsSource.py b/opengate/contrib/tps/treatmentPlanPhsSource.py index d5190252a..9f6feefc9 100644 --- a/opengate/contrib/tps/treatmentPlanPhsSource.py +++ b/opengate/contrib/tps/treatmentPlanPhsSource.py @@ -121,12 +121,12 @@ def initialize_tpPhssource(self): source.batch_size = 30000 # POSITION: - source.override_position = True + source.translate_position = True source.position.translation = self._get_pbs_position(spot) print("source.position.translation: ", source.position.translation) # ROTATION: - source.override_direction = True + source.rotate_direction = True source.position.rotation = self._get_pbs_rotation(spot) # add weight diff --git a/opengate/sources/phspsources.py b/opengate/sources/phspsources.py index 5a5353da6..46777e97c 100644 --- a/opengate/sources/phspsources.py +++ b/opengate/sources/phspsources.py @@ -153,9 +153,9 @@ def generate(self, source): f"in the phsp file and no source.particle" ) - # if override_position is set to True, the position + # if translate_position is set to True, the position # supplied will be added to the phsp file position - if ui.override_position: + if ui.translate_position: x = batch[ui.position_key_x] + ui.position.translation[0] y = batch[ui.position_key_y] + ui.position.translation[1] z = batch[ui.position_key_z] + ui.position.translation[2] @@ -168,9 +168,9 @@ def generate(self, source): source.SetPositionZBatch(batch[ui.position_key_z]) # direction is a rotation of the stored direction - # if override_direction is set to True, the direction + # if rotate_direction is set to True, the direction # in the root file will be rotated based on the supplied rotation matrix - if ui.override_direction: + if ui.rotate_direction: # create point vectors self.points = np.column_stack( ( @@ -260,8 +260,8 @@ def set_default_user_info(user_info): # change position and direction of the source # position is relative to the stored coordinates # direction is a rotation of the stored direction - user_info.override_position = False - user_info.override_direction = False + user_info.translate_position = False + user_info.rotate_direction = False user_info.position = Box() user_info.position.translation = [0, 0, 0] user_info.position.rotation = Rotation.identity().as_matrix() diff --git a/opengate/tests/src/test060_phsp_source_helpers.py b/opengate/tests/src/test060_phsp_source_helpers.py index 4112e1a8c..506ba9c64 100644 --- a/opengate/tests/src/test060_phsp_source_helpers.py +++ b/opengate/tests/src/test060_phsp_source_helpers.py @@ -334,7 +334,7 @@ def test_source_translation( source.particle = "proton" source.batch_size = 3000 source.n = number_of_particles - source.override_position = True + source.translate_position = True source.position.translation = [3 * cm, 0 * cm, 0 * cm] print(source) @@ -362,9 +362,9 @@ def test_source_rotation( source.particle = "proton" source.batch_size = 3000 source.n = number_of_particles - # source.override_position = True + # source.translate_position = True # source.position.translation = [3 * cm, 1 * cm, 0 * cm] - source.override_direction = True + source.rotate_direction = True # rotation = Rotation.from_euler("zyx", [30, 20, 10], degrees=True) rotation = Rotation.from_euler("x", [30], degrees=True) source.position.rotation = rotation.as_matrix() diff --git a/opengate/tests/src/test060_phsp_source_particlename_direct.py b/opengate/tests/src/test060_phsp_source_particletype_direct.py similarity index 100% rename from opengate/tests/src/test060_phsp_source_particlename_direct.py rename to opengate/tests/src/test060_phsp_source_particletype_direct.py diff --git a/opengate/tests/src/test060_phsp_source_particlename_fromphs_particlename_wip.py b/opengate/tests/src/test060_phsp_source_particletype_fromphs_particlename_wip.py similarity index 100% rename from opengate/tests/src/test060_phsp_source_particlename_fromphs_particlename_wip.py rename to opengate/tests/src/test060_phsp_source_particletype_fromphs_particlename_wip.py diff --git a/opengate/tests/src/test060_phsp_source_particlename_fromphs_pdgcode.py b/opengate/tests/src/test060_phsp_source_particletype_fromphs_pdgcode.py similarity index 100% rename from opengate/tests/src/test060_phsp_source_particlename_fromphs_pdgcode.py rename to opengate/tests/src/test060_phsp_source_particletype_fromphs_pdgcode.py From b52563579ad85c8368fa95d5a307d7a85f88506e Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Tue, 28 Nov 2023 09:16:03 +0100 Subject: [PATCH 21/46] added documentation to new PhaseSpaceSource features --- docs/source/user_guide_2_2_sources.md | 47 +++++++++++++++++++++++++-- 1 file changed, 44 insertions(+), 3 deletions(-) diff --git a/docs/source/user_guide_2_2_sources.md b/docs/source/user_guide_2_2_sources.md index 899680513..c8fb7bd9d 100644 --- a/docs/source/user_guide_2_2_sources.md +++ b/docs/source/user_guide_2_2_sources.md @@ -133,7 +133,7 @@ Like all objects, by default, the source is located according to the coordinate ### Phase-Space sources -A phase-space source read particles properties (position, direction, energy, etc.) from a root file and use them as events. Here is an example: +A phase-space source reads particles properties (position, direction, energy, etc.) from a root file and use them as events. Here is an example: ```python source = sim.add_source("PhaseSpaceSource", "phsp_source") @@ -149,15 +149,56 @@ source.n = 20000 In that case, the key "PrePositionLocal" in the root tree file will be used to define the position of all generated particles. The flag "global_flag" is False so the position will be relative to the mother volume (the plane here) ; otherwise, position is considered as global (in the world coordinate system). + Limitation: the particle timestamps is NOT read from the phsp and not considered (yet) -The particle type can be set by ```source.particle = "proton"``` option (all generated particles will be for example proton), or read in the phsp file by using the PDGCode: +The particle type can be set by ```source.particle = "proton"``` option. Using this option all generated particles will be for example protons, overriding the particle type specified in the phase space. +Alternatively, by setting ```source.particle = None``` the particle type is read from the phase space file using the PDGCode. ```source.PDGCode_key = PDGCode``` specifies the name of the entry in the phase space file. +Full listing: ```python source.PDGCode_key = "PDGCode" source.particle = None ``` +The PDGCode is defined by the particle data group (see https://pdg.lbl.gov/2023/web/viewer.html?file=../reviews/rpp2022-rev-monte-carlo-numbering.pdf). +Here is a short overview of common particle types and its corresponding PDG Code +``` +proton: 2212 +neutron: 2211 +electron: 11 +gamma: 22 +carbon ion C12: 1000060120 +``` + +The naming of the phsp file entries generated by e.g. a GATE phase space actor changed over time, most notably from GATE v9 to GATE v10. +Setting ```source.position_key = "PrePositionLocal"``` will cause the phsp source to look for particle positions in ```PrePositionLocal_X, PrePositionLocal_Y, PrePositionLocal_Z```. +```source.direction_key = "PreDirectionLocal"``` will do the corresponding for the particle direction vector components in ```PreDirectionLocal_X, PreDirectionLocal_y, PreDirectionLocal_Z```. + +Alternatively, it is possible to directly set the individual components: +``` +source.position_key = None"PrePositionLocal" +source.position_key_x = Position_X" +source.position_key_y = Position_X +source.position_key_z = Position_X +source.direction_key = None +source.direction_key_x = Direction_X +source.direction_key_y = Direction_X +source.direction_key_z = Direction_X +``` + +In case of simulating a realistic particle mix, for example the output after a linac, a phsp file could contain a mixture of particles. +Typically, one would be interested in simulating a given number of primary particles (e.g. protons), simulating, but not counting as secondary particles (e.g. secondary electrons) in the number of particles to simulate. +This can be acchieved by setting ```generate_until_next_primary = True```. Furthermore, the PDG code of the primary particle needs to be specified, as well as a +lower energy threshold (in order to identify secondary particles of the same type as the primary particle). +The example below will consider protons above 90 MeV as primary particles. Every primary particle found in the phsp file will increase the number of particles simulated counter, while secondary particles (e.g. all other particles in the phsp file) +will be simulated, but not be considered in (e.g. not increasing) the number of particles simulated. +``` +source.generate_until_next_primary = True +source.primary_lower_energy_threshold = 90.0 * MeV +source.primary_PDGCode = 2212 +``` + For multithread: you need to indicate the ```entry_start``` for all threads, as an array, so that each thread starts in the phsp file at a different position. This done for example as follows (see ```test019_linac_phsp_source_MT.py```). Warning, if the phsp reach its end, it will cycle and start back at the beginning. ```python @@ -166,7 +207,7 @@ nb_of_threads = 4 source.entry_start = [total_nb_of_particle * p for p in range(nb_of_threads)] ``` -See all test019 as examples. +See all test019 and test060 as examples. ### GAN sources (Generative Adversarial Network) From e1c61038ad32a472047701a3900d70318946905c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 4 Dec 2023 20:23:06 +0000 Subject: [PATCH 22/46] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/pre-commit/mirrors-clang-format: v17.0.5 → v17.0.6](https://github.com/pre-commit/mirrors-clang-format/compare/v17.0.5...v17.0.6) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c43f13afd..cd9239084 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -9,7 +9,7 @@ repos: hooks: - id: black - repo: https://github.com/pre-commit/mirrors-clang-format - rev: v17.0.5 + rev: v17.0.6 hooks: - id: clang-format ci: From b972cb16df3895e175470ba71295a378e2f73f4f Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Fri, 8 Dec 2023 15:39:24 +0100 Subject: [PATCH 23/46] Update tests from PR#205 I merged the PR because the tests were ok, but they were old and did not take into account some, modifications --- opengate/tests/src/test060_phsp_source_helpers.py | 3 ++- opengate/tests/src/test061_TPPhsSource_helpers.py | 15 ++++++--------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/opengate/tests/src/test060_phsp_source_helpers.py b/opengate/tests/src/test060_phsp_source_helpers.py index a0b24a7d0..98a9dd651 100644 --- a/opengate/tests/src/test060_phsp_source_helpers.py +++ b/opengate/tests/src/test060_phsp_source_helpers.py @@ -397,7 +397,8 @@ def test_source_untilPrimary( source.primary_PDGCode = 2212 print(source) - output = sim.start() + sim.run() + output = sim.output def get_first_entry_of_key( diff --git a/opengate/tests/src/test061_TPPhsSource_helpers.py b/opengate/tests/src/test061_TPPhsSource_helpers.py index ae468d78c..d99e6d8d4 100644 --- a/opengate/tests/src/test061_TPPhsSource_helpers.py +++ b/opengate/tests/src/test061_TPPhsSource_helpers.py @@ -94,9 +94,7 @@ def create_test_Phs( f.particle = particle ta1.filters.append(f) - phys = sim.get_physics_user_info() - # ~ phys.physics_list_name = "FTFP_BERT" - phys.physics_list_name = "QGSP_BIC_EMZ" + sim.physics_manager.physics_list_name = "QGSP_BIC_EMZ" ########################################################################################## # Source @@ -113,8 +111,8 @@ def create_test_Phs( source.energy.mono = 150 * MeV source.n = number_of_particles - # output = sim.run() - output = sim.start(start_new_process=True) + sim.run() + output = sim.output def create_PhS_withoutSource( @@ -175,9 +173,7 @@ def create_PhS_withoutSource( ta1.output = phs_name ta1.debug = False - phys = sim.get_physics_user_info() - # ~ phys.physics_list_name = "FTFP_BERT" - phys.physics_list_name = "QGSP_BIC_EMZ" + sim.physics_manager.physics_list_name = "QGSP_BIC_EMZ" # ########################################################################################## # # Source @@ -229,7 +225,8 @@ def test_source_rotation_A( # depending on the rotation of the gantry, the rotation of the phase space to catch the particles is different plane.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() - output = sim.start() + sim.run() + output = sim.output def get_first_entry_of_key( From 86729fb924c3bf705785b0efbb2aae1f2c11511a Mon Sep 17 00:00:00 2001 From: majacquet Date: Mon, 4 Sep 2023 14:59:51 +0200 Subject: [PATCH 24/46] Update of Vsource and PhaseSpaceSource to include the activity parameter --- .../opengate_lib/GateGenericSource.cpp | 23 +++++++------- .../opengate_lib/GateGenericSource.h | 12 ++++---- .../opengate_lib/GatePhaseSpaceSource.cpp | 21 ++++++++++--- .../opengate_lib/GatePhaseSpaceSource.h | 4 ++- .../opengate_lib/GateVSource.cpp | 30 +++++++++++++++++++ core/opengate_core/opengate_lib/GateVSource.h | 20 +++++++++---- opengate/sources/phspsources.py | 2 ++ 7 files changed, 86 insertions(+), 26 deletions(-) diff --git a/core/opengate_core/opengate_lib/GateGenericSource.cpp b/core/opengate_core/opengate_lib/GateGenericSource.cpp index 374950c1b..43c1db60e 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.cpp +++ b/core/opengate_core/opengate_lib/GateGenericSource.cpp @@ -13,9 +13,12 @@ #include GateGenericSource::GateGenericSource() : GateVSource() { - fNumberOfGeneratedEvents = 0; + /*fNumberOfGeneratedEvents = 0; fMaxN = 0; fActivity = 0; + fHalfLife = -1; + fLambda = -1; + */ fInitGenericIon = false; fA = 0; fZ = 0; @@ -23,8 +26,6 @@ GateGenericSource::GateGenericSource() : GateVSource() { fInitConfine = false; fWeight = -1; fWeightSigma = -1; - fHalfLife = -1; - fDecayConstant = -1; fTotalSkippedEvents = 0; fCurrentSkippedEvents = 0; fTotalZeroEvents = 0; @@ -75,14 +76,15 @@ void GateGenericSource::InitializeUserInfo(py::dict &user_info) { CreateSPS(); // get user info about activity or nb of events + /* fMaxN = DictGetInt(user_info, "n"); fActivity = DictGetDouble(user_info, "activity"); fInitialActivity = fActivity; // half life ? fHalfLife = DictGetDouble(user_info, "half_life"); - fDecayConstant = log(2) / fHalfLife; - fUserParticleLifeTime = DictGetDouble(user_info, "user_particle_life_time"); + fLambda = log(2) / fHalfLife; + */ // weight fWeight = DictGetDouble(user_info, "weight"); @@ -108,9 +110,7 @@ void GateGenericSource::InitializeUserInfo(py::dict &user_info) { void GateGenericSource::UpdateActivity(double time) { if (!fTAC_Times.empty()) return UpdateActivityWithTAC(time); - if (fHalfLife <= 0) - return; - fActivity = fInitialActivity * exp(-fDecayConstant * (time - fStartTime)); + GateVSource::UpdateActivity(time); } void GateGenericSource::UpdateActivityWithTAC(double time) { @@ -150,7 +150,7 @@ double GateGenericSource::PrepareNextTime(double current_simulation_time) { fCurrentZeroEvents = 0; auto cse = fCurrentSkippedEvents; fCurrentSkippedEvents = 0; - + // if MaxN is below zero, we check the time if (fMaxN <= 0) { if (fEffectiveEventTime < fStartTime) @@ -159,14 +159,15 @@ double GateGenericSource::PrepareNextTime(double current_simulation_time) { return -1; // get next time according to current fActivity - double next_time = - fEffectiveEventTime - log(G4UniformRand()) * (1.0 / fActivity); + double next_time = CalcNextTime(fEffectiveEventTime); if (next_time >= fEndTime) return -1; return next_time; } // check according to t MaxN + //std::cout<= fMaxN) { return -1; } diff --git a/core/opengate_core/opengate_lib/GateGenericSource.h b/core/opengate_core/opengate_lib/GateGenericSource.h index 3d1e65b72..3c7ff0234 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.h +++ b/core/opengate_core/opengate_lib/GateGenericSource.h @@ -26,7 +26,7 @@ class GateGenericSource : public GateVSource { void InitializeUserInfo(py::dict &user_info) override; - double PrepareNextTime(double current_simulation_time) override; + double PrepareNextTime(double current_simulation_time); void PrepareNextRun() override; @@ -34,7 +34,7 @@ class GateGenericSource : public GateVSource { /// Current number of simulated events in this source /// (do not include skipped events) - unsigned long fNumberOfGeneratedEvents; + //unsigned long fNumberOfGeneratedEvents; /// Count the number of skipped events /// (e.g. Acceptance Angle or in GANSource) @@ -51,15 +51,17 @@ class GateGenericSource : public GateVSource { const std::vector &activities); protected: - unsigned long fMaxN; + //unsigned long fMaxN; // We cannot not use a std::unique_ptr // (or maybe by controlling the deletion during the CleanWorkerThread ?) GateSingleParticleSource *fSPS; - + /* double fActivity; double fInitialActivity; double fHalfLife; - double fDecayConstant; + double fLambda; + */ + G4ParticleDefinition *fParticleDefinition; double fEffectiveEventTime; double fUserParticleLifeTime; diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp index b888601fa..683809dba 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp @@ -14,7 +14,8 @@ GatePhaseSpaceSource::GatePhaseSpaceSource() : GateVSource() { fCharge = 0; fMass = 0; - fMaxN = 0; + fCurrentBatchSize = 0; + //fMaxN = 0; fGlobalFag = false; } @@ -31,7 +32,7 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) { GateVSource::InitializeUserInfo(user_info); // Number of events to generate - fMaxN = DictGetInt(user_info, "n"); + //fMaxN = DictGetInt(user_info, "n"); // global (world) or local (mother volume) coordinate system fGlobalFag = DictGetBool(user_info, "global_flag"); @@ -69,8 +70,20 @@ void GatePhaseSpaceSource::PrepareNextRun() { double GatePhaseSpaceSource::PrepareNextTime(double current_simulation_time) { // check according to t MaxN - auto &l = fThreadLocalDataPhsp.Get(); - if (l.fNumberOfGeneratedEvents >= fMaxN) { + + UpdateActivity(current_simulation_time); + if (fMaxN <= 0) { + if (current_simulation_time < fStartTime) + return fStartTime; + if (current_simulation_time >= fEndTime) + return -1; + + double next_time = CalcNextTime(current_simulation_time); + if (next_time >= fEndTime) + return -1; + return next_time; + } + if (fNumberOfGeneratedEvents >= fMaxN) { return -1; } return fStartTime; // FIXME timing ? diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h index b5998b90e..0b61d3e99 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h @@ -58,7 +58,9 @@ class GatePhaseSpaceSource : public GateVSource { bool fGlobalFag; bool fUseParticleTypeFromFile; - unsigned long fMaxN; + //unsigned long fMaxN; + long fNumberOfGeneratedEvents; + size_t fCurrentBatchSize; void SetPDGCodeBatch(const py::array_t &fPDGCode) const; diff --git a/core/opengate_core/opengate_lib/GateVSource.cpp b/core/opengate_core/opengate_lib/GateVSource.cpp index ba139cb29..3101740cd 100644 --- a/core/opengate_core/opengate_lib/GateVSource.cpp +++ b/core/opengate_core/opengate_lib/GateVSource.cpp @@ -10,6 +10,7 @@ #include "GateHelpers.h" #include "GateHelpersDict.h" #include "GateHelpersGeometry.h" +#include "G4RandomTools.hh" G4Mutex SourceOrientationMutex = G4MUTEX_INITIALIZER; @@ -20,6 +21,11 @@ GateVSource::GateVSource() { fMother = ""; fLocalTranslation = G4ThreeVector(); fLocalRotation = G4RotationMatrix(); + fNumberOfGeneratedEvents = 0; + fMaxN = 0; + fActivity = 0; + fHalfLife = -1; + fLambda = -1; } GateVSource::~GateVSource() {} @@ -30,7 +36,31 @@ void GateVSource::InitializeUserInfo(py::dict &user_info) { fStartTime = DictGetDouble(user_info, "start_time"); fEndTime = DictGetDouble(user_info, "end_time"); fMother = DictGetStr(user_info, "mother"); + + // get user info about activity or nb of events + fMaxN = DictGetInt(user_info, "n"); + fActivity = DictGetDouble(user_info, "activity"); + fInitialActivity = fActivity; + + // half life ? + fHalfLife = DictGetDouble(user_info, "half_life"); + fLambda = log(2) / fHalfLife; + +} + + +void GateVSource::UpdateActivity(double time) { + if (fHalfLife <= 0) + return; + fActivity = fInitialActivity * exp(-fLambda * (time - fStartTime)); } + +double GateVSource::CalcNextTime(double current_simulation_time) { + double next_time = + current_simulation_time - log(G4UniformRand()) * (1.0 / fActivity); + return next_time; +} + void GateVSource::PrepareNextRun() { SetOrientationAccordingToMotherVolume(); } diff --git a/core/opengate_core/opengate_lib/GateVSource.h b/core/opengate_core/opengate_lib/GateVSource.h index ae2631963..dc0d9414b 100644 --- a/core/opengate_core/opengate_lib/GateVSource.h +++ b/core/opengate_core/opengate_lib/GateVSource.h @@ -28,6 +28,10 @@ class GateVSource { // Called at initialisation to set the source properties from a single dict virtual void InitializeUserInfo(py::dict &user_info); + + virtual void UpdateActivity(double time); + + double CalcNextTime(double current_simulation_time); virtual void PrepareNextRun(); @@ -40,6 +44,8 @@ class GateVSource { std::string fName; double fStartTime; double fEndTime; + unsigned long fNumberOfGeneratedEvents; + std::string fMother; std::vector fTranslations; std::vector fRotations; @@ -47,11 +53,15 @@ class GateVSource { G4ThreeVector fLocalTranslation; G4RotationMatrix fLocalRotation; - struct threadLocalT { - G4ThreeVector fGlobalTranslation; - G4RotationMatrix fGlobalRotation; - }; - G4Cache fThreadLocalData; + G4ThreeVector fGlobalTranslation; + G4RotationMatrix fGlobalRotation; + +protected: + unsigned long fMaxN; + double fActivity; + double fInitialActivity; + double fHalfLife; + double fLambda; }; #endif // GateVSource_h diff --git a/opengate/sources/phspsources.py b/opengate/sources/phspsources.py index b4fdc628a..a249b7d9d 100644 --- a/opengate/sources/phspsources.py +++ b/opengate/sources/phspsources.py @@ -234,6 +234,8 @@ def set_default_user_info(user_info): # initial user info user_info.phsp_file = None user_info.n = 1 + user_info.activity = 0 + user_info.half_life = -1 # negative value is not half_life user_info.particle = "" # FIXME later as key user_info.entry_start = 0 # if a particle name is supplied, the particle type is set to it From 9e89aebbf6a54f1b7c24bba1bcb704c84e2cfe93 Mon Sep 17 00:00:00 2001 From: majacquet Date: Fri, 22 Sep 2023 13:25:40 +0200 Subject: [PATCH 25/46] Implementation of activity for phaseSpace and merge with master --- core/opengate_core/opengate_lib/GateGenericSource.cpp | 1 + core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp | 3 ++- core/opengate_core/opengate_lib/GateVSource.cpp | 9 ++++++--- core/opengate_core/opengate_lib/GateVSource.h | 9 ++++++++- 4 files changed, 17 insertions(+), 5 deletions(-) diff --git a/core/opengate_core/opengate_lib/GateGenericSource.cpp b/core/opengate_core/opengate_lib/GateGenericSource.cpp index 43c1db60e..fc6281c75 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.cpp +++ b/core/opengate_core/opengate_lib/GateGenericSource.cpp @@ -89,6 +89,7 @@ void GateGenericSource::InitializeUserInfo(py::dict &user_info) { // weight fWeight = DictGetDouble(user_info, "weight"); fWeightSigma = DictGetDouble(user_info, "weight_sigma"); + fUserParticleLifeTime = DictGetDouble(user_info, "user_particle_life_time"); // get the user info for the particle InitializeParticle(user_info); diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp index 683809dba..c3b5cf8d8 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp @@ -83,7 +83,8 @@ double GatePhaseSpaceSource::PrepareNextTime(double current_simulation_time) { return -1; return next_time; } - if (fNumberOfGeneratedEvents >= fMaxN) { + auto &l = fThreadLocalDataPhsp.Get(); + if (l.fNumberOfGeneratedEvents >= fMaxN) { return -1; } return fStartTime; // FIXME timing ? diff --git a/core/opengate_core/opengate_lib/GateVSource.cpp b/core/opengate_core/opengate_lib/GateVSource.cpp index 3101740cd..00d41e457 100644 --- a/core/opengate_core/opengate_lib/GateVSource.cpp +++ b/core/opengate_core/opengate_lib/GateVSource.cpp @@ -25,7 +25,9 @@ GateVSource::GateVSource() { fMaxN = 0; fActivity = 0; fHalfLife = -1; - fLambda = -1; + fDecayConstant = -1; + + } GateVSource::~GateVSource() {} @@ -37,6 +39,7 @@ void GateVSource::InitializeUserInfo(py::dict &user_info) { fEndTime = DictGetDouble(user_info, "end_time"); fMother = DictGetStr(user_info, "mother"); + // get user info about activity or nb of events fMaxN = DictGetInt(user_info, "n"); fActivity = DictGetDouble(user_info, "activity"); @@ -44,7 +47,7 @@ void GateVSource::InitializeUserInfo(py::dict &user_info) { // half life ? fHalfLife = DictGetDouble(user_info, "half_life"); - fLambda = log(2) / fHalfLife; + fDecayConstant = log(2) / fHalfLife; } @@ -52,7 +55,7 @@ void GateVSource::InitializeUserInfo(py::dict &user_info) { void GateVSource::UpdateActivity(double time) { if (fHalfLife <= 0) return; - fActivity = fInitialActivity * exp(-fLambda * (time - fStartTime)); + fActivity = fInitialActivity * exp(-fDecayConstant * (time - fStartTime)); } double GateVSource::CalcNextTime(double current_simulation_time) { diff --git a/core/opengate_core/opengate_lib/GateVSource.h b/core/opengate_core/opengate_lib/GateVSource.h index dc0d9414b..a24b5fe05 100644 --- a/core/opengate_core/opengate_lib/GateVSource.h +++ b/core/opengate_core/opengate_lib/GateVSource.h @@ -61,7 +61,14 @@ class GateVSource { double fActivity; double fInitialActivity; double fHalfLife; - double fLambda; + double fDecayConstant; + + struct threadLocalT { + G4ThreeVector fGlobalTranslation; + G4RotationMatrix fGlobalRotation; + }; + G4Cache fThreadLocalData; + }; #endif // GateVSource_h From d99f9691551eab9741c1f43c40fde676a710c2a3 Mon Sep 17 00:00:00 2001 From: majacquet Date: Wed, 15 Nov 2023 17:46:36 +0100 Subject: [PATCH 26/46] Update of source activity for PhaseSpace and creation of a test --- opengate/sources/phspsources.py | 12 +- .../src/test065_PhaseSpaceSource_Activity.py | 118 ++++++++++++++++++ 2 files changed, 124 insertions(+), 6 deletions(-) create mode 100644 opengate/tests/src/test065_PhaseSpaceSource_Activity.py diff --git a/opengate/sources/phspsources.py b/opengate/sources/phspsources.py index a249b7d9d..e8de0061c 100644 --- a/opengate/sources/phspsources.py +++ b/opengate/sources/phspsources.py @@ -156,12 +156,12 @@ def generate(self, source): # if translate_position is set to True, the position # supplied will be added to the phsp file position if ui.translate_position: - x = batch[ui.position_key_x] + ui.position.translation[0] - y = batch[ui.position_key_y] + ui.position.translation[1] - z = batch[ui.position_key_z] + ui.position.translation[2] - source.SetPositionXBatch(x) - source.SetPositionYBatch(y) - source.SetPositionZBatch(z) + batch[ui.position_key_x] += ui.position.translation[0] + batch[ui.position_key_y] += ui.position.translation[1] + batch[ui.position_key_z] += ui.position.translation[2] + source.SetPositionXBatch(batch[ui.position_key_x]) + source.SetPositionYBatch(batch[ui.position_key_y]) + source.SetPositionZBatch(batch[ui.position_key_z]) else: source.SetPositionXBatch(batch[ui.position_key_x]) source.SetPositionYBatch(batch[ui.position_key_y]) diff --git a/opengate/tests/src/test065_PhaseSpaceSource_Activity.py b/opengate/tests/src/test065_PhaseSpaceSource_Activity.py new file mode 100644 index 000000000..16e09fdc2 --- /dev/null +++ b/opengate/tests/src/test065_PhaseSpaceSource_Activity.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import uproot +import numpy as np + + + + +def validation_test(n,n_measured): + n_part_minus = n - 3*np.sqrt(n) + n_part_max = n + 3 * np.sqrt(n) + if n_part_minus <=n_measured <= n_part_max: + return True + else: + return False + + +if __name__ == "__main__": + paths = gate.get_default_test_paths(__file__) + output_path = "/home/mjacquet/Software/gatePython/opengate/opengate/tests/output/test019/" + + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + # ui.visu = True + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + # ui.running_verbose_level = gate.EVENT + ui.number_of_threads = 1 + ui.random_seed = "auto" + + # units + m = gate.g4_units("m") + km = gate.g4_units("km") + mm = gate.g4_units("mm") + cm = gate.g4_units("cm") + nm = gate.g4_units("nm") + Bq = gate.g4_units("Bq") + MeV = gate.g4_units("MeV") + keV = gate.g4_units("keV") + gcm3 = gate.g4_units("g/cm3") + + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 2 * m] + world.material = "G4_Galactic" + + # source_1 + nb_part = 1000 + plane_1 = sim.add_volume("Box", "plan_1") + plane_1.material = "G4_Galactic" + plane_1.mother = world.name + plane_1.size = [1 * m, 1 * m, 1 * nm] + plane_1.color = [0.8, 0.2, 0.1, 1] + + source_1 = sim.add_source("PhaseSpaceSource", "phsp_source_global_1") + source_1.mother = plane_1.name + source_1.phsp_file = paths.output / "test019/test019_hits_phsp_source_global.root" + source_1.position_key = "PrePosition" + source_1.direction_key = "PreDirection" + source_1.weight_key = "Weight" + source_1.global_flag = False + source_1.particle = "gamma" + source_1.batch_size = 100 + source_1.override_position = True + source_1.position.translation = [0, 0, 300] + source_1.activity = nb_part*Bq + + + # add phase space plan 1 + + phsp_plane_1 = sim.add_volume("Box", "phase_space_plane_1") + phsp_plane_1.mother = world.name + phsp_plane_1.material = "G4_Galactic" + phsp_plane_1.size = [1 * m, 1 * m, 1 * nm] + phsp_plane_1.translation = [0, 0, - 1000 *nm] + phsp_plane_1.color = [1, 0, 0, 1] # red + + phsp_actor = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp_actor.mother = phsp_plane_1.name + phsp_actor.attributes = [ + "KineticEnergy", + ] + + + phsp_actor.output = paths.output / "test065_output_data.root" + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + p = sim.get_physics_user_info() + p.physics_list_name = "G4EmStandardPhysics_option3" + p.enable_decay = False + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1 * mm + sim.physics_manager.global_production_cuts.positron = 1 * mm + # + # # go ! + output = sim.start() + + + # + # # print results + stats = sim.output.get_actor("Stats") + h = output.get_actor("PhaseSpace") + print(stats) + # + f_phsp = uproot.open(paths.output / "test065_output_data.root") + arr = f_phsp["PhaseSpace"].arrays() + print("Number of detected events :", len(arr)) + print("Number of expected events :", nb_part, "+- 3*" + str(int(np.round(np.sqrt(nb_part))))) + + is_ok = validation_test(nb_part,len(arr)) + gate.test_ok(is_ok) From f36bb2b9b4af84f512f165520fddacd8a18f7152 Mon Sep 17 00:00:00 2001 From: majacquet Date: Wed, 15 Nov 2023 17:47:22 +0100 Subject: [PATCH 27/46] Update --- opengate/tests/src/test065_PhaseSpaceSource_Activity.py | 1 - 1 file changed, 1 deletion(-) diff --git a/opengate/tests/src/test065_PhaseSpaceSource_Activity.py b/opengate/tests/src/test065_PhaseSpaceSource_Activity.py index 16e09fdc2..e81de212a 100644 --- a/opengate/tests/src/test065_PhaseSpaceSource_Activity.py +++ b/opengate/tests/src/test065_PhaseSpaceSource_Activity.py @@ -19,7 +19,6 @@ def validation_test(n,n_measured): if __name__ == "__main__": paths = gate.get_default_test_paths(__file__) - output_path = "/home/mjacquet/Software/gatePython/opengate/opengate/tests/output/test019/" # create the simulation sim = gate.Simulation() From 8c087dfddd8375787d5676bbdb849a10dee4c6e1 Mon Sep 17 00:00:00 2001 From: majacquet Date: Mon, 20 Nov 2023 14:17:11 +0100 Subject: [PATCH 28/46] Correction of the close geometry boolean on the MotionVolumeActor --- core/opengate_core/opengate_lib/GateMotionVolumeActor.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/opengate_core/opengate_lib/GateMotionVolumeActor.cpp b/core/opengate_core/opengate_lib/GateMotionVolumeActor.cpp index 220f1c48a..bdd5e2363 100644 --- a/core/opengate_core/opengate_lib/GateMotionVolumeActor.cpp +++ b/core/opengate_core/opengate_lib/GateMotionVolumeActor.cpp @@ -63,7 +63,7 @@ void GateMotionVolumeActor::MoveGeometry(int run_id) { rot->set(r.rep3x3()); // close the geometry manager - gm->CloseGeometry(false, false, pv); + gm->CloseGeometry(true, false, pv); // G4RunManager::GetRunManager()->GeometryHasBeenModified(true); } From a69e78793d777cc286c08560c447ce18aa7ac3a2 Mon Sep 17 00:00:00 2001 From: majacquet Date: Mon, 20 Nov 2023 16:18:09 +0100 Subject: [PATCH 29/46] Update phspSource_Activity after rebasing Rename test Forgot to change phspsource.ui.n = 0 by default during the rebasing lint --- opengate/sources/phspsources.py | 4 +- ...vity.py => test066_phspsource_activity.py} | 58 ++++++++++--------- 2 files changed, 32 insertions(+), 30 deletions(-) rename opengate/tests/src/{test065_PhaseSpaceSource_Activity.py => test066_phspsource_activity.py} (70%) diff --git a/opengate/sources/phspsources.py b/opengate/sources/phspsources.py index e8de0061c..2c1e49fde 100644 --- a/opengate/sources/phspsources.py +++ b/opengate/sources/phspsources.py @@ -233,9 +233,9 @@ def set_default_user_info(user_info): SourceBase.set_default_user_info(user_info) # initial user info user_info.phsp_file = None - user_info.n = 1 + user_info.n = 0 user_info.activity = 0 - user_info.half_life = -1 # negative value is not half_life + user_info.half_life = -1 # negative value is not half_life user_info.particle = "" # FIXME later as key user_info.entry_start = 0 # if a particle name is supplied, the particle type is set to it diff --git a/opengate/tests/src/test065_PhaseSpaceSource_Activity.py b/opengate/tests/src/test066_phspsource_activity.py similarity index 70% rename from opengate/tests/src/test065_PhaseSpaceSource_Activity.py rename to opengate/tests/src/test066_phspsource_activity.py index e81de212a..979c7c012 100644 --- a/opengate/tests/src/test065_PhaseSpaceSource_Activity.py +++ b/opengate/tests/src/test066_phspsource_activity.py @@ -4,21 +4,22 @@ import opengate as gate import uproot import numpy as np +from opengate.tests import utility - - -def validation_test(n,n_measured): - n_part_minus = n - 3*np.sqrt(n) +def validation_test(n, n_measured): + n_part_minus = n - 3 * np.sqrt(n) n_part_max = n + 3 * np.sqrt(n) - if n_part_minus <=n_measured <= n_part_max: + if n_part_minus <= n_measured <= n_part_max: return True else: return False if __name__ == "__main__": - paths = gate.get_default_test_paths(__file__) + paths = utility.get_default_test_paths( + __file__, "test066_PhsSource_Activity", output_folder="test066" + ) # create the simulation sim = gate.Simulation() @@ -34,15 +35,15 @@ def validation_test(n,n_measured): ui.random_seed = "auto" # units - m = gate.g4_units("m") - km = gate.g4_units("km") - mm = gate.g4_units("mm") - cm = gate.g4_units("cm") - nm = gate.g4_units("nm") - Bq = gate.g4_units("Bq") - MeV = gate.g4_units("MeV") - keV = gate.g4_units("keV") - gcm3 = gate.g4_units("g/cm3") + m = gate.g4_units.m + km = gate.g4_units.km + mm = gate.g4_units.mm + cm = gate.g4_units.cm + nm = gate.g4_units.nm + Bq = gate.g4_units.Bq + MeV = gate.g4_units.MeV + keV = gate.g4_units.keV + gcm3 = gate.g4_units.g / gate.g4_units.cm3 # adapt world size world = sim.world @@ -59,7 +60,7 @@ def validation_test(n,n_measured): source_1 = sim.add_source("PhaseSpaceSource", "phsp_source_global_1") source_1.mother = plane_1.name - source_1.phsp_file = paths.output / "test019/test019_hits_phsp_source_global.root" + source_1.phsp_file = paths.output_ref / ".." / "test019" / "test019_hits.root" source_1.position_key = "PrePosition" source_1.direction_key = "PreDirection" source_1.weight_key = "Weight" @@ -68,8 +69,7 @@ def validation_test(n,n_measured): source_1.batch_size = 100 source_1.override_position = True source_1.position.translation = [0, 0, 300] - source_1.activity = nb_part*Bq - + source_1.activity = nb_part * Bq # add phase space plan 1 @@ -77,7 +77,7 @@ def validation_test(n,n_measured): phsp_plane_1.mother = world.name phsp_plane_1.material = "G4_Galactic" phsp_plane_1.size = [1 * m, 1 * m, 1 * nm] - phsp_plane_1.translation = [0, 0, - 1000 *nm] + phsp_plane_1.translation = [0, 0, -1000 * nm] phsp_plane_1.color = [1, 0, 0, 1] # red phsp_actor = sim.add_actor("PhaseSpaceActor", "PhaseSpace") @@ -85,9 +85,8 @@ def validation_test(n,n_measured): phsp_actor.attributes = [ "KineticEnergy", ] - - - phsp_actor.output = paths.output / "test065_output_data.root" + + phsp_actor.output = paths.output / "test066_output_data.root" s = sim.add_actor("SimulationStatisticsActor", "Stats") s.track_types_flag = True @@ -101,17 +100,20 @@ def validation_test(n,n_measured): # # go ! output = sim.start() - # # # print results stats = sim.output.get_actor("Stats") h = output.get_actor("PhaseSpace") print(stats) # - f_phsp = uproot.open(paths.output / "test065_output_data.root") + f_phsp = uproot.open(paths.output / "test066_output_data.root") arr = f_phsp["PhaseSpace"].arrays() print("Number of detected events :", len(arr)) - print("Number of expected events :", nb_part, "+- 3*" + str(int(np.round(np.sqrt(nb_part))))) - - is_ok = validation_test(nb_part,len(arr)) - gate.test_ok(is_ok) + print( + "Number of expected events :", + nb_part, + "+- 3*" + str(int(np.round(np.sqrt(nb_part)))), + ) + + is_ok = validation_test(nb_part, len(arr)) + utility.test_ok(is_ok) From 3c8b5dae34a082988ffe1ca3613aafb6b4de92a9 Mon Sep 17 00:00:00 2001 From: majacquet Date: Tue, 21 Nov 2023 11:17:39 +0100 Subject: [PATCH 30/46] Move the parameters n, activity and half_life to SourceBase The parameters were moved previously to VSource.h but were not moved inside the python part --- opengate/sources/generic.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/opengate/sources/generic.py b/opengate/sources/generic.py index 0d384cbfd..d1d7d0b8b 100644 --- a/opengate/sources/generic.py +++ b/opengate/sources/generic.py @@ -219,6 +219,9 @@ def set_default_user_info(user_info): user_info.mother = __world_name__ user_info.start_time = None user_info.end_time = None + user_info.n = 0 + user_info.activity = 0 + user_info.half_life = -1 # negative value is no half_life def __init__(self, user_info): # type_name MUST be defined in class that inherit from SourceBase @@ -305,11 +308,8 @@ def set_default_user_info(user_info): # initial user info user_info.particle = "gamma" user_info.ion = Box() - user_info.n = 0 - user_info.activity = 0 user_info.weight = -1 user_info.weight_sigma = -1 - user_info.half_life = -1 # negative value is no half_life user_info.user_particle_life_time = -1 # negative means : by default user_info.tac_times = None user_info.tac_activities = None @@ -497,7 +497,6 @@ class TemplateSource(SourceBase): def set_default_user_info(user_info): SourceBase.set_default_user_info(user_info) # initial user info - user_info.n = 0 user_info.float_value = None user_info.vector_value = None From 3a9c3f124f1ae05ee1e151025b880e2232e554b3 Mon Sep 17 00:00:00 2001 From: majacquet Date: Thu, 7 Dec 2023 16:29:38 +0100 Subject: [PATCH 31/46] Correction of a matrix rotation error to position the source of a multi run simulation. Implementation, for generic sources, of the initial momentum rotation in accordance with the motion actor rotation --- .../opengate_lib/GateGANSource.cpp | 15 +- .../opengate_lib/GateGANSource.h | 2 +- .../opengate_lib/GateGenericSource.cpp | 29 ++- .../opengate_lib/GateGenericSource.h | 4 + .../opengate_lib/GateVoxelsSource.cpp | 16 +- opengate/tests/src/test068_rotation_source.py | 171 ++++++++++++++++++ 6 files changed, 231 insertions(+), 6 deletions(-) create mode 100644 opengate/tests/src/test068_rotation_source.py diff --git a/core/opengate_core/opengate_lib/GateGANSource.cpp b/core/opengate_core/opengate_lib/GateGANSource.cpp index 8e4e193b2..3786e504b 100644 --- a/core/opengate_core/opengate_lib/GateGANSource.cpp +++ b/core/opengate_core/opengate_lib/GateGANSource.cpp @@ -88,7 +88,20 @@ void GateGANSource::PrepareNextRun() { // (no need to update th fSPS pos in GateGenericSource) // GateVSource::PrepareNextRun(); // FIXME remove this function ? - GateGenericSource::PrepareNextRun(); + // GateGenericSource::PrepareNextRun(); + GateVSource::PrepareNextRun(); + // This global transformation is given to the SPS that will + // generate particles in the correct coordinate system + auto &l = fThreadLocalData.Get(); + auto *pos = fSPS->GetPosDist(); + pos->SetCentreCoords(l.fGlobalTranslation); + + // orientation according to mother volume + auto rotation = l.fGlobalRotation; + G4ThreeVector r1(rotation(0, 0), rotation(1, 0), rotation(2, 0)); + G4ThreeVector r2(rotation(0, 1), rotation(1, 1), rotation(2, 1)); + pos->SetPosRot1(r1); + pos->SetPosRot2(r2); } void GateGANSource::SetGeneratorFunction(ParticleGeneratorType &f) { diff --git a/core/opengate_core/opengate_lib/GateGANSource.h b/core/opengate_core/opengate_lib/GateGANSource.h index fd0fb768a..bb1be58a4 100644 --- a/core/opengate_core/opengate_lib/GateGANSource.h +++ b/core/opengate_core/opengate_lib/GateGANSource.h @@ -28,7 +28,7 @@ class GateGANSource : public GateGenericSource { void InitializeUserInfo(py::dict &user_info) override; - void PrepareNextRun() override; + void PrepareNextRun(); void GeneratePrimaries(G4Event *event, double current_simulation_time) override; diff --git a/core/opengate_core/opengate_lib/GateGenericSource.cpp b/core/opengate_core/opengate_lib/GateGenericSource.cpp index fc6281c75..4c8009dae 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.cpp +++ b/core/opengate_core/opengate_lib/GateGenericSource.cpp @@ -35,6 +35,9 @@ GateGenericSource::GateGenericSource() : GateVSource() { fParticleDefinition = nullptr; fEffectiveEventTime = -1; fEffectiveEventTime = -1; + + + } GateGenericSource::~GateGenericSource() { @@ -182,15 +185,32 @@ void GateGenericSource::PrepareNextRun() { // This global transformation is given to the SPS that will // generate particles in the correct coordinate system auto &l = fThreadLocalData.Get(); + //auto user_info_pos = py::dict(puser_info["position"]); + //auto pos_init = DictGetG4ThreeVector(user_info_pos, "translation"); auto *pos = fSPS->GetPosDist(); pos->SetCentreCoords(l.fGlobalTranslation); // orientation according to mother volume auto rotation = l.fGlobalRotation; - G4ThreeVector r1(rotation(0, 0), rotation(0, 1), rotation(0, 2)); - G4ThreeVector r2(rotation(1, 0), rotation(1, 1), rotation(1, 2)); + G4ThreeVector r1(rotation(0, 0), rotation(1, 0), rotation(2, 0)); + G4ThreeVector r2(rotation(0, 1), rotation(1, 1), rotation(2, 1)); pos->SetPosRot1(r1); pos->SetPosRot2(r2); + + + auto* ang = fSPS->GetAngDist(); + + if (fangType == "momentum") { + auto new_d = rotation * fInitializeMomentum; + ang->SetParticleMomentumDirection(new_d); + } + if (fangType == "focused") { + auto vec_f = fInitiliazeFocusPoint - fInitTranslation; + auto rot_f = rotation*vec_f; + auto new_f = rot_f + l.fGlobalTranslation; + ang->SetFocusPoint(new_f); + } + } void GateGenericSource::UpdateEffectiveEventTime( @@ -228,7 +248,6 @@ void GateGenericSource::GeneratePrimaries(G4Event *event, // (acceptance angle is included) fSPS->SetParticleTime(current_simulation_time); fSPS->GeneratePrimaryVertex(event); - // update the time according to skipped events fEffectiveEventTime = current_simulation_time; auto &l = fThreadLocalDataAA.Get(); @@ -302,6 +321,7 @@ void GateGenericSource::InitializePosition(py::dict puser_info) { std::vector l = {"sphere", "point", "box", "disc", "cylinder"}; CheckIsIn(pos_type, l); auto translation = DictGetG4ThreeVector(user_info, "translation"); + fInitTranslation = translation; if (pos_type == "point") { pos->SetPosDisType("Point"); } @@ -367,6 +387,7 @@ void GateGenericSource::InitializeDirection(py::dict puser_info) { auto user_info = py::dict(puser_info["direction"]); auto *ang = fSPS->GetAngDist(); auto ang_type = DictGetStr(user_info, "type"); + fangType = ang_type; std::vector ll = {"iso", "momentum", "focused", "beam2d"}; // FIXME check on py side ? CheckIsIn(ang_type, ll); @@ -384,11 +405,13 @@ void GateGenericSource::InitializeDirection(py::dict puser_info) { if (ang_type == "momentum") { ang->SetAngDistType("planar"); // FIXME really ?? auto d = DictGetG4ThreeVector(user_info, "momentum"); + fInitializeMomentum = d; ang->SetParticleMomentumDirection(d); } if (ang_type == "focused") { ang->SetAngDistType("focused"); auto f = DictGetG4ThreeVector(user_info, "focus_point"); + fInitiliazeFocusPoint = f; ang->SetFocusPoint(f); } if (ang_type == "beam2d") { diff --git a/core/opengate_core/opengate_lib/GateGenericSource.h b/core/opengate_core/opengate_lib/GateGenericSource.h index 3c7ff0234..9bdb3c3dd 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.h +++ b/core/opengate_core/opengate_lib/GateGenericSource.h @@ -63,6 +63,10 @@ class GateGenericSource : public GateVSource { */ G4ParticleDefinition *fParticleDefinition; + G4ThreeVector fInitializeMomentum; + G4ThreeVector fInitiliazeFocusPoint; + G4ThreeVector fInitTranslation; + G4String fangType; double fEffectiveEventTime; double fUserParticleLifeTime; diff --git a/core/opengate_core/opengate_lib/GateVoxelsSource.cpp b/core/opengate_core/opengate_lib/GateVoxelsSource.cpp index 3551d27d2..5e93acfc0 100644 --- a/core/opengate_core/opengate_lib/GateVoxelsSource.cpp +++ b/core/opengate_core/opengate_lib/GateVoxelsSource.cpp @@ -17,9 +17,23 @@ GateVoxelsSource::GateVoxelsSource() : GateGenericSource() { GateVoxelsSource::~GateVoxelsSource() {} void GateVoxelsSource::PrepareNextRun() { - GateGenericSource::PrepareNextRun(); + //GateGenericSource::PrepareNextRun(); // rotation and translation to apply, according to mother volume + GateVSource::PrepareNextRun(); + // This global transformation is given to the SPS that will + // generate particles in the correct coordinate system auto &l = fThreadLocalData.Get(); + auto *pos = fSPS->GetPosDist(); + pos->SetCentreCoords(l.fGlobalTranslation); + + // orientation according to mother volume + auto rotation = l.fGlobalRotation; + G4ThreeVector r1(rotation(0, 0), rotation(1, 0), rotation(2, 0)); + G4ThreeVector r2(rotation(0, 1), rotation(1, 1), rotation(2, 1)); + pos->SetPosRot1(r1); + pos->SetPosRot2(r2); + + //auto &l = fThreadLocalData.Get(); fVoxelPositionGenerator->fGlobalRotation = l.fGlobalRotation; fVoxelPositionGenerator->fGlobalTranslation = l.fGlobalTranslation; // the direction is 'isotropic' so we don't care about rotating the direction. diff --git a/opengate/tests/src/test068_rotation_source.py b/opengate/tests/src/test068_rotation_source.py new file mode 100644 index 000000000..a246000da --- /dev/null +++ b/opengate/tests/src/test068_rotation_source.py @@ -0,0 +1,171 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +from scipy.spatial.transform import Rotation +import numpy as np +import uproot + + + +#The test generates two differents generic sources, momentum and focused, which is attached to a plan. +#The plan is randomly rotated and we verify that the generated particles have a direction which is in accordance +# with the applied rotations and the transmissions. + + +def test068(tab,nb_run,nb_part,nb_source): + nb_event = len(tab) + nb_event_theo = nb_run * nb_part*nb_source + err = np.sqrt(nb_event_theo) + if nb_event_theo - 3*err < nb_event < nb_event_theo + 3*err : + return True + else : + return False + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__) + output_path = paths.output + + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + # ui.visu = True + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + # ui.running_verbose_level = gate.logger.EVENT + ui.number_of_threads = 1 + ui.random_seed = "auto" + + # units + m = gate.g4_units.m + km = gate.g4_units.km + mm = gate.g4_units.mm + cm = gate.g4_units.cm + nm = gate.g4_units.nm + Bq = gate.g4_units.Bq + MeV = gate.g4_units.MeV + keV = gate.g4_units.keV + sec = gate.g4_units.s + gcm3 = gate.g4_units["g/cm3"] + + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + world.material = "G4_Galactic" + + # nb of particles + n = 100 + + # Plan to attach the source + nb_source = 2 + plan = sim.add_volume("Box","source_plan") + plan.mother = world.name + plan.material = "G4_Galactic" + plan.size = [10*cm,10*cm,1*nm] + plan.translation = np.array([0*cm,0*cm,0*cm]) + + source1 = sim.add_source("GenericSource", "photon_source") + source1.particle = "gamma" + source1.position.type = "box" + source1.mother = plan.name + source1.position.size = [10 * cm, 10 * cm, 1 * nm] + source1.direction.type = "momentum" + # source1.direction.focus_point = [0*cm, 0*cm, -5 *cm] + source1.direction.momentum = [0,0,-1] + source1.energy.type = "mono" + source1.energy.mono = 1 * MeV + source1.activity = n*Bq / ui.number_of_threads + + source2 = sim.add_source("GenericSource", "photon_source_2") + source2.particle = "gamma" + source2.position.type = "box" + source2.mother = plan.name + source2.position.size = [10 * cm, 10 * cm, 1 * nm] + source2.direction.type = "focused" + source2.direction.focus_point = [0*cm, 0*cm, -5 *cm] + # source1.direction.momentum = [0, 0, -1] + source2.energy.type = "mono" + source2.energy.mono = 1 * MeV + source2.activity = n * Bq / ui.number_of_threads + + + # Phase Space + + phsp_plan = sim.add_volume("Box","phsp_plan") + phsp_plan.mother = world.name + phsp_plan.material = "G4_Galactic" + phsp_plan.size = [10*cm,10*cm,1*nm] + phsp_plan.translation = np.array([0*cm,0*cm,-5*cm]) + + phsp = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp.mother = phsp_plan.name + phsp.attributes = ["EventID"] + phsp.output = output_path / "test068.root" + + + #MOTION ACTOR for the source and actors + motion_source = sim.add_actor("MotionVolumeActor", "Move_source") + motion_source.mother = plan.name + motion_source.rotations = [] + motion_source.translations = [] + + motion_phsp = sim.add_actor("MotionVolumeActor", "Move_phsp") + motion_phsp.mother = phsp_plan.name + motion_phsp.rotations = [] + motion_phsp.translations = [] + sim.run_timing_intervals = [] + + for i in range(10): + x_translation = 20 * cm * np.random.random() + y_translation = 20 * cm * np.random.random() + z_translation = 20 * cm * np.random.random() + + x_rot = 360 * np.random.random() + y_rot = 360 * np.random.random() + z_rot = 360 * np.random.random() + + motion_source.translations.append([x_translation,y_translation,z_translation]) + rot = Rotation.from_euler("xyz",[x_rot,y_rot,z_rot], degrees=True).as_matrix() + vec_source_phsp = phsp_plan.translation - plan.translation + # vec_source_phsp = vec_source_phsp.reshape((len(vec_source_phsp),1)) + # print(vec_source_phsp) + translation_phsp = np.dot(rot,vec_source_phsp) + np.array([x_translation,y_translation,z_translation]) + motion_source.rotations.append(rot) + motion_phsp.rotations.append(rot) + motion_phsp.translations.append(translation_phsp) + sim.run_timing_intervals.append([i* sec, (i+1)* sec]) + # stat actor + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + + + + # Physic list and cuts + p = sim.get_physics_user_info() + p.physics_list_name = "G4EmStandardPhysics_option3" + p.enable_decay = False + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1 * mm + sim.physics_manager.global_production_cuts.positron = 1 * mm + + # go ! + + output = sim.run() + + + #Validation test + f_phsp = uproot.open(phsp.output) + arr = f_phsp["PhaseSpace"].arrays() + print("Number of detected events :", len(arr)) + nb_event = len(arr) + nb_part = n + nb_run = len(sim.run_timing_intervals) + err = np.sqrt(nb_part*nb_run*nb_source) + print("Number of expected events : [" +str(int(-3*err + nb_run*nb_part*nb_source)) + "," + str(int(+3*err + nb_run*nb_part*nb_source))+ "]") + is_ok = test068(arr,len(sim.run_timing_intervals),nb_part,2) + utility.test_ok(is_ok) + From 96e1f156976c81a9dd11c5ff720349eea58bf47e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 8 Dec 2023 12:56:49 +0000 Subject: [PATCH 32/46] [pre-commit.ci] Automatic python and c++ formatting --- .../opengate_lib/GateGANSource.cpp | 4 +- .../opengate_lib/GateGenericSource.cpp | 21 +++--- .../opengate_lib/GateGenericSource.h | 8 +-- .../opengate_lib/GatePhaseSpaceSource.cpp | 4 +- .../opengate_lib/GatePhaseSpaceSource.h | 2 +- .../opengate_lib/GateVSource.cpp | 12 +--- core/opengate_core/opengate_lib/GateVSource.h | 5 +- .../opengate_lib/GateVoxelsSource.cpp | 10 +-- opengate/tests/src/test068_rotation_source.py | 66 ++++++++++--------- 9 files changed, 62 insertions(+), 70 deletions(-) diff --git a/core/opengate_core/opengate_lib/GateGANSource.cpp b/core/opengate_core/opengate_lib/GateGANSource.cpp index 3786e504b..fafe749e9 100644 --- a/core/opengate_core/opengate_lib/GateGANSource.cpp +++ b/core/opengate_core/opengate_lib/GateGANSource.cpp @@ -88,8 +88,8 @@ void GateGANSource::PrepareNextRun() { // (no need to update th fSPS pos in GateGenericSource) // GateVSource::PrepareNextRun(); // FIXME remove this function ? - // GateGenericSource::PrepareNextRun(); - GateVSource::PrepareNextRun(); + // GateGenericSource::PrepareNextRun(); + GateVSource::PrepareNextRun(); // This global transformation is given to the SPS that will // generate particles in the correct coordinate system auto &l = fThreadLocalData.Get(); diff --git a/core/opengate_core/opengate_lib/GateGenericSource.cpp b/core/opengate_core/opengate_lib/GateGenericSource.cpp index 4c8009dae..62c8df11e 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.cpp +++ b/core/opengate_core/opengate_lib/GateGenericSource.cpp @@ -35,9 +35,6 @@ GateGenericSource::GateGenericSource() : GateVSource() { fParticleDefinition = nullptr; fEffectiveEventTime = -1; fEffectiveEventTime = -1; - - - } GateGenericSource::~GateGenericSource() { @@ -154,7 +151,7 @@ double GateGenericSource::PrepareNextTime(double current_simulation_time) { fCurrentZeroEvents = 0; auto cse = fCurrentSkippedEvents; fCurrentSkippedEvents = 0; - + // if MaxN is below zero, we check the time if (fMaxN <= 0) { if (fEffectiveEventTime < fStartTime) @@ -170,8 +167,8 @@ double GateGenericSource::PrepareNextTime(double current_simulation_time) { } // check according to t MaxN - //std::cout<= fMaxN) { return -1; } @@ -185,8 +182,8 @@ void GateGenericSource::PrepareNextRun() { // This global transformation is given to the SPS that will // generate particles in the correct coordinate system auto &l = fThreadLocalData.Get(); - //auto user_info_pos = py::dict(puser_info["position"]); - //auto pos_init = DictGetG4ThreeVector(user_info_pos, "translation"); + // auto user_info_pos = py::dict(puser_info["position"]); + // auto pos_init = DictGetG4ThreeVector(user_info_pos, "translation"); auto *pos = fSPS->GetPosDist(); pos->SetCentreCoords(l.fGlobalTranslation); @@ -196,9 +193,8 @@ void GateGenericSource::PrepareNextRun() { G4ThreeVector r2(rotation(0, 1), rotation(1, 1), rotation(2, 1)); pos->SetPosRot1(r1); pos->SetPosRot2(r2); - - - auto* ang = fSPS->GetAngDist(); + + auto *ang = fSPS->GetAngDist(); if (fangType == "momentum") { auto new_d = rotation * fInitializeMomentum; @@ -206,11 +202,10 @@ void GateGenericSource::PrepareNextRun() { } if (fangType == "focused") { auto vec_f = fInitiliazeFocusPoint - fInitTranslation; - auto rot_f = rotation*vec_f; + auto rot_f = rotation * vec_f; auto new_f = rot_f + l.fGlobalTranslation; ang->SetFocusPoint(new_f); } - } void GateGenericSource::UpdateEffectiveEventTime( diff --git a/core/opengate_core/opengate_lib/GateGenericSource.h b/core/opengate_core/opengate_lib/GateGenericSource.h index 9bdb3c3dd..f59581f4f 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.h +++ b/core/opengate_core/opengate_lib/GateGenericSource.h @@ -34,7 +34,7 @@ class GateGenericSource : public GateVSource { /// Current number of simulated events in this source /// (do not include skipped events) - //unsigned long fNumberOfGeneratedEvents; + // unsigned long fNumberOfGeneratedEvents; /// Count the number of skipped events /// (e.g. Acceptance Angle or in GANSource) @@ -51,9 +51,9 @@ class GateGenericSource : public GateVSource { const std::vector &activities); protected: - //unsigned long fMaxN; - // We cannot not use a std::unique_ptr - // (or maybe by controlling the deletion during the CleanWorkerThread ?) + // unsigned long fMaxN; + // We cannot not use a std::unique_ptr + // (or maybe by controlling the deletion during the CleanWorkerThread ?) GateSingleParticleSource *fSPS; /* double fActivity; diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp index c3b5cf8d8..de2b90931 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp @@ -15,7 +15,7 @@ GatePhaseSpaceSource::GatePhaseSpaceSource() : GateVSource() { fCharge = 0; fMass = 0; fCurrentBatchSize = 0; - //fMaxN = 0; + // fMaxN = 0; fGlobalFag = false; } @@ -32,7 +32,7 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) { GateVSource::InitializeUserInfo(user_info); // Number of events to generate - //fMaxN = DictGetInt(user_info, "n"); + // fMaxN = DictGetInt(user_info, "n"); // global (world) or local (mother volume) coordinate system fGlobalFag = DictGetBool(user_info, "global_flag"); diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h index 0b61d3e99..24e6edb22 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h @@ -58,7 +58,7 @@ class GatePhaseSpaceSource : public GateVSource { bool fGlobalFag; bool fUseParticleTypeFromFile; - //unsigned long fMaxN; + // unsigned long fMaxN; long fNumberOfGeneratedEvents; size_t fCurrentBatchSize; diff --git a/core/opengate_core/opengate_lib/GateVSource.cpp b/core/opengate_core/opengate_lib/GateVSource.cpp index 00d41e457..e62826075 100644 --- a/core/opengate_core/opengate_lib/GateVSource.cpp +++ b/core/opengate_core/opengate_lib/GateVSource.cpp @@ -7,10 +7,10 @@ #include "GateVSource.h" #include "G4PhysicalVolumeStore.hh" +#include "G4RandomTools.hh" #include "GateHelpers.h" #include "GateHelpersDict.h" #include "GateHelpersGeometry.h" -#include "G4RandomTools.hh" G4Mutex SourceOrientationMutex = G4MUTEX_INITIALIZER; @@ -26,8 +26,6 @@ GateVSource::GateVSource() { fActivity = 0; fHalfLife = -1; fDecayConstant = -1; - - } GateVSource::~GateVSource() {} @@ -38,8 +36,7 @@ void GateVSource::InitializeUserInfo(py::dict &user_info) { fStartTime = DictGetDouble(user_info, "start_time"); fEndTime = DictGetDouble(user_info, "end_time"); fMother = DictGetStr(user_info, "mother"); - - + // get user info about activity or nb of events fMaxN = DictGetInt(user_info, "n"); fActivity = DictGetDouble(user_info, "activity"); @@ -48,23 +45,20 @@ void GateVSource::InitializeUserInfo(py::dict &user_info) { // half life ? fHalfLife = DictGetDouble(user_info, "half_life"); fDecayConstant = log(2) / fHalfLife; - } - void GateVSource::UpdateActivity(double time) { if (fHalfLife <= 0) return; fActivity = fInitialActivity * exp(-fDecayConstant * (time - fStartTime)); } - + double GateVSource::CalcNextTime(double current_simulation_time) { double next_time = current_simulation_time - log(G4UniformRand()) * (1.0 / fActivity); return next_time; } - void GateVSource::PrepareNextRun() { SetOrientationAccordingToMotherVolume(); } double GateVSource::PrepareNextTime(double current_simulation_time) { diff --git a/core/opengate_core/opengate_lib/GateVSource.h b/core/opengate_core/opengate_lib/GateVSource.h index a24b5fe05..bc8a4bc92 100644 --- a/core/opengate_core/opengate_lib/GateVSource.h +++ b/core/opengate_core/opengate_lib/GateVSource.h @@ -28,9 +28,9 @@ class GateVSource { // Called at initialisation to set the source properties from a single dict virtual void InitializeUserInfo(py::dict &user_info); - + virtual void UpdateActivity(double time); - + double CalcNextTime(double current_simulation_time); virtual void PrepareNextRun(); @@ -68,7 +68,6 @@ class GateVSource { G4RotationMatrix fGlobalRotation; }; G4Cache fThreadLocalData; - }; #endif // GateVSource_h diff --git a/core/opengate_core/opengate_lib/GateVoxelsSource.cpp b/core/opengate_core/opengate_lib/GateVoxelsSource.cpp index 5e93acfc0..8922e43f7 100644 --- a/core/opengate_core/opengate_lib/GateVoxelsSource.cpp +++ b/core/opengate_core/opengate_lib/GateVoxelsSource.cpp @@ -17,9 +17,9 @@ GateVoxelsSource::GateVoxelsSource() : GateGenericSource() { GateVoxelsSource::~GateVoxelsSource() {} void GateVoxelsSource::PrepareNextRun() { - //GateGenericSource::PrepareNextRun(); - // rotation and translation to apply, according to mother volume - GateVSource::PrepareNextRun(); + // GateGenericSource::PrepareNextRun(); + // rotation and translation to apply, according to mother volume + GateVSource::PrepareNextRun(); // This global transformation is given to the SPS that will // generate particles in the correct coordinate system auto &l = fThreadLocalData.Get(); @@ -32,8 +32,8 @@ void GateVoxelsSource::PrepareNextRun() { G4ThreeVector r2(rotation(0, 1), rotation(1, 1), rotation(2, 1)); pos->SetPosRot1(r1); pos->SetPosRot2(r2); - - //auto &l = fThreadLocalData.Get(); + + // auto &l = fThreadLocalData.Get(); fVoxelPositionGenerator->fGlobalRotation = l.fGlobalRotation; fVoxelPositionGenerator->fGlobalTranslation = l.fGlobalTranslation; // the direction is 'isotropic' so we don't care about rotating the direction. diff --git a/opengate/tests/src/test068_rotation_source.py b/opengate/tests/src/test068_rotation_source.py index a246000da..8967921f6 100644 --- a/opengate/tests/src/test068_rotation_source.py +++ b/opengate/tests/src/test068_rotation_source.py @@ -8,21 +8,21 @@ import uproot - -#The test generates two differents generic sources, momentum and focused, which is attached to a plan. -#The plan is randomly rotated and we verify that the generated particles have a direction which is in accordance +# The test generates two differents generic sources, momentum and focused, which is attached to a plan. +# The plan is randomly rotated and we verify that the generated particles have a direction which is in accordance # with the applied rotations and the transmissions. -def test068(tab,nb_run,nb_part,nb_source): +def test068(tab, nb_run, nb_part, nb_source): nb_event = len(tab) - nb_event_theo = nb_run * nb_part*nb_source + nb_event_theo = nb_run * nb_part * nb_source err = np.sqrt(nb_event_theo) - if nb_event_theo - 3*err < nb_event < nb_event_theo + 3*err : + if nb_event_theo - 3 * err < nb_event < nb_event_theo + 3 * err: return True - else : + else: return False + if __name__ == "__main__": paths = utility.get_default_test_paths(__file__) output_path = paths.output @@ -62,11 +62,11 @@ def test068(tab,nb_run,nb_part,nb_source): # Plan to attach the source nb_source = 2 - plan = sim.add_volume("Box","source_plan") + plan = sim.add_volume("Box", "source_plan") plan.mother = world.name plan.material = "G4_Galactic" - plan.size = [10*cm,10*cm,1*nm] - plan.translation = np.array([0*cm,0*cm,0*cm]) + plan.size = [10 * cm, 10 * cm, 1 * nm] + plan.translation = np.array([0 * cm, 0 * cm, 0 * cm]) source1 = sim.add_source("GenericSource", "photon_source") source1.particle = "gamma" @@ -75,10 +75,10 @@ def test068(tab,nb_run,nb_part,nb_source): source1.position.size = [10 * cm, 10 * cm, 1 * nm] source1.direction.type = "momentum" # source1.direction.focus_point = [0*cm, 0*cm, -5 *cm] - source1.direction.momentum = [0,0,-1] + source1.direction.momentum = [0, 0, -1] source1.energy.type = "mono" source1.energy.mono = 1 * MeV - source1.activity = n*Bq / ui.number_of_threads + source1.activity = n * Bq / ui.number_of_threads source2 = sim.add_source("GenericSource", "photon_source_2") source2.particle = "gamma" @@ -86,28 +86,26 @@ def test068(tab,nb_run,nb_part,nb_source): source2.mother = plan.name source2.position.size = [10 * cm, 10 * cm, 1 * nm] source2.direction.type = "focused" - source2.direction.focus_point = [0*cm, 0*cm, -5 *cm] + source2.direction.focus_point = [0 * cm, 0 * cm, -5 * cm] # source1.direction.momentum = [0, 0, -1] source2.energy.type = "mono" source2.energy.mono = 1 * MeV source2.activity = n * Bq / ui.number_of_threads - # Phase Space - phsp_plan = sim.add_volume("Box","phsp_plan") + phsp_plan = sim.add_volume("Box", "phsp_plan") phsp_plan.mother = world.name phsp_plan.material = "G4_Galactic" - phsp_plan.size = [10*cm,10*cm,1*nm] - phsp_plan.translation = np.array([0*cm,0*cm,-5*cm]) + phsp_plan.size = [10 * cm, 10 * cm, 1 * nm] + phsp_plan.translation = np.array([0 * cm, 0 * cm, -5 * cm]) phsp = sim.add_actor("PhaseSpaceActor", "PhaseSpace") phsp.mother = phsp_plan.name phsp.attributes = ["EventID"] phsp.output = output_path / "test068.root" - - #MOTION ACTOR for the source and actors + # MOTION ACTOR for the source and actors motion_source = sim.add_actor("MotionVolumeActor", "Move_source") motion_source.mother = plan.name motion_source.rotations = [] @@ -128,22 +126,24 @@ def test068(tab,nb_run,nb_part,nb_source): y_rot = 360 * np.random.random() z_rot = 360 * np.random.random() - motion_source.translations.append([x_translation,y_translation,z_translation]) - rot = Rotation.from_euler("xyz",[x_rot,y_rot,z_rot], degrees=True).as_matrix() + motion_source.translations.append([x_translation, y_translation, z_translation]) + rot = Rotation.from_euler( + "xyz", [x_rot, y_rot, z_rot], degrees=True + ).as_matrix() vec_source_phsp = phsp_plan.translation - plan.translation # vec_source_phsp = vec_source_phsp.reshape((len(vec_source_phsp),1)) # print(vec_source_phsp) - translation_phsp = np.dot(rot,vec_source_phsp) + np.array([x_translation,y_translation,z_translation]) + translation_phsp = np.dot(rot, vec_source_phsp) + np.array( + [x_translation, y_translation, z_translation] + ) motion_source.rotations.append(rot) motion_phsp.rotations.append(rot) motion_phsp.translations.append(translation_phsp) - sim.run_timing_intervals.append([i* sec, (i+1)* sec]) + sim.run_timing_intervals.append([i * sec, (i + 1) * sec]) # stat actor s = sim.add_actor("SimulationStatisticsActor", "Stats") s.track_types_flag = True - - # Physic list and cuts p = sim.get_physics_user_info() p.physics_list_name = "G4EmStandardPhysics_option3" @@ -156,16 +156,20 @@ def test068(tab,nb_run,nb_part,nb_source): output = sim.run() - - #Validation test + # Validation test f_phsp = uproot.open(phsp.output) arr = f_phsp["PhaseSpace"].arrays() print("Number of detected events :", len(arr)) nb_event = len(arr) nb_part = n nb_run = len(sim.run_timing_intervals) - err = np.sqrt(nb_part*nb_run*nb_source) - print("Number of expected events : [" +str(int(-3*err + nb_run*nb_part*nb_source)) + "," + str(int(+3*err + nb_run*nb_part*nb_source))+ "]") - is_ok = test068(arr,len(sim.run_timing_intervals),nb_part,2) + err = np.sqrt(nb_part * nb_run * nb_source) + print( + "Number of expected events : [" + + str(int(-3 * err + nb_run * nb_part * nb_source)) + + "," + + str(int(+3 * err + nb_run * nb_part * nb_source)) + + "]" + ) + is_ok = test068(arr, len(sim.run_timing_intervals), nb_part, 2) utility.test_ok(is_ok) - From 98e0997fe8ca4b84419347481f76ed04988dcbe0 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 11 Dec 2023 11:55:06 +0100 Subject: [PATCH 33/46] In test060_phsp_source_untilNextPrimary remove a useless test related to another test --- opengate/tests/src/test060_phsp_source_untilNextPrimary.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/opengate/tests/src/test060_phsp_source_untilNextPrimary.py b/opengate/tests/src/test060_phsp_source_untilNextPrimary.py index de675743c..c90d014c0 100644 --- a/opengate/tests/src/test060_phsp_source_untilNextPrimary.py +++ b/opengate/tests/src/test060_phsp_source_untilNextPrimary.py @@ -70,11 +70,6 @@ def main(): phs_file_name_out=paths.output / "test_source_untilPrimary.root", ) - is_ok = t.check_value_from_root_file( - file_name_root=paths.output / "test_source_translation.root", - key="PrePositionLocal_X", - ref_value=30 * mm, - ) # load data from root file data_ref, keys_ref, m_ref = phsp.load( paths.output / "test_source_untilPrimary.root" From 01c3e98a30ade020d5231610b642bbe426721efc Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 11 Dec 2023 11:55:31 +0100 Subject: [PATCH 34/46] Update tests 66 and 68 with new merged commits --- opengate/tests/src/test066_phspsource_activity.py | 10 +++++----- opengate/tests/src/test068_rotation_source.py | 9 ++++----- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/opengate/tests/src/test066_phspsource_activity.py b/opengate/tests/src/test066_phspsource_activity.py index 979c7c012..c25495cd9 100644 --- a/opengate/tests/src/test066_phspsource_activity.py +++ b/opengate/tests/src/test066_phspsource_activity.py @@ -67,7 +67,7 @@ def validation_test(n, n_measured): source_1.global_flag = False source_1.particle = "gamma" source_1.batch_size = 100 - source_1.override_position = True + source_1.translate_position = True source_1.position.translation = [0, 0, 300] source_1.activity = nb_part * Bq @@ -90,15 +90,15 @@ def validation_test(n, n_measured): s = sim.add_actor("SimulationStatisticsActor", "Stats") s.track_types_flag = True - p = sim.get_physics_user_info() - p.physics_list_name = "G4EmStandardPhysics_option3" - p.enable_decay = False + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.enable_decay = False sim.physics_manager.global_production_cuts.gamma = 1 * mm sim.physics_manager.global_production_cuts.electron = 1 * mm sim.physics_manager.global_production_cuts.positron = 1 * mm # # # go ! - output = sim.start() + sim.run() + output = sim.output # # # print results diff --git a/opengate/tests/src/test068_rotation_source.py b/opengate/tests/src/test068_rotation_source.py index 8967921f6..709bd4d61 100644 --- a/opengate/tests/src/test068_rotation_source.py +++ b/opengate/tests/src/test068_rotation_source.py @@ -145,16 +145,15 @@ def test068(tab, nb_run, nb_part, nb_source): s.track_types_flag = True # Physic list and cuts - p = sim.get_physics_user_info() - p.physics_list_name = "G4EmStandardPhysics_option3" - p.enable_decay = False + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.enable_decay = False sim.physics_manager.global_production_cuts.gamma = 1 * mm sim.physics_manager.global_production_cuts.electron = 1 * mm sim.physics_manager.global_production_cuts.positron = 1 * mm # go ! - - output = sim.run() + sim.run() + output = sim.output # Validation test f_phsp = uproot.open(phsp.output) From 4ce3b0d421e2eee31b0209cac9d0b0ece5aa5f99 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 11 Dec 2023 16:42:13 +0100 Subject: [PATCH 35/46] When rotate the momentum of the source, maybe sometimes we do not want (eg: test042). Add a flag to activate the option --- core/opengate_core/opengate_lib/GateGANSource.h | 2 +- core/opengate_core/opengate_lib/GateGenericSource.cpp | 7 +++++-- core/opengate_core/opengate_lib/GateGenericSource.h | 8 ++++++-- opengate/sources/generic.py | 1 + opengate/tests/src/test068_rotation_source.py | 2 ++ 5 files changed, 15 insertions(+), 5 deletions(-) diff --git a/core/opengate_core/opengate_lib/GateGANSource.h b/core/opengate_core/opengate_lib/GateGANSource.h index bb1be58a4..fd0fb768a 100644 --- a/core/opengate_core/opengate_lib/GateGANSource.h +++ b/core/opengate_core/opengate_lib/GateGANSource.h @@ -28,7 +28,7 @@ class GateGANSource : public GateGenericSource { void InitializeUserInfo(py::dict &user_info) override; - void PrepareNextRun(); + void PrepareNextRun() override; void GeneratePrimaries(G4Event *event, double current_simulation_time) override; diff --git a/core/opengate_core/opengate_lib/GateGenericSource.cpp b/core/opengate_core/opengate_lib/GateGenericSource.cpp index 62c8df11e..6cd829c9d 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.cpp +++ b/core/opengate_core/opengate_lib/GateGenericSource.cpp @@ -35,6 +35,7 @@ GateGenericSource::GateGenericSource() : GateVSource() { fParticleDefinition = nullptr; fEffectiveEventTime = -1; fEffectiveEventTime = -1; + fforceRotation = false; } GateGenericSource::~GateGenericSource() { @@ -106,6 +107,8 @@ void GateGenericSource::InitializeUserInfo(py::dict &user_info) { fCurrentSkippedEvents = 0; fTotalSkippedEvents = 0; fEffectiveEventTime = -1; + + fforceRotation = DictGetDouble(user_info, "force_rotation"); } void GateGenericSource::UpdateActivity(double time) { @@ -196,11 +199,11 @@ void GateGenericSource::PrepareNextRun() { auto *ang = fSPS->GetAngDist(); - if (fangType == "momentum") { + if (fangType == "momentum" && fforceRotation) { auto new_d = rotation * fInitializeMomentum; ang->SetParticleMomentumDirection(new_d); } - if (fangType == "focused") { + if (fangType == "focused" && fforceRotation) { auto vec_f = fInitiliazeFocusPoint - fInitTranslation; auto rot_f = rotation * vec_f; auto new_f = rot_f + l.fGlobalTranslation; diff --git a/core/opengate_core/opengate_lib/GateGenericSource.h b/core/opengate_core/opengate_lib/GateGenericSource.h index f59581f4f..5034481d6 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.h +++ b/core/opengate_core/opengate_lib/GateGenericSource.h @@ -26,7 +26,7 @@ class GateGenericSource : public GateVSource { void InitializeUserInfo(py::dict &user_info) override; - double PrepareNextTime(double current_simulation_time); + double PrepareNextTime(double current_simulation_time) override; void PrepareNextRun() override; @@ -84,6 +84,10 @@ class GateGenericSource : public GateVSource { double fWeight; double fWeightSigma; + // Force the rotation of momentum and focal point to follow rotation of the + // source, eg: needed for motion actor + bool fforceRotation; + // angular acceptance management struct threadLocalT { GateAcceptanceAngleTesterManager *fAAManager; @@ -112,7 +116,7 @@ class GateGenericSource : public GateVSource { virtual void InitializeEnergy(py::dict user_info); - virtual void UpdateActivity(double time); + virtual void UpdateActivity(double time) override; void UpdateEffectiveEventTime(double current_simulation_time, unsigned long skipped_particle); diff --git a/opengate/sources/generic.py b/opengate/sources/generic.py index d1d7d0b8b..71cbea957 100644 --- a/opengate/sources/generic.py +++ b/opengate/sources/generic.py @@ -313,6 +313,7 @@ def set_default_user_info(user_info): user_info.user_particle_life_time = -1 # negative means : by default user_info.tac_times = None user_info.tac_activities = None + user_info.force_rotation = False # ion user_info.ion = Box() user_info.ion.Z = 0 # Z: Atomic Number diff --git a/opengate/tests/src/test068_rotation_source.py b/opengate/tests/src/test068_rotation_source.py index 709bd4d61..7c2fe7761 100644 --- a/opengate/tests/src/test068_rotation_source.py +++ b/opengate/tests/src/test068_rotation_source.py @@ -74,6 +74,7 @@ def test068(tab, nb_run, nb_part, nb_source): source1.mother = plan.name source1.position.size = [10 * cm, 10 * cm, 1 * nm] source1.direction.type = "momentum" + source1.force_rotation = True # source1.direction.focus_point = [0*cm, 0*cm, -5 *cm] source1.direction.momentum = [0, 0, -1] source1.energy.type = "mono" @@ -86,6 +87,7 @@ def test068(tab, nb_run, nb_part, nb_source): source2.mother = plan.name source2.position.size = [10 * cm, 10 * cm, 1 * nm] source2.direction.type = "focused" + source2.force_rotation = True source2.direction.focus_point = [0 * cm, 0 * cm, -5 * cm] # source1.direction.momentum = [0, 0, -1] source2.energy.type = "mono" From a15ce4f51dc10b84b1e205538335804169c19a18 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 11 Dec 2023 16:51:55 +0100 Subject: [PATCH 36/46] Remove print name of the volume during the simulation --- opengate/geometry/utility.py | 1 - 1 file changed, 1 deletion(-) diff --git a/opengate/geometry/utility.py b/opengate/geometry/utility.py index 24b70fed4..3bec06d18 100644 --- a/opengate/geometry/utility.py +++ b/opengate/geometry/utility.py @@ -188,7 +188,6 @@ def get_transform_world_to_local(volume): ctr = volume.translation_list[i] crot = volume.rotation_list[i] for vol in volume.ancestor_volumes[::-1]: - print(vol.name) crot = np.matmul(vol.rotation_list[0], crot) ctr = vol.rotation_list[0].dot(ctr) + vol.translation_list[0] cumulative_translation.append(ctr) From 7bf29999b301b2a8916a8f8951c697454ed62711 Mon Sep 17 00:00:00 2001 From: Hermann Fuchs Date: Tue, 12 Dec 2023 15:39:21 +0100 Subject: [PATCH 37/46] Added developer documentation for segmentation faults. Based on feedback from D. Sarrut and N. Krah --- docs/source/developer_guide.md | 51 ++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/docs/source/developer_guide.md b/docs/source/developer_guide.md index 64e8ad65e..1afe9716b 100644 --- a/docs/source/developer_guide.md +++ b/docs/source/developer_guide.md @@ -475,6 +475,57 @@ Below are a list of hints (compared to boost-python). - Pure virtual need a trampoline class - Python debug: `python -q -X faulthandler` + +### Memory managment between Python and Geant4 - Why is there a segmentation fault? + +Code and memory (objects) can be on the python and/or C++(Geant4) side. +To avoid memory leakage, the G4 objects need to be deleted at some point, either by the C++ side, meaning G4, or by the garbage collector on the python side. This requires some understanding of how life time is managed in G4. GATE 10 uses the G4RunManager to initialize, run, and deconstruct a simulation. The G4RunManager's destructor triggers a nested sequence of calls to the destructors many objects the run manager handles, e.g. geometry and physics. Those objects, if they are created on the python-side, e.g. because the Construct() method is implemented in python, should never be deleted on the python side because the G4RunManager is responsible for deletion and segfaults if the objects exist now longer at the time of supposed destruction. In this case, the no py::nodelete option in pybind11 is the correct way. + +​ +​There is also the python side concerning objects which are deleted by the G4RunManager: If python stores a reference to the object, and it generally should (see below), this reference will not find the object anymore once the G4RunManager has been deleted. Therefore, we have implement a mechanism which makes sure that references to G4 objects are unset before the G4RunManager is garbage collected (and thus its destructor is called). This is done via the close() and release_g4_references() methods combined with the “with SimulationEngine as se” context clause. +​ + +​The situation is somewhat different for objects that are not automatically deleted by the Geant4 session. This concerns objects which a G4 user would manually destruct and delete. The "user" from a G4 perspective is the python-side of G4, specifically the SimulationEngine, which creates and controls the G4RunManager. The objects in question should be destroyed on the python side, and in fact they (usually) are via garbage collection. To this end, they should not be bound with the py::nodelete option - that is what the "usually" in the previous sentence referred to. ​ + +In any case, G4 objects created on the python side should not be destroyed (garbage collected) too early, i.e. before the end of the G4 simulation. It is important to note here that garbage collection in python is strongly tied to reference counting, meaning, objects may be garbage collected as soon as there is no reference anymore pointing to them. Or in other words: you can prevent garbage collection by keeping a reference. In practice: If you create a G4 object that is required to persist for the duration of the simulation and you create this object in a local scope, e.g. in a function, you need to make sure to store a reference somewhere that lives beyond the local scope. Otherwise, when the function goes out of scope, the local reference no longer exists and the G4 object may be garbage collected. This is the reason why GATE 10 stores references to G4 objects in class attributes, either as plane reference, or via lists or dictionaries. A nasty detail: if a G4 object is only referenced locally (implementation error), the moment when the segfault occurs is not always the same because garbage collection in python is scheduled (for performance reasons), meaning objects are collected any time after its last reference is released, but not necessarily at that moment. This can cause a seeming randomness in the segfaults. + +​ +​So, after a long explanation, the key points are: + +#. check the nature of your G4 object, for instance G4TriangularFacet + +​#. if it is designed to be deleted by the run manager, add the py::no_delete option in the binding + +​#. In any case, make sure to store a non-local persistent reference on the python side, ideally as a class attribute + +#. add a “release statement” in the release_g4_references() method for each G4 reference (strictly only to RunManager-handled objects, but does not hurt for the others) in the class, and make sure the release_g4_references() method is called by the class’s close(). + +An example may be the implementation of the Tesselated Volume (STL). + +Pybindings in pyG4TriangularFacet.cpp: + +``py::class_(m, "G4TriangularFacet")`` +will cause a segmentation fault + +``py::class_>(m, "G4TriangularFacet")`` +will work as python is instructed not to delete the object. + +On the python side in geometry/solids.py: + +```python +def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + self.g4_solid = None + self.facetArray = None + self.tessellated_solid = None +``` + +All data which is sent to Geant4 is included as a variable to the class. As such it is only deleted at the end of the simulation, not when the function ends. +Which, would cause a segmentation fault. + + + ### ITK - The following [issue](https://github.com/OpenGATE/opengate/issues/216) occured in the past and was solved by updating ITK: From a3ba9758030af1d8fff86e0a45caab75aba626cb Mon Sep 17 00:00:00 2001 From: nkrah <34276065+nkrah@users.noreply.github.com> Date: Wed, 13 Dec 2023 01:27:07 +0100 Subject: [PATCH 38/46] Correct typos and proof-read section in developer guide --- docs/source/developer_guide.md | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/docs/source/developer_guide.md b/docs/source/developer_guide.md index 1afe9716b..18b37dee2 100644 --- a/docs/source/developer_guide.md +++ b/docs/source/developer_guide.md @@ -478,27 +478,27 @@ Below are a list of hints (compared to boost-python). ### Memory managment between Python and Geant4 - Why is there a segmentation fault? -Code and memory (objects) can be on the python and/or C++(Geant4) side. -To avoid memory leakage, the G4 objects need to be deleted at some point, either by the C++ side, meaning G4, or by the garbage collector on the python side. This requires some understanding of how life time is managed in G4. GATE 10 uses the G4RunManager to initialize, run, and deconstruct a simulation. The G4RunManager's destructor triggers a nested sequence of calls to the destructors many objects the run manager handles, e.g. geometry and physics. Those objects, if they are created on the python-side, e.g. because the Construct() method is implemented in python, should never be deleted on the python side because the G4RunManager is responsible for deletion and segfaults if the objects exist now longer at the time of supposed destruction. In this case, the no py::nodelete option in pybind11 is the correct way. +Code and memory (objects) can be on the python and/or C++(Geant4)-side. +To avoid memory leakage, the G4 objects need to be deleted at some point, either by the C++-side, meaning G4, or by the garbage collector on the python-side. This requires some understanding of how lifetime is managed in G4. GATE 10 uses the G4RunManager to initialize, run, and deconstruct a simulation. The G4RunManager's destructor triggers a nested sequence of calls to the destructors of many objects the run manager handles, e.g. geometry and physics. These objects, if they are created on the python-side, e.g. because the Construct() method is implemented in python, should never be deleted on the python-side because the G4RunManager is responsible for deletion and segfaults occur if the objects exist now longer at the time of supposed destruction. In this case, the no `py::nodelete` option in pybind11 is the correct way. ​ -​There is also the python side concerning objects which are deleted by the G4RunManager: If python stores a reference to the object, and it generally should (see below), this reference will not find the object anymore once the G4RunManager has been deleted. Therefore, we have implement a mechanism which makes sure that references to G4 objects are unset before the G4RunManager is garbage collected (and thus its destructor is called). This is done via the close() and release_g4_references() methods combined with the “with SimulationEngine as se” context clause. +​There is also the python-side concerning objects which are deleted by the G4RunManager: If python stores a reference to the object, and it generally should (see below), this reference will not find the object anymore once the G4RunManager has been deleted. Therefore, we have implemented a mechanism which makes sure that references to G4 objects are unset before the G4RunManager is garbage collected (and thus its destructor is called). This is done via the close() and release_g4_references() methods combined with a `with SimulationEngine as se` context clause. ​ -​The situation is somewhat different for objects that are not automatically deleted by the Geant4 session. This concerns objects which a G4 user would manually destruct and delete. The "user" from a G4 perspective is the python-side of G4, specifically the SimulationEngine, which creates and controls the G4RunManager. The objects in question should be destroyed on the python side, and in fact they (usually) are via garbage collection. To this end, they should not be bound with the py::nodelete option - that is what the "usually" in the previous sentence referred to. ​ +​The situation is somewhat different for objects that are not automatically deleted by the Geant4 session. This concerns objects which a G4 user would manually destruct and delete. The "user", from a G4 perspective, is the python-side of GATE, specifically the SimulationEngine, which creates and controls the G4RunManager. The objects in question should be destroyed on the python side, and in fact they (normally) are via garbage collection. To this end, they should **not** be bound with the `py::nodelete` option - which would prevent garbage collection. ​ -In any case, G4 objects created on the python side should not be destroyed (garbage collected) too early, i.e. before the end of the G4 simulation. It is important to note here that garbage collection in python is strongly tied to reference counting, meaning, objects may be garbage collected as soon as there is no reference anymore pointing to them. Or in other words: you can prevent garbage collection by keeping a reference. In practice: If you create a G4 object that is required to persist for the duration of the simulation and you create this object in a local scope, e.g. in a function, you need to make sure to store a reference somewhere that lives beyond the local scope. Otherwise, when the function goes out of scope, the local reference no longer exists and the G4 object may be garbage collected. This is the reason why GATE 10 stores references to G4 objects in class attributes, either as plane reference, or via lists or dictionaries. A nasty detail: if a G4 object is only referenced locally (implementation error), the moment when the segfault occurs is not always the same because garbage collection in python is scheduled (for performance reasons), meaning objects are collected any time after its last reference is released, but not necessarily at that moment. This can cause a seeming randomness in the segfaults. +In any case, G4 objects created on the python-side should not be destroyed (garbage collected) too early, i.e. not before the end of the G4 simulation, to avoid segfaults when G4 tries to a find the object to which it holds a pointer. It is important to note here that garbage collection in python is strongly tied to reference counting, meaning, objects may be garbage collected as soon as there are no more references pointing to them. Or in other words: you can prevent garbage collection by keeping a reference. In practice: If you create a G4 object that is required to persist for the duration of the simulation and you create this object in a local scope, e.g. in a function, you need to make sure to store a reference somewhere so that it lives beyond the local scope. Otherwise, when the function goes out of scope, the local reference no longer exists and the G4 object may be garbage collected. This is the reason why GATE 10 stores references to G4 objects in class attributes, either as plane reference, or via lists or dictionaries. A nasty detail: in case of a G4 object that is only referenced locally (implementation error), the moment when the segfault occurs might vary because garbage collection in python is typically scheduled (for performance reasons), meaning objects are collected any time after their last reference is released, but not necessarily at that moment. This can cause a seeming randomness in the segfaults. ​ ​So, after a long explanation, the key points are: -#. check the nature of your G4 object, for instance G4TriangularFacet +#. check the nature of your G4 object -​#. if it is designed to be deleted by the run manager, add the py::no_delete option in the binding +​#. if it is designed to be deleted by the G4RunManager, add the `py::no_delete` option in the pybind11 binding -​#. In any case, make sure to store a non-local persistent reference on the python side, ideally as a class attribute +​#. In any case, make sure to store a non-local persistent reference on the python-side, ideally as a class attribute -#. add a “release statement” in the release_g4_references() method for each G4 reference (strictly only to RunManager-handled objects, but does not hurt for the others) in the class, and make sure the release_g4_references() method is called by the class’s close(). +#. add a “release statement” in the release_g4_references() method for each G4 reference (strictly only to RunManager-handled objects, but it does not hurt for others) in the class, and make sure the release_g4_references() method is called by the class’s close() method. An example may be the implementation of the Tesselated Volume (STL). @@ -510,18 +510,18 @@ will cause a segmentation fault ``py::class_>(m, "G4TriangularFacet")`` will work as python is instructed not to delete the object. -On the python side in geometry/solids.py: +On the python-side in geometry/solids.py: ```python def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) self.g4_solid = None + self.g4_tessellated_solid = None self.facetArray = None - self.tessellated_solid = None ``` -All data which is sent to Geant4 is included as a variable to the class. As such it is only deleted at the end of the simulation, not when the function ends. +All data which is sent to Geant4 is included as a variable to the class. As such, it is only deleted at the end of the simulation, not when the function ends. Which, would cause a segmentation fault. From 54ae79df46fb02ddc49b8f2b47fcf079e5c7e714 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 18 Dec 2023 09:45:59 +0100 Subject: [PATCH 39/46] During Github actions, remove the correct folder for brew --- .github/workflows/main.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 856cf8f93..b9d3b6146 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -120,7 +120,7 @@ jobs: shell: bash -l {0} run: | brew update - rm -rf /usr/local/bin/python3.11-config /usr/local/bin/2to3-3.11 /usr/local/bin/idle3.11 /usr/local/bin/pydoc3.11 /usr/local/bin/python3.11 + rm -rf /usr/local/bin/python3.1*-config /usr/local/bin/2to3-3.1* /usr/local/bin/idle3.1* /usr/local/bin/pydoc3.1* /usr/local/bin/python3.1* rm -rf /usr/local/bin/python3-config /usr/local/bin/2to3 /usr/local/bin/idle3 /usr/local/bin/pydoc3 /usr/local/bin/python3 brew install --force --verbose --overwrite \ ccache \ From 5f18424bb5f65b9fd3167480446884e5dde0dc85 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 17:55:08 +0100 Subject: [PATCH 40/46] first version --- core/external/pybind11 | 2 +- core/opengate_core/opengate_core.cpp | 3 + .../opengate_lib/GateFluenceActor.cpp | 91 ++++++++++++++ .../opengate_lib/GateFluenceActor.h | 58 +++++++++ .../opengate_lib/GateGANSource.cpp | 1 + .../opengate_lib/pyGateFluenceActor.cpp | 22 ++++ opengate/actors/actorbuilders.py | 3 +- opengate/actors/doseactors.py | 111 +++++++++++++++++- opengate/tests/data | 2 +- 9 files changed, 289 insertions(+), 4 deletions(-) create mode 100644 core/opengate_core/opengate_lib/GateFluenceActor.cpp create mode 100644 core/opengate_core/opengate_lib/GateFluenceActor.h create mode 100644 core/opengate_core/opengate_lib/pyGateFluenceActor.cpp diff --git a/core/external/pybind11 b/core/external/pybind11 index a34596bfe..a67d78657 160000 --- a/core/external/pybind11 +++ b/core/external/pybind11 @@ -1 +1 @@ -Subproject commit a34596bfe1947b4a6b0bcc4218e1f72d0c2e9b4c +Subproject commit a67d786571cc6d2c2ecec33985e00e6d22dc68e5 diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index 23d4a137c..86e314219 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -249,6 +249,8 @@ void init_GateKineticEnergyFilter(py::module &); void init_GateDoseActor(py::module &m); +void init_GateFluenceActor(py::module &m); + void init_GateLETActor(py::module &m); void init_GateARFActor(py::module &m); @@ -475,6 +477,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateEventAction(m); init_GateTrackingAction(m); init_GateDoseActor(m); + init_GateFluenceActor(m); init_GateLETActor(m); init_GateSimulationStatisticsActor(m); init_GatePhaseSpaceActor(m); diff --git a/core/opengate_core/opengate_lib/GateFluenceActor.cpp b/core/opengate_core/opengate_lib/GateFluenceActor.cpp new file mode 100644 index 000000000..bc1938646 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateFluenceActor.cpp @@ -0,0 +1,91 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ------------------------------------ -------------- */ + +#include "G4Navigator.hh" +#include "G4RunManager.hh" +#include "G4Threading.hh" + +#include "GateFluenceActor.h" +#include "GateHelpersDict.h" +#include "GateHelpersImage.h" + +#include +#include +#include + +// Mutex that will be used by thread to write the output image +G4Mutex SetPixelFluenceMutex = G4MUTEX_INITIALIZER; + +GateFluenceActor::GateFluenceActor(py::dict &user_info) + : GateVActor(user_info, true) { + // Create the image pointer + // (the size and allocation will be performed on the py side) + cpp_fluence_image = Image3DType::New(); + // Action for this actor: during stepping + fActions.insert("SteppingAction"); + fActions.insert("BeginOfRunAction"); + fActions.insert("BeginOfEventAction"); + fActions.insert("EndSimulationAction"); + // translation + fInitialTranslation = DictGetG4ThreeVector(user_info, "translation"); +} + +void GateFluenceActor::ActorInitialize() {} + +void GateFluenceActor::BeginOfRunAction(const G4Run *) { + Image3DType::RegionType region = + cpp_fluence_image->GetLargestPossibleRegion(); + size_fluence = region.GetSize(); + + // Important ! The volume may have moved, so we re-attach each run + AttachImageToVolume(cpp_fluence_image, fPhysicalVolumeName, + fInitialTranslation); + // compute volume of a dose voxel + auto sp = cpp_fluence_image->GetSpacing(); + fVoxelVolume = sp[0] * sp[1] * sp[2]; +} + +void GateFluenceActor::BeginOfEventAction(const G4Event *event) {} + +void GateFluenceActor::SteppingAction(G4Step *step) { + // same method to consider only entering tracks + if (step->GetPreStepPoint()->GetStepStatus() == fGeomBoundary) { + // the pre-position is at the edge + auto preGlobal = step->GetPreStepPoint()->GetPosition(); + auto dir = step->GetPreStepPoint()->GetMomentumDirection(); + auto touchable = step->GetPreStepPoint()->GetTouchable(); + + // consider position in the local volume, slightly shifted by 0.1 nm because + // otherwise, it can be considered as outside the volume by isInside. + auto position = preGlobal + 0.1 * CLHEP::nm * dir; + auto localPosition = + touchable->GetHistory()->GetTransform(0).TransformPoint(position); + + // convert G4ThreeVector to itk PointType + Image3DType::PointType point; + point[0] = localPosition[0]; + point[1] = localPosition[1]; + point[2] = localPosition[2]; + + // get weight + auto w = step->GetTrack()->GetWeight(); + + // get pixel index + Image3DType::IndexType index; + bool isInside = + cpp_fluence_image->TransformPhysicalPointToIndex(point, index); + + // set value + if (isInside) { + G4AutoLock FluenceMutex(&SetPixelFluenceMutex); + // add hit + ImageAddValue(cpp_fluence_image, index, w); + } // else : outside the image + } +} + +void GateFluenceActor::EndSimulationAction() {} diff --git a/core/opengate_core/opengate_lib/GateFluenceActor.h b/core/opengate_core/opengate_lib/GateFluenceActor.h new file mode 100644 index 000000000..d296b6d89 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateFluenceActor.h @@ -0,0 +1,58 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateFluenceActor_h +#define GateFluenceActor_h + +#include "G4Cache.hh" +#include "G4VPrimitiveScorer.hh" +#include "GateVActor.h" +#include "itkImage.h" +#include +#include + +namespace py = pybind11; + +class GateFluenceActor : public GateVActor { + +public: + // Constructor + GateFluenceActor(py::dict &user_info); + + virtual void ActorInitialize(); + + // Main function called every step in attached volume + virtual void SteppingAction(G4Step *); + + // Called every time a Run starts (all threads) + virtual void BeginOfRunAction(const G4Run *run); + + virtual void BeginOfEventAction(const G4Event *event); + + virtual void EndSimulationAction(); + + // Image type is 3D float by default + typedef itk::Image Image3DType; + + typedef itk::Image Image4DType; + typedef itk::Image ImageInt4DType; + using Size4DType = Image4DType::SizeType; + Size4DType size_4D; + + // The image is accessible on py side (shared by all threads) + Image3DType::Pointer cpp_fluence_image; + Image3DType::SizeType size_fluence; + double fVoxelVolume; + int NbOfThreads = 0; + + std::string fPhysicalVolumeName; + + G4ThreeVector fInitialTranslation; + std::string fHitType; +}; + +#endif // GateFluenceActor_h diff --git a/core/opengate_core/opengate_lib/GateGANSource.cpp b/core/opengate_core/opengate_lib/GateGANSource.cpp index fafe749e9..12c82a0d8 100644 --- a/core/opengate_core/opengate_lib/GateGANSource.cpp +++ b/core/opengate_core/opengate_lib/GateGANSource.cpp @@ -188,6 +188,7 @@ void GateGANSource::GenerateOnePrimary(G4Event *event, // (if it is not set by the GAN, we may avoid to sample at each iteration) if (fPosition_is_set_by_GAN || fCurrentSkippedEvents == 0) position = GeneratePrimariesPosition(); + // FIXME change position is not set by GAN // direction // (if it is not set by the GAN, we may avoid to sample at each iteration) diff --git a/core/opengate_core/opengate_lib/pyGateFluenceActor.cpp b/core/opengate_core/opengate_lib/pyGateFluenceActor.cpp new file mode 100644 index 000000000..c21fe9e3e --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateFluenceActor.cpp @@ -0,0 +1,22 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include + +namespace py = pybind11; + +#include "GateFluenceActor.h" + +void init_GateFluenceActor(py::module &m) { + py::class_, + GateVActor>(m, "GateFluenceActor") + .def(py::init()) + .def_readwrite("cpp_fluence_image", &GateFluenceActor::cpp_fluence_image) + .def_readwrite("fPhysicalVolumeName", + &GateFluenceActor::fPhysicalVolumeName); +} diff --git a/opengate/actors/actorbuilders.py b/opengate/actors/actorbuilders.py index d597ea94e..284d89a6e 100644 --- a/opengate/actors/actorbuilders.py +++ b/opengate/actors/actorbuilders.py @@ -1,5 +1,5 @@ from .arfactors import ARFActor, ARFTrainingDatasetActor -from .doseactors import DoseActor, LETActor +from .doseactors import DoseActor, LETActor, FluenceActor from .digitizers import ( DigitizerAdderActor, DigitizerReadoutActor, @@ -24,6 +24,7 @@ actor_type_names = { SimulationStatisticsActor, DoseActor, + FluenceActor, LETActor, SourceInfoActor, PhaseSpaceActor, diff --git a/opengate/actors/doseactors.py b/opengate/actors/doseactors.py index 292f63c7b..7c83a73d1 100644 --- a/opengate/actors/doseactors.py +++ b/opengate/actors/doseactors.py @@ -200,7 +200,6 @@ def StartSimulationAction(self): self.output_origin = self.user_info.output_origin def EndSimulationAction(self): - # print(lol) g4.GateDoseActor.EndSimulationAction(self) # Get the itk image from the cpp side @@ -515,3 +514,113 @@ def EndSimulationAction(self): itk.imwrite(self.py_numerator_image, ensure_filename_is_str(fPath)) fPath = fPath.replace("_numerator", "_denominator") itk.imwrite(self.py_denominator_image, ensure_filename_is_str(fPath)) + + +class FluenceActor(g4.GateFluenceActor, ActorBase): + """ + FluenceActor: compute a 3D map of fluence + + FIXME: add scatter order and uncertainty + """ + + type_name = "FluenceActor" + + def set_default_user_info(user_info): + ActorBase.set_default_user_info(user_info) + # required user info, default values + mm = g4_units.mm + user_info.size = [10, 10, 10] + user_info.spacing = [1 * mm, 1 * mm, 1 * mm] + user_info.output = "fluence.mhd" + user_info.translation = [0, 0, 0] + user_info.physical_volume_index = None + user_info.uncertainty = False + user_info.scatter = False + + def __init__(self, user_info): + ActorBase.__init__(self, user_info) + g4.GateFluenceActor.__init__(self, user_info.__dict__) + # attached physical volume (at init) + self.g4_phys_vol = None + # default image (py side) + self.py_fluence_image = None + + def __str__(self): + u = self.user_info + s = f'FluenceActor "{u.name}": dim={u.size} spacing={u.spacing} {u.output} tr={u.translation}' + return s + + def __getstate__(self): + # superclass getstate + DoseActor.__getstate__(self) + return self.__dict__ + + def initialize(self, volume_engine=None): + super().initialize(volume_engine) + # create itk image (py side) + size = np.array(self.user_info.size) + spacing = np.array(self.user_info.spacing) + self.py_fluence_image = create_3d_image(size, spacing) + # compute the center, using translation and half pixel spacing + self.img_origin_during_run = ( + -size * spacing / 2.0 + spacing / 2.0 + self.user_info.translation + ) + # for initialization during the first run + self.first_run = True + + def StartSimulationAction(self): + # init the origin and direction according to the physical volume + # (will be updated in the BeginOfRun) + attached_to_volume = self.volume_engine.get_volume(self.user_info.mother) + if self.user_info.physical_volume_index is None: + physical_volume_index = 0 + else: + physical_volume_index = self.user_info.physical_volume_index + try: + self.g4_phys_vol = attached_to_volume.g4_physical_volumes[ + physical_volume_index + ] + except IndexError: + fatal( + f"Error in the FluenceActor {self.user_info.name}. " + f"Could not find the physical volume with index {physical_volume_index} " + f"in volume '{self.user_info.mother}' to which this actor is attached. " + ) + align_image_with_physical_volume( + attached_to_volume, + self.py_fluence_image, + initial_translation=self.user_info.translation, + ) + + # Set the real physical volume name + self.fPhysicalVolumeName = str(self.g4_phys_vol.GetName()) + + # FIXME for multiple run and motion + if not self.first_run: + warning(f"Not implemented yet: FluenceActor with several runs") + # send itk image to cpp side, copy data only the first run. + update_image_py_to_cpp( + self.py_fluence_image, self.cpp_fluence_image, self.first_run + ) + + # now, indicate the next run will not be the first + self.first_run = False + + def EndSimulationAction(self): + g4.GateFluenceActor.EndSimulationAction(self) + + # Get the itk image from the cpp side + # Currently a copy. Maybe later as_pyarray ? + self.py_fluence_image = get_cpp_image(self.cpp_fluence_image) + + # set the property of the output image: + origin = self.img_origin_during_run + self.py_fluence_image.SetOrigin(origin) + + # write the image at the end of the run + # FIXME : maybe different for several runs + if self.user_info.output: + out_p = ensure_filename_is_str( + self.simulation.get_output_path(self.user_info.output) + ) + itk.imwrite(self.py_fluence_image, out_p) diff --git a/opengate/tests/data b/opengate/tests/data index 29f9f6523..a0e299e40 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit 29f9f6523f0d75dd878420f6e3924d38ed474adb +Subproject commit a0e299e40bb878b7bdeed76d0f9f1a6bfa993557 From a93482d6c288c7e3fe12bf365106e034a6814392 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 17:57:53 +0100 Subject: [PATCH 41/46] no options yet --- opengate/actors/doseactors.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/opengate/actors/doseactors.py b/opengate/actors/doseactors.py index 7c83a73d1..e021b8b31 100644 --- a/opengate/actors/doseactors.py +++ b/opengate/actors/doseactors.py @@ -567,6 +567,9 @@ def initialize(self, volume_engine=None): ) # for initialization during the first run self.first_run = True + # no options yet + if self.user_info.uncertainty or self.user_info.scatter: + fatal(f"FluenceActor : uncertainty and scatter not implemented yet") def StartSimulationAction(self): # init the origin and direction according to the physical volume From 1303e33a4497df06117895a37da063f6c4c61c47 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 18:13:29 +0100 Subject: [PATCH 42/46] remove debug print --- opengate/engines.py | 7 ++++--- opengate/managers.py | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/opengate/engines.py b/opengate/engines.py index e6f591f08..7fbea3432 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -1245,9 +1245,10 @@ def initialize_g4_verbose(self): else: # no Geant4 output ui = UIsessionSilent() - self.simulation.add_g4_command_after_init( - f"/tracking/verbose {self.simulation.g4_verbose_level_tracking}" - ) + if self.simulation.g4_verbose_level_tracking >= 0: + self.simulation.add_g4_command_after_init( + f"/tracking/verbose {self.simulation.g4_verbose_level_tracking}" + ) # it is also possible to set ui=None for 'default' output # we must keep a ref to ui_session self.ui_session = ui diff --git a/opengate/managers.py b/opengate/managers.py index b5b0bf165..9d32a1076 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -918,7 +918,7 @@ class Simulation(GateObject): ), "g4_verbose": (False, {"doc": "Switch on/off Geant4's verbose output."}), "g4_verbose_level_tracking": ( - 0, + -1, { "doc": "Activate verbose tracking in Geant4 via G4 command '/tracking/verbose g4_verbose_level_tracking'." }, From 51ff3c0084f8662f4f0f82682a861713307ece8c Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 19:10:59 +0100 Subject: [PATCH 43/46] remove scatter option --- .../src/test067_cbct_fluence_actor_mt.py | 102 ++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100755 opengate/tests/src/test067_cbct_fluence_actor_mt.py diff --git a/opengate/tests/src/test067_cbct_fluence_actor_mt.py b/opengate/tests/src/test067_cbct_fluence_actor_mt.py new file mode 100755 index 000000000..0d5b5bfa4 --- /dev/null +++ b/opengate/tests/src/test067_cbct_fluence_actor_mt.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "gate_test067", "test067") + + sim = gate.Simulation() + + ui = sim.user_info + # ui.visu = True + ui.visu_type = "vrml" + ui.random_engine = "MersenneTwister" + ui.random_seed = "auto" + ui.number_of_threads = 2 + + # units + m = gate.g4_units.m + cm = gate.g4_units.cm + nm = gate.g4_units.nm + cm3 = gate.g4_units.cm3 + keV = gate.g4_units.keV + MeV = gate.g4_units.MeV + mm = gate.g4_units.mm + Bq = gate.g4_units.Bq + gcm3 = gate.g4_units.g / cm3 + + # world + world = sim.world + world.size = [3 * m, 3 * m, 3 * m] + world.material = "G4_AIR" + + # waterbox + waterbox = sim.add_volume("Box", "waterbox") + waterbox.size = [20 * cm, 20 * cm, 20 * cm] + waterbox.translation = [0 * cm, 0 * cm, 15 * cm] + waterbox.material = "G4_WATER" + + # CBCT gantry source + gantry = sim.add_volume("Box", "CBCT_gantry") + gantry.size = [0.2 * m, 0.2 * m, 0.2 * m] + gantry.translation = [1000 * mm, 0, 0] + gantry.material = "G4_AIR" + gantry.color = [0, 1, 1, 1] + + # CBCT detector plane + detector_plane = sim.add_volume("Box", "CBCT_detector_plane") + detector_plane.size = [10 * mm, 1409.6 * mm, 1409.6 * mm] + detector_plane.translation = [-536 * mm, 0, 0] + detector_plane.material = "G4_AIR" + detector_plane.color = [1, 0, 0, 1] + + # physics + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option1" + sim.physics_manager.set_production_cut("world", "all", 10 * mm) + + # actor + detector_actor = sim.add_actor("FluenceActor", "detector_actor") + detector_actor.mother = detector_plane.name + detector_actor.output = paths.output / "fluence.mhd" + detector_actor.spacing = [10 * mm, 5 * mm, 5 * mm] + detector_actor.size = [1, 100, 100] + + # source + source = sim.add_source("GenericSource", "mysource") + source.mother = gantry.name + source.particle = "gamma" + source.energy.mono = 60 * keV + source.position.type = "box" + source.position.size = [1 * nm, 2 * 8 * mm, 2 * 8 * mm] + source.direction.type = "focused" + # FIXME warning in world coord ! Should be in mother coord system + source.direction.focus_point = [gantry.translation[0] - 60 * mm, 0, 0] + source.n = 50000 / ui.number_of_threads + + if ui.visu: + source.n = 100 + + # statistics + stats = sim.add_actor("SimulationStatisticsActor", "stats") + stats.track_types_flag = True + stats.output = paths.output / "stats.txt" + + # run + output = sim.run() + + # print output statistics + stats = output.get_actor("stats") + print(stats) + + # check images + is_ok = utility.assert_images( + paths.gate_output / "detector.mhd", + paths.output / "fluence.mhd", + stats, + tolerance=44, + axis="y", + ) + + utility.test_ok(is_ok) From c0e63676be19b63ff1006954f1c7ab8d73df74f7 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 18 Dec 2023 16:30:22 +0100 Subject: [PATCH 44/46] Debug according master branch --- opengate/tests/src/test067_cbct_fluence_actor_mt.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opengate/tests/src/test067_cbct_fluence_actor_mt.py b/opengate/tests/src/test067_cbct_fluence_actor_mt.py index 0d5b5bfa4..d2161e999 100755 --- a/opengate/tests/src/test067_cbct_fluence_actor_mt.py +++ b/opengate/tests/src/test067_cbct_fluence_actor_mt.py @@ -84,7 +84,8 @@ stats.output = paths.output / "stats.txt" # run - output = sim.run() + sim.run() + output = sim.output # print output statistics stats = output.get_actor("stats") From 83b937ded26d74471961b698627b10a4780876bf Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 18 Dec 2023 20:53:38 +0000 Subject: [PATCH 45/46] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/psf/black: 23.11.0 → 23.12.0](https://github.com/psf/black/compare/23.11.0...23.12.0) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index cd9239084..6abf66756 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -5,7 +5,7 @@ repos: - id: end-of-file-fixer - id: trailing-whitespace - repo: https://github.com/psf/black - rev: 23.11.0 + rev: 23.12.0 hooks: - id: black - repo: https://github.com/pre-commit/mirrors-clang-format From 619b150c7da0245280f27d9d0c52c2e043682ea9 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Tue, 19 Dec 2023 17:40:58 +0100 Subject: [PATCH 46/46] There are a lot of brew failure with macos. According to https://github.com/actions/runner-images/issues/8642 https://github.com/actions/runner-images/issues/8649 the macos runners are evolving leading to error of VM It seems there is no problem with macos-11 --- .github/workflows/main.yml | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index b9d3b6146..bff6f8a1f 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -17,10 +17,10 @@ jobs: strategy: fail-fast: false matrix: - os: [ubuntu-latest, macos-latest, windows-latest] + os: [ubuntu-latest, macos-11, windows-latest] python-version: [3.8, 3.9, '3.10', '3.11'] exclude: - - os: macos-latest + - os: macos-11 python-version: '3.11' env: @@ -36,7 +36,7 @@ jobs: if [[ ${{ matrix.os }} == "ubuntu-latest" ]]; then export GIT_SSL_NO_VERIFY=1 git submodule update --init --recursive - elif [[ ${{ matrix.os }} == "macos-latest" ]]; then + elif [[ ${{ matrix.os }} == "macos-11" ]]; then export GIT_SSL_NO_VERIFY=1 git submodule update --init --recursive else @@ -56,7 +56,7 @@ jobs: varOS=`cat /etc/os-release | grep "VERSION=" | grep -oP '(?<=\").*?(?=\")'` varOS=($varOS) echo "release=${varOS[0]}" >> $GITHUB_OUTPUT - elif [[ ${{ matrix.os }} == "macos-latest" ]]; then + elif [[ ${{ matrix.os }} == "macos-11" ]]; then varOS=`sw_vers | grep "ProductVersion:"` varOS="${varOS#*:}" echo "release=${varOS:1}" >> $GITHUB_OUTPUT @@ -106,17 +106,17 @@ jobs: mv dist_opengate/* dist/ fi - uses: conda-incubator/setup-miniconda@v2 - if: (matrix.os == 'macos-latest') || (matrix.os == 'windows-latest') + if: (matrix.os == 'macos-11') || (matrix.os == 'windows-latest') with: auto-update-conda: true activate-environment: opengate_core python-version: ${{ matrix.python-version }} - name: Set up Homebrew - if: matrix.os == 'macos-latest' + if: matrix.os == 'macos-11' id: set-up-homebrew uses: Homebrew/actions/setup-homebrew@master - name: Create opengate_core Wheel Mac - if: matrix.os == 'macos-latest' + if: matrix.os == 'macos-11' shell: bash -l {0} run: | brew update