| 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