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
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069 #include "G4UniversalFluctuation.hh"
00070 #include "G4PhysicalConstants.hh"
00071 #include "G4SystemOfUnits.hh"
00072 #include "Randomize.hh"
00073 #include "G4Poisson.hh"
00074 #include "G4Step.hh"
00075 #include "G4Material.hh"
00076 #include "G4DynamicParticle.hh"
00077 #include "G4ParticleDefinition.hh"
00078
00079
00080
00081 using namespace std;
00082
00083 G4UniversalFluctuation::G4UniversalFluctuation(const G4String& nam)
00084 :G4VEmFluctuationModel(nam),
00085 particle(0),
00086 minNumberInteractionsBohr(10.0),
00087 theBohrBeta2(50.0*keV/proton_mass_c2),
00088 minLoss(10.*eV),
00089 nmaxCont(16.),
00090 rate(0.55),
00091 fw(4.)
00092 {
00093 lastMaterial = 0;
00094
00095 particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct
00096 = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 = esmall
00097 = e1 = e2 = 0;
00098
00099 }
00100
00101
00102
00103 G4UniversalFluctuation::~G4UniversalFluctuation()
00104 {}
00105
00106
00107
00108 void G4UniversalFluctuation::InitialiseMe(const G4ParticleDefinition* part)
00109 {
00110 particle = part;
00111 particleMass = part->GetPDGMass();
00112 G4double q = part->GetPDGCharge()/eplus;
00113 chargeSquare = q*q;
00114 }
00115
00116
00117
00118 G4double G4UniversalFluctuation::SampleFluctuations(const G4Material* material,
00119 const G4DynamicParticle* dp,
00120 G4double& tmax,
00121 G4double& length,
00122 G4double& averageLoss)
00123 {
00124
00125
00126
00127
00128
00129
00130
00131
00132 G4double meanLoss = averageLoss;
00133 G4double tkin = dp->GetKineticEnergy();
00134
00135 if (meanLoss < minLoss) { return meanLoss; }
00136
00137 if(!particle) { InitialiseMe(dp->GetDefinition()); }
00138
00139 G4double tau = tkin/particleMass;
00140 G4double gam = tau + 1.0;
00141 G4double gam2 = gam*gam;
00142 G4double beta2 = tau*(tau + 2.0)/gam2;
00143
00144 G4double loss(0.), siga(0.);
00145
00146
00147
00148
00149
00150 if ((particleMass > electron_mass_c2) &&
00151 (meanLoss >= minNumberInteractionsBohr*tmax))
00152 {
00153 G4double massrate = electron_mass_c2/particleMass ;
00154 G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
00155 (1.+massrate*(2.*gam+massrate)) ;
00156 if (tmaxkine <= 2.*tmax)
00157 {
00158 electronDensity = material->GetElectronDensity();
00159 siga = sqrt((1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
00160 * electronDensity * chargeSquare);
00161
00162
00163 G4double sn = meanLoss/siga;
00164
00165
00166 if (sn >= 2.0) {
00167
00168 G4double twomeanLoss = meanLoss + meanLoss;
00169 do {
00170 loss = G4RandGauss::shoot(meanLoss,siga);
00171 } while (0.0 > loss || twomeanLoss < loss);
00172
00173
00174 } else {
00175
00176 G4double neff = sn*sn;
00177 loss = meanLoss*CLHEP::RandGamma::shoot(neff,1.0)/neff;
00178 }
00179
00180 return loss;
00181 }
00182 }
00183
00184
00185
00186 if (material != lastMaterial) {
00187 f1Fluct = material->GetIonisation()->GetF1fluct();
00188 f2Fluct = material->GetIonisation()->GetF2fluct();
00189 e1Fluct = material->GetIonisation()->GetEnergy1fluct();
00190 e2Fluct = material->GetIonisation()->GetEnergy2fluct();
00191 e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
00192 e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
00193 ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy();
00194 ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
00195 e0 = material->GetIonisation()->GetEnergy0fluct();
00196 esmall = 0.5*sqrt(e0*ipotFluct);
00197 lastMaterial = material;
00198 }
00199
00200
00201 if(tmax <= e0) { return meanLoss; }
00202
00203 G4double losstot = 0.;
00204 G4int nstep = 1;
00205 if(meanLoss < 25.*ipotFluct)
00206 {
00207 if(G4UniformRand() < 0.04*meanLoss/ipotFluct)
00208 { nstep = 1; }
00209 else
00210 {
00211 nstep = 2;
00212 meanLoss *= 0.5;
00213 }
00214 }
00215
00216 for (G4int istep=0; istep < nstep; ++istep) {
00217
00218 loss = 0.;
00219
00220 G4double a1 = 0. , a2 = 0., a3 = 0. ;
00221
00222 if(tmax > ipotFluct) {
00223 G4double w2 = log(2.*electron_mass_c2*beta2*gam2)-beta2;
00224
00225 if(w2 > ipotLogFluct) {
00226 G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
00227 a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
00228 if(w2 > e2LogFluct) {
00229 a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
00230 }
00231 if(a1 < nmaxCont) {
00232
00233 G4double sa1 = sqrt(a1);
00234 if(G4UniformRand() < exp(-sa1))
00235 {
00236 e1 = esmall;
00237 a1 = meanLoss*(1.-rate)/e1;
00238 a2 = 0.;
00239 e2 = e2Fluct;
00240 }
00241 else
00242 {
00243 a1 = sa1 ;
00244 e1 = sa1*e1Fluct;
00245 e2 = e2Fluct;
00246 }
00247
00248 } else {
00249
00250
00251 a1 /= fw;
00252 e1 = fw*e1Fluct;
00253 e2 = e2Fluct;
00254 }
00255 }
00256 }
00257
00258 G4double w1 = tmax/e0;
00259 if(tmax > e0) {
00260 a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*log(w1));
00261 }
00262
00263 G4double emean = 0.;
00264 G4double sig2e = 0., sige = 0.;
00265 G4double p1 = 0., p2 = 0., p3 = 0.;
00266
00267
00268 if(a1 > nmaxCont)
00269 {
00270 emean += a1*e1;
00271 sig2e += a1*e1*e1;
00272 }
00273 else if(a1 > 0.)
00274 {
00275 p1 = G4double(G4Poisson(a1));
00276 loss += p1*e1;
00277 if(p1 > 0.) {
00278 loss += (1.-2.*G4UniformRand())*e1;
00279 }
00280 }
00281
00282
00283
00284 if(a2 > nmaxCont)
00285 {
00286 emean += a2*e2;
00287 sig2e += a2*e2*e2;
00288 }
00289 else if(a2 > 0.)
00290 {
00291 p2 = G4double(G4Poisson(a2));
00292 loss += p2*e2;
00293 if(p2 > 0.)
00294 loss += (1.-2.*G4UniformRand())*e2;
00295 }
00296
00297 if(emean > 0.)
00298 {
00299 sige = sqrt(sig2e);
00300 loss += max(0.,G4RandGauss::shoot(emean,sige));
00301 }
00302
00303
00304 G4double lossc = 0.;
00305 if(a3 > 0.) {
00306 emean = 0.;
00307 sig2e = 0.;
00308 sige = 0.;
00309 p3 = a3;
00310 G4double alfa = 1.;
00311 if(a3 > nmaxCont)
00312 {
00313 alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
00314 G4double alfa1 = alfa*log(alfa)/(alfa-1.);
00315 G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
00316 emean += namean*e0*alfa1;
00317 sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
00318 p3 = a3-namean;
00319 }
00320
00321 G4double w2 = alfa*e0;
00322 G4double w = (tmax-w2)/tmax;
00323 G4int nb = G4Poisson(p3);
00324 if(nb > 0) {
00325 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
00326 }
00327 }
00328
00329 if(emean > 0.)
00330 {
00331 sige = sqrt(sig2e);
00332 lossc += max(0.,G4RandGauss::shoot(emean,sige));
00333 }
00334
00335 loss += lossc;
00336
00337 losstot += loss;
00338 }
00339
00340
00341 return losstot;
00342
00343 }
00344
00345
00346
00347
00348 G4double G4UniversalFluctuation::Dispersion(
00349 const G4Material* material,
00350 const G4DynamicParticle* dp,
00351 G4double& tmax,
00352 G4double& length)
00353 {
00354 if(!particle) { InitialiseMe(dp->GetDefinition()); }
00355
00356 electronDensity = material->GetElectronDensity();
00357
00358 G4double gam = (dp->GetKineticEnergy())/particleMass + 1.0;
00359 G4double beta2 = 1.0 - 1.0/(gam*gam);
00360
00361 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
00362 * electronDensity * chargeSquare;
00363
00364 return siga;
00365 }
00366
00367
00368
00369 void
00370 G4UniversalFluctuation::SetParticleAndCharge(const G4ParticleDefinition* part,
00371 G4double q2)
00372 {
00373 if(part != particle) {
00374 particle = part;
00375 particleMass = part->GetPDGMass();
00376 }
00377 chargeSquare = q2;
00378 }
00379
00380