pctweplfit update

.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

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,
     )
 

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}$.

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"

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)