diff --git a/opengate/tests/src/source/testXXX_source_confine_in_the_detector_volume.py b/opengate/tests/src/source/testXXX_source_confine_in_the_detector_volume.py new file mode 100644 index 0000000000..9281edf06d --- /dev/null +++ b/opengate/tests/src/source/testXXX_source_confine_in_the_detector_volume.py @@ -0,0 +1,171 @@ +#!/usr/bin/env python3 +""" +Test and Guide: Intrinsic Activity in Detectors +----------------------------------------------- +This script demonstrates how to: +1. Define a PET scanner geometry with hierarchical volumes (Module -> Unit -> Block -> Crystal). +2. Define intrinsic activity (e.g., Lu176 background) confined specifically to the crystal volumes. +3. Calculate the total activity based on the volume of the crystals and a known concentration. +""" + +import numpy as np +import opengate as gate +from opengate.geometry.utility import get_grid_repetition, get_circular_repetition + +# ---------- units ---------- +u = gate.g4_units +mm, cm = u.mm, u.cm +keV, Bq, sec = u.keV, u.Bq, u.s + +# ---------- simulation/world and material ---------- +sim = gate.Simulation() + +sim.volume_manager.material_database.add_material_weights( + "Vacuum", ["H"], [1], 1e-12 * u.g_cm3 +) + +sim.volume_manager.material_database.add_material_weights( + "LYSO", + ["Lu", "Y", "Si", "O"], # Elements + [ + 0.714467891, + 0.04033805, + 0.063714272, + 0.181479788, + ], # Mass Fractions (must sum to ~1.0) + 7.36 * u.g_cm3, # Density +) + +sim.world.size = [90 * cm, 90 * cm, 90 * cm] +sim.world.material = "Vacuum" + +sim.world.color = [1, 1, 1, 0] + +# ---------- module container ---------- +# Module: The main sector of the ring. Contains 5 Units in Z. +module = sim.add_volume("Box", "module") +module.mother = "world" +module.size = [64 * mm, 30 * mm, 180.0 * mm] +module.material = "Vacuum" +module.color = [1, 1, 1, 0] +translations_ring, rotations_ring = get_circular_repetition( + 20, [0 * mm, 300 * mm, 0 * mm], start_angle_deg=90, axis=[0, 0, 1] +) +module.translation = translations_ring +module.rotation = rotations_ring + +# # ---------- Unit (1×5) inside Module ---------- +Unit_det = sim.add_volume("Box", "Unit") +Unit_det.mother = module.name +Unit_det.size = [64 * mm, 30.0 * mm, 36 * mm] +Unit_det.material = "Vacuum" +Unit_det.color = [1, 0, 0, 0] +Unit_det.translation = get_grid_repetition([1, 1, 5], [0 * mm, 0 * mm, 36 * mm]) + +# # ---------- Block (4×1) inside Unit ---------- +Block_det = sim.add_volume("Box", "Block") +Block_det.mother = Unit_det.name +Block_det.size = [16 * mm, 30.0 * mm, 36 * mm] +Block_det.material = "Vacuum" +Block_det.color = [0, 1, 0, 1] +Block_det.translation = get_grid_repetition([4, 1, 1], [16 * mm, 0.0 * mm, 0.0 * mm]) + +# # ---------- Crystal (4×9) inside Block ---------- +crystal_det = sim.add_volume("Box", "crystal") +crystal_det.mother = Block_det.name +crystal_det.size = [4 * mm, 30.0 * mm, 4 * mm] +crystal_det.material = "Vacuum" +crystal_det.color = [0, 0, 1, 0] +crystal_det.translation = get_grid_repetition([4, 1, 9], [4 * mm, 0.0 * mm, 4 * mm]) + +# # ---------- Crystal Layer (Active Material) same size as crystal_det ---------- +# NOTE: This is the final layer of the geometry definition. +# It is defined as a single volume (no repetition here) inside the 'crystal' volume. +# We use THIS volume ('crystal_layer') for source confinement. +crystal_det_layer = sim.add_volume("Box", "crystal_layer") +crystal_det_layer.mother = crystal_det.name +crystal_det_layer.size = [4 * mm, 30.0 * mm, 4 * mm] +crystal_det_layer.material = "LYSO" +crystal_det_layer.color = [0, 1, 1, 0] + +# # ---------------------- Degitizer defination ---------------------- +crystal_det_layer = sim.volume_manager.volumes["crystal_layer"] +Block_det = sim.volume_manager.volumes["Block"] + +hc = sim.add_actor("DigitizerHitsCollectionActor", "Hits") +hc.attached_to = crystal_det_layer.name +hc.authorize_repeated_volumes = True +hc.output_filename = None + +hc.attributes = [ + "EventID", + "TrackID", + "TotalEnergyDeposit", + "GlobalTime", + "EventPosition", + "Position", + "ProcessDefinedStep", + "PostPosition", + "PreStepUniqueVolumeID", +] + +sc = sim.add_actor("DigitizerReadoutActor", "Singles") +sc.input_digi_collection = "Hits" +sc.group_volume = Block_det.name +sc.discretize_volume = crystal_det_layer.name +sc.policy = "EnergyWeightedCentroidPosition" +sc.output_filename = r"singles.root" + +# ---------------------- Change physics ---------------------- +sim.physics_manager.physics_list_name = "G4EmPenelopePhysics" +sim.physics_manager.enable_decay = True + +# ------------------ Intrinsic Activity Source --------------------- +# To simulate intrinsic activity (e.g., Lu176 in LYSO), we confine the source +# to the detector volume. + +source_radius = 380.0 # mm +source_hight = 255.0 # mm +Lu176_concentration = 0.270 # Bq/mm3 + +# Calculate total volume of all crystals to determine total activity +# Volume = (Crystal Size) * (Crystals per Block) * (Blocks per Unit) * (Units per Module) * (Modules) +det_total_vol = (4 * 30.0 * 4) * (4 * 9) * 5 * 4 * 20 # mm3 +total_activity = det_total_vol * Lu176_concentration + +source = sim.add_source("GenericSource", "myConfSource") +# 'ion Z A' -> Lu is Z=71, Mass=176 +source.particle = "ion 71 176" # Lu176 +source.position.type = "cylinder" +source.position.radius = source_radius * mm +source.position.dz = (source_hight / 2) * mm +source.position.translation = [0.0 * mm, 0.0 * mm, 0.0 * mm] +source.half_life = 3.7e10 * u.year +source.position.confine = crystal_det_layer.name +# Note: We confine to "crystal_layer". Since "crystal_layer" is inside repeated volumes +# (Block -> Unit -> Module), the source is automatically distributed across all instances. +# Activity is set for the *entire* confined volume name. +source.activity = total_activity * Bq + +stats = sim.add_actor("SimulationStatisticsActor", "Stats") +stats.track_types_flag = True +stats.output_filename = "stats.txt" + + +# ---------- Setup scan and viewer ---------- +sim.g4_verbose = True +sim.g4_verbose_level = 1 +sim.visu = False # True +# sim.visu_type = "vrml_file_only" # "vrml" +# sim.visu_filename = "geometry.wrl" + +sim.number_of_threads = 1 + +sim.random_engine = "MersenneTwister" +sim.random_seed = "auto" + +sim.run_timing_intervals = [[0 * sec, 0.01 * sec]] + +# ---------- Run ---------- +# sim.run(start_new_process=True) # not available on windows +sim.run()