Positron annihilation w/o two gammas and even annihilating?

furutaka · April 12, 2023, 5:52am

Hi!

I’m calculating responses of a germanium crystal to gamma rays (up to a few MeV at present).
I’m still using version 11.0.2 of Geant4. To do the simulation, I used examples/basic/B1 as its base and modified. The physcs lists used are

  // Physics list
  G4VModularPhysicsList* physicsList = new FTFP_BERT;
  physicsList->ReplacePhysics(new G4EmLowEPPhysics());
  physicsList->SetVerboseLevel(1);
  runManager->SetUserInitialization(physicsList);

In the simulation, all the tracks and reactions are recorded in ROOT Ntuples and analyzed afterwards.

While analyzing the results I found some strange positron annihilation events in which two gammas are NOT emitted, and the positron itself is NOT annihilated but continues to exist.
The following is an example observed when a 1980-keV gamma ray is incident on the germanium crystal (the name of each process is gotten via G4Process::GetPrcessName():

evtID = 35372891
  (track,parent,particle) = (1,0,22), 3steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,        -150,      1980,         0,    Envelope,      initStep
               0,           0,           0,      1980,         0,    Envelope,Transportation
               0,           0,     58.7526,         0,         0,      Shape2,          conv
  (track,parent,particle) = (4,1,1000320700), 2steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,     58.7526, 0.0108039,         0,      Shape2,      initStep
    -3.99833e-07, -9.3874e-07,     58.7526,         0, 0.0108039,      Shape2,       ionIoni
  (track,parent,particle) = (3,1,-11), 17steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,     58.7526,   753.531,         0,      Shape2,      initStep
       0.0166704,     0.13674,     58.7944,   663.003,   90.5285,      Shape2,           msc
        0.094517,    0.183877,      58.834,   593.345,   69.6576,      Shape2,           msc
        0.151473,    0.240856,     58.8636,   530.803,   62.5418,      Shape2,           msc
        0.210505,    0.212396,     58.9062,   496.691,   34.1126,      Shape2,           msc
        0.241111,    0.240432,     58.9425,   460.527,   36.1642,      Shape2,           msc
        0.242372,    0.273407,      58.955,   437.386,   23.1408,      Shape2,           msc
        0.236917,    0.304803,     58.9687,   346.861,   90.5251,      Shape2,           msc
        0.223654,    0.334319,     58.9763,   317.055,   29.8053,      Shape2,           msc
        0.203938,    0.353037,     58.9925,   295.514,   21.5409,      Shape2,           msc
         0.20379,    0.351938,          59,   290.127,    5.3873,      Shape2,Transportation
         2.31919,    -5.45019,     76.4093,   287.013,   2.02656,    Envelope,         eIoni
         8.32559,     -16.418,     104.733,    282.01,   5.00226,    Envelope,       annihil
         18.5236,    -30.5561,     142.525,   272.934,   9.07652,    Envelope,           msc
         18.6629,    -30.6433,     142.936,   271.474, 0.0248335,    Envelope,         eIoni
         21.8451,    -32.3272,         150,   270.417,   1.05696,    Envelope,Transportation
         35.2265,     -39.382,         180,    264.86,   5.55707,       World,    OutOfWorld
  (track,parent,particle) = (6,3,11), 2steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
         18.6629,    -30.6433,     142.936,   1.43527,         0,    Envelope,      initStep
         18.6329,     -30.622,     142.937,         0,   1.43527,    Envelope,         eIoni
  (track,parent,particle) = (5,3,11), 2steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
         2.31919,    -5.45019,     76.4093,   1.08812,         0,    Envelope,      initStep
         2.33994,    -5.45897,     76.4025,         0,   1.08812,    Envelope,         eIoni
  (track,parent,particle) = (2,1,11), 3steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,     58.7526,    204.46,         0,      Shape2,      initStep
     -0.00072669,  -0.0103068,     58.7736,   59.6119,   144.848,      Shape2,         eIoni
      0.00131827, -0.00806186,     58.7733,         0,   59.6119,      Shape2,         eIoni

(the volume names “Shape2” and “Envelope” stand for germanium crystal and its surrounding air, respectively)

Such events occured for several times out of 10⁸ primaries.
Are these really positron annihilation evens?
Any thoughts?

Thanks in advance.
Yours,
Kazuyoshi

John_McFee · April 13, 2023, 9:27pm

Hi,
It looks like the positron escapes the world with ~264 keV of kinetic energy before it has a chance to annihilate with an electron:
“35.2265, -39.382, 180, 264.86, 5.55707, World, OutOfWorld”

So there should not be any annihilation gamma rays, as is the case.

I do not know what process “annihil” refers to in:
“8.32559, -16.418, 104.733, 282.01, 5.00226, Envelope, annihil”

It cannot be the annhilation of the positron with an electron because the positron continues to exist after the process occurs. It is at far too low an energy for any exotic processes to be happening. So I am puzzled.

All the other entries seem to make sense.

furutaka · April 13, 2023, 10:41pm

Hi @John_McFee. Thanks for the response.

looks like the positron escapes the world with ~264 keV of kinetic energy before it has a chance to annihilate with an electron:
“35.2265, -39.382, 180, 264.86, 5.55707, World, OutOfWorld”

So there should not be any annihilation gamma rays, as is the case.

Each line of the data shown was written out using G4UserSteppingAction::UserSteppingAction() and filled into a ROOT Ntuple in that order. So I believe the positron said it annihil before it was gone OutOfWorld.

I do not know what process “annihil” refers to in:
“8.32559, -16.418, 104.733, 282.01, 5.00226, Envelope, annihil”

It cannot be the annhilation of the positron with an electron because the positron continues to exist after the process occurs. It is at far too low an energy for any exotic processes to be happening. So I am puzzled.

So am I

Here’s another example in which two 511-keV gammas are emitted and which I think quite normal:

evtID = 30007
  (track,parent,particle) = (1,0,22), 3steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,        -150,      1980,         0,    Envelope,      initStep
               0,           0,           0,      1980,         0,    Envelope,Transportation
               0,           0,     47.1229,         0,         0,      Shape2,          conv
  (track,parent,particle) = (4,1,1000320730), 2steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,     47.1229,0.00387778,         0,      Shape2,      initStep
     1.21651e-07,-7.11514e-07,     47.1229,         0,0.00387778,      Shape2,       ionIoni
  (track,parent,particle) = (3,1,-11), 4steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,     47.1229,   284.514,         0,      Shape2,      initStep
      0.00792904, -0.00046371,     47.1615,   71.9131,   212.601,      Shape2,         eIoni
      0.00781958,  0.00359035,     47.1631,         0,   71.9131,      Shape2,         eIoni
      0.00781958,  0.00359035,     47.1631,         0,         0,      Shape2,       annihil
  (track,parent,particle) = (6,3,22), 5steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
      0.00781958,  0.00359035,     47.1631,   510.999,         0,      Shape2,      initStep
         2.04043,     1.64927,     47.8839,   437.072,   0.03819,      Shape2,         compt
         6.69162,     7.47388,     55.7186,   257.412,   0.03819,      Shape2,         compt
         21.8675,     2.36956,     52.0903,   229.715,   0.03819,      Shape2,         compt
         23.6451,   -0.222614,     52.4412,         0,     1.413,      Shape2,          phot
  (track,parent,particle) = (10,6,11), 3steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
         23.6451,   -0.222614,     52.4412,   228.302,         0,      Shape2,      initStep
         23.6543,   -0.240996,     52.4224,   18.4572,   209.845,      Shape2,         eIoni
         23.6543,   -0.240712,     52.4227,         0,   18.4572,      Shape2,         eIoni
  (track,parent,particle) = (9,6,11), 2steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
         21.8675,     2.36956,     52.0903,   27.6593,         0,      Shape2,      initStep
         21.8681,     2.37003,     52.0899,         0,   27.6593,      Shape2,         eIoni
  (track,parent,particle) = (8,6,11), 2steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
         6.69162,     7.47388,     55.7186,   179.621,         0,      Shape2,      initStep
         6.68905,     7.48609,     55.7329,         0,   179.621,      Shape2,         eIoni
  (track,parent,particle) = (7,6,11), 2steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
         2.04043,     1.64927,     47.8839,   73.8891,         0,      Shape2,      initStep
          2.0437,     1.65012,     47.8811,         0,   73.8891,      Shape2,         eIoni
  (track,parent,particle) = (5,3,22), 4steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
      0.00781958,  0.00359035,     47.1631,   510.999,         0,      Shape2,      initStep
        -8.24535,    -6.67852,     44.2363,   334.372,    1.2554,      Shape2,         compt
        -18.8727,    -7.67742,     53.0292,    164.51,   0.12938,      Shape2,         compt
        -20.3195,     -7.8774,     45.1189,         0,    11.107,      Shape2,          phot
  (track,parent,particle) = (13,5,11), 2steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
        -20.3195,     -7.8774,     45.1189,   153.403,         0,      Shape2,      initStep
        -20.3262,    -7.88866,     45.1122,         0,   153.403,      Shape2,         eIoni
  (track,parent,particle) = (12,5,11), 3steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
        -18.8727,    -7.67742,     53.0292,   169.733,         0,      Shape2,      initStep
        -18.8818,    -7.67827,     53.0439,   35.5967,   134.137,      Shape2,         eIoni
        -18.8828,    -7.67812,     53.0448,         0,   35.5967,      Shape2,         eIoni
  (track,parent,particle) = (11,5,11), 3steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
        -8.24535,    -6.67852,     44.2363,   175.371,         0,      Shape2,      initStep
        -8.25184,    -6.68912,      44.223,   4.47536,   170.896,      Shape2,         eIoni
        -8.25186,    -6.68916,      44.223,         0,   4.47536,      Shape2,         eIoni
  (track,parent,particle) = (2,1,11), 5steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,     47.1229,   673.484,         0,      Shape2,      initStep
      -0.0263752,    0.066541,     47.2337,   453.158,   220.326,      Shape2,         eIoni
      -0.0704182,   0.0367982,      47.291,   207.836,   245.322,      Shape2,         eIoni
       -0.092582,   0.0392047,     47.2998,    48.973,   158.863,      Shape2,         eIoni
      -0.0932645,      0.0375,      47.301,         0,    48.973,      Shape2,         eIoni

The positron in track 3 (3,1,-11) annihilated at “0.00781958, 0.00359035, 47.1631” in “Shape2” with its kinetic energy=0 and ceased to exist:

      0.00781958,  0.00359035,     47.1631,         0,         0,      Shape2,       annihil

Then, at that point, two tracks of 510.999-keV gammas were created:

  (track,parent,particle) = (6,3,22), 5steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
      0.00781958,  0.00359035,     47.1631,   510.999,         0,      Shape2,      initStep

and

  (track,parent,particle) = (5,3,22), 4steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
      0.00781958,  0.00359035,     47.1631,   510.999,         0,      Shape2,      initStep

A total of 3,068,108 e⁺ annihilation events occured for the 10⁸ cases of 1980-keV primary gamma, and all these positrons annihilated and emitted two gammas (some were boosted because the positron had finite kinetic energy just before annihilating).

Only 7 positrons did not annihilate (as shown in my first post) although they said they annihil.
Are these real or physical events? or caused by some bugs?

Yours,
Kazuyoshi

John_McFee · April 14, 2023, 5:16am

Good question and I do not know whether or not it is a bug.

It is possible to have an electron-positron annihilation at low energy with only neutrinos emitted. That would explain the lack of gammas but then there should be neutrino tracks, which you do not see.The probability for that process is ~10000 times less likely than annihilation with photon emission, so you would expect ~300 annihilations with neutrino emission (since you say you got ~3000000 annihilations with gamma emission). The other problem, of course, is that the positron seems to continue on and exit the world boundary AFTER it is annihilated. That makes no sense.

Have you tried to duplicate the problem with any other physics lists?

furutaka · April 14, 2023, 6:41am

No, not yet, and I’ll do that now with QBBC as with the original exmapleB1… But geant4-v11.0.2 is almost one year old and at first it has to be recompiled because the OS on which it is used (Fedora x86_64) has been updated/upgraded since the release of v11.0.2.

Kazuyoshi

furutaka · April 14, 2023, 1:05pm

With QBBC, such an event also occurred but only once in 10⁸ primaries, although strictly speaking it’s not run in exactly the same configuration; the previous one, FTFP_BERT version with G4EmLowEPPhysics, was run around July of 2022 and on the previous version of Fedora linux (36), while this time the QBBC version was built with the same version of geant4 as before but recompiled on Fedora-37 (on one of my computer at home )

evtID = 51427561
  (track,parent,particle) = (1,0,22), 3steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,        -150,      1980,         0,    Envelope,      initStep
               0,           0,           0,      1980,         0,    Envelope,Transportation
               0,           0,   0.0203307,         0,         0,      Shape2,          conv
  (track,parent,particle) = (3,1,-11), 12steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,   0.0203307,   873.044,         0,      Shape2,      initStep
       0.0172057, -0.00246648,   0.0622744,   850.641,   22.4026,      Shape2,           msc
       0.0451361,   -0.025171,   0.0867381,   818.321,   32.3203,      Shape2,           msc
        0.067114,  -0.0142757,    0.130091,   790.173,   28.1481,      Shape2,       annihil
        0.131491,   0.0295229,    0.114927,    738.09,   52.0827,      Shape2,           msc
        0.178125,   0.0864511,    0.110663,   701.676,   36.4139,      Shape2,           msc
        0.225396,    0.114309,   0.0953141,   650.855,   50.8212,      Shape2,           msc
        0.238526,    0.125564,    0.040598,   617.976,    32.879,      Shape2,           msc
        0.243811,    0.138146,  0.00683397,   601.803,   16.1732,      Shape2,           msc
        0.240937,    0.142102, 3.21791e-16,   593.785,    8.0179,      Shape2,Transportation
        -30.6691,     64.5022,        -150,   568.554,   25.2305,    Envelope,Transportation
        -40.5761,     89.0953,        -180,   562.284,   6.27037,       World,    OutOfWorld
  (track,parent,particle) = (2,1,11), 3steps
               X,           Y,           Z,      Ekin,      Edep,      volume,       process
               0,           0,   0.0203307,   84.9585,         0,      Shape2,      initStep
     0.000539715,  0.00788838,   0.0175252,   69.9967,   14.9618,      Shape2,           msc
      0.00369103,  0.00589879,   0.0164618,         0,   69.9967,      Shape2,         eIoni

Kazuyoshi

furutaka · April 14, 2023, 1:33pm

Sorry, I made a mistake not to analyze all the simulation data… Such events occurred twice out of 10⁸ primaries with QBBC!

Kazuyoshi

John_McFee · April 15, 2023, 5:52am

There are a couple of more things to check.

Have you inactivated any processes either in code or in macros?
Are you able to run the code on the latest build of Geant4 (11.1.1) and see if the problem persists?

It could be that there is something subtle about the track printout that we are missing, although I doubt that. However, the problem is rare enough (7 out of 3,068,108 annihilations) and weird enough that it might be a bug, so filing a bug report is a good idea.

Does anybody else have any ideas about what we might be observing?

furutaka · April 15, 2023, 7:53am

No, AFAIK…

I just rebuilt the original FTFP_BERT & G4EmLowEPPhysics version w/ v11.1.1,

**************************************************************
 Geant4 version Name: geant4-11-01-patch-01 [MT]   (10-February-2023)
  << in Multi-threaded mode >> 
                       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...

and made it run for 10⁸ primaries.

The result again contains seven non-annihilating annihil events. (There are other 3,068,096 annihil events in which a positron was annihilated and two gammas were emitted.)

Kazuyoshi

John_McFee · April 15, 2023, 3:17pm

I have one last suggestion and then I am tapped out.

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.

furutaka · April 16, 2023, 9:08am

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 · April 16, 2023, 6:27pm

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 · April 17, 2023, 9:13pm

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?

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 · April 17, 2023, 10:45pm

Hi @dswkey, thanks for your interest.

Not yet done; I will give it a try…

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.

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();

No.

Kazuyoshi

John_McFee · April 18, 2023, 12:11am

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 · April 18, 2023, 5:19am

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 · April 18, 2023, 3:20pm

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 · April 18, 2023, 8:13pm

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 · April 18, 2023, 8:29pm

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 · April 18, 2023, 8:32pm

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();
    }
  }