Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4WentzelOKandVIxSection Class Reference

#include <G4WentzelOKandVIxSection.hh>

Public Member Functions

 G4WentzelOKandVIxSection (G4bool combined=true)
 
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 ComputeSecondTransportMoment (G4double CosThetaMax)
 
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 72 of file G4WentzelOKandVIxSection.hh.

Constructor & Destructor Documentation

G4WentzelOKandVIxSection::G4WentzelOKandVIxSection ( G4bool  combined = true)
explicit

Definition at line 66 of file G4WentzelOKandVIxSection.cc.

66  :
67  temp(0.,0.,0.),
68  numlimit(0.1),
69  nwarnings(0),
70  nwarnlimit(50),
71  isCombined(combined),
72  cosThetaMax(-1.0),
73  alpha2(fine_structure_const*fine_structure_const)
74 {
75  fNistManager = G4NistManager::Instance();
76  fG4pow = G4Pow::GetInstance();
77  theElectron = G4Electron::Electron();
78  thePositron = G4Positron::Positron();
79  theProton = G4Proton::Proton();
80  lowEnergyLimit = 1.0*eV;
81  G4double p0 = electron_mass_c2*classic_electr_radius;
82  coeff = twopi*p0*p0;
83  particle = nullptr;
84 
85  fNucFormfactor = fExponentialNF;
86 
87  // Thomas-Fermi screening radii
88  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
89 
90  if(0.0 == ScreenRSquare[0]) {
91  G4double a0 = electron_mass_c2/0.88534;
92  G4double constn = 6.937e-6/(MeV*MeV);
93 
94  ScreenRSquare[0] = alpha2*a0*a0;
95  ScreenRSquareElec[0] = ScreenRSquare[0];
96  for(G4int j=1; j<100; ++j) {
97  G4double x = a0*fG4pow->Z13(j);
98  if(1 == j) { ScreenRSquare[j] = 0.5*alpha2*a0*a0; }
99  else {
100  ScreenRSquare[j] = 0.5*(1 + G4Exp(-j*j*0.001))*alpha2*x*x;
101  ScreenRSquareElec[j] = 0.5*alpha2*x*x;
102  }
103  x = fNistManager->GetA27(j);
104  FormFactor[j] = constn*x*x;
105  }
106  }
107  currentMaterial = 0;
108  factB = factD = formfactA = screenZ = 0.0;
109  cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;
110 
111  factB1= 0.5*CLHEP::pi*fine_structure_const;
112 
113  tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
114  ecut = etag = DBL_MAX;
115  targetZ = 0;
116  targetMass = proton_mass_c2;
117  particle = nullptr;
118 }
a0
const G4double a0
Definition: G4IonCoulombCrossSection.cc:68
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
CLHEP::proton_mass_c2
static constexpr double proton_mass_c2
Definition: PhysicalConstants.h:80
G4NistManager::GetA27
G4double GetA27(G4int Z) const
Definition: G4NistManager.hh:586
G4int
int G4int
Definition: G4Types.hh:78
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4Pow::Z13
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
CLHEP::classic_electr_radius
static constexpr double classic_electr_radius
Definition: PhysicalConstants.h:97
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::fine_structure_const
static constexpr double fine_structure_const
Definition: PhysicalConstants.h:96
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
fExponentialNF
Definition: G4NuclearFormfactorType.hh:50
CLHEP::pi
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection ( )
virtual

Definition at line 122 of file G4WentzelOKandVIxSection.cc.

123 {}

Member Function Documentation

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

Definition at line 250 of file G4WentzelOKandVIxSection.hh.

252 {
253  G4double xsec = 0.0;
254  G4double cost1 = std::max(cosTMin,cosTetMaxElec);
255  G4double cost2 = std::max(cosTMax,cosTetMaxElec);
256  if(cost1 > cost2) {
257  xsec = kinFactor*(cost1 - cost2)/
258  ((1.0 - cost1 + screenZ)*(1.0 - cost2 + screenZ));
259  }
260  return xsec;
261 }
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 236 of file G4WentzelOKandVIxSection.hh.

238 {
239  G4double xsec = 0.0;
240  if(cosTMax < cosTMin) {
241  xsec = targetZ*kinFactor*(cosTMin - cosTMax)/
242  ((1.0 - cosTMin + screenZ)*(1.0 - cosTMax + screenZ));
243  }
244  return xsec;
245 }
G4double
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4WentzelOKandVIxSection::ComputeSecondTransportMoment ( G4double  CosThetaMax)

Definition at line 398 of file G4WentzelOKandVIxSection.cc.

399 {
400  return 0.0;
401 }
G4double G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom ( G4double  CosThetaMax)

Definition at line 210 of file G4WentzelOKandVIxSection.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4WentzelOKandVIxSection::GetCosThetaElec ( ) const
inline

Definition at line 228 of file G4WentzelOKandVIxSection.hh.

229 {
230  return cosTetMaxElec;
231 }
G4double G4WentzelOKandVIxSection::GetCosThetaNuc ( ) const
inline

Definition at line 221 of file G4WentzelOKandVIxSection.hh.

222 {
223  return cosTetMaxNuc;
224 }
G4double G4WentzelOKandVIxSection::GetMomentumSquare ( ) const
inline

Definition at line 214 of file G4WentzelOKandVIxSection.hh.

215 {
216  return mom2;
217 }

Here is the caller graph for this function:

void G4WentzelOKandVIxSection::Initialise ( const G4ParticleDefinition *  p,
G4double  CosThetaLim 
)

Definition at line 127 of file G4WentzelOKandVIxSection.cc.

129 {
130  SetupParticle(p);
131  tkin = mom2 = momCM2 = 0.0;
132  ecut = etag = DBL_MAX;
133  targetZ = 0;
134 
135  // cosThetaMax is below 1.0 only when MSC is combined with SS
136  if(isCombined) { cosThetaMax = cosThetaLim; }
137 
138  G4double a = G4EmParameters::Instance()->FactorForAngleLimit()
139  *CLHEP::hbarc/CLHEP::fermi;
140  factorA2 = 0.5*a*a;
141  currentMaterial = nullptr;
142 
143  fNucFormfactor = G4EmParameters::Instance()->NuclearFormfactorType();
144 
145  //G4cout << "G4WentzelOKandVIxSection::Initialise mass= " << mass
146  // << " " << p->GetParticleName()
147  // << " cosThetaMax= " << cosThetaMax << G4endl;
148 
149 }
G4WentzelOKandVIxSection::SetupParticle
void SetupParticle(const G4ParticleDefinition *)
Definition: G4WentzelOKandVIxSection.cc:153
CLHEP::hbarc
static constexpr double hbarc
Definition: PhysicalConstants.h:65
G4EmParameters::NuclearFormfactorType
G4NuclearFormfactorType NuclearFormfactorType() const
Definition: G4EmParameters.cc:776
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition: G4EmParameters.cc:67
CLHEP::fermi
static constexpr double fermi
Definition: SystemOfUnits.h:83
G4double
double G4double
Definition: G4Types.hh:76
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4EmParameters::FactorForAngleLimit
G4double FactorForAngleLimit() const
Definition: G4EmParameters.cc:563

Here is the call graph for this function:

Here is the caller graph for this function:

G4ThreeVector & G4WentzelOKandVIxSection::SampleSingleScattering ( G4double  CosThetaMin,
G4double  CosThetaMax,
G4double  elecRatio = 0.0 
)

Definition at line 309 of file G4WentzelOKandVIxSection.cc.

312 {
313  temp.set(0.0,0.0,1.0);
314  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
315 
316  G4double formf = formfactA;
317  G4double cost1 = cosTMin;
318  G4double cost2 = cosTMax;
319  if(elecRatio > 0.0) {
320  if(rndmEngineMod->flat() <= elecRatio) {
321  formf = 0.0;
322  cost1 = std::max(cost1,cosTetMaxElec);
323  cost2 = std::max(cost2,cosTetMaxElec);
324  }
325  }
326  if(cost1 > cost2) {
327  G4double w1 = 1. - cost1 + screenZ;
328  G4double w2 = 1. - cost2 + screenZ;
329  G4double z1 = w1*w2/(w1 + rndmEngineMod->flat()*(w2 - w1)) - screenZ;
330 
331  G4double fm = 1.0;
332  if(fNucFormfactor == fExponentialNF) {
333  fm += formf*z1;
334  fm = 1.0/(fm*fm);
335  } else if(fNucFormfactor == fGaussianNF) {
336  fm = G4Exp(-2*formf*z1);
337  } else if(fNucFormfactor == fFlatNF) {
338  static const G4double ccoef = 0.00508/MeV;
339  G4double x = std::sqrt(2.*mom2*z1)*ccoef*2.;
340  fm = FlatFormfactor(x);
341  fm *= FlatFormfactor(x*0.6
342  *fG4pow->A13(fNistManager->GetAtomicMassAmu(targetZ) ) );
343  }
344 
345  G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2. - z1))
346  *fm*fm/(1.0 + z1*factD);
347 
348  if(rndmEngineMod->flat() <= grej) {
349  // exclude "false" scattering due to formfactor and spin effect
350  G4double cost = 1.0 - z1;
351  if(cost > 1.0) { cost = 1.0; }
352  else if(cost < -1.0) { cost =-1.0; }
353  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
354  //G4cout << "sint= " << sint << G4endl;
355  G4double phi = twopi*rndmEngineMod->flat();
356  temp.set(sint*cos(phi),sint*sin(phi),cost);
357  }
358  }
359  return temp;
360 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
fFlatNF
Definition: G4NuclearFormfactorType.hh:52
CLHEP::HepRandomEngine::flat
virtual double flat()=0
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
fGaussianNF
Definition: G4NuclearFormfactorType.hh:51
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4Pow::A13
G4double A13(G4double A) const
Definition: G4Pow.hh:132
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:359
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
fExponentialNF
Definition: G4NuclearFormfactorType.hh:50

Here is the call graph for this function:

Here is the caller graph for this function:

void G4WentzelOKandVIxSection::SetTargetMass ( G4double  value)
inline

Definition at line 206 of file G4WentzelOKandVIxSection.hh.

207 {
208  targetMass = value;
209  factD = std::sqrt(mom2)/value;
210 }
value
const XML_Char int const XML_Char * value
Definition: expat.h:331

Here is the caller graph for this function:

G4double G4WentzelOKandVIxSection::SetupKinematic ( G4double  kinEnergy,
const G4Material *  mat 
)
inline

Definition at line 187 of file G4WentzelOKandVIxSection.hh.

188 {
189  if(ekin != tkin || mat != currentMaterial) {
190  currentMaterial = mat;
191  tkin = ekin;
192  mom2 = tkin*(tkin + 2.0*mass);
193  invbeta2 = 1.0 + mass*mass/mom2;
194  factB = spin/invbeta2;
195  cosTetMaxNuc = cosThetaMax;
196  if(isCombined) {
197  cosTetMaxNuc = std::max(cosTetMaxNuc,
198  1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2);
199  }
200  }
201  return cosTetMaxNuc;
202 }
G4Material::GetIonisation
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
G4IonisParamMat::GetInvA23
G4double GetInvA23() const
Definition: G4IonisParamMat.hh:129
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112

Here is the call graph for this function:

Here is the caller graph for this function:

void G4WentzelOKandVIxSection::SetupParticle ( const G4ParticleDefinition *  p)

Definition at line 153 of file G4WentzelOKandVIxSection.cc.

154 {
155  particle = p;
156  mass = particle->GetPDGMass();
157  spin = particle->GetPDGSpin();
158  if(0.0 != spin) { spin = 0.5; }
159  G4double q = std::fabs(particle->GetPDGCharge()/eplus);
160  chargeSquare = q*q;
161  charge3 = chargeSquare*q;
162  tkin = 0.0;
163  currentMaterial = 0;
164  targetZ = 0;
165 }
p
const char * p
Definition: xmltok.h:285
eplus
static constexpr double eplus
Definition: G4SIunits.hh:199
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
G4ParticleDefinition::GetPDGSpin
G4double GetPDGSpin() const
Definition: G4ParticleDefinition.hh:126
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4WentzelOKandVIxSection::SetupTarget ( G4int  Z,
G4double  cut = DBL_MAX 
)

Definition at line 170 of file G4WentzelOKandVIxSection.cc.

171 {
172  G4double cosTetMaxNuc2 = cosTetMaxNuc;
173  if(Z != targetZ || tkin != etag) {
174  etag = tkin;
175  targetZ = Z;
176  if(targetZ > 99) { targetZ = 99; }
177  G4double massT = proton_mass_c2;
178  if(targetZ > 1) {
179  massT = fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2;
180  }
181  SetTargetMass(massT);
182 
183  kinFactor = coeff*Z*chargeSquare*invbeta2/mom2;
184 
185  if(1 == Z) {
186  screenZ = ScreenRSquare[targetZ]/mom2;
187  } else if(mass > MeV) {
188  screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
189  ScreenRSquare[targetZ]/mom2;
190  } else {
191  G4double tau = tkin/mass;
192  screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z
193  *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))*
194  ScreenRSquareElec[targetZ]/mom2;
195  }
196  if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
197  cosTetMaxNuc2 = 0.0;
198  }
199  formfactA = FormFactor[targetZ]*mom2;
200 
201  cosTetMaxElec = 1.0;
202  ComputeMaxElectronScattering(cut);
203  }
204  return cosTetMaxNuc2;
205 }
CLHEP::proton_mass_c2
static constexpr double proton_mass_c2
Definition: PhysicalConstants.h:80
CLHEP::amu_c2
static constexpr double amu_c2
Definition: PhysicalConstants.h:82
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:359
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4Pow::Z23
G4double Z23(G4int Z) const
Definition: G4Pow.hh:154
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

Here is the caller graph for this function:

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