From 3cd6352d9ab91c239113e956bac7c3b902af2c3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Fri, 24 Apr 2026 19:41:38 +0200 Subject: [PATCH 1/6] ENH: Make pctweplfit handle opengate>=10.1.0 --- applications/pctweplfit/pctweplfit.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/applications/pctweplfit/pctweplfit.py b/applications/pctweplfit/pctweplfit.py index f5b77e9..3cc20ea 100644 --- a/applications/pctweplfit/pctweplfit.py +++ b/applications/pctweplfit/pctweplfit.py @@ -110,10 +110,12 @@ def add_detector(name, translation, attach_to_phantom=False): "PreGlobalTime", "KineticEnergy", ] - particle_filter = sim.add_filter("ParticleFilter", "Filter" + name) - particle_filter.particle = "proton" - - phase_space.filters.append(particle_filter) + if int(gate.utility.version("opengate").split(".")[1]) > 0: + F = gate.actors.filters.GateFilterBuilder() + phase_space.filter = F.ParticleName == "proton" + else: + particle_filter = sim.add_filter("ParticleFilter", "Filter" + name) + particle_filter.particle = "proton" add_detector("In", [0 * mm, 0 * mm, (-detector_distance_mm / 2 - epsilon_mm) * mm]) add_detector("Out", [0 * mm, 0 * mm, (detector_distance_mm / 2 + epsilon_mm) * mm]) From 6628ea370e1d98c423a7396694eaa27f0cf7f6c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Fri, 24 Apr 2026 19:45:18 +0200 Subject: [PATCH 2/6] ENH: Define an upper bound for phantom thickness in pctweplfit The previous behavior allowed to have values for phantom thickness larger than the proton range in water, which caused issues because no protons could traverse the phantom. Now, the maximum thickness is given by the approximate proton range in water, which is of about 26 cm. The parameter --phantom-length-samples defines the (evenly spaced) number of thicknesses to compute between 0 and the minimum between the proton range and the spacing between the detectors. --- applications/pctweplfit/pctweplfit.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/applications/pctweplfit/pctweplfit.py b/applications/pctweplfit/pctweplfit.py index 3cc20ea..4e9059f 100644 --- a/applications/pctweplfit/pctweplfit.py +++ b/applications/pctweplfit/pctweplfit.py @@ -10,6 +10,8 @@ epsilon_mm = 1e-5 +PROTON_RANGE = 260 # mm + def pv(verbose, *args, **kwargs): if verbose: @@ -124,7 +126,7 @@ def add_detector(name, translation, attach_to_phantom=False): for x in np.linspace( -phantom_length_mm / 2, phantom_length_mm / 2, number_of_detectors ): - add_detector(str(int(x)), [0 * mm, 0 * mm, x * mm], True) + add_detector(str(x), [0 * mm, 0 * mm, x * mm], True) # Particle stats stat = sim.add_actor("SimulationStatisticsActor", "stat") @@ -346,7 +348,9 @@ def pctweplfit( savefig, verbose, ): - phantom_lengths = np.linspace(0.0, detector_distance, phantom_length_samples) + phantom_lengths = np.linspace( + 0.0, min(PROTON_RANGE, detector_distance), phantom_length_samples + ) results = [] with Pool(maxtasksperchild=1) as pool: From 202a0e8df312e4c9eb22c7f786f876680a17e2ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Mon, 27 Apr 2026 09:39:40 +0200 Subject: [PATCH 3/6] ENH: Add test for pctweplfit --- .github/workflows/build-test-package.yml | 5 +- applications/pctweplfit/pctweplfit.py | 61 +++++++++++++++--------- pyproject.toml | 2 +- test/pct_application_test.py | 30 ++++++++++++ 4 files changed, 74 insertions(+), 24 deletions(-) diff --git a/.github/workflows/build-test-package.yml b/.github/workflows/build-test-package.yml index 4c76fe4..1baa555 100644 --- a/.github/workflows/build-test-package.yml +++ b/.github/workflows/build-test-package.yml @@ -51,7 +51,10 @@ jobs: echo "Installing wheel: $wheel" pip install $wheel - pip install pytest uproot + pip install pytest uproot opengate + + # Force the installation of Geant4 data, required by opengate + opengate_info # Run tests with warning detection. # SWIG < 4.4 triggers the Py_LIMITED_API "builtin type swig…" warning when diff --git a/applications/pctweplfit/pctweplfit.py b/applications/pctweplfit/pctweplfit.py index 4e9059f..af6f8a4 100644 --- a/applications/pctweplfit/pctweplfit.py +++ b/applications/pctweplfit/pctweplfit.py @@ -3,6 +3,7 @@ import json import sys from multiprocessing import Pool, Manager +import warnings import numpy as np import itk @@ -28,6 +29,7 @@ def tof_fit_mc( path_type, number_of_detectors, visu, + seed, verbose, ): import opengate as gate @@ -49,6 +51,7 @@ def tof_fit_mc( sim.g4_verbose = False sim.progress_bar = verbose sim.number_of_threads = 1 + sim.random_seed = seed # Misc yellow = [1, 1, 0, 1] @@ -89,6 +92,7 @@ def tof_fit_mc( # Physics list sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option4" + sim.physics_manager.set_user_limits_particles(["proton"]) # Phase spaces def add_detector(name, translation, attach_to_phantom=False): @@ -346,35 +350,46 @@ def pctweplfit( visu, display, savefig, + seed, verbose, ): phantom_lengths = np.linspace( 0.0, min(PROTON_RANGE, detector_distance), phantom_length_samples ) + seed_seq = np.random.SeedSequence(seed) + seeds = seed_seq.generate_state(len(phantom_lengths)) + results = [] - with Pool(maxtasksperchild=1) as pool: - for phantom_length in phantom_lengths: - result = pool.apply_async( - tof_fit_mc, - ( - phantom_length, - output, - phantom_width, - detector_distance, - number_of_particles, - initial_energy, - path_type, - number_of_detectors, - visu, - verbose, - ), - ) - results.append(result) - pool.close() - pool.join() - for result in results: - result.get() + + # Catch DeprecationWarnings coming from opengate and making PCT tests fail + # To be removed once the warnings are gone from opengate + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=DeprecationWarning) + + with Pool(maxtasksperchild=1) as pool: + for phantom_length, seed_phantom_length in zip(phantom_lengths, seeds): + result = pool.apply_async( + tof_fit_mc, + ( + phantom_length, + output, + phantom_width, + detector_distance, + number_of_particles, + initial_energy, + path_type, + number_of_detectors, + visu, + seed_phantom_length, + verbose, + ), + ) + results.append(result) + pool.close() + pool.join() + for result in results: + result.get() tof_fit( phantom_lengths, @@ -468,6 +483,7 @@ def build_parser(): help="Write polynomial fit plot to disk", action="store_true", ) + parser.add_argument("--seed", help="Seed for random number generator", type=int) parser.add_argument( "--verbose", "-v", @@ -493,6 +509,7 @@ def process(args_info: argparse.Namespace): visu=args_info.visu, display=args_info.display, savefig=args_info.savefig, + seed=args_info.seed, verbose=args_info.verbose, ) diff --git a/pyproject.toml b/pyproject.toml index 68ecc1c..086950d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,7 +38,7 @@ dependencies = [ "itk-rtk == 2.7.*" ] [dependency-groups] -test = ["pytest", "uproot"] +test = ["pytest", "uproot", "opengate"] [project.scripts] pctbinning = "itk.pctbinning:main" diff --git a/test/pct_application_test.py b/test/pct_application_test.py index 0314b0a..75536ef 100644 --- a/test/pct_application_test.py +++ b/test/pct_application_test.py @@ -1,3 +1,4 @@ +import json import pytest import itk import urllib.request @@ -43,3 +44,32 @@ def test_pairprotons_application( pairs_test = itk.array_from_image(itk.imread(output0000)) pairs_baseline = itk.array_from_image(itk.imread(baseline_pairs_mhd)) assert np.array_equal(pairs_test, pairs_baseline) + + +def test_weplfit_application(tmp_path): + output = tmp_path / "weplfit" + pct.pctweplfit( + f"-o {output} --path-type phantom_length -d 220 -e 200 --seed 1234 -v") + + with open(output / "tof_to_wepl_fit_deg3.json", encoding="utf-8") as f: + tof_to_wepl_fit = np.array(json.load(f)) + reference = np.array( + [ + 8507.18830081492, + -38342.09213481464, + 58268.25904916112, + -29633.875319259325, + ] + ) + assert np.allclose(tof_to_wepl_fit, reference) + with open(output / "eloss_to_wepl_fit_deg3.json", encoding="utf-8") as f: + eloss_to_wepl_fit = np.array(json.load(f)) + reference = np.array( + [ + -5.4569381663242e-06, + -0.003455118013641362, + 2.236667157315751, + -0.1768424416075149, + ] + ) + assert np.allclose(eloss_to_wepl_fit, reference) From 6b4a62570579efa357ae8fec8004d362df2a6578 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Mon, 27 Apr 2026 09:44:40 +0200 Subject: [PATCH 4/6] DOC: Make pctweplfit documentation consistent The text claimed that the examples uses energy-loss fit, whereas the command used TOF fit. --- documentation/docs/wepl_calibration.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/documentation/docs/wepl_calibration.md b/documentation/docs/wepl_calibration.md index 16300ad..512ad27 100644 --- a/documentation/docs/wepl_calibration.md +++ b/documentation/docs/wepl_calibration.md @@ -20,8 +20,8 @@ pctpairprotons \ -o pairs.mhd \ --plane-in -110 \ --plane-out 110 \ - --fit output/tof_to_wepl_fit_deg3.json \ - --fit-kind tof + --fit output/eloss_to_wepl_fit_deg3.json \ + --fit-kind energy ``` The resulting pairs are directly associated to the corresponding WEPL following the convention explained in the [PCT data format](pct_format.md), that is that $e_\text{in}=0$ and $e_\text{out}=\text{WEPL}$. From ddc545fe585aea35d348da37b680ada24c9dff76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Mon, 27 Apr 2026 10:44:39 +0200 Subject: [PATCH 5/6] ENH: Cleaner output for pctweplfit --- applications/pctweplfit/pctweplfit.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/applications/pctweplfit/pctweplfit.py b/applications/pctweplfit/pctweplfit.py index af6f8a4..1d2e37c 100644 --- a/applications/pctweplfit/pctweplfit.py +++ b/applications/pctweplfit/pctweplfit.py @@ -107,7 +107,7 @@ def add_detector(name, translation, attach_to_phantom=False): phase_space = sim.add_actor("PhaseSpaceActor", "PhaseSpace" + name) phase_space.attached_to = plane.name phase_space.output_filename = ( - f"{output}/l{int(phantom_length_mm)}_ps{name}.root" + f"{output}/output/l{int(phantom_length_mm)}_ps{name}.root" ) phase_space.attributes = [ "EventID", @@ -132,10 +132,6 @@ def add_detector(name, translation, attach_to_phantom=False): ): add_detector(str(x), [0 * mm, 0 * mm, x * mm], True) - # Particle stats - stat = sim.add_actor("SimulationStatisticsActor", "stat") - stat.output_filename = f"{output}/wepl{int(phantom_length_mm)}_stats.txt" - sim.run() @@ -155,7 +151,9 @@ def process_phantom_length( wepls_phantom_length = [] elosses_phantom_length = [] - data = uproot.concatenate(f"{output}/l{int(phantom_length)}_*.root", library="np") + data = uproot.concatenate( + f"{output}/output/l{int(phantom_length)}_*.root", library="np" + ) pv( verbose, "Loaded", From 00d351ccfd578ce5a1d89bbe3d4df92348257e3a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Thu, 30 Apr 2026 11:27:10 +0200 Subject: [PATCH 6/6] ENH: Add parameter to specify maximum phantom length in pctweplfit --- applications/pctweplfit/pctweplfit.py | 18 +++++++++++++++--- test/pct_application_test.py | 3 ++- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/applications/pctweplfit/pctweplfit.py b/applications/pctweplfit/pctweplfit.py index 1d2e37c..c4b8595 100644 --- a/applications/pctweplfit/pctweplfit.py +++ b/applications/pctweplfit/pctweplfit.py @@ -11,8 +11,6 @@ epsilon_mm = 1e-5 -PROTON_RANGE = 260 # mm - def pv(verbose, *args, **kwargs): if verbose: @@ -340,6 +338,7 @@ def pctweplfit( path_type, phantom_length_samples, phantom_width, + max_phantom_length, detector_distance, number_of_detectors, initial_energy, @@ -351,8 +350,13 @@ def pctweplfit( seed, verbose, ): + if detector_distance < max_phantom_length: + print( + f"Warning, detector distance of {detector_distance} mm is smaller than the maximum phantom length of {max_phantom_length} mm. A maximum phantom length of {max_phantom_length} will be used.", + file=sys.stderr, + ) phantom_lengths = np.linspace( - 0.0, min(PROTON_RANGE, detector_distance), phantom_length_samples + 0.0, min(max_phantom_length, detector_distance), phantom_length_samples ) seed_seq = np.random.SeedSequence(seed) @@ -434,6 +438,13 @@ def build_parser(): default=40.0, type=float, ) + parser.add_argument( + "-l", + "--max-phantom-length", + help="Maximum phantom length to consider (defaults to proton range in water)", + type=float, + default=260.0, + ) parser.add_argument( "-d", "--detector-distance", @@ -499,6 +510,7 @@ def process(args_info: argparse.Namespace): path_type=args_info.path_type, phantom_length_samples=args_info.phantom_length_samples, phantom_width=args_info.phantom_width, + max_phantom_length=args_info.max_phantom_length, detector_distance=args_info.detector_distance, number_of_detectors=args_info.number_of_detectors, initial_energy=args_info.initial_energy, diff --git a/test/pct_application_test.py b/test/pct_application_test.py index 75536ef..750ea8f 100644 --- a/test/pct_application_test.py +++ b/test/pct_application_test.py @@ -49,7 +49,8 @@ def test_pairprotons_application( def test_weplfit_application(tmp_path): output = tmp_path / "weplfit" pct.pctweplfit( - f"-o {output} --path-type phantom_length -d 220 -e 200 --seed 1234 -v") + f"-o {output} --path-type phantom_length -d 220 -e 200 -l 220 --seed 1234 -v" + ) with open(output / "tof_to_wepl_fit_deg3.json", encoding="utf-8") as f: tof_to_wepl_fit = np.array(json.load(f))