Dear Experts,
During the simulation of 50 GeV pi- and kaon- incident on G4_H material in Hadr00, Hadr03 examples with the G4ChargeExchangePhysics enabled the following exceptions are occurring:
G4WT0 > G4IonTable::GetIon() : illegal atomic number/mass Z =0 A = 2 E = 0
G4WT0 >
-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : had012
issued by : G4HadronicProcess:CheckResult()
Warning: Bad energy non-conservation detected, will re-sample the interaction
Process / Model: pi-ChargeEx / ChargeExchange
Primary: pi- (-211), E= 50139.6, target nucleus (1, 2)
E(initial - final) = 50001 MeV.
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------
that finally leads to this exception error message:
-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : had006
issued by : G4HadronicProcess::PostStepDoIt
Call for ChargeExchange
Target element H Z= 1 A= 2
Unrecoverable error in the method ApplyYourself of pi-ChargeEx
TrackID= 1 ParentID= 0 pi-
Ekin(GeV)= 50; direction= (1,0,0)
Position(mm)= (-5e+09,-3.54574e+09,1.42566e+09); material G4_H
PhysicalVolume <G4_H>
ApplyYourself does not completed after 100 attempts
*** Fatal Exception *** core dump ***
G4WT0 > G4Track (0x7fa88c29f120) - track ID = 1, parent ID = 0
G4WT0 > Particle type : pi- - creator process : not available
G4WT0 > Kinetic energy : 50 GeV - Momentum direction : (1,0,0)
G4WT0 > Step length : 4336 km - total energy deposit : 0 eV
G4WT0 > Pre-step point : (-5e+09,-3.54574e+09,1.42566e+09) - Physical volume : G4_H (G4_H)
G4WT0 > - defined by : not available
G4WT0 > Post-step point : (-5e+09,-3.54574e+09,1.42566e+09) - Physical volume : G4_H (G4_H)
G4WT0 > - defined by : pi-ChargeEx - step status : 4
G4WT0 > *** Note: Step information might not be properly updated.
G4WT0 >
-------- EEEE -------- G4Exception-END --------- EEEE -------
It works fine if I set the isotope material of H1, or with the heavier nuclei. It seems that the charge exchange physics library does not handle charge exchange processes of pions & kaons on deuterons/tritons. There is no separate definition for the recoil nucleus in case of (Z==0 && A==2) and (Z==0 && A==3) in G4ChargeExchange.cc, and in common case the ApplyYourself(…) function cannot find the recoil nucleus in Geant4 particle table/ion table, leading to energy non-conservation exception.
I have added charge exchange physics as follows in Hadr03.cc:
PhysicsList* phys = new PhysicsList;
G4ChargeExchangePhysics* ceph = new G4ChargeExchangePhysics();
phys->RegisterPhysics(ceph);
runManager->SetUserInitialization(phys);
and this is my working script:
# Macro file for "Hadr03.cc"
/control/verbose 2
/run/verbose 2
#/testhadr/det/setIsotopeMat H1 1 1 8.89e-5 g/cm3
#/testhadr/det/setIsotopeMat Fe56 26 56 7.87 g/cm3
/testhadr/det/setMat G4_H
/run/numberOfThreads 1
/testhadr/det/setSize 10000000 m
/run/initialize
/gun/particle pi-
/gun/energy 50. GeV
/process/list
/process/inactivate hadElastic
/process/inactivate pi-Inelastic
/process/inactivate Transportation
/process/inactivate kaon-Inelastic
/run/printProgress 100000
/run/beamOn 10000000
The other issue I would like to solve is that for the 50 GeV pi- charge exchange reaction on nuclei other than protons I cannot get any other reaction products but pi0. For pi- on protons I get significant amount of pi0, eta, eta’, omega and f2 with the total charge exchange cross section about 9 ub which is quite close to one that can be found in literature, e.g. arXiv:2312.01703v1 [hep-ph], Table II.
For pi- on Al, Fe and Pb I get not only one channel of reaction (with pi0 as a product), but also the surprisingly low cross section (around 0.1-0.3 ub). This behavior reproduces in Hadr00 and Hadr03 with the different values of charge exchange cross section factors (1, 100, 1e6).
For kaons- on protons I get 4.2 ub and on Al I get 23 ub, which seems to be adequate.
Here is the output of Hadr03 for pi- on Fe56:
The run is 1000000 pi- of 50 GeV through 10000 km of Fe56 (density: 7.87 g/cm3 )
Process calls frequency:
pi-ChargeEx= 1000000
MeanFreePath: 701.64 km +- 701.79 km massic: 5.5219e+05 kg/cm2
CrossSection: 1.4252e-08 cm^-1 massic: 0.0001811 um2/mg
crossSection per atom: 168.21 nbarn
Verification: crossSections from G4HadronicProcessStore
pi-ChargeEx = 629.83 um2/mg 585 mbarn
total = 629.83 um2/mg 585 mbarn
List of nuclear reactions:
pi- + Fe56 --> Mn56 + pi0: 834151 Q = 390.05 keV
pi- + Fe56 --> pi- + : 165849 Q = 0 eV
List of generated particles:
Mn56: 834151 Emean = 168.76 eV ( 168.59 eV --> 178.32 eV )
pi0: 834151 Emean = 50 GeV ( 50 GeV --> 50 GeV)
Here is the output of Hadr03 for pi- on protons:
The run is 1000000 pi- of 50 GeV through 10000 km of H1 (density: 0.0889 mg/cm3)
Process calls frequency:
pi-ChargeEx= 1000000
MeanFreePath: 14817 km +- 14819 km massic: 131.72 kg/cm2
CrossSection: 6.7489e-10 cm^-1 massic: 0.75916 um2/mg
crossSection per atom: 12.705 mubarn
Verification: crossSections from G4HadronicProcessStore
pi-ChargeEx = 1.2483 mm2/g 20.89 mbarn
total = 1.2483 mm2/g 20.89 mbarn
List of nuclear reactions:
pi- + H1 --> eta + neutron: 78036 Q = -409.59 MeV
pi- + H1 --> eta_prime + neutron: 2286 Q = -819.5 MeV
pi- + H1 --> neutron + f2(1270): 45166 Q = -1.1367 GeV
pi- + H1 --> neutron + omega: 13354 Q = -644.37 MeV
pi- + H1 --> pi- + : 258612 Q = 0 eV
pi- + H1 --> pi0 + neutron: 602545 Q = 3.3002 MeV
pi- + XXXX --> pi- + : 1 Q = 0 eV
Here I would also like to note that firstly, cross sections from G4HadronicProcessStore for charge exchange reactions in fact return the value for inelastic processes of pi-, because if I will not inactivate pi-Inelastic, there will be new line for pi-Inelastic cross section of 20.89 mbarn.
As for cross sections of charge exchange reaction, ones calculated from MeanFreePath are close to ones in literature, but there is a discrepancy. As you can see, in a list of nuclear reactions the significant part takes pi- + H1 → pi- + : which does not seem to be any kind of real interaction. For the calculation of cross section, I have used the proportion in which I find the cross section for the real reactions (with pi0, eta, eta’, omega, f2) as a part of calculated one (12.705 ub) like sigma_cexch = 12.705 * N_real/N_entries, and it gives about 9.4 ub (this value is within the error of one given in literature).
These “blank” reactions are taken into account when the cross section is calculated in Hadr03 example, although it is possible to exclude them in SteppingAction. With this, I would like to note that these reactions with the “blank” products list are possibly treated as the physical ones. These “blank” reactions also appear for pi- on heavier nuclei, and for kaon-.
Is there any bug? Thank you in advance.
_Geant4 Version: geant4-v11.2.2 multithreaded
_Operating System: Linux Mint 22 (on Win10 WSL)
_Compiler/Version: gcc 13.2.0
_CMake Version: 3.28.3