Physics list Shielding does not seem to work about opengate

[dsarrut]

Thanks Herman ! We will investigate. Note that phys list management is currently under work. We should create an additional tests to check Shielding

[nkrah]

Thanks Hermann for reporting this and for testing Gate10!

I looked into the issue you described. First, the correct physics list names are indeed "Shielding" and "Shielding_EMZ". The "_EMZ" suffix instructs Geant4 to replace the EM physics constructor inside the list by G4EmStandardPhysics_option4 (default is G4EmStandardPhysics). I guess you know that, but I write this here for whoever stumbles across the issue.

To know which physics lists are available, you get use these 3 lines of python code, either in a script, or from an interactive terminal like ipython or jupyter notebook, or similar:

import opengate sim = opengate.Simulation() print(sim.physics_manager.dump_available_physics_lists())

GATE 10 relies on Geant4 to build the physics list, namely via the G4PhysListFactory. To check if this works correctly, I ran the following simulation, in analogy to the test013_XXX scripts in the opengate/tests/src folder:

import opengate as gate
from test013_phys_lists_helpers import create_pl_sim

sim = create_pl_sim()
sim.user_info.g4_verbose = True
sim.source_manager.user_info_sources.pop("ion1")
sim.source_manager.user_info_sources.pop("ion2")
sim.physics_manager.physics_list_name = "Shielding"

sim.run()

stats = sim.output.get_actor("Stats")
print("STATS:")
print(stats)

Based on the verbose G4 output, I judge that the correct phys list has been selected by G4. On my MacBook, the simulation takes 8-9 sec. I get

Events    9938
Tracks    65433
Step      242800

and

Track types: {'e+': 1515, 'e-': 45496, 'gamma': 18422}

You can also use
sim.physics_manager.physics_list_name = "G4EmStandardPhysics",
in which case the simulation takes about 1 sec.

@GitFuchs : you mention that "no data is generated". Can you specify? Or better: provide an example simulation to reproduce this?

[nkrah]

Quick update: I have been looking at another example where the Shielding list does not yield the expected results. The underlying issues are probably similar if not the same. Therefore, do not bother putting together a minimal example for this issue at the moment. I will post back once I understood the other example.

[nkrah]

Here is a minimal working example which reproduces the error:

import opengate as gate
from scipy.spatial.transform import Rotation

def start_simulation(phlist, source_particle='carbon'):
    # units
    cm = gate.g4_units("cm")
    mm = gate.g4_units("mm")
    MeV = gate.g4_units("MeV")
    
    # create the simulation
    sim = gate.Simulation()
    
    # main options
    ui = sim.user_info
    #ui.start_new_process = True
    ui.g4_verbose = True
    ui.g4_verbose_level = 2
    ui.visu = False
    ui.random_seed = 12365478910
    ui.random_engine = "MersenneTwister"

    # physics
    sim.physics_manager.physics_list_name = phlist
    
    #  change world size
    world = sim.world
    world.size = [600 * cm, 500 * cm, 500 * cm]
    
    # target
    phantom = sim.add_volume("Box", f"phantom_{phlist}")
    phantom.size = [200 * mm, 200 * mm, 200 * mm]
    phantom.translation = [-100* mm, 0, 0]
    phantom.material = "G4_WATER"
    phantom.color = [0, 0, 1, 1]
    
    # source
    source = sim.add_source("GenericSource", "mysource")
    if source_particle == 'carbon':
        source.energy.mono = 1440 * MeV
        source.particle = "ion 6 12"
    elif source_particle == 'proton':
        source.energy.mono = 40 * MeV
        source.particle = "proton"
    else: 
        raise ValueError(f"Unknown source name: {source_particle}")

    source.position.type = "disc" 
    source.position.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix()
    source.position.sigma_x = 8 * mm
    source.position.sigma_y = 8 * mm
    source.direction.type = "momentum"
    source.direction.momentum = [-1, 0, 0]
    source.n = 10
    # ***

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

    # run simulation
    sim.run(start_new_process=True)
    
    # print results at the end
    stat = sim.output.get_actor("Stats")
    print(stat)


# start_simulation("QGSP_BIC_EMZ")
# start_simulation("QGSP_BIC_HP_EMZ", source_particle='proton')
start_simulation("Shielding", source_particle='carbon')

[nkrah]

What I observe:
With the lists Shielding and QGSP_BIC_HP_EMZ and carbon ions as source particles, the simulation does effectively not transport any particles. In other words: The step, track, and event count is 10, i.e. the number of primaries, and only C12 appears as track type. With protons, many different secondaries are produced, as expected.

Both lists have in common that they use G4HadronElasticPhysicsHP as hadronic component of the physics list (by default).

One can investigate further changing the function GateSimulationStatisticsActor::SteppingAction (in opengate/core/opengate_core/opengate_lib/) into:

void GateSimulationStatisticsActor::SteppingAction(G4Step *step) {
  // Called every step
  threadLocalData.Get().fStepCount++;
  std::cout << "DEBUG: SteppingAction:" << std::endl;
  std::cout << "DEBUG: Name of step defining process: " << step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() << std::endl;
}

This prints out the process which defines the next step to be taken.
For carbon ions, I get
"DEBUG: Name of step defining process: RadioactiveDecay"
for all events. In other words, the primary particles decay immediately after the source generates them and no real step is ever taken.

This observation is in line with other reports we had related to PR138, where manually removing "RadioActiveDecay" from the physics list seemed stop the issue from appearing.

I have hacked some verbose output into the G4 function G4SteppingManager::DefinePhysicalStepLength which has shown me that the physical interaction length proposed by the process "RadioActiveDecay" is of the order E-308 for carbons but E+308 for protons. Therefore, in the former case, this process is always chosen, while it is always ignored in the latter case. At least I think. Maybe some numerical issues?

This needs further debugging on the Geant4 side. I will post more insights once I have any.

[nkrah]

Found the bug:
In core/opengate_core/opengate_lib/GateGenericSource.cpp,
GateGenericSource::InitializeHalfTime sets the source particle's life time to 0, while it should honor the value retrieved from the user info, which defaults to -1 (= stable ion).

Should be p->SetPDGLifeTime(fHalfLife); instead of p->SetPDGLifeTime(0);

I'll make a PR to fix this ASAP.

[nkrah]

Just to add an explanation:
The physics list QGSP_BIC (and thus QGSP_BIC_EMZ) do not include G4RadioactiveDecayPhysics, while the physics lists QGSP_BIC_HP and Shielding do. Therefore, the function G4RadioactiveDecayPhysics::GetMeanFreePath does not get called in the former case and the above described bug does not have any effect.

[nkrah]

In fact, any time the physics list includes G4RadioactiveDecayPhysics and uses a generic ion as source, such as "ion 6 12", the bug will show up. This can be reproduced by using the minimal example above with physics list QGSP_BIC_EMZ (which by itself works correctly) and adding

sim.physics_manager.special_physics_constructors['G4RadioactiveDecayPhysics'] = True

In any case, the issues should be understood and will be closed once the PR is done.

[nkrah]

See PR #218