Strange results of G4FermiBreakUp in v10.4

Dear Colleagues,

my apologies for a long posting, but I have to explain the whole story. In 2008-2009 I tested multifragmentation models of G4: G4StatMF, G4FermiBreakUp along with G4ExcitationHandler. At that time I proposed several fixes because in G4 v9.1 there were obvious deviations of the results of G4StatMF and G4FermiBreakUp from those of famous Fortran prototypes written by Alexander Botvina. I have to stress that these two implementations were based on exactly same physics. In v9.2 after all the corrections implemented we managed to obtain same results from Fortran and C++ versions of G4StatMF and G4FermiBreakUp. I would like to post several plots to illustrate this, but as a new user of this forum I’m allowed only one figure per post. As can be seen from these plots, the higher excitation energy, the wider list of open decay channels, and, hence, the larger fragment multiplicity.

However, in v10.4 we found recently a very strange results of G4FermiBreakUpVI:

Decay modes of C-12, N-12, C-13, N-13 do not depend on their excitation energy. C-12 always decays exactly into three fragments, while C-12 and N-13 - into a nucleon and a residue. Even more strange behavor found for N-12. Until the excitation of 6 MeV/nucleon is reached, this nucleus remain intact, and only above this threshold decays into p + C-11 become possible. This can be seen from the following output Decays_of_N12_products_v10_4.txt (10.8 KB)

From a more verbose output Decays_of_N12_verbose_v10_4.txt (6.6 KB) it becomes clear that the list of possible decay channels may contain, for example, 8 entries but they are all represented by the p + C-11 configuration.

We also tried v10.3, and even v9.6. The results of these two versions are identical, somehow more realistic, but still very different from 9.2 and the Fortran prototype of the Fermi break-up model.

We would appreciate very much any comments, suggestions from developers to fix G4FermiBreakUpVI. We believe that the reliability of this model is crucial for proper modeling of carbon-ion therapy, where a proper description of beam fragmentation is important.

Kind regards,
Igor Pshenichnov
Alexander Svetlichnyi

Dear Colleagues,
in addition to our previous post let us show a figure made with G4 v. 9.2.


It clearly demonstrates an excellent agreement between the results of a Fortran prototype of Fermi break-up model and its G4 implementation at that time in v. 9.2.
Kind regards,
Igor Pshenichnov

We can also demonstrate that some problems appeared already in v. 9.6. Please, see the next plot.

Hello Igor,

thank you for this post, in time between 9.1 and 10.6 there were several complete upgrade of the de-excitation module. It may be too long to explain why we started these migrations.

For today, we resolve all technical problems and what is left is physics, in particular, we know that probability of reaction channels may be improved/tuned.

VI

Hello Igor,

can you, please, explain what value your plot and how you get it, namely initial state and final state.

VI

Dear Vladimir,
thank you for your comments. To produce these plots we take a specific nucleus, e.g., C-12 or N-13, and simulate its decays at a given excitation energy. Then we score the total numbers of nucleons and nuclear fragments in each event and plot the average versus the excitation energy.

I also can show plots of another kind.


Here we take a specific nucleus, C-12, simulate its decays with v9.2. with scoring the numbers of certain elements as decay products. Again, we finally plot the average numbers of certain elements vs the initial excitation energy. A peak of decays C-12 --> 3 alphas is clearly seen as expected.
The results of v9.6 apparently differ from those of Fortran Fermi break-up.

Same with v10.3:

As I explained above, it is meaningless to make same kind of plots with v10.4 as in all decays of N-12, C-13, N-13 we have nothing, but a nucleon and a residue, or 3 alphas in decays of C-12. It looks like that the list of possible decay channels is not filled properly, just one channel is open at any excitation.
Kind regards,
Igor