Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4WentzelOKandVIxSection.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 // $Id: G4WentzelOKandVIxSection.cc 109683 2018-05-08 10:36:32Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4WentzelOKandVIxSection
34 //
35 // Author: V.Ivanchenko
36 //
37 // Creation date: 09.04.2008 from G4MuMscModel
38 //
39 // Modifications:
40 //
41 // -------------------------------------------------------------------
42 //
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "Randomize.hh"
52 #include "G4Electron.hh"
53 #include "G4Positron.hh"
54 #include "G4Proton.hh"
55 #include "G4EmParameters.hh"
56 #include "G4Log.hh"
57 #include "G4Exp.hh"
58 
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 
64 
65 #ifdef G4MULTITHREADED
66 G4Mutex G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex = G4MUTEX_INITIALIZER;
67 #endif
68 
69 using namespace std;
70 
72  temp(0.,0.,0.),
73  numlimit(0.1),
74  nwarnings(0),
75  nwarnlimit(50),
76  isCombined(comb),
77  cosThetaMax(-1.0),
79 {
82  fMottXSection = nullptr;
83 
87 
88  lowEnergyLimit = 1.0*eV;
90  coeff = twopi*p0*p0;
91  particle = nullptr;
92 
94 
95  currentMaterial = nullptr;
96  factB = factD = formfactA = screenZ = 0.0;
98  = gam0pcmp = pcmp2 = 1.0;
99 
101 
102  tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
103  ecut = etag = DBL_MAX;
104  targetZ = 0;
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109 
111 {}
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114 
116  G4double cosThetaLim)
117 {
118  SetupParticle(p);
119  tkin = mom2 = momCM2 = 0.0;
120  ecut = etag = DBL_MAX;
121  targetZ = 0;
122 
123  // cosThetaMax is below 1.0 only when MSC is combined with SS
124  if(isCombined) { cosThetaMax = cosThetaLim; }
125 
128  factorA2 = 0.5*a*a;
129  currentMaterial = nullptr;
130 
132  if(0.0 == ScreenRSquare[0]) { InitialiseA(); }
133 
134  // Mott corrections
135  if((p == theElectron || p == thePositron) && !fMottXSection) {
137  fMottXSection->Initialise(p, 1.0);
138  }
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 
144 {
145  // Thomas-Fermi screening radii
146  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
147 #ifdef G4MULTITHREADED
148  G4MUTEXLOCK(&G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex);
149 #endif
150  if(0.0 == ScreenRSquare[0]) {
151  G4double a0 = electron_mass_c2/0.88534;
152  G4double constn = 6.937e-6/(MeV*MeV);
154 
155  G4double afact = fct*0.5*alpha2*a0*a0;
156  ScreenRSquare[0] = afact;
157  ScreenRSquare[1] = afact;
158  ScreenRSquareElec[1] = afact;
159  FormFactor[1] = 3.097e-6/(MeV*MeV);
160 
161  for(G4int j=2; j<100; ++j) {
162  G4double x = fG4pow->Z13(j);
163  ScreenRSquare[j] = afact*(1 + G4Exp(-j*j*0.001))*x*x;
164  ScreenRSquareElec[j] = afact*x*x;
165  x = fNistManager->GetA27(j);
166  FormFactor[j] = constn*x*x;
167  }
168  }
169 #ifdef G4MULTITHREADED
170  G4MUTEXUNLOCK(&G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex);
171 #endif
172 
173  //G4cout << "G4WentzelOKandVIxSection::Initialise mass= " << mass
174  // << " " << p->GetParticleName()
175  // << " cosThetaMax= " << cosThetaMax << G4endl;
176 
177 }
178 
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180 
182 {
183  /*
184  G4cout << "G4WentzelOKandVIxSection::SetupParticle " << p
185  << " " << particle << " " << this << G4endl;
186  G4cout << this << " " << p->GetParticleName() << G4endl;
187  */
188  particle = p;
189  mass = particle->GetPDGMass();
190  spin = particle->GetPDGSpin();
191  if(0.0 != spin) { spin = 0.5; }
192  G4double q = std::abs(particle->GetPDGCharge()/eplus);
193  chargeSquare = q*q;
194  charge3 = chargeSquare*q;
195  tkin = 0.0;
196  currentMaterial = nullptr;
197  targetZ = 0;
198 }
199 
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
201 
202 G4double
204 {
205  if(ekin != tkin || mat != currentMaterial) {
207  tkin = ekin;
208  mom2 = tkin*(tkin + 2.0*mass);
209  invbeta2 = 1.0 + mass*mass/mom2;
210  factB = spin/invbeta2;
213  : cosThetaMax;
214  }
215  return cosTetMaxNuc;
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
219 
220 G4double
222 {
223  G4double cosTetMaxNuc2 = cosTetMaxNuc;
224  if(Z != targetZ || tkin != etag) {
225  etag = tkin;
226  targetZ = std::min(Z, 99);
227  G4double massT = (1 == Z) ? CLHEP::proton_mass_c2 :
229  SetTargetMass(massT);
230 
233  fMottFactor = (1.0 + 2.0e-4*Z*Z);
234  }
235 
236  if(1 == Z) {
238  } else if(mass > MeV) {
239  screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
241  } else {
242  G4double tau = tkin/mass;
243  screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z
244  *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))*
246  }
247  if(targetZ == 1 && particle == theProton && cosTetMaxNuc2 < 0.0) {
248  cosTetMaxNuc2 = 0.0;
249  }
251 
252  cosTetMaxElec = 1.0;
254  }
255  return cosTetMaxNuc2;
256 }
257 
258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
259 
260 G4double
262 {
263  G4double xSection = 0.0;
264  if(cosTMax >= 1.0) { return xSection; }
265 
266  G4double costm = std::max(cosTMax,cosTetMaxElec);
268 
269  // scattering off electrons
270  if(costm < 1.0) {
271  G4double x = (1.0 - costm)/screenZ;
272  if(x < numlimit) {
273  G4double x2 = 0.5*x*x;
274  xSection = x2*((1.0 - 1.3333333*x + 3*x2) - fb*x*(0.6666667 - x));
275  } else {
276  G4double x1= x/(1 + x);
277  G4double xlog = G4Log(1.0 + x);
278  xSection = xlog - x1 - fb*(x + x1 - 2*xlog);
279  }
280 
281  if(xSection < 0.0) {
282  ++nwarnings;
283  if(nwarnings < nwarnlimit) {
284  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom"
285  << " scattering on e- <0"
286  << G4endl;
287  G4cout << "cross= " << xSection
288  << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
289  << " Z= " << targetZ << " "
291  G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
292  << " x= " << x << G4endl;
293  }
294  xSection = 0.0;
295  }
296  }
297  /*
298  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom: \n"
299  << " Z= " << targetZ
300  << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
301  << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
302  << " zmaxN= " << (1.0 - cosThetaMax)/screenZ
303  << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
304  */
305  // scattering off nucleus
306  if(cosTMax < 1.0) {
307  G4double x = (1.0 - cosTMax)/screenZ;
308  G4double y;
309  if(x < numlimit) {
310  G4double x2 = 0.5*x*x;
311  y = x2*((1.0 - 1.3333333*x + 3*x2) - fb*x*(0.6666667 - x));
312  } else {
313  G4double x1= x/(1 + x);
314  G4double xlog = G4Log(1.0 + x);
315  y = xlog - x1 - fb*(x + x1 - 2*xlog);
316  }
317 
318  if(y < 0.0) {
319  ++nwarnings;
320  if(nwarnings < nwarnlimit) {
321  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom"
322  << " scattering on nucleus <0"
323  << G4endl;
324  G4cout << "y= " << y
325  << " e(MeV)= " << tkin << " Z= " << targetZ << " "
327  G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
328  << " x= " << x <<G4endl;
329  }
330  y = 0.0;
331  }
332  xSection += y*targetZ;
333  }
334  xSection *= kinFactor;
335 
336  /*
337  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
338  << " screenZ= " << screenZ << " formF= " << formfactA
339  << " for " << particle->GetParticleName()
340  << " m= " << mass << " 1/v= " << sqrt(invbeta2)
341  << " p= " << sqrt(mom2)
342  << " x= " << x << G4endl;
343  */
344  return xSection;
345 }
346 
347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
348 
351  G4double cosTMax,
352  G4double elecRatio)
353 {
354  temp.set(0.0,0.0,1.0);
355  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
356 
357  G4double formf = formfactA;
358  G4double cost1 = cosTMin;
359  G4double cost2 = cosTMax;
360  if(elecRatio > 0.0) {
361  if(rndmEngineMod->flat() <= elecRatio) {
362  formf = 0.0;
363  cost1 = std::max(cost1,cosTetMaxElec);
364  cost2 = std::max(cost2,cosTetMaxElec);
365  }
366  }
367  if(cost1 > cost2) {
368 
369  G4double w1 = 1. - cost1 + screenZ;
370  G4double w2 = 1. - cost2 + screenZ;
371  G4double z1 = w1*w2/(w1 + rndmEngineMod->flat()*(w2 - w1)) - screenZ;
372 
373  G4double fm = 1.0;
375  fm += formf*z1;
376  fm = 1.0/(fm*fm);
377  } else if(fNucFormfactor == fGaussianNF) {
378  fm = G4Exp(-2*formf*z1);
379  } else if(fNucFormfactor == fFlatNF) {
380  static const G4double ccoef = 0.00508/MeV;
381  G4double x = std::sqrt(2.*mom2*z1)*ccoef*2.;
382  fm = FlatFormfactor(x);
383  fm *= FlatFormfactor(x*0.6
385  }
386  G4double grej;
387  if(fMottXSection) {
389  grej = fMottXSection->RatioMottRutherfordCosT(std::sqrt(z1))*fm*fm;
390  } else {
391  grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2. - z1))
392  *fm*fm/(1.0 + z1*factD);
393  }
394  //G4cout << "SampleSingleScattering: E= " << tkin << " z1= "
395  // << z1 << " grej= " << grej << G4endl;
396  if(fMottFactor*rndmEngineMod->flat() <= grej ) {
397  // exclude "false" scattering due to formfactor and spin effect
398  G4double cost = 1.0 - z1;
399  if(cost > 1.0) { cost = 1.0; }
400  else if(cost < -1.0) { cost =-1.0; }
401  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
402  //G4cout << "sint= " << sint << G4endl;
403  G4double phi = twopi*rndmEngineMod->flat();
404  temp.set(sint*cos(phi),sint*sin(phi),cost);
405  }
406  }
407  return temp;
408 }
409 
410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
411 
412 void
414 {
415  if(mass > MeV) {
417  G4double tau = tkin/mass;
418  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
419  (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
420  cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
421  } else {
422 
423  G4double tmax = (particle == theElectron) ? 0.5*tkin : tkin;
424  G4double t = std::min(cutEnergy, tmax);
425  G4double mom21 = t*(t + 2.0*electron_mass_c2);
426  G4double t1 = tkin - t;
427  //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
428  //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
429  if(t1 > 0.0) {
430  G4double mom22 = t1*(t1 + 2.0*mass);
431  G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
432  if(ctm < 1.0) { cosTetMaxElec = ctm; }
433  if(particle == theElectron && cosTetMaxElec < 0.0) {
434  cosTetMaxElec = 0.0;
435  }
436  }
437  }
438 }
439 
440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
441 
442 G4double
444 {
445  return 0.0;
446 }
447 
448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Float_t x
Definition: compare.C:6
void set(double x, double y, double z)
G4double ComputeSecondTransportMoment(G4double CosThetaMax)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
const G4ParticleDefinition * thePositron
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ThreeVector & SampleSingleScattering(G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4double Z13(G4int Z) const
Definition: G4Pow.hh:132
G4double A13(G4double A) const
Definition: G4Pow.cc:138
TTree * t1
Definition: plottest35.C:26
static constexpr double MeV
Definition: G4SIunits.hh:214
const G4ParticleDefinition * theProton
Float_t x1[n_points_granero]
Definition: compare.C:5
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
const char * p
Definition: xmltok.h:285
static constexpr double hbarc
const G4String & GetParticleName() const
G4double GetPDGCharge() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double GetAtomicMassAmu(const G4String &symb) const
G4double GetPDGMass() const
static constexpr double classic_electr_radius
G4double G4Log(G4double x)
Definition: G4Log.hh:230
static constexpr double proton_mass_c2
Float_t Z
static constexpr double fermi
Definition: SystemOfUnits.h:83
G4double GetPDGSpin() const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static constexpr double amu_c2
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
const G4ParticleDefinition * particle
virtual double flat()=0
void ComputeMaxElectronScattering(G4double cut)
void Initialise(const G4ParticleDefinition *, G4double cosThetaLim)
G4NuclearFormfactorType fNucFormfactor
static constexpr double electron_mass_c2
static constexpr double twopi
Definition: G4SIunits.hh:76
static G4Positron * Positron()
Definition: G4Positron.cc:94
static constexpr double eV
Definition: G4SIunits.hh:215
Float_t mat
void SetupParticle(const G4ParticleDefinition *)
G4WentzelOKandVIxSection(G4bool comb=true)
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:234
G4double ScreeningFactor() const
TFile fb("Li6.root")
G4double SetupTarget(G4int Z, G4double cut)
static constexpr double eplus
Definition: G4SIunits.hh:199
const G4double a0
void Initialise(const G4ParticleDefinition *, G4double CosThetaLim)
G4double GetInvA23() const
int G4int
Definition: G4Types.hh:78
G4NuclearFormfactorType NuclearFormfactorType() const
const G4ParticleDefinition * theElectron
G4double GetA27(G4int Z) const
static G4double ScreenRSquareElec[100]
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:227
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:233
G4ScreeningMottCrossSection * fMottXSection
G4GLOB_DLL std::ostream G4cout
G4double FactorForAngleLimit() const
void SetupKinematic(G4double kinEnergy, G4double Z)
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
static constexpr double fine_structure_const
Float_t x2[n_points_geant4]
Definition: compare.C:26
#define DBL_MAX
Definition: templates.hh:83
static G4EmParameters * Instance()
static constexpr double pi
Definition: SystemOfUnits.h:54
G4double Z23(G4int Z) const
Definition: G4Pow.hh:137
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4NistManager * Instance()
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
std::mutex G4Mutex
Definition: G4Threading.hh:84