Problem 2240

Summary: unstable particles sometimes decay prematurely in-flight, triggered by boundary crossing
Product: Geant4 Reporter: Richard Jones <richard.t.jones>
Component: trackingAssignee: asai
Status: RESOLVED FIXED    
Severity: major    
Priority: P4    
Version: 10.6   
Hardware: All   
OS: All   
Attachments: Example of in-flight decay position distribution, showing spike which should not be there.

Description Richard Jones 2020-04-29 09:46:35 CEST
Created attachment 618 [details]
Example of in-flight decay position distribution, showing spike which should not be there.

The GlueX experiment simulation uses the layered-mass geometry feature of G4. This means that we are using CoupledTransportation and ParallelWorld transportation processes. Under certain circumstances, these two act together in a way that triggers early in-flight decays of unstable particles (kaons in my test) at boundary crossings. These particles still have plenty of energy left, so they are not stopping, and I know these decays are premature because the particle still has plenty of decay distance left (see below).

I added debug statements in G4SteppingManager::DefinePhysicalStepLength() and this is a summary of what I see during one of these premature decays. 

post step process loop:
   ParallelWorld1 returns PIL=1.7976931348623e+308
   Cerenkov returns PIL=1.7976931348623e+308
   kaon-Inelastic returns PIL=2284.2840915216
   hadElastic returns PIL=3492.0118816984
   G4Decay returns PIL=35.023131390599
   CoulombScat returns PIL=53424.598349673
   hPairProd returns PIL=1.7976931348623e+308
   hBrems returns PIL=1.7976931348623e+308
   hIoni returns PIL=42.937647919513
   CoupledTransportation returns PIL=1.7976931348623e+308
along step process loop:
   hIoni returns PIL=209.13485484738
   msc returns PIL=35.021051113073, NotCandidateForSelection
   ParallelWorld1 returns PIL=2.527785383601, NotCandidateForSelection
   CoupledTransportation returns PIL=2.527785383601

At exit, proc selected is Decay, PIL=2.527785383601, and fStepStatus=fPostStepDoItProc. This triggers a step of 2.628mm to the edge of the volume, where it immediately decays. It should go another 32.4mm before decaying.

Clearly Decay is not the step limiter here. There is no decay pending at step distance 2.528mm. The winner was ParallelWorld1. But that process is not selected because it returns from AlongStepGPIL() with parameter selection == NotCandidateForSelection. The reason for that is probably because it wants to defer to CoupledTransportation, which also has a boundary at this same location. But because ParallelWorld1 has already lowered the PhysicalStep to 2.585mm, the CoupledTransportation step limit is never seen.

Normally this condition is protected against in the AlongStepProcess loop by the following assignment to fStepStatus that happens whether the process returns NotCandidateForSelection or not.

In G4ProcessManager2.cc, line 261:

         // Transportation is assumed to be the last process in the vector
      if(kp == MAXofAlongStepLoops-1)
        fStepStatus = fGeomBoundary;

The bug here is that ParallelWorld1 is also a transportation process, even though it is not the last process in the AlongStepProcs vector. If both ParallelWorld1 and CoupledTransportation share a common volume boundary then ParallelWorld1 will report it first and become the new step limiter, but fGeomBoundary will not be set. In that case, the previous value fPostStepDoItProc remains in place, even though the step size has now been reduced to that set by the geometry, not the decay process. Afteward when CoupledTransportation reports the same along-step length as was already set by ParallelWorld1, it is ignored, and fStepStatus is never actually set to fGeomBoundary. I propose replacing the code snippet above with the following.

In G4ProcessManager2.cc, line 261 (in place of the above code)
        G4ProcessType ptype = fCurrentProcess->GetProcessType();
       if (ptype == fTransportation || ptype == fParallel)
         fStepStatus = fGeomBoundary;

I tested this code, and it eliminated the premature decays. DefinePhysicalStepLength still returns with Decay as the step-limiting process, but the fStepStatus is now fGeomBoundary which is enough for to make the framework recognize that Decay is not actually happening at the end of this step.
Comment 1 asai 2020-05-14 22:22:20 CEST
Hi,

Thanks for reporting this. I'm trying to reproduce the issue with simplified test.

If both CoupledTransportation and ParallelWorldProcess are limiting the step, ParallelWorldProcess will return slightly larger distance to avoid from being selected (and returns NotCandidateForSelection as you see).

https://geant4.kek.jp/lxr/source/processes/scoring/src/G4ParallelWorldProcess.cc#L331

So, we need to understand why ParallelWorldProcess returned exactly the same distance as CoupledTransportation in your case.

I will update you once I succeed to reproduce the issue.

Kind regards,
Makoto
Comment 2 asai 2021-08-17 19:20:01 CEST
Thank you once again for your detailed report. Though I could not still succeed to reproduce the symptom, I believe I could spot the issue. The fix will be included in the next public release.