The kinematics is not correctly calculated in G4LElastic.cc:

line 115: G4double p = aParticle->GetTotalMomentum()/GeV;

the momentum is converted into GeV, while all other variables are in MeV,
especially the masses.

Apart from that, the kinematics formula is wrong:

lines 201-213:

   G4double etot = sqrt(m1*m1+p*p)+m2;
   a = etot*etot-p*p*cost*cost;
   b = 2*p*p*(m1*cost*cost-etot);
   c = p*p*p*p*sint*sint;

   G4double de = (-b-sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a);
   G4double e1 = sqrt(p*p+m1*m1)-de;
   G4double p12=e1*e1-m1*m1;
   p1 = sqrt(std::max(0.0,p12))/1000;
   p = p/1000;
   px = p1*sint*sin(phi);
   py = p1*sint*cos(phi);
   pz = p1*cost;

Problem Fix: replacing lines 201-213 with the following code, keeping line 115:

   G4double pcheck=p*1000.; // GeV->MeV
   G4double ev1 = sqrt(m1*m1+pcheck*pcheck);
   G4double etot = ev1+m2;
   a = etot*etot-pcheck*pcheck*cost*cost;
   b = -2.*(m2*m2*ev1+m2*ev1*ev1+m2*m1*m1+ev1*m1*m1);
   c = 2.*m2*ev1*m1*m1+m2*m2*ev1*ev1+m1*m1*m1*m1+m1*m1*pcheck*pcheck*cost*cost;
   G4double e1 = (-b+sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a);
   G4double p12=e1*e1-m1*m1;

   p1 = sqrt(std::max(0.0,p12))/1000.; //MeV->GeV
   px = p1*sint*sin(phi);
   py = p1*sint*cos(phi);
   pz = p1*cost;

Perhaps I am using an old an forgotten file, but I have the impression that it
is still frequently used in the physics lists for the elastic processes, or is
there already a replacement?

Cheers

Markus Hauger