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 param.mode == "gamma_tle":
+        # 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,8 +1,9 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
-import opengate as gate
 from box import Box
+
+import opengate as gate
 from opengate.contrib.dose.doserate import create_simulation
 from opengate.tests import utility
 
@@ -24,6 +25,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_e-.py

@@ -0,0 +1,72 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+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")
+    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 = 1e6
+    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 = "e-"
+
+    # 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.txt"
+
+    # Change source to alpha to get quick high local dose
+    source = sim.get_source_user_info("vox")
+    source.particle = "e-"
+    MeV = gate.g4_units.MeV
+    source.energy.mono = 1 * MeV
+
+    print("Phys list cuts:")
+    print(sim.physics_manager.dump_production_cuts())
+
+    # run
+    sim.run(start_new_process=True)
+
+    # print results
+    print()
+    print(stats)
+    stats_ref = utility.read_stats_file(paths.output_ref / "stats.txt")
+    is_ok = utility.assert_stats(stats, stats_ref, 0.10)
+
+    # dose comparison
+    print()
+    gate.exception.warning(f"Check dose")
+    h = sim.get_actor("dose")
+    print(h)
+    is_ok = (
+        utility.assert_images(
+            paths.output_ref / "edep.mhd",
+            h.edep.get_output_path(),
+            stats,
+            tolerance=15,
+            ignore_value_data2=0,
+        )
+        and is_ok
+    )
+
+    utility.test_ok(is_ok)

opengate/tests/src/source/test035c_dose_rate_vrt_gamma.py

@@ -0,0 +1,72 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+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")
+    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 = 1e6
+    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 = "gamma_tle"
+
+    # 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.txt"
+
+    # Change source to alpha to get quick high local dose
+    source = sim.get_source_user_info("vox")
+    source.particle = "gamma"
+    MeV = gate.g4_units.MeV
+    source.energy.mono = 1 * MeV
+
+    print("Phys list cuts:")
+    print(sim.physics_manager.dump_production_cuts())
+
+    # run
+    sim.run(start_new_process=True)
+
+    # print results
+    print()
+    print(stats)
+    stats_ref = utility.read_stats_file(paths.output_ref / "stats.txt")
+    is_ok = utility.assert_stats(stats, stats_ref, 0.10)
+
+    # dose comparison
+    print()
+    gate.exception.warning(f"Check dose")
+    h = sim.get_actor("dose")
+    print(h)
+    is_ok = (
+        utility.assert_images(
+            paths.output_ref / "edep.mhd",
+            h.edep.get_output_path(),
+            stats,
+            tolerance=15,
+            ignore_value_data2=0,
+        )
+        and is_ok
+    )
+
+    utility.test_ok(is_ok)