Multiple full-energy de-excitations from single nucleus?

[mkelsey]

Geant4 Version: 10.7.p02
Operating System: CentOS7
Compiler/Version: G++ 8.3
CMake Version: 3.17.3

We run our simulations using the Shielding physics list with some additions (tritium decay, optical physics); that means we’re using NeutronHP.

One of my collaborators reported seeing odd energy non-conservation in some of our simulated nuclear recoil events, and was able to isolate things to a specific problem.

It seems that nuclear de-excitation is producing multiple iterations of the full excitation energy for a single nucleus in a single G4Step. Here’s a specific example for neutron capture on C-12, producing C-13* and de-excitation:

G4WT0 > *********************************************************************************************************
G4WT0 > * G4Track Information:   Particle = neutron,   Track ID = 1,   Parent ID = 0
G4WT0 > *********************************************************************************************************
G4WT0 > 
G4WT0 > Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
G4WT0 >     0    -81.9 -1.18e+03     -221     0.001        0        0         0 Shield/InnerPoly initStep
...
G4WT0 >    35    -65.5 -1.2e+03     -152         0        0    0.666       163 Shield/InnerPoly nCapture
G4WT0 >     :----- List of 2ndaries - #SpawnInStep=  5(Rest= 0,Along= 0,Post= 5), #SpawnTotal= 39 ---------------
G4WT0 >     :     -65.5  -1.2e+03      -152      4.95              gamma
G4WT0 >     :     -65.5  -1.2e+03      -152      4.95              gamma
G4WT0 >     :     -65.5  -1.2e+03      -152      3.68              gamma
G4WT0 >     :     -65.5  -1.2e+03      -152      1.26              gamma
G4WT0 >     :     -65.5  -1.2e+03      -152  0.000321                C13
G4WT0 >     :----------------------------------------------------------------- EndOf2ndaries Info ---------------

Our user wrote:

Looking up the gamma energies from the 12C(n,gamma)13C process in the NNDC database for C-12 (note the level scheme available as pdf download, also attached here), I expect to see from one neutron capture of C12 either:

• exactly 1 gamma with 4945 keV, or
• 1 gamma with 1261 keV and 1 gamma with 3683 keV, or
• the lower probability stuff you see in the level scheme.

What I do not expect at all, is to see from one single neutron capture 2 gammas of 4945 keV, and also not 2 gammas of 4945 keV and additional gammas. This doesn’t make sense to me. Where is this energy coming from? Does Geant4 ignore energy conservation here? Btw, this doesn’t happen every time. Often enough, there is just 1 gamma of 4945 keV as it should be.

There are a couple of forum topics from 2021 which report the same issue:

Question in neutron capture (Gd155(n, gamma)Gd156) reaction

Hello everyone! I’m simulating neutron capture reaction with Gd material. The physics list is QGSP_BIC_AllHP. The data I could find is Gd156* will release several gammas with total energy about 8 MeV, which is correspond to the highest level of Gd156. However, when I do the simulation, here comes: [gd] The total energy of gammas is far more than 8 MeV. The same situation happens with Gd157. It really confuses me. What happened here?

Too many gamma-rays from neutron capture by nickel-58?

Dear all, Now we are checking the energy of gamma-rays generated from neutron capture of nickel using Geant4.10.03.p02. I defined the material using Nistmanager G4Element* elNi = NistManager->FindOrBuildElement(“Ni”); and the information of all steps is saved. When I checked each step, I confirmed that 3 gamma-rays of ~9 MeV were generated from 1 neutron capture by nickel-58. If my understanding is correct, total energy of gamma-rays generated from 1 neutron capture by nickel-58 is ~9 MeV. …

It seems as thought this problem arises from having G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION=0 (the default). Users in the above topics report that setting that flag to 1 eliminates the spurious multiple energies problem.

But that doesn’t seem entirely right to me. Generating multiple instances of a full de-excitation cascade can’t simply be a “feature” of the code. Is anyone sufficiently familiar with the non-photoevaporation side of de-excitation to understand what’s going on under the hood?

[John_McFee]

I am not intimately familiar with the non-photoevaporation de-excitation models. However, my experience has shown me that thermal neutron capture gamma ray results have to be treated with much scepticism and simulation results must be carefully verified against experimental data on an isotope by isotope basis. This opinion seems to be supported by the Geant4 documentation.

In Book For Application Developers, Release 11.2 High-precision neutron interactions (NeutronHP) p.239-240, it says:

It is highly recommended or mandatory to set the following UNIX environment variable and UI commands when running a GEANT4 application:
…
/process/had/particle_hp/do_not_adjust_final_state true
Without this UI command, a GEANT4 model that attempts to satisfy the energy and momentum conservation in some nuclear reactions, by generating artificial gamma rays…

Indeed, I have observed this on many occasions.

Further on it says:

The results concerning the production of gamma-rays have to be questioned always. For example, the information on the number and energies of gamma-rays emitted in the neutron capture process is incomplete for almost all the nuclei and is frequently also uncorrelated. When the information is available, it will be used, but one can obtain results which are quite far from reality on an event by event basis: the total energy of the cascade won’t be correct in many cases and only some specific gamma-rays which are stored in the neutron databases will be emitted. If there isn’t any information concerning these gamma-rays, GEANT4 will use a simple a model instead which is generally missing the relevant spectroscopic information.

It has been my experience that some combinations of the HP particle parameters give good (n,gamma) agreement with NNDC for gamma energies and intensities for some isotopes and poor results for others. A different combination may give reverse results. I have mentioned specifics in some of the threads that you have referenced. In particular, the two main HP particle parameters which affect the results are DoNotAdjustFinalState and UseOnlyPhotoEvaporation.

In the manual they suggest that this variability with parameter choice is a necessary consequence of the incompleteness of the low energy data set, given that the low energy reaction processes are data-driven.

[mkelsey]

Thank you very much for the detailed insight, @John_McFee! I’ll pass this back and try to figure out what parameter settings we need. Even if G4 is doing something like continuum emission, the idea that it generate multiple copies of tabulated gammas has got to be a bug somewhere.

[John_McFee]

the idea that it generate multiple copies of tabulated gammas has got to be a bug somewhere

Yes @mkelsey, I definitely agree with that. If you find parameter settings or any other tricks that work, I’d be grateful if you’d let me know.

Do you think it warrants a bug report or do you think the designers have accepted it as a necessary evil?

[mkelsey]

I would definitely like to file a bug report on it, but as you know the barrier is set high for that I’ll need to construct a simple standalone example (rather than using our CDMS simulation framework), and probably find random seed values to set so that the problem happens on the first event.

As it is, its an “occasional” problem we see, which causes a dribble of weird outliers in our energy distributions.

If we can identify the NeutronHP settings that make it go away, I will definitely follow up here. I’ve linked this post in our internal JIRA ticket for it

[ajbiffl219]

Ok, so after looking into the capture results in carbon, G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION=1 does indeed maintain energy conservation, at the cost of slightly distorting the gamma spectrum. If you want the correct gamma distribution but don’t care about energy conservation, G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION=0 with /process/had/particle_hp/do_not_adjust_final_state true will give the correct gamma yields.

The problem seems to happen for nuclei with discrete gamma multiplicity data. When not using the PhotonEvaporation classes, each “type” of gamma is sampled as an independent Poisson process, and the lists of gammas are concatenated together - so you get the correct yields, but no limits on the energy released in any individual event.