F19 metastable state doesn't have decay time

Physics Processes, Models and Cross Sections
1

Hi all,

I’m simulating simulating neutron inelastic scattering on Fluor 19, using 2.5MeV neutrons with the QGSP_BIC_HP physics list, in Geant4 version 10.5.
This results in a gamma line at 197KeV, amongst others, which I expect.

However, from literature I found that the excited state in Fluor that produces this 197KeV gamma is “meta-stable”, with a decay time in the order of 100-200ns.
When I look at the timing information of my 197KeV gammas, this decay time seems to be instantaneous, with every gamma being created within 3ns.
When looking at the creation process for these gammas, it’s neutronInelastic, without any reference to an excited state.

Is there a way to implement this decay time in my simulations?

Thanks in advance for your help,
Koen

2

@ribon : metastable state F19[197 keV] do exist in Geant4 tables. It can be observed with rdecay01; its mean life is ~130 ns.
What is the threshold for explicit generation of G4 ions : millisecond ? microsecond ?
And what is the name of the function and/or UI command which allows users to play with this value ?

3

The time threshold for isomer production depends on the Geant4 version and whether Radioactive Decay is used or not. In the case of QGSP_BIC_HP, Radioactive Decay is always included, so that threshold is 1 microsec for Geant4 10.5 and 10.6, while it has been reduced to 1 nsec in 10.7.p02,. 10.7.p03, and 11.0.

There is no way yet to change this value via UI command.
It is possible to change it via C++ interface, as follows, for instance to set it to 1 nsec :

G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
deex->SetMaxLifeTime( 1.0*nsec/std::log(2.0) );

4

@ribon Thank you for proposed solution.
I’ve set the parameters you’ve proposed in the C++ interface.
I set it during initialization, and called it again during the finalizing, and it was set to 1ns in both cases.

I still don’t see the creation of any additional F19 excited states or a delay in the 197KeV gamma creation times though.
Do I need to pass the G4DeexPrecoParemeters to some variable?

5

There is no need for you to pass the G4DeexPrecoParameter anywhere.
So, in principle, it should work. The fact that it does not for you could be due to the fact that the setting is overwritten (internally, by the Geant4 code) and re-set it to 1 microsec.
I would suggest you to try out G4 10.7.p03, where the default is 1 nanosec, so it should work fine out of the box.

6

Few conclusions from Hadr03, with G4 11.0 @ribon, please, correct me if I am wrong.

1- Both ParticleHP, Bertini, INCL (?) are unable to generate excited states ; only ground state.
so, the code suggested by Alberto is inefficient for these packages.

2- Only Binary Cascade can do. Unfortunately, in this particular case, the only channel generated by BIC is neutron(2.5 MeV) + F19 → F20 + gammas.
And F20 has not metastable excited states.

7

HI Alberto, few additional questions :

  • what exactly the following two functions do ?
    G4NuclideTable::SetThresholdOfHalfLife( )
    G4DeexPrecoParameters::SetMaxLifeTime( )

  • Are they independent or must be used together ?

  • In which circumstances they can be in conflict ?

  • when they must be set : before or after initialization ?

8

Hi Michel,
thanks for the questions! Looking at the code with more attention, to try to answer to your questions, I have realized why the solution I have suggested in my first post does not work, and hopefully I have found another solution.

The reason why my first solution does not work is because, when Radioactive Decay is used (as in QGSP_BIC_HP physics list), then the time threshold in G4NuclideTable - 1000.0ns - overwrites the time threshold in G4DeexPrecoParameters. In other words, whatever value we try to set in G4DeexPrecoParameters - e.g. 1.0ns - it is the value in G4NuclideTable which is used!

Regarding the new solution I would like to propose, this consists of setting 1.0*ns the time threshold in G4NuclideTable, which can be done via the following UI command (it should work both before initialization, but also after, i.e. before /run/beamOn):

/particle/nuclideTable/min_halflife 1.0 ns

Regarding your questions, the two time thresholds are independent, when Radioactive Decay is not used.
The time threshold in G4NuclideTable is used to create the set of isomers that can be used in a Geant4 application. As I said, the default is 1000 ns, and can be changed via UI command both before and after initialization.
The time threshold in G4DeexPrecoParameters, instead, is related to the isomers that can be created by the nuclear de-excitation. Note that the intranuclear cascade models of Geant4 - BERT, BIC, INCL - cannot create isomers. Only the de-excitation part - which can be invoked by BIC and INCL (but not by BERT, unless the variant BERP is considered) - can create isomers. The default value of this time threshold in G4DeexPrecoParameters has changed with the Geant4 version: in G4 10.4 (and any patch), 10.5 (and any patch), 10.6 and 10.6.p01 it was set to 1000 seconds; in G4 10.6.p02, 10.6.p03, 10.7 and 10.7.p01 it was set to 1000 ns; then, since 10.7.p02 (including therefore 10.7.p03, 11.0 and 11.0.p01) it is 1 ns. The reason why the time threshold in G4DeexPrecoParameters has been introduced (in G4 10.4), and its original high default value (1000 s), was simply to save CPU time for cases, like HEP, in which isomers are not important. In fact, only isomers with a mean lifetime above this threshold are created by de-excitation. Then, some unphysical results reported by NA61/SHINE, which uses FTFP_BERT physics list (so, without Radioactive Decay), has led us to decrease this time threshold, first to 1000 ns, then 1 ns. Note that this threshold can be changed only via C++ interface, and only before initialization. What happens when this time threshold is smaller than the time threshold of G4NuclideTable? It happens that the candidate isomer is not found in the list of existing isomers, and therefore the de-excitation considers the corresponding ground state (and, to conserve energy, it includes in the kinetic energy the excitation energy of the would-be isomer state).
What happens instead when the candidate isomer has an average lifetime which is above the threshold in G4NuclideTable, but below the time threshold in G4DeexPrecoParameters (as it happened when the latter was set to 1000 seconds)? In this case the isomer is not created (in the sense that it does not become a G4Track object, and therefore it is not transported by Geant4), but its decay is forced, producing secondaries at the same spatial position where the candidate isomer has been created, but with a time which is sampled from the exponential distribution. The unphysical results that NA61/SHINE reported was due to the fact that this decay happens in-the-flight, i.e. with a relatively high Lorentz boost, whereas in reality the isomer would loose energy by ionization very quickly, and would decay at rest, so its decay products would not receive any Lorentz boost…

9

Hi Alberto, Michel,

Thank you for helping me understand the problem.
I’ve changed the parameters as you suggested. I did (I only run in the C++ interface):

    G4NuclideTable* NuclTable = G4NuclideTable::GetInstance();
    NuclTable->SetThresholdOfHalfLife(1);
    G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
    deex->SetMaxLifeTime(1/std::log(2.0));

With these parameters I ran a simulation running a beam of 2.5 into a pure block of Fluorine. For each event I output the particleID of all particles in the event. Also, for gammas I output the lifetime, creation process info and time (global, proper and local).
This was the result for the events with 197KeV gammas. The numbers show the particle ID, and the additional information is outputted for the gamma.

image

As you can see, it’s inelastic scattering of neutron → gamma, neutron and Fluor ion.
The global time still shows only a time of 1.7ns, instead of the expected 132+ns time after the decay for the 197KeV gamma.
Do you have any idea if there’s something else that could be happening?

10

I can only repeat what I have said.
The adjustment of time life threshold is applicable to BinaryCascade only.
With your physics list, neutrons 2.5 MeV are managed by ParticleHP package, which is unable to create explicitly excited states whatever their time life.
Here, a macro for example Hadr03.
koen.mac.txt (330 Bytes)

11

@ribon , In summary, the time threshold is 1 nanosecond without radioactiveRecay, and 1 microsecond with radioactiveDecay, right ? Of course without any documentation or warning …
This looks like a bug, or at least an inconstancy in Geant4.

12

First, regarding the issue of still getting prompt gammas from 2.5 MeV neutrons on F19 with QGSP_BIC_HP, the problem, as pointed out by Michel, is not the setting of parameters - which are set correctly to 1 nanosec - but the model which is used. Isomers are created only by de-excitation, which is called by BIC (Binary Cascade model) but not for HP (High-Precision neutron transport). In fact, for neutrons below 20 MeV, it is HP, not BIC, which is used in the case of the QGSP_BIC_HP physics list.

If you really need HP - this should be the case if you need a precise treatment of low-energy neutrons - I would suggest to try out the following, but I don’t think, unfortunately, that it will solve the problem: set the environmental variable:

export G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION=1

(for Geant4 >= 11, this environmental variable has been replaced by the UI command:
/process/had/particle_hp/use_photo_evaporation true
).
If this does not work (as I expect), then there is no solution with the current versions of Geant4, but I will discuss with my colleagues whether it is possible to find a solution for future versions of Geant4.

If instead you do not need a precise treatment of low-energy neutrons, I would suggest to try out QGSP_BIC (note, no HP), with RadioactiveDecay on top of it, e.g.

physList->RegisterPhysics( new G4RadioactiveDecayPhysics );

Then you need to set to 1 nanosecond the isomer time threshold (as you have done in your latest attempt).
In this case, de-excitation is called through BIC, so you should be able to see the gammas with a proper time.

Regarding Michel’s question the answer is no: regardless of Radioactive Decay, the default isomers that are created in Geant4 have a half-life of 1 microsec.
The difference - when taking the default - between Radioactive on/off is what happens to the isomers with shorter mean lifetime than ~1.44 microsec (i.e. 1 microsec/ln(2)).
In the case Radioactive Decay is ON, then their gamma de-excitation are prompt;
in the case Radioactive Decay is OFF, then their gamma de-excitation is also prompt if their mean lifetime is below 1 nanosec, while if it is in between\ 1 nanosec and ~1.44 microsec, then there is no gamma emission at all, and the isomer is replaced by the corresponding group state, assigning
by hand the excitation energy of the isomer as kinetic energy of the ground state.
And yes, unfortunately due to several reasons (lazy developers who do not write documentation, different people putting their hands in the code, frequent changes to fix problems, etc.), a clear documentation of the treatment of isomers in Geant4 is still missing…

13

HI Alberto,
Here, printout de /process/had/verbose, without and with G4RadioactiveDecayPhysics registered.

withoutRDM.out.txt (1.5 KB)
withRDM.out.txt (2.4 KB)

14

Here, an output of example Hadr03.
In PhysicsList.cc (attached) I have added the initialization of G4NuclideTable time threshold, and selected INCLXX physics constructor.
As you can see, 4262 states F19[197 keV] are produced, and can be passed to RadioactiveDecay process for deexcitation.
This is not done here. For a ’ complete ’ simulation, have a look to Hadr06 or Hadr07.

@ribon, the command /particle/nuclideTable/min_halflife seems to do nothing …

koen.mac.txt (258 Bytes)
koen.out.txt (5.5 KB)

PhysicsList.cc.txt (5.3 KB)

15

Hi Michel,
it seems to me that, from your log file:
G4NuclideTable::GetInstance()->SetThresholdOfHalfLife( 1.0*picosecond );
has worked correctly:
Time limit for long lived isomeres (ns) 0.0014427
Max life time 1.4427 ps
Am I missing something?

16

Yes, it works ! and F19 [197 keV] is created …
It is the corresponding UI command which makes problem …