OpenGATE / opengate

Gate 10 (beta)
http://www.opengatecollaboration.org
GNU Lesser General Public License v3.0
42 stars 38 forks source link

Edep distribution between F18 and Ga68 #342

Closed jizhang02 closed 7 months ago

jizhang02 commented 7 months ago

Hello,

When I ran the pet simulation, I found a strange thing. When I use F18 as source energy, the distribution of edep in root file is like this: image However, when I use Ga68 as source energy, the distribution of edep in root file is like this, the peak energy is not in 511Kv. image

I would like to ask if is it normal. Thank you.

The code I use is as below:

import opengate as gate
import opengate.contrib.pet.siemensbiograph as pet_biograph # pet system
import opengate.contrib.pet.philipsvereos as pet_vereos
from opengate.userhooks import check_production_cuts
import numpy as np

'''
General
'''
sim = gate.Simulation()
sim.number_of_threads = 1
ui = sim.user_info
ui.g4_verbose = False
ui.g4_verbose_level = 2
ui.visu = False
ui.visu_type = "vrml_file_only"
ui.visu_filename = "VisuFile.wrl"
ui.visu_verbose = False
sim.random_engine = "MersenneTwister"
sim.random_seed = 123456
m = gate.g4_units.m
mm = gate.g4_units.mm
cm3 = gate.g4_units.cm3
Bq = gate.g4_units.Bq
kBq = Bq * 1e3
BqmL = Bq / cm3
sec = gate.g4_units.second
ps = gate.g4_units.ps
keV = gate.g4_units.keV
MeV = gate.g4_units.MeV

'''
Geometry
'''
# world
sim.world.size = [2 * m, 2 * m, 2 * m]
sim.world.material = "G4_AIR"
sim.volume_manager.add_material_database("../MC/GateMaterials_pet.db") # add specified materials
# add a PET Biograph system
pet = pet_biograph.add_pet(sim, "pet")

# CT image
ct = sim.add_volume('Image', 'ct')
ct.mother = 'world'
ct.image = '../data/ct_cropped.mhd'
ct_info = gate.image.read_image_info(ct.image)
ct.voxel_materials = [
    [0,   4,  "G4_AIR"],
    [5,   6,  "G4_EYE_LENS_ICRP"],
    [7,   7,  "G4_SKIN_ICRP"],
    [8,  43,  "G4_BRAIN_ICRP"],
    [44, 51,  "G4_BLOOD_ICRP"],
    [52, 54,  "G4_LUNG_ICRP"],
    [55, 61,  "G4_MUSCLE_SKELETAL_ICRP"],
    [62, 63,  "G4_B-100_BONE"],
    [64, 66,  "G4_BONE_COMPACT_ICRU"],
    [67, 111, "G4_TISSUE_SOFT_ICRP"]
]

'''
Physics
'''
sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3"
sim.physics_manager.enable_decay = True
sim.physics_manager.set_production_cut("world", "all", 1 * m)
sim.physics_manager.set_production_cut(f"{pet.name}_crystal", "all", 0.1 * mm)

'''
Sources
'''
source = sim.add_source("VoxelsSource", "source")
source.mother = ct.name
source.image = '../data/activity_ga_lesion_norm_cropped.mhd'
total_yield = gate.sources.generic.get_rad_yield("Ga68")
source.activity = 158.6 * kBq * total_yield / sim.number_of_threads
source.half_life = 4062.6 * sec
source.particle = "e+"
source.energy.type = "F18"
#source.energy.type = "Ga68"
source.direction.type = "iso"
'''
Actors
'''
# hits collection
crystal = sim.volume_manager.volumes[f"{pet.name}_crystal"]
hc = sim.add_actor("DigitizerHitsCollectionActor", "Hits")
hc.mother = crystal.name
hc.output = f"pet1s.root"
hc.attributes = ["PostPosition", "TotalEnergyDeposit", "PreStepUniqueVolumeID", "GlobalTime",]

# singles collection
sc = sim.add_actor("DigitizerReadoutActor", "Singles1")
sc.output = f"pet1s.root"
sc.input_digi_collection = hc.name
sc.group_volume = crystal.name
sc.discretize_volume = crystal.name
sc.policy = "EnergyWeightedCentroidPosition"

# Detection Efficiency
ea = sim.add_actor("DigitizerEfficiencyActor", "Singles2")
ea.output = f"pet1s.root"
ea.input_digi_collection = "Singles1"
ea.efficiency = 0.86481

# energy blurring
eb = sim.add_actor("DigitizerBlurringActor", "Singles3")
eb.output = f"pet1s.root"
eb.input_digi_collection = "Singles2"
eb.blur_attribute = "TotalEnergyDeposit"
eb.blur_method = "InverseSquare"
eb.blur_resolution = 0.112
eb.blur_reference_value = 511 * keV

# time blurring
tb = sim.add_actor("DigitizerBlurringActor", "Singles4")
tb.output = f"pet1s.root"
tb.input_digi_collection = "Singles3"
tb.blur_attribute = "GlobalTime"
tb.blur_method = "Gaussian"
tb.blur_fwhm = 220.0 * ps

# EnergyWindows
ew = sim.add_actor("DigitizerEnergyWindowsActor", "Singles5")
ew.mother = crystal.name
ew.input_digi_collection = "Singles4"
ew.channels = [{"name": ew.name, "min": 425 * keV, "max": 650 * keV}] # 511KeV
ew.output = f"pet1s.root"

# projection
proj = sim.add_actor("DigitizerProjectionActor", "Projection")
proj.input_digi_collections = [ew.name]
proj.spacing = [ct_info.spacing[0] * mm, ct_info.spacing[1] * mm]
proj.size = [ct_info.size[0], ct_info.size[1]]
proj.origin_as_image_center = False
# proj.plane = 'XY' # not implemented yet
proj.output = "projection.mhd"

# add stat actor
stats = sim.add_actor("SimulationStatisticsActor", "Stats")
stats.track_types_flag = True

'''
Simulation
'''
sim.run_timing_intervals = [[0, 2.0 * sec]]
sim.user_hook_after_init = check_production_cuts # set user hook to dump production cuts from G4
sim.run()
stats = sim.output.get_actor("Stats")
print(stats)
jizhang02 commented 7 months ago

It seems to be related to the energy blurring. When I change the resolution to 0, the bin bar of 511Kv stands out again in the case of Ga68.

# energy blurring
eb = sim.add_actor("DigitizerBlurringActor", "Singles3")
eb.output = f"pet1s.root"
eb.input_digi_collection = "Singles2"
eb.blur_attribute = "TotalEnergyDeposit"
eb.blur_method = "InverseSquare"
eb.blur_resolution = 0 #0.112
eb.blur_reference_value = 511 * keV