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 #include "G4IonFluctuations.hh"
00059 #include "G4PhysicalConstants.hh"
00060 #include "G4SystemOfUnits.hh"
00061 #include "Randomize.hh"
00062 #include "G4Poisson.hh"
00063 #include "G4Material.hh"
00064 #include "G4DynamicParticle.hh"
00065
00066
00067
00068 using namespace std;
00069
00070 G4IonFluctuations::G4IonFluctuations(const G4String& nam)
00071 : G4VEmFluctuationModel(nam),
00072 particle(0),
00073 particleMass(proton_mass_c2),
00074 charge(1.0),
00075 chargeSquare(1.0),
00076 effChargeSquare(1.0),
00077 parameter(10.0*CLHEP::MeV/CLHEP::proton_mass_c2),
00078 minNumberInteractionsBohr(0.0),
00079 theBohrBeta2(50.0*keV/CLHEP::proton_mass_c2),
00080 minFraction(0.2),
00081 xmin(0.2),
00082 minLoss(0.001*eV)
00083 {
00084 kineticEnergy = 0.0;
00085 beta2 = 0.0;
00086 }
00087
00088
00089
00090 G4IonFluctuations::~G4IonFluctuations()
00091 {}
00092
00093
00094
00095 void G4IonFluctuations::InitialiseMe(const G4ParticleDefinition* part)
00096 {
00097 particle = part;
00098 particleMass = part->GetPDGMass();
00099 charge = part->GetPDGCharge()/eplus;
00100 chargeSquare = charge*charge;
00101 effChargeSquare= chargeSquare;
00102 uniFluct.InitialiseMe(part);
00103 }
00104
00105
00106
00107 G4double G4IonFluctuations::SampleFluctuations(const G4Material* material,
00108 const G4DynamicParticle* dp,
00109 G4double& tmax,
00110 G4double& length,
00111 G4double& meanLoss)
00112 {
00113
00114 if(meanLoss <= minLoss) return meanLoss;
00115
00116
00117
00118
00119
00120 if(dp->GetKineticEnergy() > parameter*charge*particleMass) {
00121 return uniFluct.SampleFluctuations(material,dp,tmax,length,meanLoss);
00122 }
00123
00124 G4double siga = Dispersion(material,dp,tmax,length);
00125 G4double loss = meanLoss;
00126
00127
00128
00129
00130
00131
00132
00133 if ( meanLoss > minFraction*kineticEnergy ) {
00134 G4double gam = (kineticEnergy - meanLoss)/particleMass + 1.0;
00135 G4double b2 = 1.0 - 1.0/(gam*gam);
00136 if(b2 < xmin*beta2) b2 = xmin*beta2;
00137 G4double x = b2/beta2;
00138 G4double x3 = 1.0/(x*x*x);
00139 siga *= 0.25*(1.0 + x)*(x3 + (1.0/b2 - 0.5)/(1.0/beta2 - 0.5) );
00140 }
00141 siga = sqrt(siga);
00142 G4double sn = meanLoss/siga;
00143 G4double twomeanLoss = meanLoss + meanLoss;
00144
00145
00146
00147 if (sn >= 2.0) {
00148
00149 do {
00150 loss = G4RandGauss::shoot(meanLoss,siga);
00151 } while (0.0 > loss || twomeanLoss < loss);
00152
00153
00154 } else if(sn > 0.1) {
00155
00156 G4double neff = sn*sn;
00157 loss = meanLoss*CLHEP::RandGamma::shoot(neff,1.0)/neff;
00158
00159
00160 } else {
00161 loss = twomeanLoss*G4UniformRand();
00162 }
00163
00164
00165 return loss;
00166 }
00167
00168
00169
00170 G4double G4IonFluctuations::Dispersion(const G4Material* material,
00171 const G4DynamicParticle* dp,
00172 G4double& tmax,
00173 G4double& length)
00174 {
00175 kineticEnergy = dp->GetKineticEnergy();
00176 G4double etot = kineticEnergy + particleMass;
00177 beta2 = kineticEnergy*(kineticEnergy + 2.*particleMass)/(etot*etot);
00178
00179 G4double electronDensity = material->GetElectronDensity();
00180
00181
00182
00183
00184
00185
00186 G4double siga = (1. - beta2*0.5)*tmax*length*electronDensity*
00187 twopi_mc2_rcl2*chargeSquare/beta2;
00188
00189
00190
00191
00192
00193 G4double Z = electronDensity/material->GetTotNbOfAtomsPerVolume();
00194
00195 G4double fac = Factor(material, Z);
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 G4double fac_cut = 1.0 + (fac - 1.0)*2.0*electron_mass_c2*beta2/(tmax*(1.0 - beta2));
00206 if(fac_cut > 0.01 && fac > 0.01) {
00207 siga *= fac_cut;
00208 }
00209
00210
00211
00212
00213 return siga;
00214 }
00215
00216
00217
00218 G4double G4IonFluctuations::Factor(const G4Material* material, G4double Z)
00219 {
00220
00221
00222
00223
00224 G4double energy = kineticEnergy *amu_c2/(particleMass*MeV) ;
00225
00226
00227 G4double s1 = RelativisticFactor(material, Z);
00228
00229
00230 if( beta2 < 3.0*theBohrBeta2*Z ) {
00231
00232 static G4double a[96][4] = {
00233 {-0.3291, -0.8312, 0.2460, -1.0220},
00234 {-0.5615, -0.5898, 0.5205, -0.7258},
00235 {-0.5280, -0.4981, 0.5519, -0.5865},
00236 {-0.5125, -0.4625, 0.5660, -0.5190},
00237 {-0.5127, -0.8595, 0.5626, -0.8721},
00238 {-0.5174, -1.1930, 0.5565, -1.1980},
00239 {-0.5179, -1.1850, 0.5560, -1.2070},
00240 {-0.5209, -0.9355, 0.5590, -1.0250},
00241 {-0.5255, -0.7766, 0.5720, -0.9412},
00242
00243 {-0.5776, -0.6665, 0.6598, -0.8484},
00244 {-0.6013, -0.6045, 0.7321, -0.7671},
00245 {-0.5781, -0.5518, 0.7605, -0.6919},
00246 {-0.5587, -0.4981, 0.7835, -0.6195},
00247 {-0.5466, -0.4656, 0.7978, -0.5771},
00248 {-0.5406, -0.4690, 0.8031, -0.5718},
00249 {-0.5391, -0.5061, 0.8024, -0.5974},
00250 {-0.5380, -0.6483, 0.7962, -0.6970},
00251 {-0.5355, -0.7722, 0.7962, -0.7839},
00252 {-0.5329, -0.7720, 0.7988, -0.7846},
00253
00254 {-0.5335, -0.7671, 0.7984, -0.7933},
00255 {-0.5324, -0.7612, 0.7998, -0.8031},
00256 {-0.5305, -0.7300, 0.8031, -0.7990},
00257 {-0.5307, -0.7178, 0.8049, -0.8216},
00258 {-0.5248, -0.6621, 0.8165, -0.7919},
00259 {-0.5180, -0.6502, 0.8266, -0.7986},
00260 {-0.5084, -0.6408, 0.8396, -0.8048},
00261 {-0.4967, -0.6331, 0.8549, -0.8093},
00262 {-0.4861, -0.6508, 0.8712, -0.8432},
00263 {-0.4700, -0.6186, 0.8961, -0.8132},
00264
00265 {-0.4545, -0.5720, 0.9227, -0.7710},
00266 {-0.4404, -0.5226, 0.9481, -0.7254},
00267 {-0.4288, -0.4778, 0.9701, -0.6850},
00268 {-0.4199, -0.4425, 0.9874, -0.6539},
00269 {-0.4131, -0.4188, 0.9998, -0.6332},
00270 {-0.4089, -0.4057, 1.0070, -0.6218},
00271 {-0.4039, -0.3913, 1.0150, -0.6107},
00272 {-0.3987, -0.3698, 1.0240, -0.5938},
00273 {-0.3977, -0.3608, 1.0260, -0.5852},
00274 {-0.3972, -0.3600, 1.0260, -0.5842},
00275
00276 {-0.3985, -0.3803, 1.0200, -0.6013},
00277 {-0.3985, -0.3979, 1.0150, -0.6168},
00278 {-0.3968, -0.3990, 1.0160, -0.6195},
00279 {-0.3971, -0.4432, 1.0050, -0.6591},
00280 {-0.3944, -0.4665, 1.0010, -0.6825},
00281 {-0.3924, -0.5109, 0.9921, -0.7235},
00282 {-0.3882, -0.5158, 0.9947, -0.7343},
00283 {-0.3838, -0.5125, 0.9999, -0.7370},
00284 {-0.3786, -0.4976, 1.0090, -0.7310},
00285 {-0.3741, -0.4738, 1.0200, -0.7155},
00286
00287 {-0.3969, -0.4496, 1.0320, -0.6982},
00288 {-0.3663, -0.4297, 1.0430, -0.6828},
00289 {-0.3630, -0.4120, 1.0530, -0.6689},
00290 {-0.3597, -0.3964, 1.0620, -0.6564},
00291 {-0.3555, -0.3809, 1.0720, -0.6454},
00292 {-0.3525, -0.3607, 1.0820, -0.6289},
00293 {-0.3505, -0.3465, 1.0900, -0.6171},
00294 {-0.3397, -0.3570, 1.1020, -0.6384},
00295 {-0.3314, -0.3552, 1.1130, -0.6441},
00296 {-0.3235, -0.3531, 1.1230, -0.6498},
00297
00298 {-0.3150, -0.3483, 1.1360, -0.6539},
00299 {-0.3060, -0.3441, 1.1490, -0.6593},
00300 {-0.2968, -0.3396, 1.1630, -0.6649},
00301 {-0.2935, -0.3225, 1.1760, -0.6527},
00302 {-0.2797, -0.3262, 1.1940, -0.6722},
00303 {-0.2704, -0.3202, 1.2100, -0.6770},
00304 {-0.2815, -0.3227, 1.2480, -0.6775},
00305 {-0.2880, -0.3245, 1.2810, -0.6801},
00306 {-0.3034, -0.3263, 1.3270, -0.6778},
00307 {-0.2936, -0.3215, 1.3430, -0.6835},
00308
00309 {-0.3282, -0.3200, 1.3980, -0.6650},
00310 {-0.3260, -0.3070, 1.4090, -0.6552},
00311 {-0.3511, -0.3074, 1.4470, -0.6442},
00312 {-0.3501, -0.3064, 1.4500, -0.6442},
00313 {-0.3490, -0.3027, 1.4550, -0.6418},
00314 {-0.3487, -0.3048, 1.4570, -0.6447},
00315 {-0.3478, -0.3074, 1.4600, -0.6483},
00316 {-0.3501, -0.3283, 1.4540, -0.6669},
00317 {-0.3494, -0.3373, 1.4550, -0.6765},
00318 {-0.3485, -0.3373, 1.4570, -0.6774},
00319
00320 {-0.3462, -0.3300, 1.4630, -0.6728},
00321 {-0.3462, -0.3225, 1.4690, -0.6662},
00322 {-0.3453, -0.3094, 1.4790, -0.6553},
00323 {-0.3844, -0.3134, 1.5240, -0.6412},
00324 {-0.3848, -0.3018, 1.5310, -0.6303},
00325 {-0.3862, -0.2955, 1.5360, -0.6237},
00326 {-0.4262, -0.2991, 1.5860, -0.6115},
00327 {-0.4278, -0.2910, 1.5900, -0.6029},
00328 {-0.4303, -0.2817, 1.5940, -0.5927},
00329 {-0.4315, -0.2719, 1.6010, -0.5829},
00330
00331 {-0.4359, -0.2914, 1.6050, -0.6010},
00332 {-0.4365, -0.2982, 1.6080, -0.6080},
00333 {-0.4253, -0.3037, 1.6120, -0.6150},
00334 {-0.4335, -0.3245, 1.6160, -0.6377},
00335 {-0.4307, -0.3292, 1.6210, -0.6447},
00336 {-0.4284, -0.3204, 1.6290, -0.6380},
00337 {-0.4227, -0.3217, 1.6360, -0.6438}
00338 } ;
00339
00340 G4int iz = G4int(Z) - 2;
00341 if( 0 > iz ) iz = 0;
00342 else if(95 < iz ) iz = 95;
00343
00344 G4double ss = 1.0 + a[iz][0]*pow(energy,a[iz][1])+
00345 + a[iz][2]*pow(energy,a[iz][3]);
00346
00347
00348 G4double slim = 0.001;
00349 if(ss < slim) s1 = 1.0/slim;
00350
00351 else if(s1*ss < 1.0) s1 = 1.0/ss;
00352 }
00353
00354 G4int i = 0 ;
00355 G4double factor = 1.0 ;
00356
00357
00358
00359
00360
00361
00362 static G4double b[5][4] = {
00363 {0.1014, 0.3700, 0.9642, 3.987},
00364 {0.1955, 0.6941, 2.522, 1.040},
00365 {0.05058, 0.08975, 0.1419, 10.80},
00366 {0.05009, 0.08660, 0.2751, 3.787},
00367 {0.01273, 0.03458, 0.3951, 3.812}
00368 } ;
00369
00370
00371 if(1.5 > charge) {
00372 if( kStateGas != material->GetState() ) i = 1 ;
00373
00374
00375 } else {
00376
00377 factor = charge * pow(charge/Z, 0.33333333);
00378
00379 if( kStateGas == material->GetState() ) {
00380 energy /= (charge * sqrt(charge)) ;
00381
00382 if(1 == (material->GetNumberOfElements())) {
00383 i = 2 ;
00384 } else {
00385 i = 3 ;
00386 }
00387
00388 } else {
00389 energy /= (charge * sqrt(charge*Z)) ;
00390 i = 4 ;
00391 }
00392 }
00393
00394 G4double x = b[i][2];
00395 G4double y = energy * b[i][3];
00396 if(y <= 0.2) x *= (y*(1.0 - 0.5*y));
00397 else x *= (1.0 - exp(-y));
00398
00399 y = energy - b[i][1];
00400
00401 G4double s2 = factor * x * b[i][0] / (y*y + x*x);
00402
00403
00404
00405
00406 return s1*effChargeSquare/chargeSquare + s2;
00407 }
00408
00409
00410
00411 G4double G4IonFluctuations::RelativisticFactor(const G4Material* mat,
00412 G4double Z)
00413 {
00414 G4double eF = mat->GetIonisation()->GetFermiEnergy();
00415 G4double I = mat->GetIonisation()->GetMeanExcitationEnergy();
00416
00417
00418 G4double bF2= 2.0*eF/electron_mass_c2;
00419 G4double f = 0.4*(1.0 - beta2)/((1.0 - 0.5*beta2)*Z);
00420 if(beta2 > bF2) f *= log(2.0*electron_mass_c2*beta2/I)*bF2/beta2;
00421 else f *= log(4.0*eF/I);
00422
00423
00424
00425
00426 return 1.0 + f;
00427 }
00428
00429
00430
00431 void G4IonFluctuations::SetParticleAndCharge(const G4ParticleDefinition* part,
00432 G4double q2)
00433 {
00434 if(part != particle) {
00435 particle = part;
00436 particleMass = part->GetPDGMass();
00437 charge = part->GetPDGCharge()/eplus;
00438 chargeSquare = charge*charge;
00439 }
00440 effChargeSquare = q2;
00441 uniFluct.SetParticleAndCharge(part, q2);
00442 }
00443
00444