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 #include "G4BohrFluctuations.hh"
00053 #include "G4PhysicalConstants.hh"
00054 #include "G4SystemOfUnits.hh"
00055 #include "Randomize.hh"
00056 #include "G4Poisson.hh"
00057 #include "G4ParticleDefinition.hh"
00058
00059
00060
00061 using namespace std;
00062
00063 G4BohrFluctuations::G4BohrFluctuations(const G4String& nam)
00064 :G4VEmFluctuationModel(nam),
00065 particle(0),
00066 minNumberInteractionsBohr(2.0),
00067 minFraction(0.2),
00068 xmin(0.2),
00069 minLoss(0.001*eV)
00070 {
00071 particleMass = proton_mass_c2;
00072 chargeSquare = 1.0;
00073 kineticEnergy = 0.0;
00074 beta2 = 0.0;
00075 }
00076
00077
00078
00079 G4BohrFluctuations::~G4BohrFluctuations()
00080 {}
00081
00082
00083
00084 void G4BohrFluctuations::InitialiseMe(const G4ParticleDefinition* part)
00085 {
00086 particle = part;
00087 particleMass = part->GetPDGMass();
00088 G4double q = part->GetPDGCharge()/eplus;
00089 chargeSquare = q*q;
00090 }
00091
00092
00093
00094 G4double G4BohrFluctuations::SampleFluctuations(const G4Material* material,
00095 const G4DynamicParticle* dp,
00096 G4double& tmax,
00097 G4double& length,
00098 G4double& meanLoss)
00099 {
00100 if(meanLoss <= minLoss) { return meanLoss; }
00101 G4double siga = Dispersion(material,dp,tmax,length);
00102 G4double loss = meanLoss;
00103
00104 G4double navr = meanLoss*meanLoss/siga;
00105
00106 if (navr >= minNumberInteractionsBohr) {
00107
00108
00109 if ( meanLoss > minFraction*kineticEnergy ) {
00110 G4double gam = (kineticEnergy - meanLoss)/particleMass + 1.0;
00111 G4double b2 = 1.0 - 1.0/(gam*gam);
00112 if(b2 < xmin*beta2) b2 = xmin*beta2;
00113 G4double x = b2/beta2;
00114 G4double x3 = 1.0/(x*x*x);
00115 siga *= 0.25*(1.0 + x)*(x3 + (1.0/b2 - 0.5)/(1.0/beta2 - 0.5) );
00116 }
00117 siga = sqrt(siga);
00118 G4double twomeanLoss = meanLoss + meanLoss;
00119
00120
00121 if(twomeanLoss < siga) {
00122 G4double x;
00123 do {
00124 loss = twomeanLoss*G4UniformRand();
00125 x = (loss - meanLoss)/siga;
00126 } while (1.0 - 0.5*x*x < G4UniformRand());
00127 } else {
00128 do {
00129 loss = G4RandGauss::shoot(meanLoss,siga);
00130 } while (0.0 > loss || loss > twomeanLoss);
00131 }
00132
00133
00134 } else {
00135 G4double n = (G4double)(G4Poisson(navr));
00136 loss = meanLoss*n/navr;
00137 }
00138
00139
00140 return loss;
00141 }
00142
00143
00144
00145 G4double G4BohrFluctuations::Dispersion(const G4Material* material,
00146 const G4DynamicParticle* dp,
00147 G4double& tmax,
00148 G4double& length)
00149 {
00150 if(!particle) { InitialiseMe(dp->GetDefinition()); }
00151
00152 G4double electronDensity = material->GetElectronDensity();
00153 kineticEnergy = dp->GetKineticEnergy();
00154 G4double etot = kineticEnergy + particleMass;
00155 beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
00156 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
00157 * electronDensity * chargeSquare;
00158
00159 return siga;
00160 }
00161
00162
00163
00164