#include <G4RayleighAngularGenerator.hh>
Inheritance diagram for G4RayleighAngularGenerator:
Public Member Functions | |
| G4RayleighAngularGenerator () | |
| virtual | ~G4RayleighAngularGenerator () |
| virtual G4ThreeVector & | SampleDirection (const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=0) |
Definition at line 61 of file G4RayleighAngularGenerator.hh.
| 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] |
| 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 }
1.4.7