Add vrt into doserate

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

opengate/tests/src/source/test035_dose_rate.py → 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

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)