diff --git a/opengate/actors/coincidences.py b/opengate/actors/coincidences.py index f1ecce5bc..ae04c23dd 100644 --- a/opengate/actors/coincidences.py +++ b/opengate/actors/coincidences.py @@ -643,7 +643,7 @@ def ccmod_ideal_singles(data): ) df = filter_pandas_tree( df, branch_name="ProcessDefinedStep", value="Rayl", accepted=False - ) # Check but I think that this process did not generate a pulse + ).copy() # Create a new branch with ideal energy info df["IdealTotalEnergyDeposit"] = df["PreKineticEnergy"] - df["PostKineticEnergy"] @@ -668,7 +668,7 @@ def ccmod_ideal_coincidences(df): # create a new attribute CoincID that groups hits from the same coincidence. We can have more that two hits/pulses in a coincidence oe nSingles = df["EventID"].value_counts() # keep only events with more than one nSingles - df = df[df["EventID"].isin(nSingles[nSingles > 1].index)] + df = df[df["EventID"].isin(nSingles[nSingles > 1].index)].copy() # Assign CoincIDs starting from 0, in order of first appearance). df["CoincID"] = pd.factorize(df["EventID"])[0] diff --git a/opengate/contrib/compton_camera/coresi_helpers.py b/opengate/contrib/compton_camera/coresi_helpers.py new file mode 100644 index 000000000..7d1f79a89 --- /dev/null +++ b/opengate/contrib/compton_camera/coresi_helpers.py @@ -0,0 +1,208 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +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 = { + "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, + }, + "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") diff --git a/opengate/contrib/compton_camera/macaco.py b/opengate/contrib/compton_camera/macaco.py new file mode 100644 index 000000000..259067182 --- /dev/null +++ b/opengate/contrib/compton_camera/macaco.py @@ -0,0 +1,246 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from pathlib import Path +from opengate.utility import g4_units + + +def add_macaco1_materials(sim): + """ + Adds the Macaco materials database to the simulation if not already present. + """ + db_folder = Path(__file__).parent.resolve() + db_filename = db_folder / "macaco_materials.db" + if not db_filename in sim.volume_manager.material_database.filenames: + sim.volume_manager.add_material_database(db_filename) + + +def add_macaco1_camera(sim, name="macaco1"): + """ + Adds a MACACO1 camera to the simulation. + - Bounding box (BB_box) + - Plane 1 (scatterer side): holder, crystal carcass, PCB, SiPM, LaBr3_VLC crystal + - Plane 2 (absorber side): holder, crystal carcass, PCB, SiPM, LaBr3_VLC crystal + """ + + # units + mm = g4_units.mm + cm = g4_units.cm + + # material + add_macaco1_materials(sim) + + # BB box (acts as camera envelope) + camera = sim.add_volume("Box", f"{name}_BB_box") + camera.mother = sim.world + camera.material = "G4_AIR" + camera.size = [16 * cm, 40 * cm, 7.6 * cm] + camera.translation = [0, 0, 0 * cm] + camera.color = [0.1, 0.1, 0.1, 0.1] + + # Scatterer + holder1 = sim.add_volume("Box", f"{name}_Holder1") + holder1.mother = camera.name + holder1.material = "Plastic" + holder1.size = [6.2 * cm, 6.2 * cm, 0.56 * cm] + holder1.translation = [0, 0, -2.74 * cm] + holder1.color = [0.1, 0.1, 0.1, 0.9] + + crys_carcase_scatt = sim.add_volume("Box", f"{name}_crysCarcaseScatt") + crys_carcase_scatt.mother = holder1.name + crys_carcase_scatt.material = "G4_Al" + crys_carcase_scatt.size = [2.72 * cm, 2.68 * cm, 0.01 * cm] + crys_carcase_scatt.translation = [0, 0, -0.265 * cm] + + pcb_scatt = sim.add_volume("Box", f"{name}_PCBScatt") + pcb_scatt.mother = camera.name + pcb_scatt.material = "PCB" + pcb_scatt.size = [10.89 * cm, 20.7 * cm, 0.4 * cm] + pcb_scatt.translation = [0, 6.25 * cm, -2.46 * cm] + pcb_scatt.color = [0.0, 0.5, 0.0, 0.9] + + sipm_scatt = sim.add_volume("Box", f"{name}_SiPMScatt") + sipm_scatt.mother = holder1.name + sipm_scatt.material = "G4_Si" + sipm_scatt.size = [2.72 * cm, 2.68 * cm, 0.04 * cm] + sipm_scatt.translation = [0, 0, 0.26 * cm] + sipm_scatt.color = [1.0, 0.5, 0.0, 0.9] + + scatterer = sim.add_volume("Box", f"{name}_scatterer") + scatterer.mother = holder1.name + scatterer.material = "LaBr3_VLC" + scatterer.size = [2.72 * cm, 2.68 * cm, 0.50 * cm] + scatterer.translation = [0, 0, -0.01 * cm] + scatterer.color = [0.4, 0.7, 1.0, 1.0] + + # Absorber + holder2 = sim.add_volume("Box", f"{name}_Holder2") + holder2.mother = camera.name + holder2.material = "Plastic" + holder2.size = [8.0 * cm, 8.0 * cm, 1.06 * cm] + holder2.translation = [0, 0, 2.51 * cm] + holder2.color = [0.1, 0.1, 0.1, 0.9] + + crys_carcase_abs = sim.add_volume("Box", f"{name}_crysCarcaseAbs") + crys_carcase_abs.mother = holder2.name + crys_carcase_abs.material = "G4_Al" + crys_carcase_abs.size = [3.24 * cm, 3.60 * cm, 0.01 * cm] + crys_carcase_abs.translation = [0, 0, -0.515 * cm] + + absorber = sim.add_volume("Box", f"{name}_absorber") + absorber.mother = holder2.name + absorber.material = "LaBr3_VLC" + absorber.size = [3.24 * cm, 3.60 * cm, 1.00 * cm] + absorber.translation = [0, 0, -0.01 * cm] + absorber.color = [0.4, 0.7, 1.0, 1.0] + + sipm_abs = sim.add_volume("Box", f"{name}_SiPMAbs") + sipm_abs.mother = holder2.name + sipm_abs.material = "G4_Si" + sipm_abs.size = [3.24 * cm, 3.60 * cm, 0.04 * cm] + sipm_abs.translation = [0, 0, 0.51 * cm] + sipm_abs.color = [1.0, 0.5, 0.0, 0.9] + + pcb_abs = sim.add_volume("Box", f"{name}_PCBAbs") + pcb_abs.mother = camera.name + pcb_abs.material = "PCB" + pcb_abs.size = [9.54 * cm, 16.0 * cm, 0.4 * cm] + pcb_abs.translation = [0, 4.50 * cm, 3.24 * cm] + pcb_abs.color = [0.0, 0.5, 0.0, 0.9] + + return { + "camera": camera, + "scatterer": scatterer, + "absorber": absorber, + "holder1": holder1, + "holder2": holder2, + } + + +def add_macaco1_camera_digitizer(sim, scatterer, absorber, output_path): + + # Units + keV = g4_units.keV + MeV = g4_units.MeV + mm = g4_units.mm + ns = g4_units.ns + + # 1) Collect step-level hits for both scatterer and absorber + hits_scatt = sim.add_actor("DigitizerHitsCollectionActor", "HitsScatt") + hits_scatt.attached_to = scatterer.name + hits_scatt.attributes = [ + "EventID", + "TrackID", + "TotalEnergyDeposit", + "GlobalTime", + "PrePosition", + "PostPosition", + "PreStepUniqueVolumeID", + ] + + hits_abs = sim.add_actor("DigitizerHitsCollectionActor", "HitsAbs") + hits_abs.attached_to = absorber.name + hits_abs.attributes = hits_scatt.attributes + scatt_collection = hits_scatt.name + abs_collection = hits_abs.name + + # 2) Process Hits into Singles + sing_scatt = sim.add_actor("DigitizerAdderActor", "SinglesScatt") + sing_scatt.input_digi_collection = scatt_collection + sing_scatt.policy = "EnergyWeightedCentroidPosition" + sing_scatt.attributes = [ + "EventID", + "GlobalTime", + "TotalEnergyDeposit", + "PostPosition", + ] + scatt_collection = sing_scatt.name + + sing_abs = sim.add_actor("DigitizerAdderActor", "SinglesAbs") + sing_abs.input_digi_collection = abs_collection + sing_abs.policy = "EnergyWeightedCentroidPosition" + sing_abs.attributes = sing_scatt.attributes + abs_collection = sing_abs.name + + # Spatial blurring + spat_scatt = sim.add_actor("DigitizerSpatialBlurringActor", "SpatScatt") + spat_scatt.attached_to = scatterer.name + spat_scatt.input_digi_collection = scatt_collection + spat_scatt.blur_attribute = "PostPosition" + spat_scatt.blur_fwhm = [2 * mm, 2 * mm, 2 * mm] + scatt_collection = spat_scatt.name + + spat_abs = sim.add_actor("DigitizerSpatialBlurringActor", "SpatAbs") + spat_abs.attached_to = absorber.name + spat_abs.input_digi_collection = abs_collection + spat_abs.blur_attribute = "PostPosition" + spat_abs.blur_fwhm = [2 * mm, 2 * mm, 2 * mm] + abs_collection = spat_abs.name + + # Energy blurring + reference_energy = 511 * keV + scatt_resolution = 0.085 # FWHM/E at 511 keV + abs_resolution = 0.125 # FWHM/E at 511 keV + + blur_scatt = sim.add_actor("DigitizerBlurringActor", "BlurScatt") + blur_scatt.attached_to = scatterer.name + blur_scatt.input_digi_collection = scatt_collection + blur_scatt.blur_attribute = "TotalEnergyDeposit" + blur_scatt.blur_method = "InverseSquare" + blur_scatt.blur_reference_value = reference_energy + blur_scatt.blur_resolution = scatt_resolution + scatt_collection = blur_scatt.name + + blur_abs = sim.add_actor("DigitizerBlurringActor", "BlurAbs") + blur_abs.attached_to = absorber.name + blur_abs.input_digi_collection = abs_collection + blur_abs.blur_attribute = "TotalEnergyDeposit" + blur_abs.blur_method = "InverseSquare" + blur_abs.blur_reference_value = reference_energy + blur_abs.blur_resolution = abs_resolution + abs_collection = blur_abs.name + + # Time blurring + time_fwhm = 10 * ns + time_sigma = time_fwhm / 2.355 + time_scatt = sim.add_actor("DigitizerBlurringActor", "TimeBlurScatt") + time_scatt.attached_to = scatterer.name + time_scatt.input_digi_collection = scatt_collection + time_scatt.blur_attribute = "GlobalTime" + time_scatt.blur_method = "Gaussian" + time_scatt.blur_sigma = time_sigma + scatt_collection = time_scatt.name + + time_abs = sim.add_actor("DigitizerBlurringActor", "TimeBlurAbs") + time_abs.attached_to = absorber.name + time_abs.input_digi_collection = abs_collection + time_abs.blur_attribute = "GlobalTime" + time_abs.blur_method = "Gaussian" + time_abs.blur_sigma = time_sigma + abs_collection = time_abs.name + + # Energy windows (thresholds) + threshold_min = 70 * keV + threshold_max = 2.0 * MeV + thr_scatt = sim.add_actor("DigitizerEnergyWindowsActor", "ThrScatt") + thr_scatt.attached_to = scatterer.name + thr_scatt.input_digi_collection = scatt_collection + thr_scatt.channels = [ + {"name": thr_scatt.name, "min": threshold_min, "max": threshold_max} + ] + scatt_collection = thr_scatt.name + + thr_abs = sim.add_actor("DigitizerEnergyWindowsActor", "ThrAbs") + thr_abs.attached_to = absorber.name + thr_abs.input_digi_collection = abs_collection + thr_abs.channels = [ + {"name": thr_abs.name, "min": threshold_min, "max": threshold_max} + ] + + # Saving root files + thr_scatt.output_filename = f"{thr_scatt.name}.root" + thr_abs.output_filename = f"{thr_abs.name}.root" + scatt_file = thr_scatt.get_output_path() + abs_file = thr_abs.get_output_path() + + return scatt_file, abs_file diff --git a/opengate/contrib/compton_camera/macaco_materials.db b/opengate/contrib/compton_camera/macaco_materials.db new file mode 100644 index 000000000..9bd8fd933 --- /dev/null +++ b/opengate/contrib/compton_camera/macaco_materials.db @@ -0,0 +1,103 @@ +[Elements] +Hydrogen: S= H ; Z= 1. ; A= 1.01 g/mole +Carbon: S= C ; Z= 6. ; A= 12.01 g/mole +Nitrogen: S= N ; Z= 7. ; A= 14.01 g/mole +Oxygen: S= O ; Z= 8. ; A= 16.00 g/mole +Aluminium: S= Al ; Z= 13. ; A= 26.98 g/mole +Argon: S= Ar ; Z= 18. ; A= 39.95 g/mole +Titanium: S= Ti ; Z= 22. ; A= 47.867 g/mole +Manganese: S= Mn ; Z= 25. ; A= 54.938 g/mole +Iron: S= Fe ; Z= 26. ; A= 55.845 g/mole +Lead: S= Pb ; Z= 82. ; A= 207.20 g/mole +Silicon: S= Si ; Z= 14. ; A= 28.09 g/mole +Cadmium: S= Cd ; Z= 48. ; A= 112.41 g/mole +Tellurium: S= Te ; Z= 52. ; A= 127.6 g/mole +Zinc: S= Zn ; Z= 30. ; A= 65.39 g/mole +Tungsten: S= W ; Z= 74. ; A= 183.84 g/mole +Gallium: S= Ga ; Z= 31. ; A= 69.723 g/mole +Germanium: S= Ge ; Z= 32. ; A= 72.61 g/mole +Bromine: S= Br ; Z= 35. ; A= 79.90 g/mole +Yttrium: S= Y ; Z= 39. ; A= 88.91 g/mole +Cesium: S= Cs ; Z= 55. ; A= 132.905 g/mole +Lanthanum: S= La ; Z= 57. ; A= 138.91 g/mole +Cerium: S= Ce ; Z= 58. ; A= 140.116 g/mole +Gadolinium: S= Gd ; Z= 64. ; A= 157.25 g/mole +Lutetium: S= Lu ; Z= 71. ; A= 174.97 g/mole +Thallium: S= Tl ; Z= 81. ; A= 204.37 g/mole +Bismuth: S= Bi ; Z= 83. ; A= 208.98 g/mole +Uranium: S= U ; Z= 92. ; A= 238.03 g/mole + +[Materials] +Vacuum: d=0.000001 mg/cm3 ; n=1 + +el: name=Hydrogen ; n=1 + +Tungsten: d=19.3 g/cm3 ; n=1 ; state=solid + +el: name=auto ; n=1 + +Lead: d=11.4 g/cm3 ; n=1 ; state=solid + +el: name=auto ; n=1 + + +Aluminium: d=2.7 g/cm3 ; n=1 ; state=solid + +el: name=auto ; n=1 + +Plastic: d=1.18 g/cm3 ; n=3; state=solid + +el: name=Carbon ; n=5 + +el: name=Hydrogen ; n=8 + +el: name=Oxygen ; n=2 + +CZT: d=5.68 g/cm3 ; n=3; state=solid + +el: name=Cadmium ; n=9 + +el: name=Zinc ; n=1 + +el: name=Tellurium ; n=10 + + +LYSO-Ce-Hilger: d=7.10 g/cm3; n=5 ; state=Solid + +el: name=Lutetium ; f=0.713838658203075 + +el: name=Yttrium ; f=0.040302477781781 + +el: name=Silicon; f=0.063721807284236 + +el: name=Oxygen; f=0.181501252152072 + +el: name=Cerium; f=0.000635804578835201 + +Polyethylene: d=0.96 g/cm3 ; n=2 + +el: name=Hydrogen ; n=2 + +el: name=Carbon ; n=1 + +PVC: d=1.65 g/cm3 ; n=3 ; state=solid + +el: name=Hydrogen ; n=3 + +el: name=Carbon ; n=2 + +el: name=Chlorine ; n=1 + +SS304: d=7.92 g/cm3 ; n=4 ; state=solid + +el: name=Iron ; f=0.695 + +el: name=Chromium ; f=0.190 + +el: name=Nickel ; f=0.095 + +el: name=Manganese ; f=0.020 + +PTFE: d= 2.18 g/cm3 ; n=2 ; state=solid + +el: name=Carbon ; n=1 + +el: name=Fluorine ; n=2 + +Plexiglass: d=1.19 g/cm3; n=3; state=solid + +el: name=Hydrogen; f=0.080538 + +el: name=Carbon; f=0.599848 + +el: name=Oxygen; f=0.319614 + +LaBr3: d=5.29 g/cm3 ; n=2 + +el: name=Bromine; n=3 + +el: name=Lanthanum; n=1 + + +CeBr3: d=5.1 g/cm3 ; n=2 + +el: name=Bromine; n=3 + +el: name=Cerium; n=1 + + +LaBr3_VLC: d=5.06 g/cm3 ; n=2 + +el: name=Bromine; n=3 + +el: name=Lanthanum; n=1 + +PCB: d=1.20 g/cm3 ; n=3 + +el: name=Hydrogen; n=8 + +el: name=Carbon; n=5 + +el: name=Oxygen; n=2 diff --git a/opengate/contrib/pet/castor_helpers.py b/opengate/contrib/pet/castor_helpers.py index d6a265a9e..d348d7c79 100644 --- a/opengate/contrib/pet/castor_helpers.py +++ b/opengate/contrib/pet/castor_helpers.py @@ -7,7 +7,7 @@ def set_hook_castor_config(sim, crystal_name, filename): """ Prepare everything to create the castor config file at the init of the simulation. - The param structure allows to retrieve the castor config at the end of the simulation. + The param structure allows retrieving the castor config at the end of the simulation. """ # create the param structure param = { @@ -21,6 +21,7 @@ def set_hook_castor_config(sim, crystal_name, filename): def create_castor_config(simulation_engine, param): + # (note: simulation_engine is not used here but must be the first param) castor_config = { "rotation": [], "size": [], diff --git a/opengate/tests/src/actors/test096_ideal_Compton_sorter.py b/opengate/tests/src/actors/test096_ideal_compton_sorter.py similarity index 92% rename from opengate/tests/src/actors/test096_ideal_Compton_sorter.py rename to opengate/tests/src/actors/test096_ideal_compton_sorter.py index a9929c1fe..cb496e41d 100755 --- a/opengate/tests/src/actors/test096_ideal_Compton_sorter.py +++ b/opengate/tests/src/actors/test096_ideal_compton_sorter.py @@ -1,13 +1,11 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -from opengate.actors.coincidences import ccmod_ideal_singles -from opengate.actors.coincidences import ccmod_ideal_coincidences -from opengate.actors.coincidences import ccmod_make_cones +from opengate.actors.coincidences import * import uproot import numpy as np import pandas as pd -from test096_ideal_Compton_sorter_helpers import * +from test096_ideal_compton_sorter_helpers import * def main(): diff --git a/opengate/tests/src/actors/test096_ideal_Compton_sorter_helpers.py b/opengate/tests/src/actors/test096_ideal_compton_sorter_helpers.py similarity index 100% rename from opengate/tests/src/actors/test096_ideal_Compton_sorter_helpers.py rename to opengate/tests/src/actors/test096_ideal_compton_sorter_helpers.py diff --git a/opengate/tests/src/actors/test097_coresi_ccmod.py b/opengate/tests/src/actors/test097_coresi_ccmod.py new file mode 100755 index 000000000..fa20ebe1b --- /dev/null +++ b/opengate/tests/src/actors/test097_coresi_ccmod.py @@ -0,0 +1,149 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.utility import g4_units +import opengate.tests.utility as utility +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__": + + # get tests paths + paths = utility.get_default_test_paths( + __file__, gate_folder="", output_folder="test097_coresi_ccmod" + ) + output_folder = paths.output + data_folder = paths.data + + # units + m = g4_units.m + mm = g4_units.mm + cm = g4_units.cm + keV = g4_units.keV + Bq = g4_units.Bq + MBq = 1e6 * g4_units.Bq + sec = g4_units.s + + # sim + sim = gate.Simulation() + sim.visu = False + sim.visu_type = "qt" + sim.random_seed = "auto" + sim.output_dir = output_folder + sim.number_of_threads = 1 + + # world + world = sim.world + world.size = [1 * m, 1 * m, 2 * m] + sim.world.material = "G4_AIR" + + # 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 = camera1 # [camera1, camera2] + phsp.attributes = [ + "TotalEnergyDeposit", + "PreKineticEnergy", + "PostKineticEnergy", + "PostPosition", + "ProcessDefinedStep", + "ParticleName", + "EventID", + "ParentID", + "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(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") + source.particle = "gamma" + source.energy.mono = 662 * keV + source.position.type = "sphere" + source.position.radius = 0.25 * mm + source.position.translation = [0, 0, 0] + 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() + + # print stats + print(stats) + + # open the root file and create coincidences + print() + root_file = uproot.open(phsp.get_output_path()) + tree = root_file["PhaseSpace"] + hits = tree.arrays(library="pd") + singles = ccmod_ideal_singles(hits) + 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 + coinc_filename = str(phsp.get_output_path()).replace(".root", "_coinc.root") + root_write_trees( + coinc_filename, ["hits", "singles", "coincidences"], [hits, singles, coinc] + ) + print(f"Output file: {coinc_filename}") + print() + + # 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)