| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27
28 #include "G4DNAScreenedRutherfordElasticModel.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4DNAMolecularMaterial.hh"
32 #include "G4Exp.hh"
33 #include "G4Log.hh"
34
35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36
37 using namespace std;
38
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40
41 G4DNAScreenedRutherfordElasticModel::
42 G4DNAScreenedRutherfordElasticModel(const G4ParticleDefinition*,
43 const G4String& nam) :
44 G4VEmModel(nam)
45 {
46 fpWaterDensity = nullptr;
47
48 lowEnergyLimit = 0 * eV;
49 intermediateEnergyLimit = 200 * eV; // Switch between two final state models
50 highEnergyLimit = 1. * MeV;
51
52 SetLowEnergyLimit(lowEnergyLimit);
53 SetHighEnergyLimit(highEnergyLimit);
54
55 verboseLevel = 0;
56 // Verbosity scale:
57 // 0 = nothing
58 // 1 = warning for energy non-conservation
59 // 2 = details of energy budget
60 // 3 = calculation of cross sections, file openings, sampling of atoms
61 // 4 = entering in methods
62
63 #ifdef SR_VERBOSE
64 if (verboseLevel > 0)
65 {
66 G4cout << "Screened Rutherford Elastic model is constructed "
67 << G4endl
68 << "Energy range: "
69 << lowEnergyLimit / eV << " eV - "
70 << highEnergyLimit / MeV << " MeV"
71 << G4endl;
72 }
73 #endif
74 fParticleChangeForGamma = nullptr;
75
76 // Selection of computation method
77 // We do not recommend "true" usage with the current cumul. proba. settings
78 fasterCode = false;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82
83 G4DNAScreenedRutherfordElasticModel::~G4DNAScreenedRutherfordElasticModel()
84 = default;
85
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87
88 void G4DNAScreenedRutherfordElasticModel::
89 Initialise(const G4ParticleDefinition* particle,
90 const G4DataVector& /*cuts*/)
91 {
92 #ifdef SR_VERBOSE
93 if (verboseLevel > 3)
94 {
95 G4cout << "Calling G4DNAScreenedRutherfordElasticModel::Initialise()"
96 << G4endl;
97 }
98 #endif
99
100 if(particle->GetParticleName() != "e-")
101 {
102 G4Exception ("*** WARNING: the G4DNAScreenedRutherfordElasticModel is not "
103 "intented to be used with another particle than the electron",
104 "",FatalException,"") ;
105 }
106
107 // Energy limits
108 if (LowEnergyLimit() < 9*eV)
109 {
110 G4Exception("*** WARNING: the G4DNAScreenedRutherfordElasticModel class is "
111 "not validated below 9 eV",
112 "",JustWarning,"") ;
113 }
114
115 if (HighEnergyLimit() > 1*MeV)
116 {
117 G4Exception("*** WARNING: the G4DNAScreenedRutherfordElasticModel class is "
118 "not validated above 1 MeV",
119 "",JustWarning,"") ;
120 }
121
122 //
123 #ifdef SR_VERBOSE
124 if( verboseLevel>0 )
125 {
126 G4cout << "Screened Rutherford elastic model is initialized " << G4endl
127 << "Energy range: "
128 << LowEnergyLimit() / eV << " eV - "
129 << HighEnergyLimit() / MeV << " MeV"
130 << G4endl;
131 }
132 #endif
133
134 if (isInitialised) { return; } // return here, prevent reinit consts + pointer
135
136 // Initialize water density pointer
137 fpWaterDensity = G4DNAMolecularMaterial::Instance()->
138 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
139
140 fParticleChangeForGamma = GetParticleChangeForGamma();
141 isInitialised = true;
142
143 // Constants for final state by Brenner & Zaider
144 // note: if called after if(isInitialised) no need for clear and resetting
145 // the values at every call
146
147 betaCoeff=
148 {
149 7.51525,
150 -0.41912,
151 7.2017E-3,
152 -4.646E-5,
153 1.02897E-7};
154
155 deltaCoeff=
156 {
157 2.9612,
158 -0.26376,
159 4.307E-3,
160 -2.6895E-5,
161 5.83505E-8};
162
163 gamma035_10Coeff =
164 {
165 -1.7013,
166 -1.48284,
167 0.6331,
168 -0.10911,
169 8.358E-3,
170 -2.388E-4};
171
172 gamma10_100Coeff =
173 {
174 -3.32517,
175 0.10996,
176 -4.5255E-3,
177 5.8372E-5,
178 -2.4659E-7};
179
180 gamma100_200Coeff =
181 {
182 2.4775E-2,
183 -2.96264E-5,
184 -1.20655E-7};
185 }
186
187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188
189 G4double G4DNAScreenedRutherfordElasticModel::
190 CrossSectionPerVolume(const G4Material* material,
191 #ifdef SR_VERBOSE
192 const G4ParticleDefinition* particleDefinition,
193 #else
194 const G4ParticleDefinition*,
195 #endif
196 G4double ekin,
197 G4double,
198 G4double)
199 {
200 #ifdef SR_VERBOSE
201 if (verboseLevel > 3)
202 {
203 G4cout << "Calling CrossSectionPerVolume() of "
204 "G4DNAScreenedRutherfordElasticModel"
205 << G4endl;
206 }
207 #endif
208
209 // Calculate total cross section for model
210
211 G4double sigma=0.;
212 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
213
214 if(ekin <= HighEnergyLimit() && ekin >= LowEnergyLimit())
215 {
216 G4double z = 10.;
217 G4double n = ScreeningFactor(ekin,z);
218 G4double crossSection = RutherfordCrossSection(ekin, z);
219 sigma = pi * crossSection / (n * (n + 1.));
220 }
221
222 #ifdef SR_VERBOSE
223 if (verboseLevel > 2)
224 {
225 G4cout << "__________________________________" << G4endl;
226 G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO START"
227 << G4endl;
228 G4cout << "=== Kinetic energy(eV)=" << ekin/eV
229 << " particle : " << particleDefinition->GetParticleName()
230 << G4endl;
231 G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm
232 << G4endl;
233 G4cout << "=== Cross section per water molecule (cm^-1)="
234 << sigma*waterDensity/(1./cm) << G4endl;
235 G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO END"
236 << G4endl;
237 }
238 #endif
239
240 return sigma*waterDensity;
241 }
242
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
244
245 G4double G4DNAScreenedRutherfordElasticModel::RutherfordCrossSection(G4double k,
246 G4double z)
247 {
248 //
249 // e^4 / K + m_e c^2 \^2
250 // sigma_Ruth(K) = Z (Z+1) -------------------- | --------------------- |
251 // (4 pi epsilon_0)^2 \ K * (K + 2 m_e c^2) /
252 //
253 // Where K is the electron non-relativistic kinetic energy
254 //
255 // NIM 155, pp. 145-156, 1978
256
257 G4double length = (e_squared * (k + electron_mass_c2))
258 / (4 * pi * epsilon0 * k * (k + 2 * electron_mass_c2));
259 G4double cross = z * (z + 1) * length * length;
260
261 return cross;
262 }
263
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
265
266 G4double G4DNAScreenedRutherfordElasticModel::ScreeningFactor(G4double k,
267 G4double z)
268 {
269 //
270 // alpha_1 + beta_1 ln(K/eV) constK Z^(2/3)
271 // n(T) = -------------------------- -----------------
272 // K/(m_e c^2) 2 + K/(m_e c^2)
273 //
274 // Where K is the electron non-relativistic kinetic energy
275 //
276 // n(T) > 0 for T < ~ 400 MeV
277 //
278 // NIM 155, pp. 145-156, 1978
279 // Formulae (2) and (5)
280
281 const G4double alpha_1(1.64);
282 const G4double beta_1(-0.0825);
283 const G4double constK(1.7E-5);
284
285 G4double numerator = (alpha_1 + beta_1 * G4Log(k / eV)) * constK
286 * std::pow(z, 2. / 3.);
287
288 k /= electron_mass_c2;
289
290 G4double denominator = k * (2 + k);
291
292 G4double value = 0.;
293 if (denominator > 0.) value = numerator / denominator;
294
295 return value;
296 }
297
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299
300 void G4DNAScreenedRutherfordElasticModel::
301 SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
302 const G4MaterialCutsCouple* /*couple*/,
303 const G4DynamicParticle* aDynamicElectron,
304 G4double,
305 G4double)
306 {
307 #ifdef SR_VERBOSE
308 if (verboseLevel > 3)
309 {
310 G4cout << "Calling SampleSecondaries() of "
311 "G4DNAScreenedRutherfordElasticModel"
312 << G4endl;
313 }
314 #endif
315
316 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
317 G4double cosTheta = 0.;
318
319 if (electronEnergy0<intermediateEnergyLimit)
320 {
321 #ifdef SR_VERBOSE
322 if (verboseLevel > 3)
323 {G4cout << "---> Using Brenner & Zaider model" << G4endl;}
324 #endif
325 cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
326 }
327
328 if (electronEnergy0>=intermediateEnergyLimit)
329 {
330 #ifdef SR_VERBOSE
331 if (verboseLevel > 3)
332 {G4cout << "---> Using Screened Rutherford model" << G4endl;}
333 #endif
334 G4double z = 10.;
335 cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
336 }
337
338 G4double phi = 2. * pi * G4UniformRand();
339
340 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
341 G4ThreeVector xVers = zVers.orthogonal();
342 G4ThreeVector yVers = zVers.cross(xVers);
343
344 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
345 G4double yDir = xDir;
346 xDir *= std::cos(phi);
347 yDir *= std::sin(phi);
348
349 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
350
351 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
352
353 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
354 }
355
356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
357
358 G4double G4DNAScreenedRutherfordElasticModel::
359 BrennerZaiderRandomizeCosTheta(G4double k)
360 {
361 // d sigma_el 1 beta(K)
362 // ------------ (K) ~ --------------------------------- + ---------------------------------
363 // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2
364 //
365 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
366 //
367 // Phys. Med. Biol. 29 N.4 (1983) 443-447
368
369 // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
370
371 k /= eV;
372
373 G4double beta = G4Exp(CalculatePolynomial(k, betaCoeff));
374 G4double delta = G4Exp(CalculatePolynomial(k, deltaCoeff));
375 G4double gamma;
376
377 if (k > 100.)
378 {
379 gamma = CalculatePolynomial(k, gamma100_200Coeff);
380 // Only in this case it is not the exponent of the polynomial
381 }
382 else
383 {
384 if (k > 10)
385 {
386 gamma = G4Exp(CalculatePolynomial(k, gamma10_100Coeff));
387 }
388 else
389 {
390 gamma = G4Exp(CalculatePolynomial(k, gamma035_10Coeff));
391 }
392 }
393
394 // ***** Original method
395
396 if (!fasterCode)
397 {
398
399 G4double oneOverMax = 1.
400 / (1. / (4. * gamma * gamma) + beta
401 / ((2. + 2. * delta) * (2. + 2. * delta)));
402
403 G4double cosTheta = 0.;
404 G4double leftDenominator = 0.;
405 G4double rightDenominator = 0.;
406 G4double fCosTheta = 0.;
407
408 do
409 {
410 cosTheta = 2. * G4UniformRand()- 1.;
411
412 leftDenominator = (1. + 2.*gamma - cosTheta);
413 rightDenominator = (1. + 2.*delta + cosTheta);
414 if ( (leftDenominator * rightDenominator) != 0. )
415 {
416 fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator)
417 + beta/(rightDenominator*rightDenominator));
418 }
419 }
420 while (fCosTheta < G4UniformRand());
421
422 return cosTheta;
423 }
424
425 // ***** Alternative method using cumulative probability
426
427 if (fasterCode)
428 {
429
430 //
431 // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
432 //
433 // An integral of differential cross-section formula shown above this member function
434 // (integral variable: cos(theta), integral interval: [-1, x]) is as follows:
435 //
436 // 1.0 + x beta * (1 + x)
437 // I = --------------------- + ---------------------- (1)
438 // (a - x) * (a + 1.0) (b + x) * (b - 1.0)
439 //
440 // where a = 1.0 + 2.0 * gamma(K), b = 1.0 + 2.0 * delta(K)
441 //
442 // Then, a cumulative probability (cp) is as follows:
443 //
444 // cp 1.0 + x beta * (1 + x)
445 // ---- = --------------------- + ---------------------- (2)
446 // S (a - x) * (a + 1.0) (b + x) * (b - 1.0)
447 //
448 // where 1/S is the integral of differnetical cross-section (1) on interval [-1, 1]
449 //
450 // 1 2.0 2.0 * beta
451 // --- = ----------------------- + ----------------------- (3)
452 // S (a - 1.0) * (a + 1.0) (b + 1.0) * (b - 1.0)
453 //
454 // x is calculated from the quadratic equation derived from (2) and (3):
455 //
456 // A * x^2 + B * x + C = 0
457 //
458 // where A, B, anc C are coefficients of the equation:
459 // A = S * {(b - 1.0) - beta * (a + 1.0)} + cp * (a + 1.0) * (b - 1.0),
460 // B = S * {(b - 1.0) * (b + 1.0) + beta * (a - 1.0) * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * (a - b)
461 // C = S * {b * (b - 1.0) + beta * a * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * ab
462 //
463
464 // sampling cumulative probability
465 G4double cp = G4UniformRand();
466
467 G4double a = 1.0 + 2.0 * gamma;
468 G4double b = 1.0 + 2.0 * delta;
469 G4double a1 = a - 1.0;
470 G4double a2 = a + 1.0;
471 G4double b1 = b - 1.0;
472 G4double b2 = b + 1.0;
473 G4double c1 = a - b;
474 G4double c2 = a * b;
475
476 G4double S = 2.0 / (a1 * a2) + 2.0 * beta / (b1 * b2); S = 1.0 / S;
477
478 // coefficients for the quadratic equation
479 G4double A = S * (b1 - beta * a2) + cp * a2 * b1;
480 G4double B = S * (b1 * b2 + beta * a1 * a2) - cp * a2 * b1 * c1;
481 G4double C = S * (b * b1 + beta * a * a2) - cp * a2 * b1 * c2;
482
483 // calculate cos(theta)
484 return (-1.0 * B + std::sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
485
486 /*
487 G4double cosTheta = -1;
488 G4double cumul = 0;
489 G4double value = 0;
490 G4double leftDenominator = 0.;
491 G4double rightDenominator = 0.;
492
493 // Number of integration steps in the -1,1 range
494 G4int iMax=200;
495
496 G4double random = G4UniformRand();
497
498 // Cumulate differential cross section
499 for (G4int i=0; i<iMax; i++)
500 {
501 cosTheta = -1 + i*2./(iMax-1);
502 leftDenominator = (1. + 2.*gamma - cosTheta);
503 rightDenominator = (1. + 2.*delta + cosTheta);
504 if ( (leftDenominator * rightDenominator) != 0. )
505 {
506 cumul = cumul + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
507 }
508 }
509
510 // Select cosTheta
511 for (G4int i=0; i<iMax; i++)
512 {
513 cosTheta = -1 + i*2./(iMax-1);
514 leftDenominator = (1. + 2.*gamma - cosTheta);
515 rightDenominator = (1. + 2.*delta + cosTheta);
516 if (cumul !=0 && (leftDenominator * rightDenominator) != 0.)
517 value = value + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)) / cumul;
518 if (random < value) break;
519 }
520
521 return cosTheta;
522 */
523 }
524
525 return 0.;
526 }
527
528 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
529
530 G4double G4DNAScreenedRutherfordElasticModel::
531 CalculatePolynomial(G4double k,
532 std::vector<G4double>& vec)
533 {
534 // Sum_{i=0}^{size-1} vector_i k^i
535 //
536 // Phys. Med. Biol. 29 N.4 (1983) 443-447
537
538 G4double result = 0.;
539 size_t size = vec.size();
540
541 while (size > 0)
542 {
543 size--;
544
545 result *= k;
546 result += vec[size];
547 }
548
549 return result;
550 }
551
552 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553
554 G4double G4DNAScreenedRutherfordElasticModel::
555 ScreenedRutherfordRandomizeCosTheta(G4double k,
556 G4double z)
557 {
558
559 // d sigma_el sigma_Ruth(K)
560 // ------------ (K) ~ -----------------------------
561 // d Omega (1 + 2 n(K) - cos(theta))^2
562 //
563 // We extract cos(theta) distributed as (1 + 2 n(K) - cos(theta))^-2
564 //
565 // Maximum is for theta=0: 1/(4 n(K)^2) (When n(K) is positive, that is always satisfied within the validity of the process)
566 //
567 // Phys. Med. Biol. 45 (2000) 3171-3194
568
569 // ***** Original method
570
571 if (!fasterCode)
572 {
573
574 G4double n = ScreeningFactor(k, z);
575
576 G4double oneOverMax = (4. * n * n);
577
578 G4double cosTheta = 0.;
579 G4double fCosTheta;
580
581 do
582 {
583 cosTheta = 2. * G4UniformRand()- 1.;
584 fCosTheta = (1 + 2.*n - cosTheta);
585 if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta);
586 }
587 while (fCosTheta < G4UniformRand());
588
589 return cosTheta;
590 }
591
592 // ***** Alternative method using cumulative probability
593
594
595 //
596 // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
597 //
598 // The cumulative probability (cp) is calculated by integrating
599 // the differential cross-section fomula with cos(theta):
600 //
601 // n(K) * (1.0 + cos(theta))
602 // cp = ---------------------------------
603 // 1.0 + 2.0 * n(K) - cos(theta)
604 //
605 // Then, cos(theta) is as follows:
606 //
607 // cp * (1.0 + 2.0 * n(K)) - n(K)
608 // cos(theta) = --------------------------------
609 // n(k) + cp
610 //
611 // where, K is kinetic energy, n(K) is screeing factor, and cp is cumulative probability
612 //
613
614 G4double n = ScreeningFactor(k, z);
615 G4double cp = G4UniformRand();
616 G4double numerator = cp * (1.0 + 2.0 * n) - n;
617 G4double denominator = n + cp;
618 return numerator / denominator;
619
620 /*
621 G4double cosTheta = -1;
622 G4double cumul = 0;
623 G4double value = 0;
624 G4double n = ScreeningFactor(k, z);
625 G4double fCosTheta;
626
627 // Number of integration steps in the -1,1 range
628 G4int iMax=200;
629
630 G4double random = G4UniformRand();
631
632 // Cumulate differential cross section
633 for (G4int i=0; i<iMax; i++)
634 {
635 cosTheta = -1 + i*2./(iMax-1);
636 fCosTheta = (1 + 2.*n - cosTheta);
637 if (fCosTheta !=0.) cumul = cumul + 1./(fCosTheta*fCosTheta);
638 }
639
640 // Select cosTheta
641 for (G4int i=0; i<iMax; i++)
642 {
643 cosTheta = -1 + i*2./(iMax-1);
644 fCosTheta = (1 + 2.*n - cosTheta);
645 if (cumul !=0.) value = value + (1./(fCosTheta*fCosTheta)) / cumul;
646 if (random < value) break;
647 }
648 return cosTheta;
649 */
650
651
652 //return 0.;
653 }
654
655
656