Positron annihilation w/o two gammas and even annihilating?

furutaka · 2023-04-16T09:08:13.154Z

John_McFee:

Were the runs all done in multithreaded mode? If yes, redo them in single threaded mode to see if the problem goes away. If the problem exists in single threaded mode, file a bug report.

The answer is yes: all runs were done in multi-thread mode, with four, eight or sixteen threads. Then, I also tried to make it run in single thread with the same binary which was built with MT-enabled Geant4, as well as a binary recompiled with geant4-v11.1.1 with multi-thread disabled,

**************************************************************
 Geant4 version Name: geant4-11-01-patch-01    (10-February-2023)
                       Copyright : Geant4 Collaboration
                      References : NIM A 506 (2003), 250-303
                                 : IEEE-TNS 53 (2006), 270-278
                                 : NIM A 835 (2016), 186-225
                             WWW : http://geant4.org/
**************************************************************

<<< Geant4 Physics List simulation engine: FTFP_BERT

G4VModularPhysicsList::ReplacePhysics: G4EmStandard with type : 2 is replaced with G4EmLowEPPhysics
Visualization Manager instantiating with verbosity "warnings (3)"...
Visualization Manager initialising...
Registering graphics systems...

Again, such non-annihilating annihil events were observed, seven times with MT-enabled g4 run in single thread, and six times with MT-disabled one.

Kazuyoshi

John_McFee · 2023-04-16T18:27:14.210Z

Just as a point of curiosity, could you provide me your code or at least the modifications you’ve made to example B1 to produce the root files? I am running Fedora 37 latest kernel with Geant4 11.1.1 MT and would like to try duplicating your results with and without multithreading.

In any case, unless any other ideas are posted, I would suggest that you file a bug report.

dsawkey · 2023-04-17T21:13:23.104Z

Are you able to reproduce this in an unmodified Geant4 example? e.g. set /tracking/verbose 1, and if the file is too large, pipe through grep.

Maybe something is wrong with the analysis. What particle is this?

furutaka:

(track,parent,particle) = (4,1,1000320700), 2steps

furutaka:

the name of each process is gotten via G4Process::GetPrcessName():

Do you use theStep->GetProcessDefinedStep() to get the process that limits the step?

Do you enable either biasing or entanglement?

The annihilation code is here:
https://geant4.kek.jp/lxr/source/processes/electromagnetic/standard/src/G4eplusAnnihilation.cc#L160
You could add a G4cout statement (or breakpoint in debug mode) if the number of secondaries equals 0. Interestingly, the number of secondaries is checked a few lines down, and code is only run if the number is greater than 0. This implies the number might be zero.

Still, simply reading the code here,
https://geant4.kek.jp/lxr/source/processes/electromagnetic/standard/src/G4eeToTwoGammaModel.cc#L181
the primary should always be killed.

Another debugging option for rare occurence, is to check for the occurence in your UserSteppingAction (say), then save the random number seed. You can then re-run that Event, say in debug mode, or with increased verbosity. This is done in example extended/optical/LXe, for example.

furutaka · 2023-04-17T22:45:33.861Z

Hi @dswkey, thanks for your interest.

dsawkey:

Are you able to reproduce this in an unmodified Geant4 example? e.g. set /tracking/verbose 1, and if the file is too large, pipe through grep.

Not yet done; I will give it a try…

dsawkey:

Maybe something is wrong with the analysis. What particle is this?

(track,parent,particle) = (4,1,1000320700), 2steps

I believe it’s a ⁷⁰Ge nucleus; in my simulations, primary gamma rays were incident into a germanium crystal. According to the Monte Carlo Particle Numbering Scheme,

Nuclear codes are given as 10-digit numbers ±10LZZZAAAI.

dsawkey:

Do you use theStep->GetProcessDefinedStep() to get the process that limits the step?

Yes:

  const G4StepPoint *postPoint = step->GetPostStepPoint();
  const G4ThreeVector postPos = postPoint->GetPosition();
  const G4VProcess* process
    = step->GetPostStepPoint()->GetProcessDefinedStep();
  G4String procName = " UserLimit";
  if (process != nullptr) procName = process->GetProcessName();

dsawkey:

Do you enable either biasing or entanglement?

No.

Kazuyoshi

John_McFee · 2023-04-18T00:11:13.914Z

I ran the unmodified B1 example code with 10^6 1980keV incident gammas and /tracking/verbose 2. The output was huge as expected but I filtered it through grep (a few passes required). I found one weird event and I have attached that output.
double_annihil_example.txt (5.0 KB)

The weird part is in Track ID=3 where it seems that the annihilation of the e+ with no gammas emitted is followed by energy loss by ionization, then annihilation of the e+ again with 511keV gammas emitted.

Maybe I’m reading it wrong.

furutaka · 2023-04-18T05:19:38.360Z

furutaka:

dsawkey:

Are you able to reproduce this in an unmodified Geant4 example? e.g. set /tracking/verbose 1, and if the file is too large, pipe through grep.

Not yet done; I will give it a try…

The answer is YES.

exampleB1 of geant4-v11.1.1 (-DGEANT4_BUILD_MULTITHREADED=OFF) was run with the following macro which is a truncated/modified version of run2.mac of the example,

# Macro file for example B1
# 
# To be run preferably in batch, without graphics:
# % exampleB1 run2.mac
#
#/run/numberOfThreads 4
/run/initialize
#
/control/verbose 2
/run/verbose 2
/tracking/verbose 1
#
# gamma 6 MeV to the direction (0.,0.,1.)
# 100000000 events
#
/gun/particle gamma
/gun/energy 6 MeV
#
/run/printProgress 100
/run/beamOn 100000000

and filtered with the following awk script, (./exampleB1 run2mod10to8.mac | awk -f findStrangeAnnihil.awk)

# findStrangeAnnihil.awk:
/^* G4Track Information:   Particle = e\+/ {  # find positron tracks
    ANNIHIL = 0;  # to flag 'annihil' process included in the present track
    #
    # store the 'annihil' track data into array "LINE[]"
    #
    LINE[2] = $0; # "G4Track Information:" line
    maxIdx = 0;   # to point the line number of the final empty line
    # 3 : *********************************** ...
    # 4 :
    # 5 : Step#       X          Y          Z ...
    # 6 : 0   1.339 cm  -6.915 cm   13.43 cm  ...  initStep
    for(idx=3; idx<=6; idx++){ # 3: ***, 4: empty, 5: Step#..., 6: 0 ...
	getline;
	LINE[idx] = $0;
    }
    LINE[1] = LINE[3]; # separator ("*****...*****")
    ### store the track data after 'initStep'
    while(NF != 0){ # find the last empty line. (already done for the "initStep" line)
	getline;
	if ($NF=="annihil") {ANNIHIL = 1;}  # 'annihil' process in this track
	LINE[idx] = $0;    # re-use the loop variable "idx"
	if(NF != 0) idx++; # advance if not the last empty line (end of track)
    }
    maxIdx = idx; # the number of stored lines including the last empty one
    #
    # check if the 'annihil' process is the last one (and annihilating)
    #
    if (ANNIHIL==1) { # 'annihil' process included
	numF = split(LINE[maxIdx-1], data); # split/separate the last line into an array "data"
	if (data[numF] != "annihil")           # 'annihil' not in the last step; that is,
	for (j=1;j<=maxIdx;j++) print LINE[j]; #  this is a non-annihilating 'annihil' track.
    }
}

a total of 193 such non-annihilating annihil tracks in 10⁸ primaries were observed.

Kazuyoshi

John_McFee · 2023-04-18T15:20:06.536Z

I tried to duplicate furutaka’s results on exampleB1 just now. I am running geant4-v11.1.1 distribution, built with -DGEANT4_BUILD_MULTITHREADED=ON. I used exampleB1 from the distribution with furutaka’s run2mod10to8.mac file (with
‘#/run/numberOfThreads 4’ changed to ‘/run/numberOfThreads 1’) to force single threading. When I run the command

./exampleB1 run2mod10to8.mac | awk -f findStrangeAnnihil.awk

I get NO strange annihilation events output. Still, as I said in an earlier post, I have found one strange event in 10^6 primaries in an earlier run. Sorry to muddy the waters.

furutaka · 2023-04-18T20:13:15.424Z

John_McFee:

I tried to duplicate furutaka’s results on exampleB1 just now. I am running geant4-v11.1.1 distribution, built with -DGEANT4_BUILD_MULTITHREADED=ON. I used exampleB1 from the distribution with furutaka’s run2mod10to8.mac file (with
‘#/run/numberOfThreads 4’ changed to ‘/run/numberOfThreads 1’) to force single threading. When I run the command

./exampleB1 run2mod10to8.mac | awk -f findStrangeAnnihil.awk

I get NO strange annihilation events output. Still, as I said in an earlier post, I have found one strange event in 10^6 primaries in an earlier run. Sorry to muddy the waters.

John, as your double_anihil_example.txt clearly indicated, in MT case, a word like G4WT0 > is prepended to the * G4Track Information: line and therefore the line can NOT be detected by the REGEXP of my awk script as it is. Was it modified to so that it can also detect such lines in MT case?

Kazuyoshi

dsawkey · 2023-04-18T20:29:38.146Z

Thanks Kazuyoshi and John for verifying that this occurs with example B1. I ran a single event (same macro as given above but setting the random number state) with verbose tracking. The output is attached.

I believe there is a bug in tracking.

The anomolous event happens at step 7 (line 3735 and following). There is:

 >>DefinePhysicalStepLength (List of proposed StepLengths):
    ++ProposedStep(PostStep ) = 5.82593e+288 pc  : ProcName = positronNuclear (No ForceCondition)
    ++ProposedStep(PostStep ) = 5.82593e+288 pc  : ProcName = CoulombScat (No ForceCondition)
    ++ProposedStep(PostStep ) =   43.3825 um  : ProcName = annihil (No ForceCondition)
    ++ProposedStep(PostStep ) =   3.02858 mm  : ProcName = eBrem (No ForceCondition)
    ++ProposedStep(PostStep ) =   2.78022 cm  : ProcName = eIoni (No ForceCondition)
    ++ProposedStep(PostStep ) = 5.82593e+288 pc  : ProcName = Transportation (Forced)
    ++ProposedStep(AlongStep) =     1.646 mm  : ProcName = eIoni (CandidateForSelection)
    ++ProposedStep(AlongStep) =   43.0397 um  : ProcName = msc (NotCandidateForSelection)
    ++ProposedStep(AlongStep) =   43.0397 um  : ProcName = Transportation (CandidateForSelection)

The shortest proposed step is for msc and Transportation, yet ‘Process defined step’ is set to annihil. However, as expected from the proposed step, annihil->PostStepDoIt() is not called, no secondaries are created, the primary e+ is not killed.

I don’t understand why this is happening.

out10648.txt (435.7 KB)

dsawkey · 2023-04-18T20:32:08.651Z

If you add the lines below to B1’s SteppingAction, it will be easy to locate these events in the verbose output. (partly copied from Kazuyoshi.)

void SteppingAction::UserSteppingAction(const G4Step* step)
{
  const G4StepPoint *postPoint = step->GetPostStepPoint();
  const G4ThreeVector postPos = postPoint->GetPosition();
  const G4VProcess* process
    = step->GetPostStepPoint()->GetProcessDefinedStep();
  G4String procName = " UserLimit";
  if (process != nullptr) procName = process->GetProcessName();
  if (procName == "annihil") {
    const auto en = step->GetTrack()->GetKineticEnergy();
    if (en>0) {
      G4cout << "Annihilation! " << en << G4endl;
      //G4RunManager::GetRunManager()->rndmSaveThisEvent();
    }
  }

furutaka · 2023-04-18T21:50:28.113Z

Well, then, what shall we do?

By the way,

dsawkey:

The anomolous event happens at step 7 (line 3735 and following).

You mean step 7 starting from line 4191 of your attached file, right?

Kazuyoshi

John_McFee · 2023-04-18T22:24:22.064Z

Oops. I missed that. I should have looked closer at the awk script. Thanks.

John_McFee · 2023-04-19T15:33:14.497Z

I modified Furutaka’s awk script to handle the different output from the multithreaded Geant4 build and re-ran exampleB1 in single thread mode. I got 164 “non-annihilating” e+ events compared to Furutaka’s 193 for a single threaded Geant4 build. The difference in numbers is likely just statistics and differing intial seeds.

dsawkey · 2023-04-19T19:08:53.792Z

furutaka:

Well, then, what shall we do?

It would be great if you filed a bug report. I think we need input from the tracking people to proceed.

furutaka:

dsawkey:

The anomolous event happens at step 7 (line 3735 and following).

You mean step 7 starting from line 4191 of your attached file, right?

Lines 3735-4191, which I believe is step 6.

furutaka · 2023-04-19T23:17:18.919Z

dsawkey:

It would be great if you filed a bug report. I think we need input from the tracking people to proceed.

Just filled a report

Thanks for your help and interest, @John_McFee and @dsawkey !

Yours,
Kazuyoshi

civanch · 2023-04-20T10:46:30.335Z

Hello,

thank you for this post. This is a copy of the comment to the bug-report. In fact, it is correct to make a bug report even if you are not 100% sure that it is the bug.

In Geant4 it is allowed to define a track step by some process but before sampling of final state PostStep some conditions are checked and a decision to continue tracking may be made. One of the reasons for that - so called “integral approach for tracking”, which allows to take into account dependence of a cross section on energy due to energy loss of a particle at the step.

So, likely what you observe is normal, with low probability this may happens.

Let us think how to double check this assumption before closing this bug report.
For example, you may use UI command: “/process/em/integral false”. If in such run this type of events does not happens anymore, then my explanation is correct.

VI

furutaka · 2023-04-20T20:01:07.799Z

Dear @civanch,

civanch:

For example, you may use UI command: “/process/em/integral false”. If in such run this type of events does not happens anymore, then my explanation is correct.

Thanks… Indeed, by adding /process/em/integral false before /run/initialize, such strange events dispappeared (no such events in 10⁸ primaries).

Where can I find the description in the G4 documentation on this command and the “integral approach for tracking”?

Kazuyoshi

furutaka · 2023-04-20T23:11:27.857Z

I could only find two references to the “integral approach for tracking” in the Physics Reference manual, here and here (and not in other manuals).

And before the release of Geant4 version 11, the setting seemed to be noted in the output of a simulation:

$ for dir in geant4{.,-}*[0-9]; do echo -n ${dir};grep "integral approach" ${dir}/source/processes/electromagnetic/utils/src/G4EmParameters.cc;done
geant4.10.06.p03  os << "Use integral approach for tracking                 " <<integral << "\n";
  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4.10.07  os << "Use integral approach for tracking                 " <<integral << "\n";
  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4.10.07.p01  os << "Use integral approach for tracking                 " <<integral << "\n";
  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4.10.07.p02  os << "Use integral approach for tracking                 " <<integral << "\n";
  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4.10.07.p03  os << "Use integral approach for tracking                 " <<integral << "\n";
  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4-v11.0.0  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4-v11.0.1  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4-v11.0.2  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4-v11.0.3  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4-v11.0.4  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4-v11.1.0  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";
geant4-v11.1.1  os << "X-section factor for integral approach             " <<lambdaFactor << "\n";

John_McFee · 2023-04-21T02:22:12.012Z

I also added the /process/em/integral false before /run/initialize and got no strange events in 10^8 primaries with my multithreaded Geant4 11.1.1 build (but running B1example in single threaded mode).

dsawkey · 2023-04-21T13:55:20.947Z

Thanks Vladimir for the explanation!