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