G4RayleighAngularGenerator Class Reference

#include <G4RayleighAngularGenerator.hh>

Inheritance diagram for G4RayleighAngularGenerator:

G4VEmAngularDistribution

Public Member Functions

 G4RayleighAngularGenerator ()
virtual ~G4RayleighAngularGenerator ()
virtual G4ThreeVectorSampleDirection (const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=0)

Detailed Description

Definition at line 61 of file G4RayleighAngularGenerator.hh.


Constructor & Destructor Documentation

G4RayleighAngularGenerator::G4RayleighAngularGenerator (  ) 

Definition at line 64 of file G4RayleighAngularGenerator.cc.

00065   : G4VEmAngularDistribution("CullenGenerator")
00066 {
00067   G4double x = cm/(h_Planck*c_light);
00068   fFactor = 0.5*x*x;
00069 }

G4RayleighAngularGenerator::~G4RayleighAngularGenerator (  )  [virtual]

Definition at line 73 of file G4RayleighAngularGenerator.cc.

00074 {}


Member Function Documentation

G4ThreeVector & G4RayleighAngularGenerator::SampleDirection ( const G4DynamicParticle dp,
G4double  out_energy,
G4int  Z,
const G4Material mat = 0 
) [virtual]

Implements G4VEmAngularDistribution.

Definition at line 79 of file G4RayleighAngularGenerator.cc.

References G4VEmAngularDistribution::fLocalDirection, G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), and CLHEP::detail::n.

00082 {
00083   G4double ekin = dp->GetKineticEnergy();
00084   G4double xx = fFactor*ekin*ekin;
00085   G4double cost;
00086 
00087   G4double n0 = PP6[Z] - 1.0;
00088   G4double n1 = PP7[Z] - 1.0;
00089   G4double n2 = PP8[Z] - 1.0;
00090   G4double b0 = PP3[Z];
00091   G4double b1 = PP4[Z];
00092   G4double b2 = PP5[Z];
00093   G4double w0 = 0.0;
00094   G4double w1 = 0.0;
00095   G4double w2 = 0.0;
00096 
00097   const G4double numlim = 0.02;
00098   G4double x  = 2*xx*b0;
00099   if(x < numlim)   { w0 = n0*x*(1 - 0.5*(n0 - 1)*x*(1 - (n0 - 2)*x/3.)); }
00100   else             { w0 = (1 - std::pow(1 + x,-n0)); }
00101 
00102   if(PP1[Z] > 0.0) {
00103     x  = 2*xx*b1;
00104     if(x < numlim) { w1 = n1*x*(1 - 0.5*(n1 - 1)*x*(1 - (n1 - 2)*x/3.)); }
00105     else           { w1 = (1 - std::pow(1 + x,-n1)); }
00106   }
00107   if(PP2[Z] > 0.0) {
00108     x  = 2*xx*b2;
00109     if(x < numlim) { w2 = n2*x*(1 - 0.5*(n2 - 1)*x*(1 - (n2 - 2)*x/3.)); }
00110     else           { w2 = (1 - std::pow(1 + x,-n2)); }
00111   }
00112 
00113   G4double x0= w0*PP0[Z]/(b0*n0);
00114   G4double x1= w1*PP1[Z]/(b1*n1);
00115   G4double x2= w2*PP2[Z]/(b2*n2);
00116 
00117   do {
00118 
00119     G4double w = w0;
00120     G4double n = n0;
00121     G4double b = b0;
00122   
00123     x = G4UniformRand()*(x0 + x1 + x2);
00124     if(x > x0) {
00125       x -= x0;
00126       if(x <= x1 ) {
00127         w = w1;
00128         n = n1;
00129         b = b1;
00130       } else { 
00131         w = w2;
00132         n = n2;
00133         b = b2;
00134       }
00135     }
00136     n = 1.0/n;
00137 
00138     // sampling of angle
00139     G4double y = w*G4UniformRand();
00140     if(y < numlim) { x = y*n*( 1 + 0.5*(n + 1)*y*(1 - (n + 2)*y/3.)); }
00141     else           { x = 1.0/std::pow(1 - y, n) - 1.0; }
00142     cost = 1.0 - x/(b*xx);
00143     //G4cout << "cost = " << cost << " w= " << w << " n= " << n 
00144     //     << " b= " << b << " x= " << x << " xx= " << xx << G4endl;  
00145   } while (2*G4UniformRand() > 1.0 + cost*cost || cost < -1.0);
00146 
00147   G4double phi  = twopi*G4UniformRand();
00148   G4double sint = sqrt((1. - cost)*(1.0 + cost));
00149   fLocalDirection.set(sint*cos(phi),sint*sin(phi),cost);
00150   fLocalDirection.rotateUz(dp->GetMomentumDirection());
00151 
00152   return fLocalDirection;
00153 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:15 2013 for Geant4 by  doxygen 1.4.7