Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4WentzelVIRelXSection Class Reference

#include <G4WentzelVIRelXSection.hh>

Public Member Functions

 G4WentzelVIRelXSection (G4bool combined=true)
 
virtual ~G4WentzelVIRelXSection ()
 
void Initialise (const G4ParticleDefinition *, G4double CosThetaLim)
 
void SetupParticle (const G4ParticleDefinition *)
 
G4double SetupTarget (G4int Z, G4double cut)
 
G4double ComputeTransportCrossSectionPerAtom (G4double CosThetaMax)
 
G4ThreeVectorSampleSingleScattering (G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio)
 
G4double ComputeNuclearCrossSection (G4double CosThetaMin, G4double CosThetaMax)
 
G4double ComputeElectronCrossSection (G4double CosThetaMin, G4double CosThetaMax)
 
G4double SetupKinematic (G4double kinEnergy, const G4Material *mat, G4double cut, G4double tmass)
 
G4double GetMomentumSquare () const
 
G4double GetCosThetaNuc () const
 
G4double GetCosThetaElec () const
 

Detailed Description

Definition at line 71 of file G4WentzelVIRelXSection.hh.

Constructor & Destructor Documentation

G4WentzelVIRelXSection::G4WentzelVIRelXSection ( G4bool  combined = true)
explicit

Definition at line 67 of file G4WentzelVIRelXSection.cc.

67  :
68  temp(0.,0.,0.),
69  numlimit(0.1),
70  nwarnings(0),
71  nwarnlimit(50),
72  isCombined(combined),
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;
82  coeff = twopi*p0*p0;
83  particle = nullptr;
84 
85  // Thomas-Fermi screening radii
86  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
87 
88  if(0.0 == ScreenRSquare[0]) {
89  G4double a0 = electron_mass_c2/0.88534;
90  G4double constn = 6.937e-6/(MeV*MeV);
91 
92  ScreenRSquare[0] = alpha2*a0*a0;
93  for(G4int j=1; j<100; ++j) {
94  G4double x = a0*fG4pow->Z13(j);
95  //ScreenRSquare[j] = 0.5*(1 + exp(-j*j*0.001))*alpha2*x*x;
96  ScreenRSquare[j] = 0.5*alpha2*x*x;
97  x = fNistManager->GetA27(j);
98  FormFactor[j] = constn*x*x;
99  }
100  }
101  currentMaterial = 0;
102  factB = factD = formfactA = screenZ = 0.0;
103  cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;
104  cosThetaMax = 1.0;
105 
106  factB1= 0.5*CLHEP::pi*fine_structure_const;
107 
108  tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
109  ecut = etag = DBL_MAX;
110  targetZ = 0;
111  targetMass = proton_mass_c2;
112 }
const G4double a0
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double GetA27(G4int Z) const
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
static G4Proton * Proton()
Definition: G4Proton.cc:93
static constexpr double eV
Definition: G4SIunits.hh:215
float proton_mass_c2
Definition: hepunit.py:275
float electron_mass_c2
Definition: hepunit.py:274
static G4Positron * Positron()
Definition: G4Positron.cc:94
static G4Electron * Electron()
Definition: G4Electron.cc:94
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

G4WentzelVIRelXSection::~G4WentzelVIRelXSection ( )
virtual

Definition at line 116 of file G4WentzelVIRelXSection.cc.

117 {}

Member Function Documentation

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

Definition at line 259 of file G4WentzelVIRelXSection.hh.

261 {
262  G4double xsec = 0.0;
263  G4double cost1 = std::max(cosTMin,cosTetMaxElec);
264  G4double cost2 = std::max(cosTMax,cosTetMaxElec);
265  if(cost1 > cost2) {
266  xsec = kinFactor*(cost1 - cost2)/
267  ((1.0 - cost1 + screenZ)*(1.0 - cost2 + screenZ));
268  }
269  return xsec;
270 }
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 245 of file G4WentzelVIRelXSection.hh.

247 {
248  G4double xsec = 0.0;
249  if(cosTMax < cosTMin) {
250  xsec = targetZ*kinFactor*(cosTMin - cosTMax)/
251  ((1.0 - cosTMin + screenZ)*(1.0 - cosTMax + screenZ));
252  }
253  return xsec;
254 }
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom ( G4double  CosThetaMax)

Definition at line 184 of file G4WentzelVIRelXSection.cc.

185 {
186  G4double xsec = 0.0;
187  if(cosTMax >= 1.0) { return xsec; }
188 
189  G4double xSection = 0.0;
190  G4double x = 0;
191  G4double y = 0;
192  G4double x1= 0;
193  G4double x2= 0;
194  G4double xlog = 0.0;
195 
196  G4double costm = std::max(cosTMax,cosTetMaxElec);
197  G4double fb = screenZ*factB;
198 
199  // scattering off electrons
200  if(costm < 1.0) {
201  x = (1.0 - costm)/screenZ;
202  if(x < numlimit) {
203  x2 = 0.5*x*x;
204  y = x2*(1.0 - 1.3333333*x + 3*x2);
205  if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
206  } else {
207  x1= x/(1 + x);
208  xlog = G4Log(1.0 + x);
209  y = xlog - x1;
210  if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
211  }
212 
213  if(y < 0.0) {
214  ++nwarnings;
215  if(nwarnings < nwarnlimit) {
216  G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
217  << G4endl;
218  G4cout << "y= " << y
219  << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
220  << " Z= " << targetZ << " "
221  << particle->GetParticleName() << G4endl;
222  G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
223  << " x= " << x << G4endl;
224  }
225  y = 0.0;
226  }
227  xSection = y;
228  }
229  /*
230  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
231  << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
232  << " cut(MeV)= " << ecut/MeV
233  << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
234  << " zmaxN= " << (1.0 - cosThetaMax)/screenZ
235  << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
236  */
237  // scattering off nucleus
238  if(cosTMax < 1.0) {
239  x = (1.0 - cosTMax)/screenZ;
240  if(x < numlimit) {
241  x2 = 0.5*x*x;
242  y = x2*(1.0 - 1.3333333*x + 3*x2);
243  if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
244  } else {
245  x1= x/(1 + x);
246  xlog = G4Log(1.0 + x);
247  y = xlog - x1;
248  if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
249  }
250 
251  if(y < 0.0) {
252  ++nwarnings;
253  if(nwarnings < nwarnlimit) {
254  G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
255  << G4endl;
256  G4cout << "y= " << y
257  << " e(MeV)= " << tkin << " Z= " << targetZ << " "
258  << particle->GetParticleName() << G4endl;
259  G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
260  << " x= " << " x1= " << x1 <<G4endl;
261  }
262  y = 0.0;
263  }
264  xSection += y*targetZ;
265  }
266  xSection *= kinFactor;
267  /*
268  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
269  << " screenZ= " << screenZ << " formF= " << formfactA
270  << " for " << particle->GetParticleName()
271  << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
272  << " x= " << x
273  << G4endl;
274  */
275  return xSection;
276 }
tuple x
Definition: test.py:50
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
G4double G4Log(G4double x)
Definition: G4Log.hh:230
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

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4WentzelVIRelXSection::GetCosThetaElec ( ) const
inline

Definition at line 237 of file G4WentzelVIRelXSection.hh.

238 {
239  return cosTetMaxElec;
240 }
G4double G4WentzelVIRelXSection::GetCosThetaNuc ( ) const
inline

Definition at line 230 of file G4WentzelVIRelXSection.hh.

231 {
232  return cosTetMaxNuc;
233 }
G4double G4WentzelVIRelXSection::GetMomentumSquare ( ) const
inline

Definition at line 223 of file G4WentzelVIRelXSection.hh.

224 {
225  return mom2;
226 }
void G4WentzelVIRelXSection::Initialise ( const G4ParticleDefinition p,
G4double  CosThetaLim 
)

Definition at line 121 of file G4WentzelVIRelXSection.cc.

123 {
124  SetupParticle(p);
125  tkin = mom2 = momCM2 = 0.0;
126  ecut = etag = DBL_MAX;
127  targetZ = 0;
128 
129  // cosThetaMax is below 1.0 only when MSC is combined with SS
130  if(isCombined) { cosThetaMax = cosThetaLim; }
131 
134  factorA2 = 0.5*a*a;
135  currentMaterial = nullptr;
136 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static constexpr double hbarc
void SetupParticle(const G4ParticleDefinition *)
static G4EmParameters * Instance()
static constexpr double fermi
Definition: SystemOfUnits.h:83
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
G4double FactorForAngleLimit() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4ThreeVector & G4WentzelVIRelXSection::SampleSingleScattering ( G4double  CosThetaMin,
G4double  CosThetaMax,
G4double  elecRatio 
)

Definition at line 281 of file G4WentzelVIRelXSection.cc.

284 {
285  temp.set(0.0,0.0,1.0);
286 
287  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
288 
289  G4double formf = formfactA;
290  G4double cost1 = cosTMin;
291  G4double cost2 = cosTMax;
292  if(elecRatio > 0.0) {
293  if(rndmEngine->flat() <= elecRatio) {
294  formf = 0.0;
295  cost1 = std::max(cost1,cosTetMaxElec);
296  cost2 = std::max(cost2,cosTetMaxElec);
297  }
298  }
299  if(cost1 < cost2) { return temp; }
300 
301  G4double w1 = 1. - cost1 + screenZ;
302  G4double w2 = 1. - cost2 + screenZ;
303  G4double z1 = w1*w2/(w1 + rndmEngine->flat()*(w2 - w1)) - screenZ;
304 
305  G4double fm = 1.0 + formf*z1;
306  G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))
307  /( (1.0 + z1*factD)*fm*fm );
308  // "false" scattering
309  if(rndmEngine->flat() > grej ) { return temp; }
310  // }
311  G4double cost = 1.0 - z1;
312 
313  if(cost > 1.0) { cost = 1.0; }
314  else if(cost < -1.0) { cost =-1.0; }
315  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
316  //G4cout << "sint= " << sint << G4endl;
317  G4double phi = twopi*rndmEngine->flat();
318  G4double vx1 = sint*cos(phi);
319  G4double vy1 = sint*sin(phi);
320 
321  // only direction is changed
322  temp.set(vx1,vy1,cost);
323  return temp;
324 }
void set(double x, double y, double z)
virtual double flat()=0
static constexpr double twopi
Definition: G4SIunits.hh:76
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4WentzelVIRelXSection::SetupKinematic ( G4double  kinEnergy,
const G4Material mat,
G4double  cut,
G4double  tmass 
)
inline

Definition at line 180 of file G4WentzelVIRelXSection.hh.

184 {
185  if(kinEnergy != tkin || mat != currentMaterial ||
186  ecut != cut || tmass != targetMass) {
187 
188  currentMaterial = mat;
189  ecut = cut;
190  tkin = kinEnergy;
191  G4double momLab2 = tkin*(tkin + 2.0*mass);
192 
193  G4double etot = tkin + mass;
194  G4double ptot = std::sqrt(momLab2);
195  G4double m12 = mass*mass;
196 
197  targetMass = tmass;
198 
199  // relativistic reduced mass from publucation
200  // A.P. Martynenko, R.N. Faustov, Teoret. mat. Fiz. 64 (1985) 179
201 
202  //incident particle & target nucleus
203  G4double Ecm = std::sqrt(m12 + targetMass*targetMass + 2.0*etot*targetMass);
204  G4double mu_rel = mass*targetMass/Ecm;
205  G4double momCM = ptot*targetMass/Ecm;
206  // relative system
207  mom2 = momCM*momCM;
208  invbeta2 = 1.0 + mu_rel*mu_rel/mom2;
209 
210  factB = spin/invbeta2;
211  factD = std::sqrt(mom2)/tmass;
212  if(isCombined) {
213  G4double cost = 1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2;
214  if(cost > cosTetMaxNuc) { cosTetMaxNuc = cost; }
215  }
216  }
217  return cosTetMaxNuc;
218 
219 }
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
G4double GetInvA23() const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4WentzelVIRelXSection::SetupParticle ( const G4ParticleDefinition p)

Definition at line 140 of file G4WentzelVIRelXSection.cc.

141 {
142  particle = p;
143  mass = particle->GetPDGMass();
144  spin = particle->GetPDGSpin();
145  if(0.0 != spin) { spin = 0.5; }
146  G4double q = std::fabs(particle->GetPDGCharge()/eplus);
147  chargeSquare = q*q;
148  charge3 = chargeSquare*q;
149  tkin = 0.0;
150  currentMaterial = nullptr;
151  targetZ = 0;
152 }
const char * p
Definition: xmltok.h:285
static constexpr double eplus
Definition: G4SIunits.hh:199
G4double GetPDGMass() const
G4double GetPDGSpin() const
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4WentzelVIRelXSection::SetupTarget ( G4int  Z,
G4double  cut 
)

Definition at line 157 of file G4WentzelVIRelXSection.cc.

158 {
159  G4double cosTetMaxNuc2 = cosTetMaxNuc;
160  if(Z != targetZ || tkin != etag) {
161  etag = tkin;
162  targetZ = Z;
163 
164  kinFactor = coeff*targetZ*chargeSquare*invbeta2/mom2;
165 
166  screenZ = ScreenRSquare[targetZ]/mom2;
167  if(Z > 1) {
168  screenZ *= std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare);
169  }
170  if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
171  cosTetMaxNuc2 = 0.0;
172  }
173  formfactA = FormFactor[targetZ]*mom2;
174 
175  cosTetMaxElec = 1.0;
176  ComputeMaxElectronScattering(cut);
177  }
178  return cosTetMaxNuc2;
179 }
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double G4double
Definition: G4Types.hh:76

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: