00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 #include "G4DipBustGenerator.hh"
00054 #include "G4PhysicalConstants.hh"
00055 #include "Randomize.hh"
00056
00057 G4DipBustGenerator::G4DipBustGenerator(const G4String&)
00058 : G4VEmAngularDistribution("DipBustGen")
00059 {}
00060
00061 G4DipBustGenerator::~G4DipBustGenerator()
00062 {}
00063
00064 G4ThreeVector&
00065 G4DipBustGenerator::SampleDirection(const G4DynamicParticle* dp,
00066 G4double, G4int, const G4Material*)
00067 {
00068 G4double c, cosTheta, delta, cofA, signc = 1., a, power = 1./3.;
00069
00070 G4double eTkin = dp->GetKineticEnergy();
00071
00072 c = 4. - 8.*G4UniformRand();
00073 a = c;
00074
00075 if( c < 0. )
00076 {
00077 signc = -1.;
00078 a = -c;
00079 }
00080 delta = std::sqrt(a*a+4.);
00081 delta += a;
00082 delta *= 0.5;
00083
00084 cofA = -signc*std::pow(delta, power);
00085
00086 cosTheta = cofA - 1./cofA;
00087
00088 G4double tau = eTkin/electron_mass_c2;
00089 G4double beta = std::sqrt(tau*(tau + 2.))/(tau + 1.);
00090
00091 cosTheta = (cosTheta + beta)/(1 + cosTheta*beta);
00092
00093 G4double sinTheta = std::sqrt((1 - cosTheta)*(1 + cosTheta));
00094 G4double phi = twopi*G4UniformRand();
00095
00096 fLocalDirection.set(sinTheta*std::cos(phi), sinTheta*std::sin(phi),cosTheta);
00097 fLocalDirection.rotateUz(dp->GetMomentumDirection());
00098
00099 return fLocalDirection;
00100
00101 }
00102
00103 G4double G4DipBustGenerator::PolarAngle(const G4double eTkin,
00104 const G4double,
00105 const G4int )
00106 {
00107 G4double c, cosTheta, delta, cofA, signc = 1., a, power = 1./3.;
00108 G4double gamma, beta, theta;
00109
00110 c = 4. - 8.*G4UniformRand();
00111 a = c;
00112
00113 if( c < 0. )
00114 {
00115 signc = -1.;
00116 a = -c;
00117 }
00118 delta = std::sqrt(a*a+4.);
00119 delta += a;
00120 delta *= 0.5;
00121
00122 cofA = -signc*std::pow(delta, power);
00123
00124 cosTheta = cofA - 1./cofA;
00125
00126 gamma = 1. + eTkin/electron_mass_c2;
00127 beta = std::sqrt(1. - 1./gamma/gamma);
00128
00129 cosTheta = (cosTheta + beta)/(1 + cosTheta*beta);
00130
00131 theta = std::acos(cosTheta);
00132
00133 if( theta < 0. ) theta = 0.;
00134 if( theta > pi ) theta = pi;
00135
00136
00137 return theta;
00138 }
00139
00140 void G4DipBustGenerator::PrintGeneratorInformation() const
00141 {
00142 G4cout << "\n" << G4endl;
00143 G4cout << "Angular Generator based on classical formula from" << G4endl;
00144 G4cout << "J.D. Jackson, Classical Electrodynamics, Wiley, New York 1975"
00145 << G4endl;
00146 }