Geant4  10.00.p02
G4WentzelOKandVIxSection.cc
Go to the documentation of this file.
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 81865 2014-06-06 11:32:58Z 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 //
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "Randomize.hh"
53 #include "G4Electron.hh"
54 #include "G4Positron.hh"
55 #include "G4Proton.hh"
56 #include "G4LossTableManager.hh"
57 #include "G4Log.hh"
58 #include "G4Exp.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
65 
66 using namespace std;
67 
69  numlimit(0.1),
70  nwarnings(0),
71  nwarnlimit(50),
72  alpha2(fine_structure_const*fine_structure_const)
73 {
79  lowEnergyLimit = 1.0*eV;
80  G4double p0 = electron_mass_c2*classic_electr_radius;
81  coeff = twopi*p0*p0;
82  particle = 0;
83 
84  // Thomas-Fermi screening radii
85  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
86 
87  if(0.0 == ScreenRSquare[0]) {
88  G4double a0 = electron_mass_c2/0.88534;
89  G4double constn = 6.937e-6/(MeV*MeV);
90 
91  ScreenRSquare[0] = alpha2*a0*a0;
93  for(G4int j=1; j<100; ++j) {
94  G4double x = a0*fG4pow->Z13(j);
95  if(1 == j) { ScreenRSquare[j] = 0.5*alpha2*a0*a0; }
96  else {
97  ScreenRSquare[j] = 0.5*(1 + G4Exp(-j*j*0.001))*alpha2*x*x;
98  ScreenRSquareElec[j] = 0.5*alpha2*x*x;
99  }
100  x = fNistManager->GetA27(j);
101  FormFactor[j] = constn*x*x;
102  }
103  }
104  currentMaterial = 0;
105  elecXSRatio = factB = factD = formfactA = screenZ = 0.0;
107 
108  factB1= 0.5*CLHEP::pi*fine_structure_const;
109 
110  tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
111  ecut = etag = DBL_MAX;
112  targetZ = 0;
113  cosThetaMax = 1.0;
114  targetMass = proton_mass_c2;
115  particle = 0;
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119 
121 {}
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 
126  G4double CosThetaLim)
127 {
128  SetupParticle(p);
129  tkin = mom2 = momCM2 = 0.0;
130  ecut = etag = DBL_MAX;
131  targetZ = 0;
132  cosThetaMax = CosThetaLim;
134  *CLHEP::hbarc/CLHEP::fermi;
135  factorA2 = 0.5*a*a;
136  currentMaterial = 0;
137 
138  //G4cout << "G4WentzelOKandVIxSection::Initialise mass= " << mass
139  // << " " << p->GetParticleName() << G4endl;
140 
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144 
146 {
147  particle = p;
148  mass = particle->GetPDGMass();
149  spin = particle->GetPDGSpin();
150  if(0.0 != spin) { spin = 0.5; }
151  G4double q = std::fabs(particle->GetPDGCharge()/eplus);
152  chargeSquare = q*q;
153  charge3 = chargeSquare*q;
154  tkin = 0.0;
155  currentMaterial = 0;
156  targetZ = 0;
157 }
158 
159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
160 
161 G4double
163 {
164  G4double cosTetMaxNuc2 = cosTetMaxNuc;
165  if(Z != targetZ || tkin != etag) {
166  etag = tkin;
167  targetZ = Z;
168  if(targetZ > 99) { targetZ = 99; }
169  G4double massT = proton_mass_c2;
170  if(targetZ > 1) {
171  massT = fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2;
172  }
173  SetTargetMass(massT);
174  //G4double tmass2 = targetMass*targetMass;
175  //G4double etot = tkin + mass;
176  //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass;
177  //momCM2 = mom2*tmass2/invmass2;
178  //gam0pcmp = (etot + targetMass)*targetMass/invmass2;
179  //pcmp2 = tmass2/invmass2;
180 
182 
183  if(1 == Z) {
185  } else if(mass > MeV) {
186  screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
188  } else {
189  G4double tau = tkin/mass;
190  screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z
191  *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))*
193  }
194  if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
195  cosTetMaxNuc2 = 0.0;
196  }
198 
199  cosTetMaxElec = 1.0;
201  }
202  return cosTetMaxNuc2;
203 }
204 
205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
206 
207 G4double
209 {
210  G4double xsec = 0.0;
211  if(cosTMax >= 1.0) { return xsec; }
212 
213  G4double xSection = 0.0;
214  G4double x = 0;
215  G4double y = 0;
216  G4double x1= 0;
217  G4double x2= 0;
218  G4double xlog = 0.0;
219 
220  G4double costm = std::max(cosTMax,cosTetMaxElec);
221  G4double fb = screenZ*factB;
222 
223  // scattering off electrons
224  if(costm < 1.0) {
225  x = (1.0 - costm)/screenZ;
226  if(x < numlimit) {
227  x2 = 0.5*x*x;
228  y = x2*(1.0 - 1.3333333*x + 3*x2);
229  if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
230  } else {
231  x1= x/(1 + x);
232  xlog = G4Log(1.0 + x);
233  y = xlog - x1;
234  if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
235  }
236 
237  if(y < 0.0) {
238  ++nwarnings;
239  if(nwarnings < nwarnlimit) {
240  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
241  << G4endl;
242  G4cout << "y= " << y
243  << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
244  << " Z= " << targetZ << " "
246  G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
247  << " x= " << x << G4endl;
248  }
249  y = 0.0;
250  }
251  xSection = y;
252  }
253  /*
254  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
255  << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
256  << " cut(MeV)= " << ecut/MeV
257  << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
258  << " zmaxN= " << (1.0 - cosThetaMax)/screenZ
259  << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
260  */
261  // scattering off nucleus
262  if(cosTMax < 1.0) {
263  x = (1.0 - cosTMax)/screenZ;
264  if(x < numlimit) {
265  x2 = 0.5*x*x;
266  y = x2*(1.0 - 1.3333333*x + 3*x2);
267  if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
268  } else {
269  x1= x/(1 + x);
270  xlog = G4Log(1.0 + x);
271  y = xlog - x1;
272  if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
273  }
274 
275  if(y < 0.0) {
276  ++nwarnings;
277  if(nwarnings < nwarnlimit) {
278  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
279  << G4endl;
280  G4cout << "y= " << y
281  << " e(MeV)= " << tkin << " Z= " << targetZ << " "
283  G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
284  << " x= " << " x1= " << x1 <<G4endl;
285  }
286  y = 0.0;
287  }
288  xSection += y*targetZ;
289  }
290  xSection *= kinFactor;
291  /*
292  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
293  << " screenZ= " << screenZ << " formF= " << formfactA
294  << " for " << particle->GetParticleName()
295  << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
296  << " x= " << x
297  << G4endl;
298  */
299  return xSection;
300 }
301 
302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303 
306  G4double cosTMax,
307  G4double elecRatio)
308 {
309  G4ThreeVector v(0.0,0.0,1.0);
310 
311  G4double formf = formfactA;
312  G4double cost1 = cosTMin;
313  G4double cost2 = cosTMax;
314  if(elecRatio > 0.0) {
315  if(G4UniformRand() <= elecRatio) {
316  formf = 0.0;
317  cost1 = std::max(cost1,cosTetMaxElec);
318  cost2 = std::max(cost2,cosTetMaxElec);
319  }
320  }
321  if(cost1 < cost2) { return v; }
322 
323  G4double w1 = 1. - cost1 + screenZ;
324  G4double w2 = 1. - cost2 + screenZ;
325  G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;
326 
327  //G4double fm = 1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass);
328  G4double fm = 1.0 + formf*z1;
329  //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm );
330  G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm );
331  // "false" scattering
332  if( G4UniformRand() > grej ) { return v; }
333  // }
334  G4double cost = 1.0 - z1;
335 
336  if(cost > 1.0) { cost = 1.0; }
337  else if(cost < -1.0) { cost =-1.0; }
338  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
339  //G4cout << "sint= " << sint << G4endl;
340  G4double phi = twopi*G4UniformRand();
341  G4double vx1 = sint*cos(phi);
342  G4double vy1 = sint*sin(phi);
343 
344  // only direction is changed
345  v.set(vx1,vy1,cost);
346  return v;
347 }
348 
349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
350 
351 void
353 {
354  if(mass > MeV) {
355  G4double ratio = electron_mass_c2/mass;
356  G4double tau = tkin/mass;
357  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
358  (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
359  //tmax = std::min(tmax, targetZ*targetZ*10*eV);
360  cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
361  } else {
362 
363  G4double tmax = tkin;
364  if(particle == theElectron) { tmax *= 0.5; }
365  //tmax = std::min(tmax, targetZ*targetZ*10*eV);
366  G4double t = std::min(cutEnergy, tmax);
367  G4double mom21 = t*(t + 2.0*electron_mass_c2);
368  G4double t1 = tkin - t;
369  //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
370  //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
371  if(t1 > 0.0) {
372  G4double mom22 = t1*(t1 + 2.0*mass);
373  G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
374  if(ctm < 1.0) { cosTetMaxElec = ctm; }
375  if(particle == theElectron && cosTetMaxElec < 0.0) { cosTetMaxElec = 0.0; }
376  }
377  }
378 }
379 
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
const G4ParticleDefinition * theElectron
const G4ParticleDefinition * theProton
const G4double a0
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
void SetupParticle(const G4ParticleDefinition *)
static const double MeV
Definition: G4SIunits.hh:193
static G4double ScreenRSquareElec[100]
static G4LossTableManager * Instance()
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector SampleSingleScattering(G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio=0.0)
const G4double pi
const G4ParticleDefinition * particle
G4double FactorForAngleLimit() const
void ComputeMaxElectronScattering(G4double cut)
G4double GetA27(G4int Z)
G4double a
Definition: TRTMaterials.hh:39
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4double SetupTarget(G4int Z, G4double cut=DBL_MAX)
const G4String & GetParticleName() const
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4double Z13(G4int Z) const
Definition: G4Pow.hh:129
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const double eV
Definition: G4SIunits.hh:194
#define fm
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4double GetPDGMass() const
const G4ParticleDefinition * thePositron
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetAtomicMassAmu(const G4String &symb) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double Z23(G4int Z) const
Definition: G4Pow.hh:153
G4double GetPDGSpin() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:178
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
static const double fermi
Definition: G4SIunits.hh:93
void Initialise(const G4ParticleDefinition *, G4double CosThetaLim)