From cea02e84c00334237f48a81fac40b61dd979ccb5 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Thu, 7 May 2026 16:29:05 +0200 Subject: [PATCH] Add vrt into doserate --- opengate/contrib/dose/doserate.py | 27 ++++++- ...035_dose_rate.py => test035a_dose_rate.py} | 9 ++- .../src/source/test035b_dose_rate_vrt.py | 75 +++++++++++++++++++ 3 files changed, 105 insertions(+), 6 deletions(-) rename opengate/tests/src/source/{test035_dose_rate.py => test035a_dose_rate.py} (94%) create mode 100644 opengate/tests/src/source/test035b_dose_rate_vrt.py diff --git a/opengate/contrib/dose/doserate.py b/opengate/contrib/dose/doserate.py index cb4776488..ab99ef580 100755 --- a/opengate/contrib/dose/doserate.py +++ b/opengate/contrib/dose/doserate.py @@ -2,11 +2,13 @@ # -*- coding: utf-8 -*- import pathlib -from opengate.managers import Simulation -from opengate.utility import g4_units, g4_best_unit + +from opengate.geometry.materials import HounsfieldUnit_to_material from opengate.image import get_translation_between_images_center, read_image_info from opengate.logger import INFO -from opengate.geometry.materials import HounsfieldUnit_to_material +from opengate.managers import Simulation +from opengate.sources.utility import set_source_energy_spectrum +from opengate.utility import g4_best_unit, g4_units def create_simulation(param): @@ -102,9 +104,26 @@ def create_simulation(param): sim.physics_manager.set_production_cut("world", "all", 1 * m) sim.physics_manager.set_production_cut("ct", "all", 2 * mm) + if param.mode == "e-": + # electron source + source.particle = "e-" + set_source_energy_spectrum(source, param.radionuclide) + sim.physics_manager.set_production_cut("ct", "all", 1 * m) + elif "gamma" in param.mode: + # electron source + source.particle = "gamma" + set_source_energy_spectrum(source, param.radionuclide) + sim.physics_manager.set_production_cut("ct", "all", 1 * m) + # add dose actor (get the same size as the source) source_info = read_image_info(param.activity_image) - dose = sim.add_actor("DoseActor", "dose") + if param.mode == "gamma_tle": + dose = sim.add_actor("TLEDoseActor", "dose") + dose.tle_threshold_type = "energy" + dose.tle_threshold = source.energy.spectrum_energies[-1] + print(f"tle_threshold = {dose.tle_threshold/keV} keV") + else: + dose = sim.add_actor("DoseActor", "dose") dose.output_filename = "edep.mhd" dose.dose_uncertainty.active = True dose.dose_squared.active = True diff --git a/opengate/tests/src/source/test035_dose_rate.py b/opengate/tests/src/source/test035a_dose_rate.py similarity index 94% rename from opengate/tests/src/source/test035_dose_rate.py rename to opengate/tests/src/source/test035a_dose_rate.py index 992a82231..896e36953 100755 --- a/opengate/tests/src/source/test035_dose_rate.py +++ b/opengate/tests/src/source/test035a_dose_rate.py @@ -1,13 +1,17 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -import opengate as gate +from pathlib import Path + from box import Box + +import opengate as gate from opengate.contrib.dose.doserate import create_simulation from opengate.tests import utility if __name__ == "__main__": - paths = utility.get_default_test_paths(__file__, "", output_folder="test035") + paths = utility.get_default_test_paths(__file__, "", output_folder="test035a") + paths.output_ref = Path(str(paths.output_ref)[:-1]) dr_data = paths.data / "dose_rate_data" # set param @@ -24,6 +28,7 @@ param.verbose = True param.density_tolerance_gcm3 = 0.05 param.output_folder = str(paths.output) + param.mode = "" # Create the simu # Note that the returned sim object can be modified to change source or cuts or whatever other parameters diff --git a/opengate/tests/src/source/test035b_dose_rate_vrt.py b/opengate/tests/src/source/test035b_dose_rate_vrt.py new file mode 100644 index 000000000..43c2a639d --- /dev/null +++ b/opengate/tests/src/source/test035b_dose_rate_vrt.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from pathlib import Path + +import itk +from box import Box + +import opengate as gate +from opengate.contrib.dose.doserate import create_simulation +from opengate.tests import utility + + +def run_simulation_vrt(vrt_mode=""): + paths = utility.get_default_test_paths( + __file__, "", output_folder="test035b_" + vrt_mode + ) + paths.output_ref = Path(str(paths.output_ref)[:-1]) + dr_data = paths.data / "dose_rate_data" + + # set param + gcm3 = gate.g4_units.g_cm3 + param = Box() + param.ct_image = str(dr_data / "29_CT_5mm_crop.mhd") + param.table_mat = str(dr_data / "Schneider2000MaterialsTable.txt") + param.table_density = str(dr_data / "Schneider2000DensitiesTable.txt") + param.activity_image = str(dr_data / "activity_test_crop_4mm.mhd") + param.radionuclide = "Lu177" + param.activity_bq = 5e5 + param.number_of_threads = 1 + param.visu = False + param.verbose = True + param.density_tolerance_gcm3 = 0.05 + param.output_folder = str(paths.output) + param.mode = vrt_mode + + # Create the simu + # Note that the returned sim object can be modified to change source or cuts or whatever other parameters + sim = create_simulation(param) + + # stats + stats = sim.get_actor("Stats") + stats.output_filename = "stats035_" + vrt_mode + ".txt" + + print("Phys list cuts:") + print(sim.physics_manager.dump_production_cuts()) + + # run + sim.run(start_new_process=True) + return paths + + +if __name__ == "__main__": + paths_original_simu = run_simulation_vrt(vrt_mode="") + paths_vrt_e_simu = run_simulation_vrt(vrt_mode="e-") + paths_vrt_gamma_simu = run_simulation_vrt(vrt_mode="gamma") + + # dose comparison between original simu and vrt simu (electron + gamma) + print() + gate.exception.warning(f"Check dose") + dose_vrt_e_simu = itk.imread(paths_vrt_e_simu.output / "edep_edep.mhd") + dose_vrt_gamma_simu = itk.imread(paths_vrt_gamma_simu.output / "edep_edep.mhd") + array_vrt_e_simu = itk.GetArrayFromImage(dose_vrt_e_simu) + array_vrt_gamma_simu = itk.GetArrayFromImage(dose_vrt_gamma_simu) + array_vrt = array_vrt_e_simu + array_vrt_gamma_simu + dose_vrt = itk.GetImageFromArray(array_vrt) + dose_vrt.CopyInformation(dose_vrt_e_simu) + itk.imwrite(dose_vrt, paths_original_simu.output / "edep_edep_vrt.mhd") + is_ok = utility.assert_images( + paths_original_simu.output / "edep_edep.mhd", + paths_original_simu.output / "edep_edep_vrt.mhd", + tolerance=30, + ignore_value_data2=0, + ) + utility.test_ok(is_ok)