Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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  intermediateEnergyLimit = 200. * eV; // Switch between two final state models
47 
49 // SetHighEnergyLimit(10.*keV);
51 
52  verboseLevel = 0;
53  // Verbosity scale:
54  // 0 = nothing
55  // 1 = warning for energy non-conservation
56  // 2 = details of energy budget
57  // 3 = calculation of cross sections, file openings, sampling of atoms
58  // 4 = entering in methods
59 
60 #ifdef UEHARA_VERBOSE
61  if (verboseLevel)
62  {
63  G4cout << "Screened Rutherford Elastic model is constructed " << G4endl
64  << "Energy range: "
65  << LowEnergyLimit() / eV << " eV - "
66  << HighEnergyLimit() / MeV << " MeV"
67  << G4endl;
68  }
69 #endif
70 
72 
73  // Selection of computation method
74  // We do not recommend "true" usage with the current cumul. proba. settings
75  fasterCode = false;
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 
82 {
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86 
87 void
90  const G4DataVector& /*cuts*/)
91 {
92 #ifdef UEHARA_VERBOSE
93  if (verboseLevel > 3)
94  {
95  G4cout << "Calling G4DNAUeharaScreenedRutherfordElasticModel::Initialise()"
96  << G4endl;
97  }
98 #endif
99 
100  if(particle->GetParticleName() != "e-")
101  {
102  G4Exception("*** WARNING: the G4DNAUeharaScreenedRutherfordElasticModel is "
103  "not intented to be used with another particle than the electron",
104  "",FatalException,"") ;
105  }
106 
107  // Energy limits
108  if(LowEnergyLimit() < 9.*CLHEP::eV)
109  {
110  G4Exception("*** WARNING : the G4DNAUeharaScreenedRutherfordElasticModel "
111  "class is not validated below 9 eV",
112  "",JustWarning,"") ;
113  }
114 
115  if (HighEnergyLimit() > 10.*CLHEP::keV)
116  {
117  G4Exception("*** WARNING: the G4DNAUeharaScreenedRutherfordElasticModel "
118  "class is used above 10 keV",
119  "",JustWarning,"") ;
120  }
121 
122 #ifdef UEHARA_VERBOSE
123  if( verboseLevel>0 )
124  {
125  G4cout << "Screened Rutherford elastic model is initialized " << G4endl
126  << "Energy range: "
127  << LowEnergyLimit() / eV << " eV - "
128  << HighEnergyLimit() / MeV << " MeV"
129  << G4endl;
130  }
131 #endif
132 
133  if (isInitialised){ return; }
134 
135  // Constants for final state by Brenner & Zaider
136  // Note: the instantiation must be placed after if (isInitialised)
137 
138  betaCoeff=
139  {
140  7.51525,
141  -0.41912,
142  7.2017E-3,
143  -4.646E-5,
144  1.02897E-7};
145 
146  deltaCoeff=
147  {
148  2.9612,
149  -0.26376,
150  4.307E-3,
151  -2.6895E-5,
152  5.83505E-8};
153 
155  {
156  -1.7013,
157  -1.48284,
158  0.6331,
159  -0.10911,
160  8.358E-3,
161  -2.388E-4};
162 
164  {
165  -3.32517,
166  0.10996,
167  -4.5255E-3,
168  5.8372E-5,
169  -2.4659E-7};
170 
172  {
173  2.4775E-2,
174  -2.96264E-5,
175  -1.20655E-7};
176 
177  // Initialize water density pointer
180  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
181 
183  isInitialised = true;
184 }
185 
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187 
188 G4double
191  const G4ParticleDefinition* /*particleDefinition*/,
192  G4double ekin,
193  G4double,
194  G4double)
195 {
196 #ifdef UEHARA_VERBOSE
197  if (verboseLevel > 3)
198  {
199  G4cout
200  << "Calling CrossSectionPerVolume() of G4DNAUeharaScreenedRutherfordElasticModel"
201  << G4endl;
202  }
203 #endif
204 
205  // Calculate total cross section for model
206 
207  G4double sigma = 0.;
208  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
209 
210  // Use activation boundaries instead
211  //if(ekin > HighEnergyLimit() || ekin < LowEnergyLimit()) return 0.;
212 
213  if(waterDensity!= 0.0)
214  {
215  G4double z = 7.42; // FROM PMB 37 (1992) 1841-1858 p1842
216  G4double n = ScreeningFactor(ekin,z);
217  G4double crossSection = RutherfordCrossSection(ekin, z);
218  sigma = pi * crossSection / (n * (n + 1.));
219 
220 #ifdef UEHARA_VERBOSE
221  if (verboseLevel > 2)
222  {
223  G4cout << "__________________________________" << G4endl;
224  G4cout << "=== G4DNAUeharaScreenedRutherfordElasticModel - XS INFO START"
225  << G4endl;
226  G4cout << "=== Kinetic energy(eV)=" << ekin/eV
227  << " particle : " << particleDefinition->GetParticleName() << G4endl;
228  G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm
229  << G4endl;
230  G4cout << "=== Cross section per water molecule (cm^-1)="
231  << sigma*waterDensity/(1./cm) << G4endl;
232  G4cout << "=== G4DNAUeharaScreenedRutherfordElasticModel - XS INFO END"
233  << G4endl;
234  }
235 #endif
236  }
237 
238  return sigma*waterDensity;
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 
243 G4double
245  G4double z)
246 {
247  //
248  // e^4 / K + m_e c^2 \^2
249  // sigma_Ruth(K) = Z (Z+1) -------------------- | --------------------- |
250  // (4 pi epsilon_0)^2 \ K * (K + 2 m_e c^2) /
251  //
252  // Where K is the electron non-relativistic kinetic energy
253  //
254  // NIM 155, pp. 145-156, 1978
255 
256  G4double length = (e_squared * (k + electron_mass_c2))
257  / (4 * pi * epsilon0 * k * (k + 2 * electron_mass_c2));
258  G4double cross = z * (z + 1) * length * length;
259 
260  return cross;
261 }
262 
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
264 
266  G4double z)
267 {
268  // From Phys Med Biol 37 (1992) 1841-1858
269  // Z=7.42 for water
270 
271  const G4double constK(1.7E-5);
272 
273  G4double beta2;
274  beta2 = 1. - 1. / ((1. + k / electron_mass_c2) * (1. + k / electron_mass_c2));
275 
276  G4double etaC;
277  if (k < 50 * keV)
278  etaC = 1.198;
279  else
280  etaC = 1.13 + 3.76 * (z * z / (137 * 137 * beta2));
281 
282  G4double numerator = etaC * constK * std::pow(z, 2. / 3.);
283 
284  k /= electron_mass_c2;
285 
286  G4double denominator = k * (2 + k);
287 
288  G4double value = 0.;
289  if (denominator > 0.)
290  value = numerator / denominator; // form 3
291 
292  return value;
293 }
294 
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296 
297 void
299 SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
300  const G4MaterialCutsCouple* /*couple*/,
301  const G4DynamicParticle* aDynamicElectron,
302  G4double,
303  G4double)
304 {
305 #ifdef UEHARA_VERBOSE
306  if (verboseLevel > 3)
307  {
308  G4cout
309  << "Calling SampleSecondaries() of G4DNAUeharaScreenedRutherfordElasticModel"
310  << G4endl;
311  }
312 #endif
313 
314  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
315 
316  G4double cosTheta = 0.;
317 
318  if (electronEnergy0<intermediateEnergyLimit)
319  {
320 #ifdef UEHARA_VERBOSE
321  if (verboseLevel > 3)
322  G4cout << "---> Using Brenner & Zaider model" << G4endl;
323 #endif
324  cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
325  }
326  else //if (electronEnergy0>=intermediateEnergyLimit)
327  {
328 #ifdef UEHARA_VERBOSE
329  if (verboseLevel > 3)
330  G4cout << "---> Using Screened Rutherford model" << G4endl;
331 #endif
332  G4double z = 7.42; // FROM PMB 37 (1992) 1841-1858 p1842
333  cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
334  }
335 
336  G4double phi = 2. * pi * G4UniformRand();
337 
338  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
339  G4ThreeVector xVers = zVers.orthogonal();
340  G4ThreeVector yVers = zVers.cross(xVers);
341 
342  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
343  G4double yDir = xDir;
344  xDir *= std::cos(phi);
345  yDir *= std::sin(phi);
346 
347  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
348 
350 
352 }
353 
354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355 
356 G4double
359 {
360  // d sigma_el 1 beta(K)
361  // ------------ (K) ~ --------------------------------- + ---------------------------------
362  // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2
363  //
364  // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
365  //
366  // Phys. Med. Biol. 29 N.4 (1983) 443-447
367 
368  // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
369 
370  k /= eV;
371 
374  G4double gamma;
375 
376  if (k > 100.)
377  {
379  // Only in this case it is not the exponent of the polynomial
380  } else
381  {
382  if (k > 10)
383  {
385  }
386  else
387  {
389  }
390  }
391 
392  // ***** Original method
393 
394  if (fasterCode == false)
395  {
396  G4double oneOverMax = 1.
397  / (1. / (4. * gamma * gamma)
398  + beta / ((2. + 2. * delta) * (2. + 2. * delta)));
399 
400  G4double cosTheta = 0.;
401  G4double leftDenominator = 0.;
402  G4double rightDenominator = 0.;
403  G4double fCosTheta = 0.;
404 
405  do
406  {
407  cosTheta = 2. * G4UniformRand()- 1.;
408 
409  leftDenominator = (1. + 2.*gamma - cosTheta);
410  rightDenominator = (1. + 2.*delta + cosTheta);
411  if ( (leftDenominator * rightDenominator) != 0. )
412  {
413  fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator)
414  + beta/(rightDenominator*rightDenominator));
415  }
416  }
417  while (fCosTheta < G4UniformRand());
418 
419  return cosTheta;
420  }
421 
422  // ***** Alternative method using cumulative probability
423 
424  else // if (fasterCode)
425  {
426 
427  //
428  // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
429  //
430  // An integral of differential cross-section formula shown above this member function
431  // (integral variable: cos(theta), integral interval: [-1, x]) is as follows:
432  //
433  // 1.0 + x beta * (1 + x)
434  // I = --------------------- + ---------------------- (1)
435  // (a - x) * (a + 1.0) (b + x) * (b - 1.0)
436  //
437  // where a = 1.0 + 2.0 * gamma(K), b = 1.0 + 2.0 * delta(K)
438  //
439  // Then, a cumulative probability (cp) is as follows:
440  //
441  // cp 1.0 + x beta * (1 + x)
442  // ---- = --------------------- + ---------------------- (2)
443  // S (a - x) * (a + 1.0) (b + x) * (b - 1.0)
444  //
445  // where 1/S is the integral of differnetical cross-section (1) on interval [-1, 1]
446  //
447  // 1 2.0 2.0 * beta
448  // --- = ----------------------- + ----------------------- (3)
449  // S (a - 1.0) * (a + 1.0) (b + 1.0) * (b - 1.0)
450  //
451  // x is calculated from the quadratic equation derived from (2) and (3):
452  //
453  // A * x^2 + B * x + C = 0
454  //
455  // where A, B, anc C are coefficients of the equation:
456  // A = S * {(b - 1.0) - beta * (a + 1.0)} + cp * (a + 1.0) * (b - 1.0),
457  // B = S * {(b - 1.0) * (b + 1.0) + beta * (a - 1.0) * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * (a - b)
458  // C = S * {b * (b - 1.0) + beta * a * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * ab
459  //
460 
461  // sampling cumulative probability
463 
464  G4double a = 1.0 + 2.0 * gamma;
465  G4double b = 1.0 + 2.0 * delta;
466  G4double a1 = a - 1.0;
467  G4double a2 = a + 1.0;
468  G4double b1 = b - 1.0;
469  G4double b2 = b + 1.0;
470  G4double c1 = a - b;
471  G4double c2 = a * b;
472 
473  G4double S = 2.0 / (a1 * a2) + 2.0 * beta / (b1 * b2); S = 1.0 / S;
474 
475  // coefficients for the quadratic equation
476  G4double A = S * (b1 - beta * a2) + cp * a2 * b1;
477  G4double B = S * (b1 * b2 + beta * a1 * a2) - cp * a2 * b1 * c1;
478  G4double C = S * (b * b1 + beta * a * a2) - cp * a2 * b1 * c2;
479 
480  // calculate cos(theta)
481  return (-1.0 * B + std::sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
482 
483  /*
484  G4double cosTheta = -1;
485  G4double cumul = 0;
486  G4double value = 0;
487  G4double leftDenominator = 0.;
488  G4double rightDenominator = 0.;
489 
490  // Number of integration steps in the -1,1 range
491  G4int iMax=200;
492 
493  G4double random = G4UniformRand();
494 
495  // Cumulate differential cross section
496  for (G4int i=0; i<iMax; i++)
497  {
498  cosTheta = -1 + i*2./(iMax-1);
499  leftDenominator = (1. + 2.*gamma - cosTheta);
500  rightDenominator = (1. + 2.*delta + cosTheta);
501  if ( (leftDenominator * rightDenominator) != 0. )
502  {
503  cumul = cumul + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
504  }
505  }
506 
507  // Select cosTheta
508  for (G4int i=0; i<iMax; i++)
509  {
510  cosTheta = -1 + i*2./(iMax-1);
511  leftDenominator = (1. + 2.*gamma - cosTheta);
512  rightDenominator = (1. + 2.*delta + cosTheta);
513  if (cumul !=0 && (leftDenominator * rightDenominator) != 0.)
514  value = value + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)) / cumul;
515  if (random < value) break;
516  }
517 
518  return cosTheta;
519  */
520  }
521 
522  return 0.;
523 
524 }
525 
526 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
527 
528 G4double
531  std::vector<G4double>& vec)
532 {
533  // Sum_{i=0}^{size-1} vector_i k^i
534  //
535  // Phys. Med. Biol. 29 N.4 (1983) 443-447
536 
537  G4double result = 0.;
538  size_t size = vec.size();
539 
540  while (size > 0)
541  {
542  size--;
543 
544  result *= k;
545  result += vec[size];
546  }
547 
548  return result;
549 }
550 
551 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
552 
553 G4double
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
566  // satisfied within the validity of the process)
567  //
568  // Phys. Med. Biol. 45 (2000) 3171-3194
569 
570  // ***** Original method
571 
572  if (fasterCode == false)
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  else // if (fasterCode)
594  {
595 
596  //
597  // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
598  //
599  // The cumulative probability (cp) is calculated by integrating
600  // the differential cross-section fomula with cos(theta):
601  //
602  // n(K) * (1.0 + cos(theta))
603  // cp = ---------------------------------
604  // 1.0 + 2.0 * n(K) - cos(theta)
605  //
606  // Then, cos(theta) is as follows:
607  //
608  // cp * (1.0 + 2.0 * n(K)) - n(K)
609  // cos(theta) = --------------------------------
610  // n(k) + cp
611  //
612  // where, K is kinetic energy, n(K) is screeing factor, and cp is cumulative probability
613  //
614 
615  G4double n = ScreeningFactor(k, z);
617  G4double numerator = cp * (1.0 + 2.0 * n) - n;
618  G4double denominator = n + cp;
619  return numerator / denominator;
620 
621  /*
622  G4double cosTheta = -1;
623  G4double cumul = 0;
624  G4double value = 0;
625  G4double n = ScreeningFactor(k, z);
626  G4double fCosTheta;
627 
628  // Number of integration steps in the -1,1 range
629  G4int iMax=200;
630 
631  G4double random = G4UniformRand();
632 
633  // Cumulate differential cross section
634  for (G4int i=0; i<iMax; i++)
635  {
636  cosTheta = -1 + i*2./(iMax-1);
637  fCosTheta = (1 + 2.*n - cosTheta);
638  if (fCosTheta !=0.) cumul = cumul + 1./(fCosTheta*fCosTheta);
639  }
640 
641  // Select cosTheta
642  for (G4int i=0; i<iMax; i++)
643  {
644  cosTheta = -1 + i*2./(iMax-1);
645  fCosTheta = (1 + 2.*n - cosTheta);
646  if (cumul !=0.) value = value + (1./(fCosTheta*fCosTheta)) / cumul;
647  if (random < value) break;
648  }
649  return cosTheta;
650  */
651  }
652 
653  return 0.;
654 }
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static constexpr double MeV
Definition: G4SIunits.hh:214
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
static constexpr double keV
Definition: G4SIunits.hh:216
size_t GetIndex() const
Definition: G4Material.hh:262
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
Double_t z
const G4String & GetParticleName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4DNAUeharaScreenedRutherfordElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAUeharaScreenedRutherfordElasticModel")
Double_t beta
static constexpr double epsilon0
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:608
double S(double temp)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:616
G4double CalculatePolynomial(G4double k, std::vector< G4double > &vec)
double G4double
Definition: G4Types.hh:76
TCanvas * c2
Definition: plot_hist.C:75
const XML_Char int const XML_Char * value
Definition: expat.h:331
static constexpr double electron_mass_c2
static G4DNAMolecularMaterial * Instance()
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
#define G4UniformRand()
Definition: Randomize.hh:53
double A(double temperature)
static constexpr double eV
Definition: G4SIunits.hh:215
static constexpr double eV
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
Hep3Vector unit() const
Hep3Vector cross(const Hep3Vector &) const
G4double G4ParticleHPJENDLHEData::G4double result
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
double C(double temp)
static const G4double cp
G4double GetKineticEnergy() const
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
Hep3Vector orthogonal() const
Char_t n[5]
static constexpr double pi
Definition: G4SIunits.hh:75
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:131
static constexpr double e_squared
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
static constexpr double keV
double B(double temperature)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)