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

Physics Processes, Models and Cross Sections
1

​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.