Problem 371 - Errors in G4LowEnergyRayleigh.cc
Summary: Errors in G4LowEnergyRayleigh.cc
Status: CLOSED FIXED
Alias: None
Product: Geant4
Classification: Unclassified
Component: config (show other problems)
Version: 4.0
Hardware: PC Linux
: P1 critical
Assignee: Maria.Grazia.Pia
URL:
Depends on:
Blocks:
 
Reported: 2002-05-09 19:37 CEST by mkippen
Modified: 2002-06-07 12:32 CEST (History)
0 users

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this problem.
Description mkippen 2002-05-09 19:37:15 CEST
I believe there is a problem in the way the Rayleigh scatter angle is
sampled in G4LowEnergyRayleigh::PostStepDoIt.  The relevant section of
code is included below for reference.

(1) The Rayleigh cross section per unit SOLID ANGLE is

   dSigma/dOmega \propto (1+cosTheta^2) * FF^2 ,

where FF is the coherent scattering form factor.  This is what is used in
the current code.  HOWEVER, we are sampling Theta, not Omega!.  The
probability density function for Theta, alone, is

   dSigma/dTheta \propto (1+cosTheta^2) * FF^2 * sin(Theta)

Ignoring the sin(Theta) term results in the incorrect angular distribution,
with too many forward-scattered photons.

(2) Whether one uses the sin(Theta) term, or not, the range of the
rejection function gReject must agree with the range of the random
deviate randomFormFactor.  This is currently not valid, since
gReject is in the range [0 --> 2 * Z^2], but randomFormfactor only goes
from [0--> 1*Z^2].  This results in a major error in the scatter angle
distribution.

Potential Solution to both problems:

G4double sqrRayl = sqrt(27/32) * (1 + cosTheta * cosTheta) * sinTheta;
gReject = sqrRayl * dataFormFactor * dataFormFactor;

This works, but is inefficient for some cases when the FF makes gReject
significantly lower than sqrRayl.  It may be better to samples the
cumulative distribution function rather than differential (see EGS4 manual
for the appropriate algorithm).

------------------------------------------------------------
  do
    {
      G4double thetaHalf = G4UniformRand() * pi / 2.;
      G4double sinThetaHalf = sin(thetaHalf);
      G4double x = sinThetaHalf / (wlPhoton/cm);
      G4double dataFormFactor = formFactorData->FindValue(x,Z-1);
      randomFormFactor = G4UniformRand() * Z * Z;
      G4double theta = thetaHalf*2;
      cosTheta = cos(theta);
      sinTheta = sin(theta);
      G4double sqrRayl = 1 + cosTheta * cosTheta;
      gReject = sqrRayl * dataFormFactor * dataFormFactor;

    } while( gReject < randomFormFactor);
------------------------------------------------------------
Comment 1 pia 2002-06-07 12:32:59 CEST
Thanks for spotting the problem with angular distributions and for proposing a
solution.
A fix, corresponding to the quick solution you proposed, has been implemented
and will be publicly available with the Geant4 4.1 release on 28 June 2002.
Further improvements to the angular distribution of G4LowEnergy Rayleigh are
planned and will be available in the next Geant4 releases.
Best regards,
Maria Grazia Pia