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 f5b77e9..c4b8595 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 @@ -26,6 +27,7 @@ def tof_fit_mc( path_type, number_of_detectors, visu, + seed, verbose, ): import opengate as gate @@ -47,6 +49,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] @@ -87,6 +90,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): @@ -101,7 +105,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", @@ -110,10 +114,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]) @@ -122,11 +128,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) - - # Particle stats - stat = sim.add_actor("SimulationStatisticsActor", "stat") - stat.output_filename = f"{output}/wepl{int(phantom_length_mm)}_stats.txt" + add_detector(str(x), [0 * mm, 0 * mm, x * mm], True) sim.run() @@ -147,7 +149,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", @@ -334,6 +338,7 @@ def pctweplfit( path_type, phantom_length_samples, phantom_width, + max_phantom_length, detector_distance, number_of_detectors, initial_energy, @@ -342,33 +347,51 @@ def pctweplfit( visu, display, savefig, + seed, verbose, ): - phantom_lengths = np.linspace(0.0, detector_distance, phantom_length_samples) + 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(max_phantom_length, 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, @@ -415,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", @@ -462,6 +492,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", @@ -479,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, @@ -487,6 +519,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/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}$. 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..750ea8f 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,33 @@ 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 -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)) + 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)