OpenGATE / opengate

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

DigitizerProjectionActor does not work on parallel world #350

Open buimhh opened 6 months ago

buimhh commented 6 months ago

I have a problem implementing a parallel world in SPECT. From my understanding, it may come from the DigitizerProjectionActor. I reused one exercise of David to demonstrate the error.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import opengate as gate
import opengate.contrib.spect.genm670 as spect_ge_nm670
from pathlib import Path

if __name__ == "__main__":
    # create the simulation
    sim = gate.Simulation()

    # main options
    sim.g4_verbose = False
    sim.visu = False
    sim.visu_type = "vrml"
    # sim.visu_type = "qt"
    sim.number_of_threads = 6
    sim.random_seed = "auto"

    # units
    m = gate.g4_units.m
    sec = gate.g4_units.second
    days = 3600 * 24 * sec
    cm = gate.g4_units.cm
    mm = gate.g4_units.mm
    nm = gate.g4_units.nm
    MeV = gate.g4_units.MeV
    keV = gate.g4_units.keV
    Bq = gate.g4_units.Bq
    kBq = 1000 * Bq
    MBq = 1000 * kBq

    # world size
    sim.world.size = [2 * m, 2 * m, 2 * m]
    sim.world.material = "G4_AIR"

    # waterbox
    wb = sim.add_volume("Box", "waterbox")
    wb.size = [30 * cm, 30 * cm, 30 * cm]
    wb.material = "G4_AIR"
    wb.color = [0, 0, 1, 1]  # blue

    spect, crystal = spect_ge_nm670.add_ge_nm67_spect_head(
        sim, "spect", collimator_type="megp", debug=sim.visu
    )

    sim.add_parallel_world("parallel_world")
    spect.mother = "parallel_world"
    spect.translation = [0, 0, -28 * cm]

    # spect digitizer
    channels = [
        {"name": f"spectrum", "min": 3 * keV, "max": 515 * keV},
        {"name": f"scatter", "min": 380 * keV, "max": 420 * keV},
        {"name": f"peak", "min": 400 * keV, "max": 450 * keV},
    ]
    # spect digitizer
    proj = spect_ge_nm670.add_digitizer(sim, "spect_crystal", channels) # issue come from this line
    # proj is the Projection image
    proj.output = Path("output") / "projection2_star.mhd"

    # Lu177 source (only the gammas)
    source = sim.add_source("GenericSource", "fake")
    source.particle = "gamma"
    source.mother = f"waterbox"

    source.energy.type = "spectrum_lines"
    source.energy.spectrum_weight = [0.114, 0.259]
    source.energy.spectrum_energy = [218.12 * keV, 440.45 * keV]
    source.position.type = "sphere"
    source.position.radius = 1 * mm
    source.position.translation = [0, 0, 0 * mm]
    source.direction.type = "iso"
    if sim.visu:
        sim.number_of_threads = 1
        source.activity = 2 * Bq
    else:
        source.activity = 4 * MBq / sim.number_of_threads
    source.half_life = 6.647 * days

    # add stat actor
    s = sim.add_actor("SimulationStatisticsActor", "stats")
    s.track_types_flag = True
    s.output = Path("output") / "stats2_star.txt"

    # phys
    sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3"
    sim.physics_manager.set_production_cut("world", "all", 1 * mm)

    # ---------------------------------------------------------------------
    # start simulation
    # sim.running_verbose_level = gate.EVENT
    time_per_slice = 1 * sec
    no_of_run = 5
    sim.run_timing_intervals = gate.runtiming.range_timing(0, time_per_slice * no_of_run, no_of_run)
    motion = sim.add_actor("MotionVolumeActor", f"Move_{spect.name}")
    motion.mother = spect.name
    motion.translations, motion.rotations, = gate.geometry.utility.volume_orbiting_transform(
        "z", 0, 180, no_of_run, spect.translation, spect.rotation)

    sim.run()

    # print stats
    stats = sim.output.get_actor("stats")
    print(stats)