Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 26 #include "G4DNAUeharaScreenedRutherfordElastic 27 #include "G4PhysicalConstants.hh" 28 #include "G4SystemOfUnits.hh" 29 #include "G4DNAMolecularMaterial.hh" 30 #include "G4Exp.hh" 31 32 //....oooOO0OOooo........oooOO0OOooo........oo 33 34 using namespace std; 35 36 // #define UEHARA_VERBOSE 37 38 //....oooOO0OOooo........oooOO0OOooo........oo 39 40 G4DNAUeharaScreenedRutherfordElasticModel:: 41 G4DNAUeharaScreenedRutherfordElasticModel(cons 42 cons 43 { 44 // Energy limits of the models 45 intermediateEnergyLimit = 200. * eV; 46 iLowEnergyLimit = 9.*eV; 47 iHighEnergyLimit = 1.*MeV; 48 49 verboseLevel = 0; 50 // Verbosity scale: 51 // 0 = nothing 52 // 1 = warning for energy non-conservation 53 // 2 = details of energy budget 54 // 3 = calculation of cross sections, file o 55 // 4 = entering in methods 56 57 fasterCode = false; 58 } 59 60 //....oooOO0OOooo........oooOO0OOooo........oo 61 62 void 63 G4DNAUeharaScreenedRutherfordElasticModel:: 64 Initialise(const G4ParticleDefinition* particl 65 const G4DataVector& /*cuts*/) 66 { 67 if(isInitialised) { return; } 68 if (verboseLevel > 3) 69 { 70 } 71 72 if(particle->GetParticleName() != "e-") 73 { 74 G4Exception("*** WARNING: the G4DNAUeharaS 75 "not intented to be used with 76 "",FatalException,"") ; 77 } 78 79 80 if( verboseLevel>1 ) 81 { 82 G4cout << "G4DNAUeharaScreenedRutherfordEl 83 << G4endl; 84 G4cout << "Energy range: " 85 << LowEnergyLimit() / eV << " eV - " 86 << HighEnergyLimit() / MeV << " MeV" 87 << G4endl; 88 } 89 // Constants for final state by Brenner & Za 90 // Note: the instantiation must be placed af 91 92 betaCoeff= 93 { 94 7.51525, 95 -0.41912, 96 7.2017E-3, 97 -4.646E-5, 98 1.02897E-7}; 99 100 deltaCoeff= 101 { 102 2.9612, 103 -0.26376, 104 4.307E-3, 105 -2.6895E-5, 106 5.83505E-8}; 107 108 gamma035_10Coeff= 109 { 110 -1.7013, 111 -1.48284, 112 0.6331, 113 -0.10911, 114 8.358E-3, 115 -2.388E-4}; 116 117 gamma10_100Coeff = 118 { 119 -3.32517, 120 0.10996, 121 -4.5255E-3, 122 5.8372E-5, 123 -2.4659E-7}; 124 125 gamma100_200Coeff= 126 { 127 2.4775E-2, 128 -2.96264E-5, 129 -1.20655E-7}; 130 131 // Initialize water density pointer 132 fpWaterDensity = 133 G4DNAMolecularMaterial::Instance()-> 134 GetNumMolPerVolTableFor(G4Material::GetMat 135 136 fParticleChangeForGamma = GetParticleChangeF 137 isInitialised = true; 138 } 139 140 //....oooOO0OOooo........oooOO0OOooo........oo 141 142 G4double 143 G4DNAUeharaScreenedRutherfordElasticModel:: 144 CrossSectionPerVolume(const G4Material* materi 145 const G4ParticleDefiniti 146 G4double ekin, 147 G4double, 148 G4double) 149 { 150 #ifdef UEHARA_VERBOSE 151 if (verboseLevel > 3) 152 { 153 G4cout 154 << "Calling CrossSectionPerVolume() of G4D 155 << G4endl; 156 } 157 #endif 158 159 // Calculate total cross section for model 160 161 G4double sigma = 0.; 162 if(ekin < iLowEnergyLimit || ekin > iHighEne 163 164 G4double waterDensity = (*fpWaterDensity)[ma 165 166 G4double z = 7.42; // FROM PMB 37 (1992) 184 167 G4double n = ScreeningFactor(ekin,z); 168 G4double crossSection = RutherfordCrossSecti 169 sigma = pi * crossSection / (n * (n + 1.)); 170 171 #ifdef UEHARA_VERBOSE 172 if (verboseLevel > 2) 173 { 174 G4cout << "_______________________________ 175 G4cout << "=== G4DNAUeharaScreenedRutherfo 176 << G4endl; 177 G4cout << "=== Kinetic energy(eV)=" << eki 178 << " particle : " << particleDefini 179 G4cout << "=== Cross section per water mol 180 << G4endl; 181 G4cout << "=== Cross section per water mol 182 << sigma*waterDensity/(1./cm) << G4 183 G4cout << "=== G4DNAUeharaScreenedRutherfo 184 << G4endl; 185 } 186 #endif 187 188 return sigma*waterDensity; 189 } 190 191 //....oooOO0OOooo........oooOO0OOooo........oo 192 193 G4double 194 G4DNAUeharaScreenedRutherfordElasticModel::Rut 195 196 { 197 // 198 // e^4 199 // sigma_Ruth(K) = Z (Z+1) ----------------- 200 // (4 pi epsilon_0) 201 // 202 // Where K is the electron non-relativistic 203 // 204 // NIM 155, pp. 145-156, 1978 205 206 G4double length = (e_squared * (k + electron 207 / (4 * pi * epsilon0 * k * (k + 2 * elec 208 G4double cross = z * (z + 1) * length * leng 209 210 return cross; 211 } 212 213 //....oooOO0OOooo........oooOO0OOooo........oo 214 215 G4double G4DNAUeharaScreenedRutherfordElasticM 216 217 { 218 // From Phys Med Biol 37 (1992) 1841-1858 219 // Z=7.42 for water 220 221 const G4double constK(1.7E-5); 222 223 G4double beta2; 224 beta2 = 1. - 1. / ((1. + k / electron_mass_c 225 226 G4double etaC; 227 if (k < 50 * keV) 228 etaC = 1.198; 229 else 230 etaC = 1.13 + 3.76 * (z * z / (137 * 137 * 231 232 G4double numerator = etaC * constK * std::po 233 234 k /= electron_mass_c2; 235 236 G4double denominator = k * (2 + k); 237 238 G4double value = 0.; 239 if (denominator > 0.) 240 value = numerator / denominator; // form 3 241 242 return value; 243 } 244 245 //....oooOO0OOooo........oooOO0OOooo........oo 246 247 void 248 G4DNAUeharaScreenedRutherfordElasticModel:: 249 SampleSecondaries(std::vector<G4DynamicParticl 250 const G4MaterialCutsCouple* 251 const G4DynamicParticle* aDy 252 G4double, 253 G4double) 254 { 255 #ifdef UEHARA_VERBOSE 256 if (verboseLevel > 3) 257 { 258 G4cout 259 << "Calling SampleSecondaries() of G4D 260 << G4endl; 261 } 262 #endif 263 264 G4double electronEnergy0 = aDynamicElectron- 265 if(electronEnergy0 < iLowEnergyLimit || elec 266 return; 267 268 G4double cosTheta = 0.; 269 270 if (electronEnergy0<intermediateEnergyLimit) 271 { 272 #ifdef UEHARA_VERBOSE 273 if (verboseLevel > 3) 274 G4cout << "---> Using Brenner & Zaider m 275 #endif 276 cosTheta = BrennerZaiderRandomizeCosTheta( 277 } 278 else //if (electronEnergy0>=intermediateEner 279 { 280 #ifdef UEHARA_VERBOSE 281 if (verboseLevel > 3) 282 G4cout << "---> Using Screened Rutherfor 283 #endif 284 G4double z = 7.42; // FROM PMB 37 (1992) 285 cosTheta = ScreenedRutherfordRandomizeCosT 286 } 287 288 G4double phi = 2. * pi * G4UniformRand(); 289 290 G4ThreeVector zVers = aDynamicElectron->GetM 291 G4ThreeVector xVers = zVers.orthogonal(); 292 G4ThreeVector yVers = zVers.cross(xVers); 293 294 G4double xDir = std::sqrt(1. - cosTheta*cosT 295 G4double yDir = xDir; 296 xDir *= std::cos(phi); 297 yDir *= std::sin(phi); 298 299 G4ThreeVector zPrimeVers((xDir*xVers + yDir* 300 301 fParticleChangeForGamma->ProposeMomentumDire 302 303 fParticleChangeForGamma->SetProposedKineticE 304 } 305 306 //....oooOO0OOooo........oooOO0OOooo........oo 307 308 G4double 309 G4DNAUeharaScreenedRutherfordElasticModel:: 310 BrennerZaiderRandomizeCosTheta(G4double k) 311 { 312 // d sigma_el 1 313 // ------------ (K) ~ ---------------------- 314 // d Omega (1 + 2 gamma(K) - cos 315 // 316 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/( 317 // 318 // Phys. Med. Biol. 29 N.4 (1983) 443-447 319 320 // gamma(K), beta(K) and delta(K) are polyno 321 322 k /= eV; 323 324 G4double beta = G4Exp(CalculatePolynomial(k, 325 G4double delta = G4Exp(CalculatePolynomial(k 326 G4double gamma; 327 328 if (k > 100.) 329 { 330 gamma = CalculatePolynomial(k, gamma100_20 331 // Only in this case it is not the exponen 332 } else 333 { 334 if (k > 10) 335 { 336 gamma = G4Exp(CalculatePolynomial(k, gam 337 } 338 else 339 { 340 gamma = G4Exp(CalculatePolynomial(k, gam 341 } 342 } 343 344 // ***** Original method 345 346 if (!fasterCode) 347 { 348 G4double oneOverMax = 1. 349 / (1. / (4. * gamma * gamma) 350 + beta / ((2. + 2. * delta) * (2. + 2. 351 352 G4double cosTheta = 0.; 353 G4double leftDenominator = 0.; 354 G4double rightDenominator = 0.; 355 G4double fCosTheta = 0.; 356 357 do 358 { 359 cosTheta = 2. * G4UniformRand()- 1.; 360 361 leftDenominator = (1. + 2.*gamma - cosTh 362 rightDenominator = (1. + 2.*delta + cosT 363 if ( (leftDenominator * rightDenominator 364 { 365 fCosTheta = oneOverMax * (1./(leftDeno 366 + beta/(righ 367 } 368 } 369 while (fCosTheta < G4UniformRand()); 370 371 return cosTheta; 372 } 373 374 // ***** Alternative method using cumulative 375 376 // if (fasterCode) 377 378 // 379 // modified by Shogo OKADA @ KEK, JP, 2016. 380 // 381 // An integral of differential cross-sectio 382 // (integral variable: cos(theta), integral 383 // 384 // 1.0 + x beta * ( 385 // I = --------------------- + ------------ 386 // (a - x) * (a + 1.0) (b + x) * 387 // 388 // where a = 1.0 + 2.0 * gamma(K), b = 1.0 389 // 390 // Then, a cumulative probability (cp) is a 391 // 392 // cp 1.0 + x beta * 393 // ---- = --------------------- + --------- 394 // S (a - x) * (a + 1.0) (b + x) 395 // 396 // where 1/S is the integral of differnetic 397 // 398 // 1 2.0 2. 399 // --- = ----------------------- + ------- 400 // S (a - 1.0) * (a + 1.0) (b + 1 401 // 402 // x is calculated from the quadratic equat 403 // 404 // A * x^2 + B * x + C = 0 405 // 406 // where A, B, anc C are coefficients of th 407 // A = S * {(b - 1.0) - beta * (a + 1.0)} 408 // B = S * {(b - 1.0) * (b + 1.0) + beta * 409 // C = S * {b * (b - 1.0) + beta * a * (a 410 // 411 412 // sampling cumulative probability 413 G4double cp = G4UniformRand(); 414 415 G4double a = 1.0 + 2.0 * gamma; 416 G4double b = 1.0 + 2.0 * delta; 417 G4double a1 = a - 1.0; 418 G4double a2 = a + 1.0; 419 G4double b1 = b - 1.0; 420 G4double b2 = b + 1.0; 421 G4double c1 = a - b; 422 G4double c2 = a * b; 423 424 G4double S = 2.0 / (a1 * a2) + 2.0 * beta / 425 426 // coefficients for the quadratic equation 427 G4double A = S * (b1 - beta * a2) + cp * a2 428 G4double B = S * (b1 * b2 + beta * a1 * a2) 429 G4double C = S * (b * b1 + beta * a * a2) - 430 431 // calculate cos(theta) 432 return (-1.0 * B + std::sqrt(B * B - 4.0 * 433 434 /* 435 G4double cosTheta = -1; 436 G4double cumul = 0; 437 G4double value = 0; 438 G4double leftDenominator = 0.; 439 G4double rightDenominator = 0.; 440 441 // Number of integration steps in the -1,1 442 G4int iMax=200; 443 444 G4double random = G4UniformRand(); 445 446 // Cumulate differential cross section 447 for (G4int i=0; i<iMax; i++) 448 { 449 cosTheta = -1 + i*2./(iMax-1); 450 leftDenominator = (1. + 2.*gamma - cosTheta 451 rightDenominator = (1. + 2.*delta + cosThet 452 if ( (leftDenominator * rightDenominator) ! 453 { 454 cumul = cumul + (1./(leftDenominator*leftDe 455 } 456 } 457 458 // Select cosTheta 459 for (G4int i=0; i<iMax; i++) 460 { 461 cosTheta = -1 + i*2./(iMax-1); 462 leftDenominator = (1. + 2.*gamma - cosTheta 463 rightDenominator = (1. + 2.*delta + cosThet 464 if (cumul !=0 && (leftDenominator * rightDe 465 value = value + (1./(leftDenominator*leftDe 466 if (random < value) break; 467 } 468 469 return cosTheta; 470 */ 471 472 473 //return 0.; 474 475 } 476 477 //....oooOO0OOooo........oooOO0OOooo........oo 478 479 G4double 480 G4DNAUeharaScreenedRutherfordElasticModel:: 481 CalculatePolynomial(G4double k, 482 std::vector<G4double>& vec 483 { 484 // Sum_{i=0}^{size-1} vector_i k^i 485 // 486 // Phys. Med. Biol. 29 N.4 (1983) 443-447 487 488 G4double result = 0.; 489 size_t size = vec.size(); 490 491 while (size > 0) 492 { 493 size--; 494 495 result *= k; 496 result += vec[size]; 497 } 498 499 return result; 500 } 501 502 //....oooOO0OOooo........oooOO0OOooo........oo 503 504 G4double 505 G4DNAUeharaScreenedRutherfordElasticModel:: 506 ScreenedRutherfordRandomizeCosTheta(G4double k 507 G4double z 508 { 509 510 // d sigma_el sigma_Ruth(K) 511 // ------------ (K) ~ ---------------------- 512 // d Omega (1 + 2 n(K) - cos(the 513 // 514 // We extract cos(theta) distributed as (1 + 515 // 516 // Maximum is for theta=0: 1/(4 n(K)^2) (Whe 517 // satisfied within the validity of the proc 518 // 519 // Phys. Med. Biol. 45 (2000) 3171-3194 520 521 // ***** Original method 522 523 if (!fasterCode) 524 { 525 G4double n = ScreeningFactor(k, z); 526 527 G4double oneOverMax = (4. * n * n); 528 529 G4double cosTheta = 0.; 530 G4double fCosTheta; 531 532 do 533 { 534 cosTheta = 2. * G4UniformRand()- 1.; 535 fCosTheta = (1 + 2.*n - cosTheta); 536 if (fCosTheta !=0.) fCosTheta = oneOverM 537 } 538 while (fCosTheta < G4UniformRand()); 539 540 return cosTheta; 541 } 542 543 // ***** Alternative method using cumulative 544 // if (fasterCode) 545 546 // 547 // modified by Shogo OKADA @ KEK, JP, 2016.2 548 // 549 // The cumulative probability (cp) is calcul 550 // the differential cross-section fomula wit 551 // 552 // n(K) * (1.0 + cos(theta)) 553 // cp = --------------------------------- 554 // 1.0 + 2.0 * n(K) - cos(theta) 555 // 556 // Then, cos(theta) is as follows: 557 // 558 // cp * (1.0 + 2.0 * n(K)) - n 559 // cos(theta) = ---------------------------- 560 // n(k) + cp 561 // 562 // where, K is kinetic energy, n(K) is scree 563 // 564 565 G4double n = ScreeningFactor(k, z); 566 G4double cp = G4UniformRand(); 567 G4double numerator = cp * (1.0 + 2.0 * n) - 568 G4double denominator = n + cp; 569 return numerator / denominator; 570 571 /* 572 G4double cosTheta = -1; 573 G4double cumul = 0; 574 G4double value = 0; 575 G4double n = ScreeningFactor(k, z); 576 G4double fCosTheta; 577 578 // Number of integration steps in the -1,1 579 G4int iMax=200; 580 581 G4double random = G4UniformRand(); 582 583 // Cumulate differential cross section 584 for (G4int i=0; i<iMax; i++) 585 { 586 cosTheta = -1 + i*2./(iMax-1); 587 fCosTheta = (1 + 2.*n - cosTheta); 588 if (fCosTheta !=0.) cumul = cumul + 1./(fCo 589 } 590 591 // Select cosTheta 592 for (G4int i=0; i<iMax; i++) 593 { 594 cosTheta = -1 + i*2./(iMax-1); 595 fCosTheta = (1 + 2.*n - cosTheta); 596 if (cumul !=0.) value = value + (1./(fCosTh 597 if (random < value) break; 598 } 599 return cosTheta; 600 */ 601 602 603 //return 0.; 604 } 605