ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAScreenedRutherfordElasticModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAScreenedRutherfordElasticModel.cc
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 
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.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 
43  const G4String& nam) :
44 G4VEmModel(nam), isInitialised(false)
45 {
46  fpWaterDensity = 0;
47 
48  lowEnergyLimit = 0 * eV;
49  intermediateEnergyLimit = 200 * eV; // Switch between two final state models
50  highEnergyLimit = 1. * MeV;
51 
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
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 
84 {
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
91  const G4DataVector& /*cuts*/)
92 {
93 #ifdef SR_VERBOSE
94  if (verboseLevel > 3)
95  {
96  G4cout << "Calling G4DNAScreenedRutherfordElasticModel::Initialise()"
97  << G4endl;
98  }
99 #endif
100 
101  if(particle->GetParticleName() != "e-")
102  {
103  G4Exception ("*** WARNING: the G4DNAScreenedRutherfordElasticModel is not "
104  "intented to be used with another particle than the electron",
105  "",FatalException,"") ;
106  }
107 
108  // Energy limits
109  if (LowEnergyLimit() < 9*eV)
110  {
111  G4Exception("*** WARNING: the G4DNAScreenedRutherfordElasticModel class is "
112  "not validated below 9 eV",
113  "",JustWarning,"") ;
114  }
115 
116  if (HighEnergyLimit() > 1*MeV)
117  {
118  G4Exception("*** WARNING: the G4DNAScreenedRutherfordElasticModel class is "
119  "not validated above 1 MeV",
120  "",JustWarning,"") ;
121  }
122 
123  //
124 #ifdef SR_VERBOSE
125  if( verboseLevel>0 )
126  {
127  G4cout << "Screened Rutherford elastic model is initialized " << G4endl
128  << "Energy range: "
129  << LowEnergyLimit() / eV << " eV - "
130  << HighEnergyLimit() / MeV << " MeV"
131  << G4endl;
132  }
133 #endif
134 
135  if (isInitialised) { return; } // return here, prevent reinit consts + pointer
136 
137  // Initialize water density pointer
139  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
140 
142  isInitialised = true;
143 
144  // Constants for final state by Brenner & Zaider
145  // note: if called after if(isInitialised) no need for clear and resetting
146  // the values at every call
147 
148  betaCoeff=
149  {
150  7.51525,
151  -0.41912,
152  7.2017E-3,
153  -4.646E-5,
154  1.02897E-7};
155 
156  deltaCoeff=
157  {
158  2.9612,
159  -0.26376,
160  4.307E-3,
161  -2.6895E-5,
162  5.83505E-8};
163 
165  {
166  -1.7013,
167  -1.48284,
168  0.6331,
169  -0.10911,
170  8.358E-3,
171  -2.388E-4};
172 
174  {
175  -3.32517,
176  0.10996,
177  -4.5255E-3,
178  5.8372E-5,
179  -2.4659E-7};
180 
182  {
183  2.4775E-2,
184  -2.96264E-5,
185  -1.20655E-7};
186 }
187 
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189 
192 #ifdef SR_VERBOSE
193  const G4ParticleDefinition* particleDefinition,
194 #else
195  const G4ParticleDefinition*,
196 #endif
197  G4double ekin,
198  G4double,
199  G4double)
200 {
201 #ifdef SR_VERBOSE
202  if (verboseLevel > 3)
203  {
204  G4cout << "Calling CrossSectionPerVolume() of "
205  "G4DNAScreenedRutherfordElasticModel"
206  << G4endl;
207  }
208 #endif
209 
210  // Calculate total cross section for model
211 
212  G4double sigma=0.;
213  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
214 
215  if(ekin <= HighEnergyLimit() && ekin >= LowEnergyLimit())
216  {
217  G4double z = 10.;
218  G4double n = ScreeningFactor(ekin,z);
219  G4double crossSection = RutherfordCrossSection(ekin, z);
220  sigma = pi * crossSection / (n * (n + 1.));
221  }
222 
223 #ifdef SR_VERBOSE
224  if (verboseLevel > 2)
225  {
226  G4cout << "__________________________________" << G4endl;
227  G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO START"
228  << G4endl;
229  G4cout << "=== Kinetic energy(eV)=" << ekin/eV
230  << " particle : " << particleDefinition->GetParticleName()
231  << G4endl;
232  G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm
233  << G4endl;
234  G4cout << "=== Cross section per water molecule (cm^-1)="
235  << sigma*waterDensity/(1./cm) << G4endl;
236  G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO END"
237  << G4endl;
238  }
239 #endif
240 
241  return sigma*waterDensity;
242 }
243 
244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
245 
247  G4double z)
248 {
249  //
250  // e^4 / K + m_e c^2 \^2
251  // sigma_Ruth(K) = Z (Z+1) -------------------- | --------------------- |
252  // (4 pi epsilon_0)^2 \ K * (K + 2 m_e c^2) /
253  //
254  // Where K is the electron non-relativistic kinetic energy
255  //
256  // NIM 155, pp. 145-156, 1978
257 
259  / (4 * pi * epsilon0 * k * (k + 2 * electron_mass_c2));
260  G4double cross = z * (z + 1) * length * length;
261 
262  return cross;
263 }
264 
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266 
268  G4double z)
269 {
270  //
271  // alpha_1 + beta_1 ln(K/eV) constK Z^(2/3)
272  // n(T) = -------------------------- -----------------
273  // K/(m_e c^2) 2 + K/(m_e c^2)
274  //
275  // Where K is the electron non-relativistic kinetic energy
276  //
277  // n(T) > 0 for T < ~ 400 MeV
278  //
279  // NIM 155, pp. 145-156, 1978
280  // Formulae (2) and (5)
281 
282  const G4double alpha_1(1.64);
283  const G4double beta_1(-0.0825);
284  const G4double constK(1.7E-5);
285 
286  G4double numerator = (alpha_1 + beta_1 * G4Log(k / eV)) * constK
287  * std::pow(z, 2. / 3.);
288 
289  k /= electron_mass_c2;
290 
291  G4double denominator = k * (2 + k);
292 
293  G4double value = 0.;
294  if (denominator > 0.) value = numerator / denominator;
295 
296  return value;
297 }
298 
299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
300 
302 SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
303  const G4MaterialCutsCouple* /*couple*/,
304  const G4DynamicParticle* aDynamicElectron,
305  G4double,
306  G4double)
307 {
308 #ifdef SR_VERBOSE
309  if (verboseLevel > 3)
310  {
311  G4cout << "Calling SampleSecondaries() of "
312  "G4DNAScreenedRutherfordElasticModel"
313  << G4endl;
314  }
315 #endif
316 
317  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
318  G4double cosTheta = 0.;
319 
320  if (electronEnergy0<intermediateEnergyLimit)
321  {
322 #ifdef SR_VERBOSE
323  if (verboseLevel > 3)
324  {G4cout << "---> Using Brenner & Zaider model" << G4endl;}
325 #endif
326  cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
327  }
328 
329  if (electronEnergy0>=intermediateEnergyLimit)
330  {
331 #ifdef SR_VERBOSE
332  if (verboseLevel > 3)
333  {G4cout << "---> Using Screened Rutherford model" << G4endl;}
334 #endif
335  G4double z = 10.;
336  cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
337  }
338 
339  G4double phi = 2. * pi * G4UniformRand();
340 
341  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
342  G4ThreeVector xVers = zVers.orthogonal();
343  G4ThreeVector yVers = zVers.cross(xVers);
344 
345  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
346  G4double yDir = xDir;
347  xDir *= std::cos(phi);
348  yDir *= std::sin(phi);
349 
350  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
351 
353 
355 }
356 
357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
358 
361 {
362  // d sigma_el 1 beta(K)
363  // ------------ (K) ~ --------------------------------- + ---------------------------------
364  // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2
365  //
366  // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
367  //
368  // Phys. Med. Biol. 29 N.4 (1983) 443-447
369 
370  // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
371 
372  k /= eV;
373 
376  G4double gamma;
377 
378  if (k > 100.)
379  {
381  // Only in this case it is not the exponent of the polynomial
382  }
383  else
384  {
385  if (k > 10)
386  {
388  }
389  else
390  {
392  }
393  }
394 
395  // ***** Original method
396 
397  if (!fasterCode)
398  {
399 
400  G4double oneOverMax = 1.
401  / (1. / (4. * gamma * gamma) + beta
402  / ((2. + 2. * delta) * (2. + 2. * delta)));
403 
404  G4double cosTheta = 0.;
405  G4double leftDenominator = 0.;
406  G4double rightDenominator = 0.;
407  G4double fCosTheta = 0.;
408 
409  do
410  {
411  cosTheta = 2. * G4UniformRand()- 1.;
412 
413  leftDenominator = (1. + 2.*gamma - cosTheta);
414  rightDenominator = (1. + 2.*delta + cosTheta);
415  if ( (leftDenominator * rightDenominator) != 0. )
416  {
417  fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator)
418  + beta/(rightDenominator*rightDenominator));
419  }
420  }
421  while (fCosTheta < G4UniformRand());
422 
423  return cosTheta;
424  }
425 
426  // ***** Alternative method using cumulative probability
427 
428  if (fasterCode)
429  {
430 
431  //
432  // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
433  //
434  // An integral of differential cross-section formula shown above this member function
435  // (integral variable: cos(theta), integral interval: [-1, x]) is as follows:
436  //
437  // 1.0 + x beta * (1 + x)
438  // I = --------------------- + ---------------------- (1)
439  // (a - x) * (a + 1.0) (b + x) * (b - 1.0)
440  //
441  // where a = 1.0 + 2.0 * gamma(K), b = 1.0 + 2.0 * delta(K)
442  //
443  // Then, a cumulative probability (cp) is as follows:
444  //
445  // cp 1.0 + x beta * (1 + x)
446  // ---- = --------------------- + ---------------------- (2)
447  // S (a - x) * (a + 1.0) (b + x) * (b - 1.0)
448  //
449  // where 1/S is the integral of differnetical cross-section (1) on interval [-1, 1]
450  //
451  // 1 2.0 2.0 * beta
452  // --- = ----------------------- + ----------------------- (3)
453  // S (a - 1.0) * (a + 1.0) (b + 1.0) * (b - 1.0)
454  //
455  // x is calculated from the quadratic equation derived from (2) and (3):
456  //
457  // A * x^2 + B * x + C = 0
458  //
459  // where A, B, anc C are coefficients of the equation:
460  // A = S * {(b - 1.0) - beta * (a + 1.0)} + cp * (a + 1.0) * (b - 1.0),
461  // B = S * {(b - 1.0) * (b + 1.0) + beta * (a - 1.0) * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * (a - b)
462  // C = S * {b * (b - 1.0) + beta * a * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * ab
463  //
464 
465  // sampling cumulative probability
467 
468  G4double a = 1.0 + 2.0 * gamma;
469  G4double b = 1.0 + 2.0 * delta;
470  G4double a1 = a - 1.0;
471  G4double a2 = a + 1.0;
472  G4double b1 = b - 1.0;
473  G4double b2 = b + 1.0;
474  G4double c1 = a - b;
475  G4double c2 = a * b;
476 
477  G4double S = 2.0 / (a1 * a2) + 2.0 * beta / (b1 * b2); S = 1.0 / S;
478 
479  // coefficients for the quadratic equation
480  G4double A = S * (b1 - beta * a2) + cp * a2 * b1;
481  G4double B = S * (b1 * b2 + beta * a1 * a2) - cp * a2 * b1 * c1;
482  G4double C = S * (b * b1 + beta * a * a2) - cp * a2 * b1 * c2;
483 
484  // calculate cos(theta)
485  return (-1.0 * B + std::sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
486 
487  /*
488  G4double cosTheta = -1;
489  G4double cumul = 0;
490  G4double value = 0;
491  G4double leftDenominator = 0.;
492  G4double rightDenominator = 0.;
493 
494  // Number of integration steps in the -1,1 range
495  G4int iMax=200;
496 
497  G4double random = G4UniformRand();
498 
499  // Cumulate differential cross section
500  for (G4int i=0; i<iMax; i++)
501  {
502  cosTheta = -1 + i*2./(iMax-1);
503  leftDenominator = (1. + 2.*gamma - cosTheta);
504  rightDenominator = (1. + 2.*delta + cosTheta);
505  if ( (leftDenominator * rightDenominator) != 0. )
506  {
507  cumul = cumul + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
508  }
509  }
510 
511  // Select cosTheta
512  for (G4int i=0; i<iMax; i++)
513  {
514  cosTheta = -1 + i*2./(iMax-1);
515  leftDenominator = (1. + 2.*gamma - cosTheta);
516  rightDenominator = (1. + 2.*delta + cosTheta);
517  if (cumul !=0 && (leftDenominator * rightDenominator) != 0.)
518  value = value + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)) / cumul;
519  if (random < value) break;
520  }
521 
522  return cosTheta;
523  */
524  }
525 
526  return 0.;
527 }
528 
529 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
530 
533  std::vector<G4double>& vec)
534 {
535  // Sum_{i=0}^{size-1} vector_i k^i
536  //
537  // Phys. Med. Biol. 29 N.4 (1983) 443-447
538 
539  G4double result = 0.;
540  size_t size = vec.size();
541 
542  while (size > 0)
543  {
544  size--;
545 
546  result *= k;
547  result += vec[size];
548  }
549 
550  return result;
551 }
552 
553 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
554 
557  G4double z)
558 {
559 
560  // d sigma_el sigma_Ruth(K)
561  // ------------ (K) ~ -----------------------------
562  // d Omega (1 + 2 n(K) - cos(theta))^2
563  //
564  // We extract cos(theta) distributed as (1 + 2 n(K) - cos(theta))^-2
565  //
566  // 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)
567  //
568  // Phys. Med. Biol. 45 (2000) 3171-3194
569 
570  // ***** Original method
571 
572  if (!fasterCode)
573  {
574 
575  G4double n = ScreeningFactor(k, z);
576 
577  G4double oneOverMax = (4. * n * n);
578 
579  G4double cosTheta = 0.;
580  G4double fCosTheta;
581 
582  do
583  {
584  cosTheta = 2. * G4UniformRand()- 1.;
585  fCosTheta = (1 + 2.*n - cosTheta);
586  if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta);
587  }
588  while (fCosTheta < G4UniformRand());
589 
590  return cosTheta;
591  }
592 
593  // ***** Alternative method using cumulative probability
594  if (fasterCode)
595  {
596 
597  //
598  // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
599  //
600  // The cumulative probability (cp) is calculated by integrating
601  // the differential cross-section fomula with cos(theta):
602  //
603  // n(K) * (1.0 + cos(theta))
604  // cp = ---------------------------------
605  // 1.0 + 2.0 * n(K) - cos(theta)
606  //
607  // Then, cos(theta) is as follows:
608  //
609  // cp * (1.0 + 2.0 * n(K)) - n(K)
610  // cos(theta) = --------------------------------
611  // n(k) + cp
612  //
613  // where, K is kinetic energy, n(K) is screeing factor, and cp is cumulative probability
614  //
615 
616  G4double n = ScreeningFactor(k, z);
618  G4double numerator = cp * (1.0 + 2.0 * n) - n;
619  G4double denominator = n + cp;
620  return numerator / denominator;
621 
622  /*
623  G4double cosTheta = -1;
624  G4double cumul = 0;
625  G4double value = 0;
626  G4double n = ScreeningFactor(k, z);
627  G4double fCosTheta;
628 
629  // Number of integration steps in the -1,1 range
630  G4int iMax=200;
631 
632  G4double random = G4UniformRand();
633 
634  // Cumulate differential cross section
635  for (G4int i=0; i<iMax; i++)
636  {
637  cosTheta = -1 + i*2./(iMax-1);
638  fCosTheta = (1 + 2.*n - cosTheta);
639  if (fCosTheta !=0.) cumul = cumul + 1./(fCosTheta*fCosTheta);
640  }
641 
642  // Select cosTheta
643  for (G4int i=0; i<iMax; i++)
644  {
645  cosTheta = -1 + i*2./(iMax-1);
646  fCosTheta = (1 + 2.*n - cosTheta);
647  if (cumul !=0.) value = value + (1./(fCosTheta*fCosTheta)) / cumul;
648  if (random < value) break;
649  }
650  return cosTheta;
651  */
652  }
653 
654  return 0.;
655 }
656 
657