Local time of a radioactive isotope is set to 0 at the end of the track with G4OpticalPhysics registered

Local time of a radioactive isotope is set to 0 at the end of the track with G4OpticalPhysics registered.

Here is some key code:

//PhysicsList in main.cc
G4VModularPhysicsList* physicsList = new QGSP_BIC_AllHP(phyVer);
physicsList->RegisterPhysics(new G4OpticalPhysics(phyVer));

//TrackingAction in TrackingAction.cc
void TrackingAction::PostUserTrackingAction(const G4Track* track)
{
    auto processDefinedStep = track->GetStep()->GetPostStepPoint()->GetProcessDefinedStep();
    G4String processName = "";
    if (processDefinedStep)
        processName = processDefinedStep->GetProcessName();

    if (processName == "Radioactivation" || processName == "RadioactiveDecay") 
    {
        G4double decayTime = track->GetLocalTime() / us;
    }
}

decayTime is 0 if G4OpticalPhysics is registered to physicsList, but track->GetGlobalTime() is not.
Unregistering G4OpticalPhysics, decayTime is a valid value and equal to track->GetGlobalTime() in a primary particle situation.
Is this a bug? Should it be reported?

Geant4 Version: 11.2.1
Operating System: Windows 11 23H2
Compiler/Version: Visual Studio 2022 17.9
CMake Version: 3.28.3


1 Like

It sounds like a bug to me. It’s more likely to be fully addressed and fixed if you can report more specifically which optical process is the offender, and also, if there is an existing example where you can reproduce the bug (assuming there’s an example where track->GetLocalTime() is saved to output).

To isolate the process, you’d do something like the following, with G4OpticalPhysics included in your physics list. First, modify your TrackingAction above with a reporting condition,

  if (decayTime <= 0.)
    G4cout << "##### Invalid local decay time " << decayTime << " us" << G4endl;

Then set up your macro with the additional commands

/run/initialize
/process/inactivate OpAbsorption
/process/inactivate OpRayleigh
/process/inactivate OpMieHG
/process/inactivate OpBoundary
/process/inactivate OpWLS
/process/inactivate OpWLS2
/process/inactivate Cerenkov
/process/inactivate Scintillation
/run/beamOn 1

/process/activate Scintillation
/run/beamOn 1

/process/inactivate Scintillation
/process/activate Cerenkov
/run/beamOn 1

and so on through each of the optical processes in turn.

I would not expect the other processes, which only operate on optical photons, to do anything. But those last two are where optical photons get created from other tracks, so they could be the culprits.

1 Like

I executed the simulation as you suggested. All processes were inactivated and then activated one by one. track->GetLocalTime() is always 0 which is abnormal.

Yup, that doesn’t sound normal. To be clear, even with all the optical processes inactivated, you’re still getting local time 0? If you don’t hear from anyone else about this (@dsawkey ?), I’d suggest creating a bug report about this.

Some thoughts. I don’t know what is happening here.

  1. Echoing Mike’s comment about giving as much detail as possible on how to reproduce this. Especially as I don’t know much about Geant4 radioactive decay.

  2. It would be great if you test with one of the radioactivity examples. However they use GetGlobalTime rather than GetLocalTime. I haven’t looked into the examples in detail. If or when the bug report comes to me, I would start with those examples, modifying to use localtime.

  3. The code in the first post is incorrect. decayTime goes out of scope before it is used. If you try to record it later, it will be nonsense.I realize the code in your first post is likely highly edited from what you actually use, but see point 1.

  4. Similarly, maybe test for particle type before recording local time.

  5. In my opinion, it’s always good to report strange results as a bug. But first please do what you can to check the problem isn’t with your code. And see again point 1–lots of bug reports get closed because the developer cannot reproduce the issue. Be prepared for some discussion.

As you suggested, I tested with the example rdecay01 (GitHub - zhikaiyici/LocalDecayTime).
As mentioned above, track->GetLocalTime() is 0 if G4OpticalPhysics is registered to the physicsList.

Result without G4OpticalPhysics:

G4WT3 > *********************************************************************************************************
G4WT3 > * G4Track Information:   Particle = Ne24,   Track ID = 1,   Parent ID = 0
G4WT3 > *********************************************************************************************************
G4WT3 >
G4WT3 > Step#       X          Y          Z         KineE         dEStep    StepLeng   TrakLeng    Volume     Process
G4WT3 >     0       0 fm       0 fm       0 fm       0 eV           0 eV       0 fm       0 fm      World   initStep
G4WT3 >     1       0 fm       0 fm       0 fm       0 eV           0 eV       0 fm       0 fm      World   Radioactivation
G4WT3 >
    :----- List of secondaries ----------------
G4WT3 >    Na24[472.207]:  energy = 22.34 eV   time = 21.92 s
G4WT3 >        anti_nu_e:  energy = 1.299 MeV  time = 21.92 s
G4WT3 >               e-:  energy = 694.4 keV  time = 21.92 s
G4WT3 >     :------------------------------------------ 

G4WT3 > ##### Global decay time: 21.9193 s
G4WT3 > ##### Local decay time: 21.9193 s

Result with G4OpticalPhysics:

G4WT2 > *********************************************************************************************************
G4WT2 > * G4Track Information:   Particle = Ne24,   Track ID = 1,   Parent ID = 0
G4WT2 > *********************************************************************************************************
G4WT2 >
G4WT2 > Step#       X          Y          Z         KineE         dEStep    StepLeng   TrakLeng    Volume     Process
G4WT2 >     0       0 fm       0 fm       0 fm       0 eV           0 eV       0 fm       0 fm      World   initStep
G4WT2 >     1       0 fm       0 fm       0 fm       0 eV           0 eV       0 fm       0 fm      World   Radioactivation
G4WT2 >
    :----- List of secondaries ----------------
G4WT2 >    Na24[472.207]:  energy = 73.33 eV   time = 5.212 s
G4WT2 >        anti_nu_e:  energy = 610.1 keV  time = 5.212 s
G4WT2 >               e-:  energy = 1.383 MeV  time = 5.212 s
G4WT2 >     :------------------------------------------

G4WT2 > ##### Global decay time: 5.21165 s
G4WT2 > ##### Local decay time: 0 ps

For more details, please see the GitHub repository (GitHub - zhikaiyici/LocalDecayTime).

Geant4 Version: 11.2.1

Thanks for running this test! It does look like something is wrong. We’ll have to investigate. Could you please file a bug report.