Geant4.10
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4WentzelOKandVIxSection Class Reference

#include <G4WentzelOKandVIxSection.hh>

Public Member Functions

 G4WentzelOKandVIxSection ()
 
virtual ~G4WentzelOKandVIxSection ()
 
void Initialise (const G4ParticleDefinition *, G4double CosThetaLim)
 
void SetupParticle (const G4ParticleDefinition *)
 
G4double SetupTarget (G4int Z, G4double cut=DBL_MAX)
 
G4double ComputeTransportCrossSectionPerAtom (G4double CosThetaMax)
 
G4ThreeVector SampleSingleScattering (G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio=0.0)
 
G4double ComputeNuclearCrossSection (G4double CosThetaMin, G4double CosThetaMax)
 
G4double ComputeElectronCrossSection (G4double CosThetaMin, G4double CosThetaMax)
 
G4double SetupKinematic (G4double kinEnergy, const G4Material *mat)
 
void SetTargetMass (G4double value)
 
G4double GetMomentumSquare () const
 
G4double GetCosThetaNuc () const
 
G4double GetCosThetaElec () const
 

Detailed Description

Definition at line 71 of file G4WentzelOKandVIxSection.hh.

Constructor & Destructor Documentation

G4WentzelOKandVIxSection::G4WentzelOKandVIxSection ( )

Definition at line 68 of file G4WentzelOKandVIxSection.cc.

References python.hepunit::classic_electr_radius, DBL_MAX, G4Electron::Electron(), python.hepunit::electron_mass_c2, python.hepunit::eV, python.hepunit::fine_structure_const, G4Exp(), G4NistManager::GetA27(), G4Pow::GetInstance(), G4NistManager::Instance(), python.hepunit::MeV, G4Positron::Positron(), G4Proton::Proton(), python.hepunit::proton_mass_c2, python.hepunit::twopi, test::x, and G4Pow::Z13().

68  :
69  numlimit(0.1),
70  nwarnings(0),
71  nwarnlimit(50),
73 {
74  fNistManager = G4NistManager::Instance();
75  fG4pow = G4Pow::GetInstance();
76  theElectron = G4Electron::Electron();
77  thePositron = G4Positron::Positron();
78  theProton = G4Proton::Proton();
79  lowEnergyLimit = 1.0*eV;
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;
92  ScreenRSquareElec[0] = ScreenRSquare[0];
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;
106  cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.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 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
G4double GetA27(G4int Z)
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4double Z13(G4int Z) const
Definition: G4Pow.hh:129
static G4Proton * Proton()
Definition: G4Proton.cc:93
float proton_mass_c2
Definition: hepunit.py:275
float electron_mass_c2
Definition: hepunit.py:274
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
static G4Positron * Positron()
Definition: G4Positron.cc:94
static G4Electron * Electron()
Definition: G4Electron.cc:94
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection ( )
virtual

Definition at line 120 of file G4WentzelOKandVIxSection.cc.

121 {}

Member Function Documentation

G4double G4WentzelOKandVIxSection::ComputeElectronCrossSection ( G4double  CosThetaMin,
G4double  CosThetaMax 
)
inline

Definition at line 237 of file G4WentzelOKandVIxSection.hh.

References G4INCL::Math::max().

Referenced by G4eCoulombScatteringModel::ComputeCrossSectionPerAtom().

239 {
240  G4double xsec = 0.0;
241  G4double cost1 = std::max(cosTMin,cosTetMaxElec);
242  G4double cost2 = std::max(cosTMax,cosTetMaxElec);
243  if(cost1 > cost2) {
244  xsec = kinFactor*(cost1 - cost2)/((1.0 - cost1 + screenZ)*(1.0 - cost2 + screenZ));
245  }
246  return xsec;
247 }
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double G4double
Definition: G4Types.hh:76
G4double G4WentzelOKandVIxSection::ComputeNuclearCrossSection ( G4double  CosThetaMin,
G4double  CosThetaMax 
)
inline

Definition at line 223 of file G4WentzelOKandVIxSection.hh.

Referenced by G4eCoulombScatteringModel::ComputeCrossSectionPerAtom().

225 {
226  G4double xsec = 0.0;
227  if(cosTMax < cosTMin) {
228  xsec = targetZ*kinFactor*(cosTMin - cosTMax)/
229  ((1.0 - cosTMin + screenZ)*(1.0 - cosTMax + screenZ));
230  }
231  return xsec;
232 }
double G4double
Definition: G4Types.hh:76
G4double G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom ( G4double  CosThetaMax)

Definition at line 204 of file G4WentzelOKandVIxSection.cc.

References G4cout, G4endl, G4Log(), G4ParticleDefinition::GetParticleName(), G4INCL::Math::max(), and test::x.

Referenced by G4WentzelVIModel::ComputeCrossSectionPerAtom().

205 {
206  G4double xsec = 0.0;
207  if(cosTMax >= 1.0) { return xsec; }
208 
209  G4double xSection = 0.0;
210  G4double x = 0;
211  G4double y = 0;
212  G4double x1= 0;
213  G4double x2= 0;
214  G4double xlog = 0.0;
215 
216  G4double costm = std::max(cosTMax,cosTetMaxElec);
217  G4double fb = screenZ*factB;
218 
219  // scattering off electrons
220  if(costm < 1.0) {
221  x = (1.0 - costm)/screenZ;
222  if(x < numlimit) {
223  x2 = 0.5*x*x;
224  y = x2*(1.0 - 1.3333333*x + 3*x2);
225  if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
226  } else {
227  x1= x/(1 + x);
228  xlog = G4Log(1.0 + x);
229  y = xlog - x1;
230  if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
231  }
232 
233  if(y < 0.0) {
234  ++nwarnings;
235  if(nwarnings < nwarnlimit) {
236  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
237  << G4endl;
238  G4cout << "y= " << y
239  << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
240  << " Z= " << targetZ << " "
241  << particle->GetParticleName() << G4endl;
242  G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
243  << " x= " << x << G4endl;
244  }
245  y = 0.0;
246  }
247  xSection = y;
248  }
249  /*
250  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
251  << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
252  << " cut(MeV)= " << ecut/MeV
253  << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
254  << " zmaxN= " << (1.0 - cosThetaMax)/screenZ
255  << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
256  */
257  // scattering off nucleus
258  if(cosTMax < 1.0) {
259  x = (1.0 - cosTMax)/screenZ;
260  if(x < numlimit) {
261  x2 = 0.5*x*x;
262  y = x2*(1.0 - 1.3333333*x + 3*x2);
263  if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
264  } else {
265  x1= x/(1 + x);
266  xlog = G4Log(1.0 + x);
267  y = xlog - x1;
268  if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
269  }
270 
271  if(y < 0.0) {
272  ++nwarnings;
273  if(nwarnings < nwarnlimit) {
274  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
275  << G4endl;
276  G4cout << "y= " << y
277  << " e(MeV)= " << tkin << " Z= " << targetZ << " "
278  << particle->GetParticleName() << G4endl;
279  G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
280  << " x= " << " x1= " << x1 <<G4endl;
281  }
282  y = 0.0;
283  }
284  xSection += y*targetZ;
285  }
286  xSection *= kinFactor;
287  /*
288  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
289  << " screenZ= " << screenZ << " formF= " << formfactA
290  << " for " << particle->GetParticleName()
291  << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
292  << " x= " << x
293  << G4endl;
294  */
295  return xSection;
296 }
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
G4double G4Log(G4double x)
Definition: G4Log.hh:227
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double G4WentzelOKandVIxSection::GetCosThetaElec ( ) const
inline

Definition at line 215 of file G4WentzelOKandVIxSection.hh.

216 {
217  return cosTetMaxElec;
218 }
G4double G4WentzelOKandVIxSection::GetCosThetaNuc ( ) const
inline

Definition at line 208 of file G4WentzelOKandVIxSection.hh.

209 {
210  return cosTetMaxNuc;
211 }
G4double G4WentzelOKandVIxSection::GetMomentumSquare ( ) const
inline

Definition at line 201 of file G4WentzelOKandVIxSection.hh.

Referenced by G4eCoulombScatteringModel::SampleSecondaries().

202 {
203  return mom2;
204 }
void G4WentzelOKandVIxSection::Initialise ( const G4ParticleDefinition p,
G4double  CosThetaLim 
)

Definition at line 125 of file G4WentzelOKandVIxSection.cc.

References test::a, DBL_MAX, G4LossTableManager::FactorForAngleLimit(), G4LossTableManager::Instance(), and SetupParticle().

Referenced by G4WentzelVIModel::Initialise(), and G4eCoulombScatteringModel::Initialise().

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 }
void SetupParticle(const G4ParticleDefinition *)
static G4LossTableManager * Instance()
G4double FactorForAngleLimit() const
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
G4ThreeVector G4WentzelOKandVIxSection::SampleSingleScattering ( G4double  CosThetaMin,
G4double  CosThetaMax,
G4double  elecRatio = 0.0 
)

Definition at line 301 of file G4WentzelOKandVIxSection.cc.

References fm, G4UniformRand, G4INCL::Math::max(), CLHEP::Hep3Vector::set(), python.hepunit::twopi, and test::v.

Referenced by G4WentzelVIModel::SampleScattering(), and G4eCoulombScatteringModel::SampleSecondaries().

304 {
305  G4ThreeVector v(0.0,0.0,1.0);
306 
307  G4double formf = formfactA;
308  G4double cost1 = cosTMin;
309  G4double cost2 = cosTMax;
310  if(elecRatio > 0.0) {
311  if(G4UniformRand() <= elecRatio) {
312  formf = 0.0;
313  cost1 = std::max(cost1,cosTetMaxElec);
314  cost2 = std::max(cost2,cosTetMaxElec);
315  }
316  }
317  if(cost1 < cost2) { return v; }
318 
319  G4double w1 = 1. - cost1 + screenZ;
320  G4double w2 = 1. - cost2 + screenZ;
321  G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;
322 
323  //G4double fm = 1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass);
324  G4double fm = 1.0 + formf*z1;
325  //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm );
326  G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm );
327  // "false" scattering
328  if( G4UniformRand() > grej ) { return v; }
329  // }
330  G4double cost = 1.0 - z1;
331 
332  if(cost > 1.0) { cost = 1.0; }
333  else if(cost < -1.0) { cost =-1.0; }
334  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
335  //G4cout << "sint= " << sint << G4endl;
336  G4double phi = twopi*G4UniformRand();
337  G4double vx1 = sint*cos(phi);
338  G4double vy1 = sint*sin(phi);
339 
340  // only direction is changed
341  v.set(vx1,vy1,cost);
342  return v;
343 }
#define G4UniformRand()
Definition: Randomize.hh:87
#define fm
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double G4double
Definition: G4Types.hh:76
void G4WentzelOKandVIxSection::SetTargetMass ( G4double  value)
inline

Definition at line 193 of file G4WentzelOKandVIxSection.hh.

Referenced by G4eCoulombScatteringModel::SampleSecondaries(), and SetupTarget().

194 {
195  targetMass = value;
196  factD = std::sqrt(mom2)/value;
197 }
const XML_Char int const XML_Char * value
G4double G4WentzelOKandVIxSection::SetupKinematic ( G4double  kinEnergy,
const G4Material mat 
)
inline

Definition at line 177 of file G4WentzelOKandVIxSection.hh.

References G4IonisParamMat::GetInvA23(), G4Material::GetIonisation(), and G4INCL::Math::max().

Referenced by G4WentzelVIModel::ComputeCrossSectionPerAtom(), G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4WentzelVIModel::ComputeGeomPathLength(), G4WentzelVIModel::ComputeTruePathLengthLimit(), and G4WentzelVIModel::ComputeTrueStepLength().

178 {
179  if(ekin != tkin || mat != currentMaterial) {
180  currentMaterial = mat;
181  tkin = ekin;
182  mom2 = tkin*(tkin + 2.0*mass);
183  invbeta2 = 1.0 + mass*mass/mom2;
184  factB = spin/invbeta2;
185  cosTetMaxNuc =
186  std::max(cosThetaMax,1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2);
187  }
188  return cosTetMaxNuc;
189 }
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:224
G4double GetInvA23() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void G4WentzelOKandVIxSection::SetupParticle ( const G4ParticleDefinition p)

Definition at line 145 of file G4WentzelOKandVIxSection.cc.

References python.hepunit::eplus, G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and G4ParticleDefinition::GetPDGSpin().

Referenced by Initialise(), and G4eCoulombScatteringModel::SetupParticle().

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 }
const char * p
Definition: xmltok.h:285
G4double GetPDGMass() const
G4double GetPDGSpin() const
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
G4double G4WentzelOKandVIxSection::SetupTarget ( G4int  Z,
G4double  cut = DBL_MAX 
)

Definition at line 162 of file G4WentzelOKandVIxSection.cc.

References G4NistManager::GetAtomicMassAmu(), python.hepunit::MeV, G4INCL::Math::min(), SetTargetMass(), and G4Pow::Z23().

Referenced by G4WentzelVIModel::ComputeCrossSectionPerAtom(), G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(), and G4WentzelVIModel::SampleScattering().

163 {
164  G4double cosTetMaxNuc2 = cosTetMaxNuc;
165  if(Z != targetZ || tkin != etag) {
166  etag = tkin;
167  targetZ = Z;
168  if(targetZ > 99) { targetZ = 99; }
169  SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2);
170  //G4double tmass2 = targetMass*targetMass;
171  //G4double etot = tkin + mass;
172  //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass;
173  //momCM2 = mom2*tmass2/invmass2;
174  //gam0pcmp = (etot + targetMass)*targetMass/invmass2;
175  //pcmp2 = tmass2/invmass2;
176 
177  kinFactor = coeff*Z*chargeSquare*invbeta2/mom2;
178 
179  if(1 == Z) {
180  screenZ = ScreenRSquare[targetZ]/mom2;
181  } else if(mass > MeV) {
182  screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
183  ScreenRSquare[targetZ]/mom2;
184  } else {
185  G4double tau = tkin/mass;
186  screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z
187  *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))*
188  ScreenRSquareElec[targetZ]/mom2;
189  }
190  if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
191  cosTetMaxNuc2 = 0.0;
192  }
193  formfactA = FormFactor[targetZ]*mom2;
194 
195  cosTetMaxElec = 1.0;
196  ComputeMaxElectronScattering(cut);
197  }
198  return cosTetMaxNuc2;
199 }
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
double G4double
Definition: G4Types.hh:76

The documentation for this class was generated from the following files: