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 #include "G4DNAScreenedRutherfordElasticModel.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4DNAMolecularMaterial.hh"
00033
00034
00035
00036 using namespace std;
00037
00038
00039
00040 G4DNAScreenedRutherfordElasticModel::G4DNAScreenedRutherfordElasticModel
00041 (const G4ParticleDefinition*, const G4String& nam)
00042 :G4VEmModel(nam),isInitialised(false)
00043 {
00044
00045 fpWaterDensity = 0;
00046
00047 killBelowEnergy = 9*eV;
00048 lowEnergyLimit = 0 * eV;
00049 intermediateEnergyLimit = 200 * eV;
00050 highEnergyLimit = 1. * MeV;
00051 SetLowEnergyLimit(lowEnergyLimit);
00052 SetHighEnergyLimit(highEnergyLimit);
00053
00054 verboseLevel= 0;
00055
00056
00057
00058
00059
00060
00061
00062 if( verboseLevel>0 )
00063 {
00064 G4cout << "Screened Rutherford Elastic model is constructed " << G4endl
00065 << "Energy range: "
00066 << lowEnergyLimit / eV << " eV - "
00067 << highEnergyLimit / MeV << " MeV"
00068 << G4endl;
00069 }
00070 fParticleChangeForGamma = 0;
00071 }
00072
00073
00074
00075 G4DNAScreenedRutherfordElasticModel::~G4DNAScreenedRutherfordElasticModel()
00076 {}
00077
00078
00079
00080 void G4DNAScreenedRutherfordElasticModel::Initialise(const G4ParticleDefinition* ,
00081 const G4DataVector& )
00082 {
00083
00084 if (verboseLevel > 3)
00085 G4cout << "Calling G4DNAScreenedRutherfordElasticModel::Initialise()" << G4endl;
00086
00087
00088
00089 if (LowEnergyLimit() < lowEnergyLimit)
00090 {
00091 G4cout << "G4DNAScreenedRutherfordElasticModel: low energy limit increased from " <<
00092 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
00093 SetLowEnergyLimit(lowEnergyLimit);
00094 }
00095
00096 if (HighEnergyLimit() > highEnergyLimit)
00097 {
00098 G4cout << "G4DNAScreenedRutherfordElasticModel: high energy limit decreased from " <<
00099 HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
00100 SetHighEnergyLimit(highEnergyLimit);
00101 }
00102
00103
00104
00105 betaCoeff.push_back(7.51525);
00106 betaCoeff.push_back(-0.41912);
00107 betaCoeff.push_back(7.2017E-3);
00108 betaCoeff.push_back(-4.646E-5);
00109 betaCoeff.push_back(1.02897E-7);
00110
00111 deltaCoeff.push_back(2.9612);
00112 deltaCoeff.push_back(-0.26376);
00113 deltaCoeff.push_back(4.307E-3);
00114 deltaCoeff.push_back(-2.6895E-5);
00115 deltaCoeff.push_back(5.83505E-8);
00116
00117 gamma035_10Coeff.push_back(-1.7013);
00118 gamma035_10Coeff.push_back(-1.48284);
00119 gamma035_10Coeff.push_back(0.6331);
00120 gamma035_10Coeff.push_back(-0.10911);
00121 gamma035_10Coeff.push_back(8.358E-3);
00122 gamma035_10Coeff.push_back(-2.388E-4);
00123
00124 gamma10_100Coeff.push_back(-3.32517);
00125 gamma10_100Coeff.push_back(0.10996);
00126 gamma10_100Coeff.push_back(-4.5255E-3);
00127 gamma10_100Coeff.push_back(5.8372E-5);
00128 gamma10_100Coeff.push_back(-2.4659E-7);
00129
00130 gamma100_200Coeff.push_back(2.4775E-2);
00131 gamma100_200Coeff.push_back(-2.96264E-5);
00132 gamma100_200Coeff.push_back(-1.20655E-7);
00133
00134
00135
00136 if( verboseLevel>0 )
00137 {
00138 G4cout << "Screened Rutherford elastic model is initialized " << G4endl
00139 << "Energy range: "
00140 << LowEnergyLimit() / eV << " eV - "
00141 << HighEnergyLimit() / MeV << " MeV"
00142 << G4endl;
00143 }
00144
00145
00146 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00147
00148 if (isInitialised) { return; }
00149 fParticleChangeForGamma = GetParticleChangeForGamma();
00150 isInitialised = true;
00151
00152 }
00153
00154
00155
00156 G4double G4DNAScreenedRutherfordElasticModel::CrossSectionPerVolume(const G4Material* material,
00157 const G4ParticleDefinition* particleDefinition,
00158 G4double ekin,
00159 G4double,
00160 G4double)
00161 {
00162 if (verboseLevel > 3)
00163 G4cout << "Calling CrossSectionPerVolume() of G4DNAScreenedRutherfordElasticModel" << G4endl;
00164
00165
00166
00167 G4double sigma=0;
00168
00169 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00170
00171 if(waterDensity!= 0.0)
00172
00173 {
00174
00175 if (ekin < highEnergyLimit)
00176 {
00177
00178 if (ekin < killBelowEnergy) return DBL_MAX;
00179
00180 G4double z = 10.;
00181 G4double n = ScreeningFactor(ekin,z);
00182 G4double crossSection = RutherfordCrossSection(ekin, z);
00183 sigma = pi * crossSection / (n * (n + 1.));
00184 }
00185
00186 if (verboseLevel > 2)
00187 {
00188 G4cout << "__________________________________" << G4endl;
00189 G4cout << "°°° G4DNAScreenedRutherfordElasticModel - XS INFO START" << G4endl;
00190 G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
00191 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00192 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00193
00194 G4cout << "°°° G4DNAScreenedRutherfordElasticModel - XS INFO END" << G4endl;
00195 }
00196
00197 }
00198
00199 return sigma*material->GetAtomicNumDensityVector()[1];
00200 }
00201
00202
00203
00204 G4double G4DNAScreenedRutherfordElasticModel::RutherfordCrossSection(G4double k, G4double z)
00205 {
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215 G4double length =(e_squared * (k + electron_mass_c2)) / (4 * pi *epsilon0 * k * ( k + 2 * electron_mass_c2));
00216 G4double cross = z * ( z + 1) * length * length;
00217
00218 return cross;
00219 }
00220
00221
00222
00223 G4double G4DNAScreenedRutherfordElasticModel::ScreeningFactor(G4double k, G4double z)
00224 {
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237 const G4double alpha_1(1.64);
00238 const G4double beta_1(-0.0825);
00239 const G4double constK(1.7E-5);
00240
00241 G4double numerator = (alpha_1 + beta_1 * std::log(k/eV)) * constK * std::pow(z, 2./3.);
00242
00243 k /= electron_mass_c2;
00244
00245 G4double denominator = k * (2 + k);
00246
00247 G4double value = 0.;
00248 if (denominator > 0.) value = numerator / denominator;
00249
00250 return value;
00251 }
00252
00253
00254
00255 void G4DNAScreenedRutherfordElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
00256 const G4MaterialCutsCouple* ,
00257 const G4DynamicParticle* aDynamicElectron,
00258 G4double,
00259 G4double)
00260 {
00261
00262 if (verboseLevel > 3)
00263 G4cout << "Calling SampleSecondaries() of G4DNAScreenedRutherfordElasticModel" << G4endl;
00264
00265 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
00266
00267 if (electronEnergy0 < killBelowEnergy)
00268 {
00269 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
00270 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
00271 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
00272 return ;
00273 }
00274
00275 G4double cosTheta = 0.;
00276
00277 if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
00278 {
00279 if (electronEnergy0<intermediateEnergyLimit)
00280 {
00281 if (verboseLevel > 3) G4cout << "---> Using Brenner & Zaider model" << G4endl;
00282 cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
00283 }
00284
00285 if (electronEnergy0>=intermediateEnergyLimit)
00286 {
00287 if (verboseLevel > 3) G4cout << "---> Using Screened Rutherford model" << G4endl;
00288 G4double z = 10.;
00289 cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
00290 }
00291
00292 G4double phi = 2. * pi * G4UniformRand();
00293
00294 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
00295 G4ThreeVector xVers = zVers.orthogonal();
00296 G4ThreeVector yVers = zVers.cross(xVers);
00297
00298 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
00299 G4double yDir = xDir;
00300 xDir *= std::cos(phi);
00301 yDir *= std::sin(phi);
00302
00303 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
00304
00305 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
00306
00307 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
00308 }
00309
00310 }
00311
00312
00313
00314 G4double G4DNAScreenedRutherfordElasticModel::BrennerZaiderRandomizeCosTheta(G4double k)
00315 {
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326 k /= eV;
00327
00328 G4double beta = std::exp(CalculatePolynomial(k,betaCoeff));
00329 G4double delta = std::exp(CalculatePolynomial(k,deltaCoeff));
00330 G4double gamma;
00331
00332 if (k > 100.)
00333 {
00334 gamma = CalculatePolynomial(k, gamma100_200Coeff);
00335
00336 }
00337 else
00338 {
00339 if (k>10)
00340 {
00341 gamma = std::exp(CalculatePolynomial(k, gamma10_100Coeff));
00342 }
00343 else
00344 {
00345 gamma = std::exp(CalculatePolynomial(k, gamma035_10Coeff));
00346 }
00347 }
00348
00349
00350
00351 G4double oneOverMax = 1. / (1./(4.*gamma*gamma) + beta/( (2.+2.*delta)*(2.+2.*delta) ));
00352
00353 G4double cosTheta = 0.;
00354 G4double leftDenominator = 0.;
00355 G4double rightDenominator = 0.;
00356 G4double fCosTheta = 0.;
00357
00358 do
00359 {
00360 cosTheta = 2. * G4UniformRand() - 1.;
00361
00362 leftDenominator = (1. + 2.*gamma - cosTheta);
00363 rightDenominator = (1. + 2.*delta + cosTheta);
00364 if ( (leftDenominator * rightDenominator) != 0. )
00365 {
00366 fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
00367 }
00368 }
00369 while (fCosTheta < G4UniformRand());
00370
00371 return cosTheta;
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412 }
00413
00414
00415
00416 G4double G4DNAScreenedRutherfordElasticModel::CalculatePolynomial(G4double k, std::vector<G4double>& vec)
00417 {
00418
00419
00420
00421
00422 G4double result = 0.;
00423 size_t size = vec.size();
00424
00425 while (size>0)
00426 {
00427 size--;
00428
00429 result *= k;
00430 result += vec[size];
00431 }
00432
00433 return result;
00434 }
00435
00436
00437
00438 G4double G4DNAScreenedRutherfordElasticModel::ScreenedRutherfordRandomizeCosTheta(G4double k, G4double z)
00439 {
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453 G4double n = ScreeningFactor(k, z);
00454
00455 G4double oneOverMax = (4.*n*n);
00456
00457 G4double cosTheta = 0.;
00458 G4double fCosTheta;
00459
00460 do
00461 {
00462 cosTheta = 2. * G4UniformRand() - 1.;
00463 fCosTheta = (1 + 2.*n - cosTheta);
00464 if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta);
00465 }
00466 while (fCosTheta < G4UniformRand());
00467
00468 return cosTheta;
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501 }
00502