GE Signa config for castor (PR#878)

opengate/tests/src/external/castor/test096_pet_castor_coinc_ge.py

@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from test096_pet_castor_helpers import *
+from opengate.actors.coincidences import CoincidenceSorter
+import os
+import sys
+
+if __name__ == "__main__":
+    paths = utility.get_default_test_paths(
+        __file__, gate_folder="", output_folder="test096_pet_castor_interface"
+    )
+
+    # check or create the root file
+    root_filename = paths.output / "output_ref.root"
+    if not os.path.exists(root_filename):
+        dependency = "test096_pet_castor.py"
+        # ignore on windows
+        if os.name == "nt":
+            utility.test_ok(True)
+            sys.exit(0)
+        cmd = "python " + str(paths.current / dependency)
+        r = os.system(cmd)
+
+    # open the root file
+    print(f"Opening {root_filename}")
+    root_file = uproot.open(root_filename)
+    root_folder = root_filename.parent
+
+    # consider the singles and hits trees
+    hits_tree = root_file["hits"]
+    n = int(hits_tree.num_entries)
+    print(f"There are {n} hits")
+
+    singles_tree = root_file["singles"]
+    n = int(singles_tree.num_entries)
+    print(f"There are {n} singles")
+
+    # time windows
+    ns = gate.g4_units.nanosecond
+    mm = gate.g4_units.mm
+
+    policy = "TakeAllGoods"
+
+    sorter = CoincidenceSorter()
+    sorter.window = 4.57 * ns
+    sorter.multiples_policy = policy
+    sorter.transaxial_plane = "XY"
+    sorter.min_transaxial_distance = 20 * mm
+    sorter.max_axial_distance = 300 * mm
+
+    coincidences = sorter.run(root_filename, "singles")
+
+    nc = len(coincidences["GlobalTime1"])
+    print(f"There are {nc} coincidences for the policy", policy)
+
+    # save to file
+    output_filename = root_folder / "coincidences.root"
+
+    hits_data = root_tree_get_branch_data(hits_tree)
+    singles_data = root_tree_get_branch_data(singles_tree, library="ak")
+    coincidences_data = root_tree_get_branch_data(coincidences)
+
+    hits_types = root_tree_get_branch_types(hits_data)
+    singles_types = root_tree_get_branch_types(singles_data)
+    coincidences_types = root_tree_get_branch_types(coincidences_data)
+
+    with uproot.recreate(output_filename) as output_file:
+        root_write_tree(output_file, "hits", hits_types, hits_data)
+        root_write_tree(output_file, "singles", singles_types, singles_data)
+        root_write_tree(
+            output_file, "coincidences", coincidences_types, coincidences_data
+        )
+
+    print(f"File {output_filename} saved")

opengate/tests/src/external/castor/test096_pet_castor_ge.py

@@ -0,0 +1,119 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import opengate.contrib.pet.ge_signa as pet_ge_signa
+import opengate.contrib.pet.castor_helpers as castor
+from scipy.spatial.transform import Rotation
+from test096_pet_castor_helpers import *
+
+if __name__ == "__main__":
+    paths = utility.get_default_test_paths(
+        __file__, gate_folder="", output_folder="test096_pet_castor_interface"
+    )
+
+    # folders
+    output_path = paths.output
+
+    # units
+    m = gate.g4_units.m
+    mm = gate.g4_units.mm
+    cm = gate.g4_units.cm
+    sec = gate.g4_units.s
+    ns = gate.g4_units.ns
+    ps = gate.g4_units.ps
+    keV = gate.g4_units.keV
+    Bq = gate.g4_units.Bq
+    MBq = 1e6 * gate.g4_units.Bq
+    gcm3 = gate.g4_units.g_cm3
+    deg = gate.g4_units.deg
+
+    # options
+    sim = gate.Simulation()
+    # sim.visu = True
+    sim.visu_type = "qt"
+    sim.random_seed = 123456
+    sim.number_of_threads = 1
+    sim.output_dir = output_path
+    # sim.progress_bar = True
+    sim.verbose_level = gate.logger.NONE
+
+    # world
+    world = sim.world
+    world.size = [1.0 * m, 1.0 * m, 0.5 * m]
+    world.material = "G4_AIR"
+
+    # create the pet
+    pet = pet_ge_signa.add_pet(sim, "pet")
+
+    # get the crystal volume
+    crystal = sim.volume_manager.get_volume("crystal")
+    optical = sim.volume_manager.get_volume("optical")
+    block = sim.volume_manager.get_volume("block")
+    unit = sim.volume_manager.get_volume("unit")
+    module = sim.volume_manager.get_volume("module")
+
+    n_crystal = len(crystal.translation)
+    n_optical = len(optical.translation)
+    n_block = len(block.translation)
+    n_unit = len(unit.translation)
+    n_module = len(module.translation)
+
+    # set a digitizer
+    output_root = output_path / "output_ref.root"
+    singles = pet_ge_signa.add_digitizer(
+        sim, pet.name, output_root, singles_name="singles"
+    )
+
+    # source and physics
+    stats = test_add_physics_and_stats(sim, "pet")
+    activity = 1e3 * Bq / sim.number_of_threads
+    if sim.visu:
+        activity = 50 * Bq
+    test_add_b2b_source(sim, activity)
+
+    # Set the function that will create the file for castor.
+    # This function is a hook because it must be run once the geometry is built by geant4.
+    # The param structure contains input parameters and will contain the output castor_config
+    # after the simulation
+    filename = paths.output / "castor_config.json"
+    param = castor.set_hook_castor_config(sim, crystal.name, filename)
+
+    # go
+    duration = 1 * sec
+    print(f"Run with {activity/Bq:.0f} Bq during {duration/sec:.03f} s")
+    print(f"Mean time between events = {((1*sec)/(activity/Bq))/ns:.3f} ns")
+    sim.run_timing_intervals = [[0, duration]]
+    sim.run()
+
+    # print info
+    events = stats.counts.events
+    print(f"Number of simulated events : {events}")
+    print(f"Number of crystals : {n_crystal}")
+    print(f"Number of opticals     : {n_optical}")
+    print(f"Number of blocks   : {n_block}")
+    print(f"Number of units   : {n_unit}")
+    print(f"Number of modules  : {n_module}")
+    total_crystals = n_crystal * n_optical * n_block * n_unit * n_module
+
+    # get the castor_config
+    castor_config = param["castor_config"]
+    n = len(castor_config["unique_volume_id"])
+    is_ok = n == total_crystals
+    utility.print_test(is_ok, f"The number of volumes is {n} vs {total_crystals}")
+
+    """
+        compare the hit's positions and the volume transformation in the castor file
+    """
+
+    # 1) read all hits/singles in the root file
+    fn = singles.get_output_path()
+    print()
+    # is_ok = assert_positions(fn, "hits", castor_config, check_pos=True)
+    # print()
+    # don't check the position for the singles as Position is replaced by weighted centroid
+    is_ok = assert_positions(fn, "singles", castor_config) and is_ok
+    # print()
+    # don't check the position for the singles as Position is replaced by blured position
+    # is_ok = assert_positions(fn, "singles_blur", castor_config) and is_ok
+
+    utility.test_ok(is_ok)

opengate/tests/src/external/castor/test096_pet_castor_ref_ge.py

@@ -0,0 +1,106 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import opengate.contrib.pet.ge_signa as pet_ge_signa
+import opengate.contrib.pet.castor_helpers as castor
+from scipy.spatial.transform import Rotation
+from test096_pet_castor_helpers import *
+import os
+
+if __name__ == "__main__":
+    paths = utility.get_default_test_paths(
+        __file__, gate_folder="", output_folder="test096_pet_castor_interface"
+    )
+
+    # folders
+    output_path = paths.output
+
+    # units
+    m = gate.g4_units.m
+    mm = gate.g4_units.mm
+    cm = gate.g4_units.cm
+    sec = gate.g4_units.s
+    ns = gate.g4_units.ns
+    ps = gate.g4_units.ps
+    keV = gate.g4_units.keV
+    Bq = gate.g4_units.Bq
+    MBq = 1e6 * Bq
+    gcm3 = gate.g4_units.g_cm3
+    deg = gate.g4_units.deg
+
+    # options
+    sim = gate.Simulation()
+    # sim.visu = True
+    sim.visu_type = "qt"
+    sim.random_seed = "auto"
+    sim.number_of_threads = 4
+    sim.output_dir = output_path
+    # sim.progress_bar = True
+    sim.verbose_level = gate.logger.NONE
+
+    # world
+    world = sim.world
+    world.size = [1.5 * m, 1.5 * m, 1.5 * m]
+    world.material = "G4_AIR"
+
+    # create the pet and move it
+    pet = pet_ge_signa.add_pet(sim, "pet")
+
+    # get the crystal volume
+    crystal = sim.volume_manager.get_volume("crystal")
+    optical = sim.volume_manager.get_volume("optical")
+    block = sim.volume_manager.get_volume("block")
+    unit = sim.volume_manager.get_volume("unit")
+    module = sim.volume_manager.get_volume("module")
+
+    n_crystal = len(crystal.translation)
+    n_optical = len(optical.translation)
+    n_block = len(block.translation)
+    n_unit = len(unit.translation)
+    n_module = len(module.translation)
+
+    # set a (fake) digitizer
+    hits_actor, blur_actor = add_test_digitizer(sim, crystal, "output_ref.root")
+    blur_actor.use_truncated_Gaussian = True
+
+    # source and physics
+    stats = test_add_physics_and_stats(sim, "pet")
+    activity = 1 * MBq
+    if sim.visu:
+        activity = 50 * Bq
+    test_add_b2b_source(sim, activity)
+
+    # Set the function that will create the file for castor.
+    # This function is a hook because it must be run once the geometry is built by geant4.
+    # The param structure contains input parameters and will contain the output castor_config
+    # after the simulation
+    filename = paths.output / "castor_config_ref.json"
+    param = castor.set_hook_castor_config(sim, crystal.name, filename)
+
+    # go
+    duration = 1 * sec
+    print(f"Run with {activity/Bq:.0f} Bq during {duration/sec:.03f} s")
+    print(f"Mean time between events = {((1*sec)/(activity/Bq))/ns:.3f} ns")
+    sim.run_timing_intervals = [[0, duration]]
+    sim.run()
+
+    # print info
+    events = stats.counts.events
+    print(f"Number of simulated events : {events}")
+
+    # open the root file
+    root_filename = hits_actor.get_output_path()
+    root_file = uproot.open(root_filename)
+    file_size = os.path.getsize(root_filename)
+    print(f"Size of the file is {file_size / 1e6:.2f} MB")
+
+    # consider the singles and hits trees
+    hits_tree = root_file["hits"]
+    n = int(hits_tree.num_entries)
+    print(f"There are {n} hits")
+
+    singles_tree = root_file["singles"]
+    n = int(singles_tree.num_entries)
+    print(f"There are {n} singles")
+
+    # utility.test_ok(is_ok)