Created
November 23, 2022 15:00
-
-
Save jizhang02/6778a73b5ba245aa404eb371584b049f to your computer and use it in GitHub Desktop.
opengate test files for simulating vivo radiation therapy
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ''' | |
| ----------------------------------------------- | |
| File Name: test000_test$ | |
| Description: test file with opengate on patient | |
| Author: Jing$ | |
| Date: 25/10/2022$ | |
| ----------------------------------------------- | |
| ''' | |
| # import libraries | |
| import opengate as gate | |
| import numpy as np | |
| import opengate_core as g4 | |
| from scipy.spatial.transform import Rotation # used for describe rotation matrix | |
| import pathlib | |
| import os | |
| ''' | |
| General | |
| ''' | |
| # path | |
| current_path = pathlib.Path(__file__).parent.resolve() | |
| data_path = current_path / '..' / 'data' | |
| output_path = current_path / '..' / 'output' | |
| # Main Simulation object | |
| sim = gate.Simulation() | |
| ui = sim.user_info | |
| ui.g4_verbose = False # set detailed mode | |
| ui.g4_verbose_level = gate.INFO# or gate.DEBUG | |
| ui.running_verbose_level = gate.RUN # or gate.EVENT | |
| ui.visu = False | |
| ui.number_of_threads = 1 | |
| ui.random_engine = 'MersenneTwister' | |
| ui.random_seed = 123654 # to make sure every run the same | |
| # add volume type database, defining the properties of the materials | |
| sim.add_material_database(current_path / '..' / 'data' / 'GateMaterials.db') | |
| # units | |
| Bq = gate.g4_units('Bq') # Becquerel radioactivity | |
| m = gate.g4_units('m') | |
| cm = gate.g4_units('cm') | |
| mm = gate.g4_units('mm') | |
| um = gate.g4_units("um") | |
| keV = gate.g4_units('keV') | |
| MeV = gate.g4_units('MeV') | |
| sec = gate.g4_units('second') | |
| gcm3 = gate.g4_units("g/cm3") | |
| ''' | |
| Geometry | |
| ''' | |
| # world | |
| world = sim.world | |
| world.size = [3 * m, 3 * m, 3 * m] | |
| world.material = "G4_AIR" | |
| patient = sim.add_volume('Image', 'patient') | |
| patient.image = current_path / '..' / 'data' / 'test_simple' / 'patient.mhd' | |
| patient.mother = 'world' # connect the geometry to a system | |
| f1 = str(current_path/'..'/'data'/'dose_rate_data'/"Schneider2000MaterialsTable.txt") | |
| f2 = str(current_path/'..'/'data'/'dose_rate_data'/"Schneider2000DensitiesTable.txt") | |
| tol = 0.05 * gcm3 | |
| patient.voxel_materials, materials = gate.HounsfieldUnit_to_material(tol, f1, f2) | |
| print(f"tol = {tol} g/cm3") | |
| print(f"mat : {len(patient.voxel_materials)} materials") | |
| # write the corresponding labeled image (every voxel value is replaced by the material label) | |
| patient.dump_label_image = output_path / "test000_volume_label.mhd" | |
| ''' | |
| Physics | |
| ''' | |
| p = sim.get_physics_user_info() | |
| p.physics_list_name = "G4EmStandardPhysics_option3" # for medical and space science applications | |
| p.enable_decay = True | |
| p.apply_cuts = True | |
| cuts = p.production_cuts # also tracking | |
| cuts.patient.gamma = 100 * um | |
| cuts.patient.electron = 100 * um | |
| cuts.patient.positron = 100 * um | |
| cuts.patient.proton = 100 * um | |
| ''' | |
| Sources | |
| ''' | |
| # ion source | |
| v_src = sim.add_source("Voxels", "voxel") | |
| v_src.mother = "patient" | |
| v_src.particle = "ion 71 177" # Lu177 | |
| v_src.image = current_path / '..' / 'data'/'test_simple' / 'dosemap.mhd' | |
| v_src.position.type = "sphere" | |
| v_src.position.type = 10 *um | |
| v_src.position.translation = gate.get_translation_between_images_center(str(patient.image), str(v_src.image)) | |
| v_src.half_life = 574560 * sec | |
| v_src.direction.type = "iso" | |
| #v_src.energy.type = "mono" | |
| #v_src.energy.mono = 0 | |
| #v_src.activity = 100000 * Bq / ui.number_of_threads | |
| v_src.energy.type = "Lu177" | |
| total_yield = gate.get_rad_yield("Lu177") # 0.997 | |
| v_src.activity = 10000 * Bq * total_yield / ui.number_of_threads | |
| #gate.set_source_rad_energy_spectrum(v_src,"Lu177") # "spectrum" mode | |
| # question: is the energy.type = "Lu177" the same as energy.type = "spectrum"? | |
| ''' | |
| Actors | |
| ''' | |
| # add DoseActor | |
| dose = sim.add_actor('DoseActor', 'dose') | |
| dose.output = output_path / 'test000-edep.mhd' # dose, uncertainty energy map | |
| dose.mother = 'patient' | |
| dose.size = [128, 128, 203] | |
| dose.spacing = [5.078 * mm, 5.078 * mm, 5.078 * mm] | |
| dose.img_coord_system = True # default is True | |
| dose.uncertainty = True | |
| dose.gray = True | |
| # add stat actor | |
| stats = sim.add_actor('SimulationStatisticsActor', 'Stats') | |
| stats.mother = 'world' | |
| stats.track_types_flag = True | |
| # run timing | |
| #sim.run_timing_intervals = [[0, 0.50 * sec]] # By default, 1 second. | |
| sim.initialize() # initialize simulation | |
| # verbose | |
| #sim.set_g4_verbose(True) | |
| sim.apply_g4_command("/run/verbose 2") | |
| sim.apply_g4_command("/tracking/verbose 2") | |
| sim.apply_g4_command("/event/verbose 2") | |
| sim.start() # start simulation | |
| # print results at the end | |
| stats = sim.get_actor('Stats') | |
| dose = sim.get_actor('dose') | |
| print(stats) | |
| stats.write(output_path/'test000stats.txt') |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment