ParticleHP does not produce any residual isotopes in Hadr03 from inelastic interactions

At the risk of potentially asking something already covered several times ([2597], [2715]), I observed that ParticleHP does not producing any residual isotopes after inelastic interactions.

I ran Hadr03 starting 20000 protons at 50 MeV in G4_lAr using QGSP_BIC_HP (w/o ParticleHP) and QGSP_BIC_AllHP (w/ ParticleHP). I get the following resulting produced particles in the output:

QGSP_BIC_HP:

 List of generated particles:
           Al26:       1  Emean =  4.2912 MeV	( 4.2912 MeV --> 4.2912 MeV)
           Al27:       1  Emean =  8.4148 MeV	( 8.4148 MeV --> 8.4148 MeV)
           Al28:       4  Emean =  5.1399 MeV	( 1.1916 MeV --> 9.21 MeV)
           Al29:       5  Emean =  3.3606 MeV	( 65.266 keV --> 7.4376 MeV)
           Ar35:       1  Emean =   55.65 keV	( 55.65 keV --> 55.65 keV)
           Ar36:      65  Emean =  554.35 keV	( 8.3926 keV --> 2.9962 MeV)
           Ar37:     317  Emean =  1.5447 MeV	( 53.216 keV --> 5.7144 MeV)
           Ar38:    1582  Emean =  1.3895 MeV	( 2.8741 keV --> 5.9685 MeV)
           Ar39:     956  Emean =  1.1432 MeV	( 9.1483 keV --> 4.5937 MeV)
           Ar40:   12722  Emean =     177 keV	( 7.5356 eV  --> 4.5244 MeV)
            B10:       3  Emean =  8.2226 MeV	( 4.9832 MeV --> 12.54 MeV)

            ... (58 nuclies and their excited states in total)

           Si31:      12  Emean =  2.6105 MeV	( 186.83 keV --> 7.1757 MeV)
           Si32:       2  Emean =  1.5243 MeV	( 511.58 keV --> 2.537 MeV)
          alpha:    1869  Emean =  8.0238 MeV	( 1.4127 MeV --> 22.428 MeV)
       deuteron:     472  Emean =  11.518 MeV	( 1.6765 MeV --> 40.168 MeV)
             e-:      18  Emean =  181.24 keV	( 10.439 keV --> 2.2209 MeV)
          gamma:   13437  Emean =   2.028 MeV	( 546.18 eV  --> 23.059 MeV)
        neutron:    9028  Emean =  7.0264 MeV	( 7.2636 keV --> 46.376 MeV)
         proton:    5131  Emean =  16.002 MeV	( 676.22 eV  --> 49.744 MeV)
         triton:     107  Emean =  9.0722 MeV	( 1.793 MeV --> 34.04 MeV)

QGSP_BIC_AllHP:

 List of generated particles:
           Ar36:      46  Emean =  99.409 keV	( 9.3016 keV --> 259.6 keV)     <-- from hadElastic
           Ar38:       4  Emean =  245.97 keV	( 14.822 keV --> 713.57 keV)    <-- from hadElastic
           Ar40:   12212  Emean =  154.45 keV	( 7.5356 eV  --> 4.5244 MeV)    <-- from hadElastic
          gamma:    6892  Emean =  3.2507 MeV	( 17.923 keV --> 48.469 MeV)
        neutron:   13778  Emean =  5.7083 MeV	( 4.843 keV --> 47.166 MeV)
         proton:    6926  Emean =  16.889 MeV	( 326.41 keV --> 49.684 MeV)

For the latter, looking also at the list of nuclear reactions, there are some references to reactions that do no specify the products correctly:

...
        proton + Ar40 --> proton + :     820   Q =       0 eV 
        proton + Ar40 --> proton + Ar40:   12212   Q = 0.0015564 meV
...

so some interactions seem to occur, but the isotope tracks might not be generated?

I also repeated the test using the other charged particles covered by ParticleHP again, without any residual isotopes.

The macro is attached below. By uncommenting either the G4HadronPhysicsQGSP_BIC_HP or G4HadronPhysicsQGSP_BIC_AllHP as well as the G4IonPhysicsXS or G4IonPhysicsPHP in src/PhysicsList.cc, one should reproduce the above results. I used G4TENDL1.4 shipped with the G4 version and the output points at the correct folder they are stored in.

I guess my questions are:

  1. Am I using ParticleHP correctly like this?
  2. Is this the expected result?

Regarding 2), I understand from [2715] that ParticleHP should only reproduce the expected residual particles statistically, but in this case, no residual isotopes are produced. I tried to read the source code to understand it better, but found inconsistencies with this statement as currently implemented:

I understand that PHP_USE_POISSON in G4ParticleHPManager [link] allows toggling between theMultiplicityMethod methods of the products (see [link]). However, PHP_USE_POISSON seems to be initialized as false, and the only branch that adjusts it can set it to false as well, making it impossible to set to true.

At a later point in G4ParticleHPProduct [link], the multiplicity for a product is defined. When G4HPMultiPoisson is selected, it samples from a Poissonian. When instead G4HPMultiBetweenInts is selected, which I understand is always the case, the multiplicity is taken as the rounded value of the mean multiplicity of the yield of a product using G4lrint. Since the yields for the residual isotopes in MF=6 files are all <<1, I assume that, therefore, the multiplicities of these are virtually always set to 0, resulting in no isotope production.

Therefore, question

  1. Is this description correct? Is this the intentional behavior?

I do not want to exclude the fact that I might be using G4TENDL and ParticleHP wrong. I want to apologize in advance for any obvious mistakes in the above. Thank you very much!

macro_lar_proton.txt (136 Bytes)

Geant4 Version: 11.4.0
Operating System: Ubuntu 24.04.3 LTS
Compiler/Version: 13.3.0
CMake Version: 3.28.3


Your problem maybe due to the cahnge of the default multiplicity method in Geant4. change the line to true:

source/processes/hadronic/models/particle_hp/include/G4ParticleHPManager.hh:254: G4bool PHP_USE_POISSON{false};

Or the enviromental variable G4PHP_MULTIPLICITY_METHOD to Poisson

The G4PHP_MULTIPLICITY_METHOD is apparently ignored. I set it to true, but the ParticleHP report still mentioned the setting to be off.

...
 Use DBRC                        0
 PHP use Poisson                 0
 PHP check                       1
...

After making the suggestion to the initialization of PHP_USE_POISSON and recompiling, it is indeed activated

...
 Use DBRC                        0
 PHP use Poisson                 1
 PHP check                       1
...

and I also get promising results:


 List of generated particles:
           Al28:       1  Emean =  1.4682 MeV   ( 1.4682 MeV --> 1.4682 MeV)
           Al29:       1  Emean =  1.4322 MeV   ( 1.4322 MeV --> 1.4322 MeV)
           Ar36:     100  Emean =  615.79 keV   ( 9.3016 keV --> 1.7113 MeV)
           Ar37:     701  Emean =  1.1722 MeV   ( 251.32 keV --> 2.3926 MeV)
           Ar38:    1846  Emean =  981.39 keV   ( 1.3704 keV --> 3.2122 MeV)
           Ar39:     990  Emean =  796.45 keV   ( 411.22 eV  --> 4.6964 MeV)
           Ar40:   12852  Emean =  153.15 keV   ( 7.5356 eV  --> 4.5244 MeV)
           Cl34:      15  Emean =  1.5072 MeV   ( 288.69 keV --> 3.3138 MeV)
           Cl35:     111  Emean =  1.6779 MeV   ( 6.5368 keV --> 4.0322 MeV)
           Cl36:     191  Emean =  1.2389 MeV   ( 4.259 keV --> 3.5596 MeV)
           Cl37:     498  Emean =  1.1888 MeV   ( 58.882 keV --> 2.4893 MeV)
           Cl38:     197  Emean =  1.0752 MeV   ( 66.955 keV --> 2.9793 MeV)
           Cl39:      88  Emean =     992 keV   ( 81.248 keV --> 2.8704 MeV)
            He3:     140  Emean =  18.974 MeV   ( 2.8485 MeV --> 33.283 MeV)
            K38:      31  Emean =  1.1755 MeV   ( 610.17 keV --> 2.0875 MeV)
            K39:     188  Emean =  811.33 keV   ( 290.82 keV --> 1.6894 MeV)
            K40:      75  Emean =  764.46 keV   ( 21.259 keV --> 2.5914 MeV)
            K41:       2  Emean =  1.5096 MeV   ( 1.4071 MeV --> 1.6122 MeV)
            P31:      23  Emean =  2.0498 MeV   ( 139.76 keV --> 3.8884 MeV)
            P32:      19  Emean =  2.6158 MeV   ( 578.4 keV --> 4.2898 MeV)
            P33:       7  Emean =  1.6046 MeV   ( 420.37 keV --> 3.1347 MeV)
            P34:      12  Emean =  1.5244 MeV   ( 1.027 MeV --> 2.1956 MeV)
            P35:       2  Emean =  2.8172 MeV   ( 2.5122 MeV --> 3.1221 MeV)
            S32:       2  Emean =  2.2872 MeV   ( 1.6938 MeV --> 2.8807 MeV)
            S33:      24  Emean =  1.4437 MeV   ( 47.539 keV --> 3.1988 MeV)
            S34:     412  Emean =  1.7176 MeV   ( 18.882 keV --> 4.4174 MeV)
            S35:     130  Emean =  1.6353 MeV   ( 34.124 keV --> 3.6414 MeV)
            S36:      37  Emean =  1.2011 MeV   ( 7.0729 keV --> 2.3247 MeV)
            S38:       3  Emean =  1.0949 MeV   ( 655.8 keV --> 1.7058 MeV)
           Si30:       4  Emean =  1.5811 MeV   ( 1.2503 MeV --> 1.9644 MeV)
           Si31:       2  Emean =  2.2495 MeV   ( 2.0487 MeV --> 2.4502 MeV)
           Si32:       2  Emean =  1.6027 MeV   ( 1.4299 MeV --> 1.7755 MeV)
          alpha:     764  Emean =  8.9939 MeV   ( 1.655 MeV --> 41.717 MeV)
       deuteron:     677  Emean =  21.003 MeV   ( 2.0967 MeV --> 40.616 MeV)
          gamma:    4117  Emean =  3.1425 MeV   ( 1.1435 keV --> 10.784 MeV)
        neutron:    9466  Emean =  5.6679 MeV   ( 7.2519 keV --> 46.254 MeV)
         proton:    5526  Emean =  17.173 MeV   ( 690.19 keV --> 49.591 MeV)
         triton:      73  Emean =  17.134 MeV   ( 1.6808 MeV --> 36.713 MeV)

I’d say it’s because of missing Fission fragment data in G4Neutronhp and G4particlehp.
Which macro are you running? Can you try using qgsp_bert_hp (which gives wrong n,* cross sections, but fission is correct via FissionVI)?