From 36d65e618035b36a63e3d775a8b5a2b6a4655ca9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Wed, 4 Feb 2026 15:51:21 +0100 Subject: [PATCH 1/6] Add proton CT runnable example --- .../user_guide_how_to_convert_example_1.rst | 2 +- opengate/bin/protonct.py | 26 +++ opengate/contrib/protonct/__init__.py | 0 opengate/contrib/protonct/protonct.py | 166 ++++++++++++++++++ pyproject.toml | 2 + 5 files changed, 195 insertions(+), 1 deletion(-) create mode 100755 opengate/bin/protonct.py create mode 100644 opengate/contrib/protonct/__init__.py create mode 100644 opengate/contrib/protonct/protonct.py diff --git a/docs/source/user_guide/user_guide_how_to_convert_example_1.rst b/docs/source/user_guide/user_guide_how_to_convert_example_1.rst index 11356499f2..6b91a311f8 100644 --- a/docs/source/user_guide/user_guide_how_to_convert_example_1.rst +++ b/docs/source/user_guide/user_guide_how_to_convert_example_1.rst @@ -2,7 +2,7 @@ From GATE 9 to 10 - Example Proton CT ************************************* -This is an example of a proton beam with a spiral phantom placed in between two proton detectors. The spiral phantom was used in the article `Filtered backprojection proton CT reconstruction along most likely paths by Rit et al `_. The data generated by this simulation can be processed by the `PCT software `_. +This is an example of a proton beam with a spiral phantom placed in between two proton detectors. The spiral phantom was used in the article `Filtered backprojection proton CT reconstruction along most likely paths by Rit et al `_. A complete example of GATE proton CT simulation can be found `here `_, which can be run using the `protonct `_ executable. The data generated by this simulation can be processed by the `PCT software `_. .. image:: ../figures/proton_ct.png diff --git a/opengate/bin/protonct.py b/opengate/bin/protonct.py new file mode 100755 index 0000000000..e137eb6be1 --- /dev/null +++ b/opengate/bin/protonct.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import click +from opengate.contrib.protonct.protonct import protonct + +CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"]) + + +@click.command(context_settings=CONTEXT_SETTINGS) +@click.option("--output", "-o", default=None, help="output folder") +@click.option("--projections", "-p", default=720, help="number of projections") +@click.option("--protons-per-projection", + "-n", + default=1000, + help="number of protons per projection") +@click.option("--seed", "-s", type=int, help="random number generator seed") +@click.option("--visu", is_flag=True, help="Enable visualization") +@click.option("--verbose", is_flag=True, help="verbose execution") +def go(output, projections, protons_per_projection, seed, visu, verbose): + protonct(output, projections, protons_per_projection, seed, visu, verbose) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/opengate/contrib/protonct/__init__.py b/opengate/contrib/protonct/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/opengate/contrib/protonct/protonct.py b/opengate/contrib/protonct/protonct.py new file mode 100644 index 0000000000..26e5ea7d5a --- /dev/null +++ b/opengate/contrib/protonct/protonct.py @@ -0,0 +1,166 @@ +import math +from scipy.spatial.transform import Rotation + +import opengate as gate + + +def protonct(output, + projections=720, + protons_per_projection=1000, + seed=None, + visu=False, + verbose=False): + + # Units + nm = gate.g4_units.nm + mm = gate.g4_units.mm + cm = gate.g4_units.cm + m = gate.g4_units.m + sec = gate.g4_units.second + MeV = gate.g4_units.MeV + Bq = gate.g4_units.Bq + + # Simulation + sim = gate.Simulation() + + sim.random_engine = "MersenneTwister" + sim.random_seed = "auto" if seed is None else seed + sim.check_volumes_overlap = False + sim.visu = visu + sim.visu_type = "vrml" + sim.g4_verbose = False + sim.progress_bar = verbose + sim.number_of_threads = 1 + + sim.run_timing_intervals = [[j * sec, (j + 1) * sec] + for j in range(projections)] + + # Misc + yellow = [1, 1, 0, 0.5] + blue = [0, 0, 1, 0.5] + + # Geometry + sim.volume_manager.add_material_database(gate.utility.get_contrib_path() / + "GateMaterials.db") + sim.world.material = "Vacuum" + sim.world.size = [4 * m, 4 * m, 4 * m] + sim.world.color = [0, 0, 0, 0] + + # Phantom + + def add_spiral(sim): + # Mother of all + spiral = sim.add_volume("Tubs", name="Spiral") + spiral.rmin = 0 * cm + spiral.rmax = 10 * cm + spiral.dz = 40 * cm + spiral.material = "Water" + spiral.color = blue + spiral.rotation = Rotation.from_euler("yz", [90, 90], + degrees=True).as_matrix() + + # Spiral rotation + tr, rot = gate.geometry.utility.volume_orbiting_transform( + "y", 0, 360, projections, spiral.translation, spiral.rotation) + spiral.add_dynamic_parametrisation(translation=tr, rotation=rot) + + # Spiral inserts + sradius = 4 + radius = list(range(0, 100 - sradius // 2, sradius)) + sangle = 139 + angles = [ + math.radians(a) for a in range(0, sangle * len(radius), sangle) + ] + posx = [radius[i] * math.cos(angles[i]) for i in range(len(radius))] + posy = [radius[i] * math.sin(angles[i]) for i in range(len(radius))] + + def add_spiral_insert( + sim, + mother, + name, + rmin=0 * mm, + rmax=1 * mm, + dz=40 * cm, + material="Aluminium", + translation=None, + color=None, + ): + + if translation is None: + translation = [0 * mm, 0 * mm, 0 * mm] + + if color is None: + color = yellow + + spiral_insert = sim.add_volume("Tubs", name=name) + spiral_insert.mother = mother.name + spiral_insert.rmin = rmin + spiral_insert.rmax = rmax + spiral_insert.dz = dz + spiral_insert.material = material + spiral_insert.translation = translation + spiral_insert.color = color + + for i in range(len(radius)): + add_spiral_insert( + sim, + spiral, + f"SpiralInsert{i:02d}", + translation=[posx[i] * mm, posy[i] * mm, 0], + ) + + add_spiral(sim) + + # Beam + source = sim.add_source("GenericSource", "mybeam") + source.particle = "proton" + source.energy.mono = 200 * MeV + source.energy.type = "mono" + source.position.type = "box" + source.position.size = [16 * mm, 1 * nm, 1 * nm] + source.position.translation = [0, 0, -1060 * mm] + source.direction.type = "focused" + source.direction.focus_point = [0, 0, -1000 * mm] + + if sim.visu: + # For visualisation speed, the number of particle is decreased + source.activity = 10 * Bq + else: + source.activity = protons_per_projection * Bq + + # Physics list + sim.physics_manager.physics_list_name = "QGSP_BIC_EMZ" + + # Phase spaces + + def add_detector(sim, name, translation): + plane = sim.add_volume("Box", "PlanePhaseSpace" + name) + plane.size = [400 * mm, 400 * mm, 1 * nm] + plane.translation = translation + plane.material = "Vacuum" + plane.color = yellow + + phase_space = sim.add_actor("PhaseSpaceActor", "PhaseSpace" + name) + phase_space.attached_to = plane.name + phase_space.attributes = [ + "RunID", + "EventID", + "TrackID", + "KineticEnergy", + "LocalTime", + "Position", + "Direction", + ] + phase_space.output_filename = f"{output}/PhaseSpace{name}.root" + + F = gate.actors.filters.GateFilterBuilder() + phase_space.filter = F.ParticleName == "proton" + + add_detector(sim, "In", [0, 0, -110 * mm]) + add_detector(sim, "Out", [0, 0, 110 * mm]) + + # Particle stats + stat = sim.add_actor("SimulationStatisticsActor", "stat") + stat.output_filename = f"{output}/protonct.txt" + + sim.run() diff --git a/pyproject.toml b/pyproject.toml index 2b22e9b3d1..62d048853a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -60,3 +60,5 @@ phid_gammas = "opengate.bin.phid_gammas:go" phid_tac = "opengate.bin.phid_tac:go" phid_atomic_relaxation = "opengate.bin.phid_atomic_relaxation:go" phid_isomeric_transition = "opengate.bin.phid_isomeric_transition:go" + +protonct = "opengate.bin.protonct:go" From 235bd1d260a7fb897249f46588b05d61be048498 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 4 Feb 2026 14:56:50 +0000 Subject: [PATCH 2/6] [pre-commit.ci] Automatic python and c++ formatting --- opengate/bin/protonct.py | 10 +++++---- opengate/contrib/protonct/protonct.py | 32 +++++++++++++-------------- 2 files changed, 22 insertions(+), 20 deletions(-) diff --git a/opengate/bin/protonct.py b/opengate/bin/protonct.py index e137eb6be1..1ac5a2eec1 100755 --- a/opengate/bin/protonct.py +++ b/opengate/bin/protonct.py @@ -10,10 +10,12 @@ @click.command(context_settings=CONTEXT_SETTINGS) @click.option("--output", "-o", default=None, help="output folder") @click.option("--projections", "-p", default=720, help="number of projections") -@click.option("--protons-per-projection", - "-n", - default=1000, - help="number of protons per projection") +@click.option( + "--protons-per-projection", + "-n", + default=1000, + help="number of protons per projection", +) @click.option("--seed", "-s", type=int, help="random number generator seed") @click.option("--visu", is_flag=True, help="Enable visualization") @click.option("--verbose", is_flag=True, help="verbose execution") diff --git a/opengate/contrib/protonct/protonct.py b/opengate/contrib/protonct/protonct.py index 26e5ea7d5a..e31e9904f0 100644 --- a/opengate/contrib/protonct/protonct.py +++ b/opengate/contrib/protonct/protonct.py @@ -4,12 +4,14 @@ import opengate as gate -def protonct(output, - projections=720, - protons_per_projection=1000, - seed=None, - visu=False, - verbose=False): +def protonct( + output, + projections=720, + protons_per_projection=1000, + seed=None, + visu=False, + verbose=False, +): # Units nm = gate.g4_units.nm @@ -32,16 +34,16 @@ def protonct(output, sim.progress_bar = verbose sim.number_of_threads = 1 - sim.run_timing_intervals = [[j * sec, (j + 1) * sec] - for j in range(projections)] + sim.run_timing_intervals = [[j * sec, (j + 1) * sec] for j in range(projections)] # Misc yellow = [1, 1, 0, 0.5] blue = [0, 0, 1, 0.5] # Geometry - sim.volume_manager.add_material_database(gate.utility.get_contrib_path() / - "GateMaterials.db") + sim.volume_manager.add_material_database( + gate.utility.get_contrib_path() / "GateMaterials.db" + ) sim.world.material = "Vacuum" sim.world.size = [4 * m, 4 * m, 4 * m] sim.world.color = [0, 0, 0, 0] @@ -56,21 +58,19 @@ def add_spiral(sim): spiral.dz = 40 * cm spiral.material = "Water" spiral.color = blue - spiral.rotation = Rotation.from_euler("yz", [90, 90], - degrees=True).as_matrix() + spiral.rotation = Rotation.from_euler("yz", [90, 90], degrees=True).as_matrix() # Spiral rotation tr, rot = gate.geometry.utility.volume_orbiting_transform( - "y", 0, 360, projections, spiral.translation, spiral.rotation) + "y", 0, 360, projections, spiral.translation, spiral.rotation + ) spiral.add_dynamic_parametrisation(translation=tr, rotation=rot) # Spiral inserts sradius = 4 radius = list(range(0, 100 - sradius // 2, sradius)) sangle = 139 - angles = [ - math.radians(a) for a in range(0, sangle * len(radius), sangle) - ] + angles = [math.radians(a) for a in range(0, sangle * len(radius), sangle)] posx = [radius[i] * math.cos(angles[i]) for i in range(len(radius))] posy = [radius[i] * math.sin(angles[i]) for i in range(len(radius))] From 09c4f158385e0b429abdf7ccbd5f0e33424b1a28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Thu, 5 Feb 2026 10:06:00 +0100 Subject: [PATCH 3/6] Point to the RTKConsortium PCT repository --- docs/source/user_guide/user_guide_how_to_convert_example_1.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/user_guide/user_guide_how_to_convert_example_1.rst b/docs/source/user_guide/user_guide_how_to_convert_example_1.rst index 6b91a311f8..971a7a80a0 100644 --- a/docs/source/user_guide/user_guide_how_to_convert_example_1.rst +++ b/docs/source/user_guide/user_guide_how_to_convert_example_1.rst @@ -2,7 +2,7 @@ From GATE 9 to 10 - Example Proton CT ************************************* -This is an example of a proton beam with a spiral phantom placed in between two proton detectors. The spiral phantom was used in the article `Filtered backprojection proton CT reconstruction along most likely paths by Rit et al `_. A complete example of GATE proton CT simulation can be found `here `_, which can be run using the `protonct `_ executable. The data generated by this simulation can be processed by the `PCT software `_. +This is an example of a proton beam with a spiral phantom placed in between two proton detectors. The spiral phantom was used in the article `Filtered backprojection proton CT reconstruction along most likely paths by Rit et al `_. A complete example of GATE proton CT simulation can be found `here `_, which can be run using the `protonct `_ executable. The data generated by this simulation can be processed by the `PCT software `_. .. image:: ../figures/proton_ct.png From bb5f0679a4d9a1d623fef2d6aed1e5babf78cba0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Mon, 1 Jun 2026 16:14:50 +0200 Subject: [PATCH 4/6] Add proton CT test --- .../src/external/pct/test101_protonct.py | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 opengate/tests/src/external/pct/test101_protonct.py diff --git a/opengate/tests/src/external/pct/test101_protonct.py b/opengate/tests/src/external/pct/test101_protonct.py new file mode 100644 index 0000000000..16c8ebfb56 --- /dev/null +++ b/opengate/tests/src/external/pct/test101_protonct.py @@ -0,0 +1,47 @@ +import uproot +from opengate.tests import utility +from opengate.contrib.protonct.protonct import protonct + + +def compare_root_files(path1, path2, tree): + tree1 = uproot.open(path1)[tree] + df1 = tree1.arrays(library="pd") + + tree2 = uproot.open(path2)[tree] + df2 = tree2.arrays(library="pd") + + return df1.equals(df2) + + +if __name__ == "__main__": + # ===================================================== + # INITIALISATION + # ===================================================== + + paths = utility.get_default_test_paths(__file__, + output_folder="test101_protonct") + + protonct(paths.output, + projections=10, + protons_per_projection=100, + seed=1234) + + # ===================================================== + # Perform test + # ===================================================== + + path_phasespace_in = paths.output / "PhaseSpaceIn.root" + path_phasespace_out = paths.output / "PhaseSpaceOut.root" + + path_reference_phasespace_in = paths.output_ref / "PhaseSpaceIn.root" + path_reference_phasespace_out = paths.output_ref / "PhaseSpaceOut.root" + + compare_root_files(path_reference_phasespace_in, path_phasespace_in, + "PhaseSpaceIn") + + is_ok = compare_root_files(path_reference_phasespace_in, + path_phasespace_in, "PhaseSpaceIn") + is_ok = is_ok and compare_root_files(path_reference_phasespace_out, + path_phasespace_out, "PhaseSpaceOut") + + utility.test_ok(is_ok) From 0bebc494eb4bc6edc88cc868bc7c11e353562e72 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 1 Jun 2026 14:23:03 +0000 Subject: [PATCH 5/6] [pre-commit.ci] Automatic python and c++ formatting --- .../src/external/pct/test101_protonct.py | 21 ++++++++----------- 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/opengate/tests/src/external/pct/test101_protonct.py b/opengate/tests/src/external/pct/test101_protonct.py index 16c8ebfb56..203901a2a7 100644 --- a/opengate/tests/src/external/pct/test101_protonct.py +++ b/opengate/tests/src/external/pct/test101_protonct.py @@ -18,13 +18,9 @@ def compare_root_files(path1, path2, tree): # INITIALISATION # ===================================================== - paths = utility.get_default_test_paths(__file__, - output_folder="test101_protonct") + paths = utility.get_default_test_paths(__file__, output_folder="test101_protonct") - protonct(paths.output, - projections=10, - protons_per_projection=100, - seed=1234) + protonct(paths.output, projections=10, protons_per_projection=100, seed=1234) # ===================================================== # Perform test @@ -36,12 +32,13 @@ def compare_root_files(path1, path2, tree): path_reference_phasespace_in = paths.output_ref / "PhaseSpaceIn.root" path_reference_phasespace_out = paths.output_ref / "PhaseSpaceOut.root" - compare_root_files(path_reference_phasespace_in, path_phasespace_in, - "PhaseSpaceIn") + compare_root_files(path_reference_phasespace_in, path_phasespace_in, "PhaseSpaceIn") - is_ok = compare_root_files(path_reference_phasespace_in, - path_phasespace_in, "PhaseSpaceIn") - is_ok = is_ok and compare_root_files(path_reference_phasespace_out, - path_phasespace_out, "PhaseSpaceOut") + is_ok = compare_root_files( + path_reference_phasespace_in, path_phasespace_in, "PhaseSpaceIn" + ) + is_ok = is_ok and compare_root_files( + path_reference_phasespace_out, path_phasespace_out, "PhaseSpaceOut" + ) utility.test_ok(is_ok) From b775f23201fae418584892ac647fa44e7c763224 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 1 Jun 2026 16:27:27 +0200 Subject: [PATCH 6/6] Update reference data for test101 protonct --- opengate/tests/data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opengate/tests/data b/opengate/tests/data index 1f89d1fa3f..291db9cd82 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit 1f89d1fa3f9fc593ca259ce1204c925be5a96a9f +Subproject commit 291db9cd82fedb8107ddf7b543a74d3c385a68de