Created attachment 628 [details]
Bundle including patches, print out text and slide showing observation of problem and proposed fix.

Hello,

My group has found a few problems with the EMD process when applying it to ion interactions at the LHC for the purpose of collimation. We've devised some patches for Geant4.10.4.p03, 10.5.p01 and 10.6.p01. We tend to use 10.4.p03 as the later versions are much much slower for our LHC accelerator model.

The aim is to simulate Pb 208 82 with 82 * 6.37TeV (~522TeV total with ~2.5TeV/n) on porous carbon, Cu and W for the collimators. For us, EMD is crucial because the resulting Pb 207, 206, 205 etc are close in charge / mass to get through the machine.

1) Energy conservation. The process doesn't conserve energy to as much as 15%. One of the parent hadronic classes detects this and resamples it. It seems to always fail for EMD of the projectile in the field of the target, but on re-sampling it seems to eventually select EMD of the target in the field of the projectile, which works. This results in EMD not producing Pb 207 (-1 proton) from Pb 208 for example at higher kinetic energies (ie LHC ion energies). Attached is a print out of all the kinematics in one example as well as some graphs of the occurrence of Pb 207 in our program. Proposed fix to the kinematics in G4EMDissociation.cc line 243:
-    pP.setE(pP.e()+Eg);
+    pP.setE(std::sqrt(pP.vect().mag2() + std::pow((mass + Eg), 2.)));

2) The cross-section calculation can result in a bad arithmetic exception (basic C++ one) that causes a program to crash. This happens only for calculating the EMD cross-section for extremely low kinetic energy ions about to come to rest. The cross-section here should be small I believe. We can circumvent this by by setting it to 0 when the kinetic energy is low. In G4EMDissociationCrossSection.cc ~line 127:

+  // 0 cross-section for particles with kinetic energy less than 2MeV to prevent
+  // possible abort signal from bad arithmetic in GetCrossSectionForProjectile
+  if (theDynamicParticle->GetKineticEnergy() < 2*CLHEP::MeV)
+    {return 0.0;}

3) A crash can occur when EMD applies to an H4 (1 proton, 3 neutrons) ion fragment. It proposes to remove 1 proton and then queries the ion table for A=3, Z=0, which then throws an exception. We can protect against this (assuming such an unstable fragment is handled elsewhere) in G4EMDissociationCrossSection.cc ~line 237:
-  if (Z < 6.0)
+  if (Z < 2.0)
+    p = 0;
+  else if (Z < 6.0)

4) EMD can create an unbelievable number of warnings from G4Fragment. In fact the cout slows our simulation and the log file is too big. Attached is a patch where we include an optional argument of G4bool warning=False to the G4Fragment constructor so we can control whether the warning is issued or not. This is so we only affect EMD and nothing else.

Credit to Andrey Abramov for the majority of this.

Hopefully some of these can be included in Geant4 or reviewed and developed into it.  We have a Geant4 model of the whole LHC and apply it regularly for proton loss and background studies and are working to achieve the same for the Pb ions.

Many thanks,
Laurie Nevay & Andrey Abramov