In the method G4NeutronCapture::ApplyYourself,
a G4Element on which capture occurs is chosen based on it's capture cross
section.
I believe the cross section here should be weighted by the
number density (per unit volume) of the given element.
Consider for example H2O. The current implementation
will choose either Hydrogen or Oxygen based only on their relative
capture cross sections, where the cross section on hydrogen should
have twice the weight.

The original code reads as follows:

for (i=0; i<n; i++)
    {
      index = theMaterial->GetElement(i)->GetIndex();
      xSec[i] = theCapture[index].GetXsec(aTrack.GetKineticEnergy());
      sum+=xSec[i];
    }

I have modified it to read:


 G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
 G4double rWeight;
 for (i=0; i<n; i++)
    {
      index = theMaterial->GetElement(i)->GetIndex();
      rWeight = NumAtomsPerVolume[i];
      xSec[i] = theCapture[index].GetXsec(aTrack.GetKineticEnergy());
      xSec[i] *= rWeight;
      sum+=xSec[i];
    }

The same applies to the ApplyYourself methods of G4NeutronHPInelastic,
G4NeutronHPElastic, and G4NeutronHPFission.

--Mark.