Problem 1815

Created attachment 375 [details]
Patched files, output files, and equation derivation

I am trying to simulate 2 MeV protons incident on lithium as a neutron source.
The new ParticleHP package seems to be appropriate at these low energies.
To get a reasonable number of interactions I used the hadronic example Hadr03 with the attached macro, protons_Li.mac.
(PhyisicList.cc was modified to use QGSP_BIC_AllHP physiclist.)
It produced the attached output file "output_proton_2_MeV".
The Q value for the p[Li7, Be7]n is -1.6443 MeV so only 0.3557 MeV is available to be share by the neutron and Be ion at 2 MeV.
As you can see more than 0.8 MeV is shared.

I thought I had found the problem in G4ParticleHpInelasticCompFs.cc at line 280 and made the following modification:

// Watts fix G4double residualZ = theBaseZ - aDefinition->GetPDGCharge();
    G4double residualZ = theBaseZ - aDefinition->GetPDGCharge()+theProjectile->GetPDGCharge();

Without the change the residual is Li7 rather than Be7.
This did not fix the problem however.

The kinetic energy, kinE, is defined near the end of G4ParticleDiscreteTwoBody.cc:

   G4double A1     =  GetTarget()->GetMass()/GetProjectileRP()->GetMass(); 
   G4double A1prim =  result->GetMass()/GetProjectileRP()->GetMass();
   G4double E1     =  (A1+1)*(A1+1)/A1/A1*anEnergy; 
   G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*E1+(1+A1)*GetQValue())

These equations look similar to the equation in http://t2.lanl.gov/nis/endf/intro19.html 
but not exactly. (Attached is a pdf file, "Two_Body_Derivation.pdf" developing the derivation.)

I think that anEnergy is the primary proton energy in the LAB system, and
E1 is should be equal to anEnergy. 
kinE is the kinetic energy of the neutron (or the Be7) in the Center of Mass system
then the equations should be:

   G4double A1     =  GetTarget()->GetMass()/GetProjectileRP()->GetMass(); 
   G4double A1prim =  result->GetMass()/GetProjectileRP()->GetMass();
   G4double E1     =  anEnergy; 
   G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*E1+(1+A1)*GetQValue())

Attached is "output_proton_2_MeV_Fix", the output result with the two patches applied, and
patched routines, G4ParticleHpInelasticCompFs.cc and G4ParticleDiscreteTwoBody.cc.
While these patches are for the ParticleHP routines, the second bug is also in the older NeutronHP package.