# The implementation of G4RayleighAngularGenerator::SampleDirection doesn't look like the theoretical formula

​When I tried to learn the implementation of rayleight scattering, I was confused that it looks different from the theoretical formula.
The formula in Physics Reference Manual:
$$\Psi(E,\theta) = [1+cos^2\theta]sin\theta \times FF^2(q)$$
and in Technical_5.pdf (in2p3.fr), the FF(q) can be discribed with PP0~PP8:
$$FF(q) = \frac{PP0^2 - PP1 - PP2}{(1+PP3\times q^2)^{PP6}} + \frac{PP1}{(1+PP4\times q^2)^{PP7}} + \frac{PP2}{(1+PP5 \times q^2)^{PP8}}$$
but in G4RayleighAngularGenerator::SampleDirection, I can’t connect formula to the code.

G4ThreeVector&
G4RayleighAngularGenerator::SampleDirection(const G4DynamicParticle* dp,
G4double, G4int Z,
const G4Material*){
G4double ekin = dp->GetKineticEnergy();
G4double xx = fFactor*ekin*ekin;
//fFactor = 0.5* cm/(h_Planck*c_light) * cm/(h_Planck*c_light); Unit is 1/eV

G4double n0 = PP6[Z] - 1.0;
G4double n1 = PP7[Z] - 1.0;
G4double n2 = PP8[Z] - 1.0;
G4double b0 = PP3[Z];
G4double b1 = PP4[Z];
G4double b2 = PP5[Z];

stactic const G4double numlim = 0.02;
G4double x = 2.*xx*b0;
G4double w0 = (x < numlim) ? n0*x*(1.0 - 0.5*(n0 - 1.0)*x*(1.0 - (n0 - 2.0) *x/3.)):1.0 - G4Exp(-n0*G4Log(1.0+x));

x = 2.*xx*b1;
G4double w1 = (x < numlim) ? n1*x*(1.0 - 0.5*(n1 - 1.0)*x*(1.0 - (n1 - 2.0)*x/3.)):1.0 - G4Exp(-n1*G4Log(1.0+x));

x = 2.*xx*b2;
G4double w2 = (x < numlim)? n1*x*(1.0 - 0.5*(n1 - 1.0)*x*(1.0 - (n1 - 2.0)*x/3.)):1.0 - G4Exp(-n2*G4Log(1.0 + x));

G4double x0 = w0*PP0[Z]/(b0*n0);
G4double x1 = w1*PP1[Z]/(b1*n1);
G4double x2 = w2*PP2[Z]/(b2*n2);

G4double cost;
do {
G4double w = w0;
G4double n = n0;
G4double b = b0;

x = G4UniformRand()*(x0 + x1 + x2);
if(x > x0){
x -= x0;
if(x <= x1){
w = w1;
n = n1;
b = b1;
} else{
w = w2;
n = n2;
b = b2;
}
n = 1.0/n;

G4double y = w*G4UniformRand();
if(y < numlim){
x = y*n*( 1. + 0.5(n + 1.)*y*(1. - (n + 2.)*y/3.));
}else{
x = G4Exp(-n*G4Log(1. - y)) - 1.0;
}
cost = 1.0 - x/(b*xx);
} while ( 2*G4UniformRand() > 1.0 + cost*cost || cost < -1.0);

G4double phi = twopi*G4UniformRand();
G4double simt = sqrt((1. - cost)*(1.0 + cost));
fLocalDirection.set(sint*cos(phi),sint*sin(phi),cost);
fLocalDirection.rotateUz(dp->GetMomentumDirection());
}
}

• what do numlim,w0,w1,w2,n0,n1,n2,b0,b1,b2,x0, x1,x2 mean.
• why do w0,w1,w2 have two representations depend on x < numlim
• I can’t find what is represent FF(q) in this code.
• why was the do {} while circulation used here.