add CORESI configuration management utilities, refactor coresi_helpers.py, and update test097 for multi-camera setup

Author: dsarrut

opengate/contrib/compton_camera/coresi_helpers.py

@@ -1,40 +1,208 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
-from box import Box
+import opengate_core as g4
+from opengate.geometry.utility import vec_g4_as_np, rot_g4_as_np
+from opengate.exception import fatal
+from opengate.utility import g4_units
+import yaml
+import uproot
+
+
+# --- Custom List for Inline YAML Formatting ---
+class FlowList(list):
+    """A custom list that will be dumped as [x, y, z] in YAML."""
+
+    pass
+
+
+def flow_list_representer(dumper, data):
+    return dumper.represent_sequence("tag:yaml.org,2002:seq", data, flow_style=True)
+
+
+yaml.add_representer(FlowList, flow_list_representer)
+
+
+def convert_to_flowlist(data):
+    """
+    Recursively traverse the dictionary.
+    Convert any list containing only simple scalars (int, float, str) to FlowList
+    so they appear as [a, b, c] in the YAML output.
+    """
+    if isinstance(data, dict):
+        return {k: convert_to_flowlist(v) for k, v in data.items()}
+    elif isinstance(data, list):
+        # Check if list is "simple" (contains only primitives, no dicts or nested lists)
+        is_simple = all(
+            isinstance(i, (int, float, str, bool)) or i is None for i in data
+        )
+
+        if is_simple:
+            return FlowList(data)
+        else:
+            # If list contains complex objects, process them recursively but keep the list as is
+            return [convert_to_flowlist(i) for i in data]
+    else:
+        return data
 
 
 def coresi_new_config():
-    config = Box(
-        {
-            "data_file": "coinc.dat",
-            "data_type": "GATE",
-            "n_events": 0,
-            "starts_at": 0,
-            "E0": [],
-            "remove_out_of_range_energies": False,
-            "energy_range": [120, 150],
-            "energy_threshold": 5,
-            "log_dir": None,
-            "cameras": {"n_cameras": 0},
-            "volume": {
-                "volume_dimensions": [10, 10, 10],  # in cm
-                "n_voxels": [50, 50, 1],
-                "volume_centre": [0, 0, 0],
+    config = {
+        "data_file": "coinc.dat",
+        "data_type": "GATE",
+        "n_events": 0,
+        "starts_at": 0,
+        "E0": [],
+        "remove_out_of_range_energies": False,
+        "energy_range": [120, 150],
+        "energy_threshold": 5,
+        "log_dir": None,
+        "cameras": {
+            "n_cameras": 0,
+            "common_attributes": {
+                "n_sca_layers": 0,
+                "sca_material": "Si",
+                "abs_material": "Si",
+                "n_absorbers": 0,
             },
-            "lm_mlem": {
-                "cone_thickness": "angular",
-                "model": "cos1rho2",
-                "last_iter": 0,
-                "first_iter": 0,
-                "n_sigma": 2,
-                "width_factor": 1,
-                "checkpoint_dir": "checkpoints",
-                "save_every": 76,
-                "sensitivity": False,
-                "sensitivity_model": "like_system_matrix",
-                "sensitivity_point_samples": 1,
+            "position_0": {
+                "frame_origin": [0, 0, 0],
+                "Ox": [1, 0, 0],  # parallel to scatterer edge
+                "Oy": [0, 1, 0],  # parallel to scatterer edge
+                "Oz": [0, 0, 1],  # orthogonal to the camera, tw the source"
             },
-        }
-    )
+        },
+        "volume": {
+            "volume_dimensions": [10, 10, 10],  # in cm?
+            "n_voxels": [50, 50, 1],  # in voxels
+            "volume_centre": [0, 0, 0],  # in cm?
+        },
+        "lm_mlem": {
+            "cone_thickness": "angular",
+            "model": "cos1rho2",
+            "last_iter": 0,
+            "first_iter": 0,
+            "n_sigma": 2,
+            "width_factor": 1,
+            "checkpoint_dir": "checkpoints",
+            "save_every": 76,
+            "sensitivity": False,
+            "sensitivity_model": "like_system_matrix",
+            "sensitivity_point_samples": 1,
+        },
+    }
+
     return config
+
+
+def set_hook_coresi_config(sim, cameras, filename):
+    """
+    Prepare everything to create the coresi config file at the init of the simulation.
+    The param structure allows retrieving the coresi config at the end of the simulation.
+    """
+    # create the param structure
+    param = {
+        "cameras": cameras,
+        "filename": filename,
+        "coresi_config": coresi_new_config(),
+    }
+    sim.user_hook_after_init = create_coresi_config
+    sim.user_hook_after_init_arg = param
+    return param
+
+
+def create_coresi_config(simulation_engine, param):
+    # (note: simulation_engine is not used here but must be the first param)
+    coresi_config = param["coresi_config"]
+    cameras = param["cameras"]
+
+    for camera in cameras.values():
+        c = coresi_config["cameras"]
+        c["n_cameras"] += 1
+        scatter_layer_names = camera["scatter_layer_names"]
+        absorber_layer_names = camera["absorber_layer_names"]
+
+        for layer_name in scatter_layer_names:
+            coresi_add_scatterer(coresi_config, layer_name)
+        for layer_name in absorber_layer_names:
+            coresi_add_absorber(coresi_config, layer_name)
+
+
+def coresi_add_scatterer(coresi_config, layer_name):
+    # find all volumes ('touchable' in Geant4 terminology)
+    touchables = g4.FindAllTouchables(layer_name)
+    if len(touchables) != 1:
+        fatal(f"Cannot find unique volume for layer {layer_name}: {touchables}")
+    touchable = touchables[0]
+
+    # current nb of scatterers
+    id = coresi_config["cameras"]["common_attributes"]["n_sca_layers"]
+    coresi_config["cameras"]["common_attributes"]["n_sca_layers"] += 1
+    layer = {
+        "center": [0, 0, 0],
+        "size": [0, 0, 0],
+    }
+    coresi_config["cameras"]["common_attributes"][f"sca_layer_{id}"] = layer
+
+    # Get the information: WARNING in cm!
+    cm = g4_units.cm
+    translation = vec_g4_as_np(touchable.GetTranslation(0)) / cm
+    solid = touchable.GetSolid(0)
+    pMin_local = g4.G4ThreeVector()
+    pMax_local = g4.G4ThreeVector()
+    solid.BoundingLimits(pMin_local, pMax_local)
+    size = [
+        (pMax_local.x - pMin_local.x) / cm,
+        (pMax_local.y - pMin_local.y) / cm,
+        (pMax_local.z - pMin_local.z) / cm,
+    ]
+    layer["center"] = translation.tolist()
+    layer["size"] = size
+
+
+def coresi_add_absorber(coresi_config, layer_name):
+    # find all volumes ('touchable' in Geant4 terminology)
+    touchables = g4.FindAllTouchables(layer_name)
+    if len(touchables) != 1:
+        fatal(f"Cannot find unique volume for layer {layer_name}: {touchables}")
+    touchable = touchables[0]
+
+    # current nb of scatterers
+    id = coresi_config["cameras"]["common_attributes"]["n_absorbers"]
+    coresi_config["cameras"]["common_attributes"]["n_absorbers"] += 1
+    layer = {
+        "center": [0, 0, 0],
+        "size": [0, 0, 0],
+    }
+    coresi_config["cameras"]["common_attributes"][f"abs_layer_{id}"] = layer
+
+    # Get the information: WARNING in cm!
+    cm = g4_units.cm
+    translation = vec_g4_as_np(touchable.GetTranslation(0)) / cm
+    solid = touchable.GetSolid(0)
+    pMin_local = g4.G4ThreeVector()
+    pMax_local = g4.G4ThreeVector()
+    solid.BoundingLimits(pMin_local, pMax_local)
+    size = [
+        (pMax_local.x - pMin_local.x) / cm,
+        (pMax_local.y - pMin_local.y) / cm,
+        (pMax_local.z - pMin_local.z) / cm,
+    ]
+    layer["center"] = translation.tolist()
+    layer["size"] = size
+
+
+def coresi_write_config(coresi_config, filename):
+    # Convert vectors to FlowList just before writing
+    formatted_config = convert_to_flowlist(coresi_config)
+
+    with open(filename, "w") as f:
+        yaml.dump(
+            formatted_config, f, default_flow_style=False, sort_keys=False, indent=2
+        )
+
+
+def coresi_convert_root_data(root_filename, branch_name, output_filename):
+    root_file = uproot.open(root_filename)
+    tree = root_file[branch_name]
+    print("todo")

opengate/tests/src/test097_coresi_ccmod.py

@@ -4,8 +4,10 @@
 import opengate as gate
 from opengate.utility import g4_units
 import opengate.tests.utility as utility
-import opengate.contrib.compton_camera.macaco as macaco
 from opengate.actors.coincidences import *
+import opengate.contrib.compton_camera.macaco as macaco
+from scipy.spatial.transform import Rotation
+import opengate.contrib.compton_camera.coresi_helpers as coresi
 
 if __name__ == "__main__":
 
@@ -38,16 +40,25 @@
     world.size = [1 * m, 1 * m, 2 * m]
     sim.world.material = "G4_AIR"
 
-    # add the camera
-    camera = macaco.add_macaco1_camera(sim, "macaco1")
+    # add two cameras
+    name1 = "macaco1"
+    camera1 = macaco.add_macaco1_camera(sim, name1)
+    camera1.translation = [0, 0, 10 * cm]
+
+    """
+    name2 = "macaco2"
+    camera2 = macaco.add_macaco1_camera(sim, name2)
+    camera2.rotation = Rotation.from_euler("x", -90, degrees=True).as_matrix()
+    camera2.translation = [0, 10 * cm, 0 * cm]
+    """
 
     # stats
     stats = sim.add_actor("SimulationStatisticsActor", "Stats")
     stats.output_filename = "stats.txt"
 
     # PhaseSpace Actor
     phsp = sim.add_actor("PhaseSpaceActor", "PhaseSpace")
-    phsp.attached_to = camera
+    phsp.attached_to = camera1  # [camera1, camera2]
     phsp.attributes = [
         "TotalEnergyDeposit",
         "PreKineticEnergy",
@@ -60,14 +71,17 @@
         "PDGCode",
         "TrackVertexKineticEnergy",
         "GlobalTime",
+        "PreStepUniqueVolumeID",
+        "PreStepUniqueVolumeIDAsInt",
     ]
     phsp.output_filename = "phsp.root"
     phsp.steps_to_store = "allsteps"
 
     # physics
     sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option2"
     sim.physics_manager.set_production_cut("world", "all", 1 * mm)
-    sim.physics_manager.set_production_cut(camera.name, "all", 0.1 * mm)
+    sim.physics_manager.set_production_cut(camera1.name, "all", 0.1 * mm)
+    # sim.physics_manager.set_production_cut(camera2.name, "all", 0.1 * mm)
 
     # source
     source = sim.add_source("GenericSource", "src")
@@ -79,11 +93,23 @@
     source.direction.type = "iso"
     source.activity = 0.847 * MBq / sim.number_of_threads
 
+    # acquisition time
     if sim.visu:
         source.activity = 1 * Bq
-
     sim.run_timing_intervals = [[0 * sec, 5 * sec]]
 
+    # special hook to prepare coresi config file.
+    # For each camera, we must find the names of all layers (scatterer and absorber).
+    # In this simple macaco1 case, there is only one of each.
+    cameras = {
+        name1: {
+            "scatter_layer_names": [f"{name1}_scatterer"],
+            "absorber_layer_names": [f"{name1}_absorber"],
+        }
+    }
+    yaml_filename = paths.output / "coresi_config.yaml"
+    param = coresi.set_hook_coresi_config(sim, cameras, yaml_filename)
+
     # go
     sim.run()
 
@@ -92,25 +118,32 @@
 
     # open the root file and create coincidences
     print()
-    print(f"Output file: {phsp.get_output_path()}")
     root_file = uproot.open(phsp.get_output_path())
     tree = root_file["PhaseSpace"]
     hits = tree.arrays(library="pd")
-    print(f"Number of hits: {len(hits)} ")
     singles = ccmod_ideal_singles(hits)
-    print(f"Found: {len(singles)} singles")
     coinc = ccmod_ideal_coincidences(singles)
+    print(f"Output file: {phsp.get_output_path()}")
+    print(f"Number of hits: {len(hits)} ")
+    print(f"Found: {len(singles)} singles")
     print(f"Found: {len(coinc)} coincidences")
 
     # write the new root with coinc
-    filename = str(phsp.get_output_path()).replace(".root", "_coinc.root")
+    coinc_filename = str(phsp.get_output_path()).replace(".root", "_coinc.root")
     root_write_trees(
-        filename, ["hits", "singles", "coincidences"], [hits, singles, coinc]
+        coinc_filename, ["hits", "singles", "coincidences"], [hits, singles, coinc]
     )
-    print(f"Output file: {filename}")
+    print(f"Output file: {coinc_filename}")
     print()
 
-    # CORESI stage1: create the configuration file
-    # coresi_config =
+    # CORESI stage1: we retrieve the coresi config built during the hook
+    coresi_config = param["coresi_config"]
+    # optional: change the volume information
+    # TODO LATER : set parameters from an image
+    coresi_config["volume"]["volume_dimensions"] = [20, 20, 0.5]
+    coresi.coresi_write_config(coresi_config, yaml_filename)
+    print(f"Coresi config file: {coinc_filename}")
 
     # CORESI stage2: convert the root file
+    data_filename = output_folder / "coincidences.dat"
+    coresi.coresi_convert_root_data(coinc_filename, "coincidences", data_filename)