| 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 //
27
28 #include "G4DNAScreenedRutherfordElasticModel.
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4DNAMolecularMaterial.hh"
32 #include "G4Exp.hh"
33 #include "G4Log.hh"
34
35 //....oooOO0OOooo........oooOO0OOooo........oo
36
37 using namespace std;
38
39 //....oooOO0OOooo........oooOO0OOooo........oo
40
41 G4DNAScreenedRutherfordElasticModel::
42 G4DNAScreenedRutherfordElasticModel(const G4Pa
43 const G4St
44 G4VEmModel(nam)
45 {
46 fpWaterDensity = nullptr;
47
48 lowEnergyLimit = 0 * eV;
49 intermediateEnergyLimit = 200 * eV; // Switc
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 o
61 // 4 = entering in methods
62
63 #ifdef SR_VERBOSE
64 if (verboseLevel > 0)
65 {
66 G4cout << "Screened Rutherford Elastic mod
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
78 fasterCode = false;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 G4DNAScreenedRutherfordElasticModel::~G4DNAScr
84 = default;
85
86 //....oooOO0OOooo........oooOO0OOooo........oo
87
88 void G4DNAScreenedRutherfordElasticModel::
89 Initialise(const G4ParticleDefinition* particl
90 const G4DataVector& /*cuts*/)
91 {
92 #ifdef SR_VERBOSE
93 if (verboseLevel > 3)
94 {
95 G4cout << "Calling G4DNAScreenedRutherford
96 << G4endl;
97 }
98 #endif
99
100 if(particle->GetParticleName() != "e-")
101 {
102 G4Exception ("*** WARNING: the G4DNAScreen
103 "intented to be used with ano
104 "",FatalException,"") ;
105 }
106
107 // Energy limits
108 if (LowEnergyLimit() < 9*eV)
109 {
110 G4Exception("*** WARNING: the G4DNAScreene
111 "not validated below 9 eV",
112 "",JustWarning,"") ;
113 }
114
115 if (HighEnergyLimit() > 1*MeV)
116 {
117 G4Exception("*** WARNING: the G4DNAScreene
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 mod
127 << "Energy range: "
128 << LowEnergyLimit() / eV << " eV -
129 << HighEnergyLimit() / MeV << " MeV
130 << G4endl;
131 }
132 #endif
133
134 if (isInitialised) { return; } // return her
135
136 // Initialize water density pointer
137 fpWaterDensity = G4DNAMolecularMaterial::Ins
138 GetNumMolPerVolTableFor(G4Ma
139
140 fParticleChangeForGamma = GetParticleChangeF
141 isInitialised = true;
142
143 // Constants for final state by Brenner & Za
144 // note: if called after if(isInitialised) n
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........oo
188
189 G4double G4DNAScreenedRutherfordElasticModel::
190 CrossSectionPerVolume(const G4Material* materi
191 #ifdef SR_VERBOSE
192 const G4ParticleDefiniti
193 #else
194 const G4ParticleDefiniti
195 #endif
196 G4double ekin,
197 G4double,
198 G4double)
199 {
200 #ifdef SR_VERBOSE
201 if (verboseLevel > 3)
202 {
203 G4cout << "Calling CrossSectionPerVolume()
204 "G4DNAScreenedRutherfordElasticMod
205 << G4endl;
206 }
207 #endif
208
209 // Calculate total cross section for model
210
211 G4double sigma=0.;
212 G4double waterDensity = (*fpWaterDensity)[ma
213
214 if(ekin <= HighEnergyLimit() && ekin >= LowE
215 {
216 G4double z = 10.;
217 G4double n = ScreeningFactor(ekin,z);
218 G4double crossSection = RutherfordCrossSec
219 sigma = pi * crossSection / (n * (n + 1.))
220 }
221
222 #ifdef SR_VERBOSE
223 if (verboseLevel > 2)
224 {
225 G4cout << "_______________________________
226 G4cout << "=== G4DNAScreenedRutherfordElas
227 << G4endl;
228 G4cout << "=== Kinetic energy(eV)=" << eki
229 << " particle : " << particleDefini
230 << G4endl;
231 G4cout << "=== Cross section per water mol
232 << G4endl;
233 G4cout << "=== Cross section per water mol
234 << sigma*waterDensity/(1./cm) << G4
235 G4cout << "=== G4DNAScreenedRutherfordElas
236 << G4endl;
237 }
238 #endif
239
240 return sigma*waterDensity;
241 }
242
243 //....oooOO0OOooo........oooOO0OOooo........oo
244
245 G4double G4DNAScreenedRutherfordElasticModel::
246
247 {
248 //
249 // e^4
250 // sigma_Ruth(K) = Z (Z+1) -----------------
251 // (4 pi epsilon_0)
252 //
253 // Where K is the electron non-relativistic
254 //
255 // NIM 155, pp. 145-156, 1978
256
257 G4double length = (e_squared * (k + electron
258 / (4 * pi * epsilon0 * k * (k + 2 * elec
259 G4double cross = z * (z + 1) * length * leng
260
261 return cross;
262 }
263
264 //....oooOO0OOooo........oooOO0OOooo........oo
265
266 G4double G4DNAScreenedRutherfordElasticModel::
267
268 {
269 //
270 // alpha_1 + beta_1 ln(K/eV) const
271 // n(T) = -------------------------- -------
272 // K/(m_e c^2) 2 + K
273 //
274 // Where K is the electron non-relativistic
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 * G4L
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 / de
294
295 return value;
296 }
297
298 //....oooOO0OOooo........oooOO0OOooo........oo
299
300 void G4DNAScreenedRutherfordElasticModel::
301 SampleSecondaries(std::vector<G4DynamicParticl
302 const G4MaterialCutsCouple*
303 const G4DynamicParticle* aDy
304 G4double,
305 G4double)
306 {
307 #ifdef SR_VERBOSE
308 if (verboseLevel > 3)
309 {
310 G4cout << "Calling SampleSecondaries() of
311 "G4DNAScreenedRutherfordElasticM
312 << G4endl;
313 }
314 #endif
315
316 G4double electronEnergy0 = aDynamicElectron-
317 G4double cosTheta = 0.;
318
319 if (electronEnergy0<intermediateEnergyLimit)
320 {
321 #ifdef SR_VERBOSE
322 if (verboseLevel > 3)
323 {G4cout << "---> Using Brenner & Zaider
324 #endif
325 cosTheta = BrennerZaiderRandomizeCosTheta(
326 }
327
328 if (electronEnergy0>=intermediateEnergyLimit
329 {
330 #ifdef SR_VERBOSE
331 if (verboseLevel > 3)
332 {G4cout << "---> Using Screened Rutherfo
333 #endif
334 G4double z = 10.;
335 cosTheta = ScreenedRutherfordRandomizeCosT
336 }
337
338 G4double phi = 2. * pi * G4UniformRand();
339
340 G4ThreeVector zVers = aDynamicElectron->GetM
341 G4ThreeVector xVers = zVers.orthogonal();
342 G4ThreeVector yVers = zVers.cross(xVers);
343
344 G4double xDir = std::sqrt(1. - cosTheta*cosT
345 G4double yDir = xDir;
346 xDir *= std::cos(phi);
347 yDir *= std::sin(phi);
348
349 G4ThreeVector zPrimeVers((xDir*xVers + yDir*
350
351 fParticleChangeForGamma->ProposeMomentumDire
352
353 fParticleChangeForGamma->SetProposedKineticE
354 }
355
356 //....oooOO0OOooo........oooOO0OOooo........oo
357
358 G4double G4DNAScreenedRutherfordElasticModel::
359 BrennerZaiderRandomizeCosTheta(G4double k)
360 {
361 // d sigma_el 1
362 // ------------ (K) ~ ----------------------
363 // d Omega (1 + 2 gamma(K) - cos
364 //
365 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/(
366 //
367 // Phys. Med. Biol. 29 N.4 (1983) 443-447
368
369 // gamma(K), beta(K) and delta(K) are polyno
370
371 k /= eV;
372
373 G4double beta = G4Exp(CalculatePolynomial(k,
374 G4double delta = G4Exp(CalculatePolynomial(k
375 G4double gamma;
376
377 if (k > 100.)
378 {
379 gamma = CalculatePolynomial(k, gamma100_20
380 // Only in this case it is not the exponen
381 }
382 else
383 {
384 if (k > 10)
385 {
386 gamma = G4Exp(CalculatePolynomial(k, gam
387 }
388 else
389 {
390 gamma = G4Exp(CalculatePolynomial(k, gam
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. * de
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 - cosThe
413 rightDenominator = (1. + 2.*delta + cosTh
414 if ( (leftDenominator * rightDenominator)
415 {
416 fCosTheta = oneOverMax * (1./(leftDenom
417 + beta/(right
418 }
419 }
420 while (fCosTheta < G4UniformRand());
421
422 return cosTheta;
423 }
424
425 // ***** Alternative method using cumulative
426
427 if (fasterCode)
428 {
429
430 //
431 // modified by Shogo OKADA @ KEK, JP, 2016.
432 //
433 // An integral of differential cross-sectio
434 // (integral variable: cos(theta), integral
435 //
436 // 1.0 + x beta * (
437 // I = --------------------- + ------------
438 // (a - x) * (a + 1.0) (b + x) *
439 //
440 // where a = 1.0 + 2.0 * gamma(K), b = 1.0
441 //
442 // Then, a cumulative probability (cp) is a
443 //
444 // cp 1.0 + x beta *
445 // ---- = --------------------- + ---------
446 // S (a - x) * (a + 1.0) (b + x)
447 //
448 // where 1/S is the integral of differnetic
449 //
450 // 1 2.0 2.
451 // --- = ----------------------- + -------
452 // S (a - 1.0) * (a + 1.0) (b + 1
453 //
454 // x is calculated from the quadratic equat
455 //
456 // A * x^2 + B * x + C = 0
457 //
458 // where A, B, anc C are coefficients of th
459 // A = S * {(b - 1.0) - beta * (a + 1.0)}
460 // B = S * {(b - 1.0) * (b + 1.0) + beta *
461 // C = S * {b * (b - 1.0) + beta * a * (a
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 /
477
478 // coefficients for the quadratic equation
479 G4double A = S * (b1 - beta * a2) + cp * a2
480 G4double B = S * (b1 * b2 + beta * a1 * a2)
481 G4double C = S * (b * b1 + beta * a * a2) -
482
483 // calculate cos(theta)
484 return (-1.0 * B + std::sqrt(B * B - 4.0 *
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
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 - cosThet
503 rightDenominator = (1. + 2.*delta + cosThe
504 if ( (leftDenominator * rightDenominator)
505 {
506 cumul = cumul + (1./(leftDenominator*lef
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 - cosThe
515 rightDenominator = (1. + 2.*delta + cosTh
516 if (cumul !=0 && (leftDenominator * right
517 value = value + (1./(leftDenominator*lef
518 if (random < value) break;
519 }
520
521 return cosTheta;
522 */
523 }
524
525 return 0.;
526 }
527
528 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
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(the
562 //
563 // We extract cos(theta) distributed as (1 +
564 //
565 // Maximum is for theta=0: 1/(4 n(K)^2) (Whe
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 = oneOverMa
586 }
587 while (fCosTheta < G4UniformRand());
588
589 return cosTheta;
590 }
591
592 // ***** Alternative method using cumulative
593
594
595 //
596 // modified by Shogo OKADA @ KEK, JP, 2016.
597 //
598 // The cumulative probability (cp) is calcu
599 // the differential cross-section fomula wi
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)) -
608 // cos(theta) = ---------------------------
609 // n(k) + cp
610 //
611 // where, K is kinetic energy, n(K) is scre
612 //
613
614 G4double n = ScreeningFactor(k, z);
615 G4double cp = G4UniformRand();
616 G4double numerator = cp * (1.0 + 2.0 * 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
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./(f
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./(fCos
646 if (random < value) break;
647 }
648 return cosTheta;
649 */
650
651
652 //return 0.;
653 }
654
655
656