ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4WentzelOKandVIxSection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4WentzelOKandVIxSection
33 //
34 // Author: V.Ivanchenko
35 //
36 // Creation date: 09.04.2008 from G4MuMscModel
37 //
38 // Modifications:
39 //
40 // -------------------------------------------------------------------
41 //
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 
46 #include "G4WentzelOKandVIxSection.hh"
47 #include "G4ScreeningMottCrossSection.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "Randomize.hh"
51 #include "G4Electron.hh"
52 #include "G4Positron.hh"
53 #include "G4Proton.hh"
54 #include "G4EmParameters.hh"
55 #include "G4Log.hh"
56 #include "G4Exp.hh"
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 
60 G4double G4WentzelOKandVIxSection::ScreenRSquareElec[] = {0.0};
61 G4double G4WentzelOKandVIxSection::ScreenRSquare[] = {0.0};
62 G4double G4WentzelOKandVIxSection::FormFactor[] = {0.0};
63 
64 #ifdef G4MULTITHREADED
65 G4Mutex G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex = G4MUTEX_INITIALIZER;
66 #endif
67 
68 using namespace std;
69 
70 G4WentzelOKandVIxSection::G4WentzelOKandVIxSection(G4bool comb) :
71  temp(0.,0.,0.),
72  numlimit(0.1),
73  nwarnings(0),
74  nwarnlimit(50),
75  isCombined(comb),
76  cosThetaMax(-1.0),
77  alpha2(fine_structure_const*fine_structure_const)
78 {
79  fNistManager = G4NistManager::Instance();
80  fG4pow = G4Pow::GetInstance();
81  fMottXSection = nullptr;
82 
83  theElectron = G4Electron::Electron();
84  thePositron = G4Positron::Positron();
85  theProton = G4Proton::Proton();
86 
87  lowEnergyLimit = 1.0*eV;
88  G4double p0 = electron_mass_c2*classic_electr_radius;
89  coeff = twopi*p0*p0;
90  particle = nullptr;
91 
92  fNucFormfactor = fExponentialNF;
93 
94  currentMaterial = nullptr;
95  factB = factD = formfactA = screenZ = 0.0;
96  cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = fMottFactor
97  = gam0pcmp = pcmp2 = 1.0;
98 
99  factB1= 0.5*CLHEP::pi*fine_structure_const;
100 
101  tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
102  ecut = etag = DBL_MAX;
103  targetZ = 0;
104  targetMass = CLHEP::proton_mass_c2;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
109 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection()
110 {
111  delete fMottXSection;
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115 
116 void G4WentzelOKandVIxSection::Initialise(const G4ParticleDefinition* p,
117  G4double cosThetaLim)
118 {
119  SetupParticle(p);
120  tkin = mom2 = momCM2 = 0.0;
121  ecut = etag = DBL_MAX;
122  targetZ = 0;
123 
124  // cosThetaMax is below 1.0 only when MSC is combined with SS
125  if(isCombined) { cosThetaMax = cosThetaLim; }
126  G4EmParameters* param = G4EmParameters::Instance();
127  G4double a = param->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
128  factorA2 = 0.5*a*a;
129  currentMaterial = nullptr;
130 
131  fNucFormfactor = param->NuclearFormfactorType();
132  if(0.0 == ScreenRSquare[0]) { InitialiseA(); }
133 
134  // Mott corrections always added
135  if((p == theElectron || p == thePositron) && !fMottXSection) {
136  fMottXSection = new G4ScreeningMottCrossSection();
137  fMottXSection->Initialise(p, 1.0);
138  }
139  /*
140  G4cout << "G4WentzelOKandVIxSection::Initialise for "
141  << p->GetParticleName() << " cosThetaMax= " << cosThetaMax
142  << " " << ScreenRSquare[0] << " coeff= " << coeff << G4endl;
143  */
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147 
148 void G4WentzelOKandVIxSection::InitialiseA()
149 {
150  // Thomas-Fermi screening radii
151  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
152 #ifdef G4MULTITHREADED
153  G4MUTEXLOCK(&G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex);
154 #endif
155  if(0.0 == ScreenRSquare[0]) {
156  G4double a0 = electron_mass_c2/0.88534;
157  G4double constn = 6.937e-6/(MeV*MeV);
158  G4double fct = G4EmParameters::Instance()->ScreeningFactor();
159 
160  G4double afact = 0.5*fct*alpha2*a0*a0;
161  ScreenRSquare[0] = afact;
162  ScreenRSquare[1] = afact;
163  ScreenRSquareElec[1] = afact;
164  FormFactor[1] = 3.097e-6/(MeV*MeV);
165 
166  for(G4int j=2; j<100; ++j) {
167  G4double x = fG4pow->Z13(j);
168  ScreenRSquare[j] = afact*(1 + G4Exp(-j*j*0.001))*x*x;
169  ScreenRSquareElec[j] = afact*x*x;
170  x = fNistManager->GetA27(j);
171  FormFactor[j] = constn*x*x;
172  }
173  }
174 #ifdef G4MULTITHREADED
175  G4MUTEXUNLOCK(&G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex);
176 #endif
177 
178  //G4cout << "G4WentzelOKandVIxSection::Initialise mass= " << mass
179  // << " " << p->GetParticleName()
180  // << " cosThetaMax= " << cosThetaMax << G4endl;
181 
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 
186 void G4WentzelOKandVIxSection::SetupParticle(const G4ParticleDefinition* p)
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
203 G4WentzelOKandVIxSection::SetupKinematic(G4double ekin, const G4Material* mat)
204 {
205  if(ekin != tkin || mat != currentMaterial) {
206  currentMaterial = mat;
207  tkin = ekin;
208  mom2 = tkin*(tkin + 2.0*mass);
209  invbeta2 = 1.0 + mass*mass/mom2;
210  factB = spin/invbeta2;
211  cosTetMaxNuc = isCombined ?
212  std::max(cosThetaMax, 1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2)
213  : cosThetaMax;
214  }
215  return cosTetMaxNuc;
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
219 
220 G4double
221 G4WentzelOKandVIxSection::SetupTarget(G4int Z, G4double cut)
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 :
228  fNistManager->GetAtomicMassAmu(Z)*CLHEP::amu_c2;
229  SetTargetMass(massT);
230 
231  kinFactor = coeff*Z*chargeSquare*invbeta2/mom2;
232  if(particle == theElectron && fMottXSection) {
233  fMottFactor = (1.0 + 2.0e-4*Z*Z);
234  }
235 
236  if(1 == Z) {
237  screenZ = ScreenRSquare[targetZ]/mom2;
238  } else if(mass > MeV) {
239  screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
240  ScreenRSquare[targetZ]/mom2;
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)))))*
245  ScreenRSquareElec[targetZ]/mom2;
246  }
247  if(targetZ == 1 && particle == theProton && cosTetMaxNuc2 < 0.0) {
248  cosTetMaxNuc2 = 0.0;
249  }
250  formfactA = FormFactor[targetZ]*mom2;
251 
252  cosTetMaxElec = 1.0;
253  ComputeMaxElectronScattering(cut);
254  }
255  //G4cout << "SetupTarget: Z= " << targetZ << " kinFactor= " << kinFactor
256  // << " fMottFactor= " << fMottFactor << " screenZ= " << screenZ <<G4endl;
257  return cosTetMaxNuc2;
258 }
259 
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261 
262 G4double
263 G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax)
264 {
265  G4double xSection = 0.0;
266  if(cosTMax >= 1.0) { return xSection; }
267 
268  G4double costm = std::max(cosTMax,cosTetMaxElec);
269  G4double fb = screenZ*factB;
270 
271  // scattering off electrons
272  if(costm < 1.0) {
273  G4double x = (1.0 - costm)/screenZ;
274  if(x < numlimit) {
275  G4double x2 = 0.5*x*x;
276  xSection = x2*((1.0 - 1.3333333*x + 3*x2) - fb*x*(0.6666667 - x));
277  } else {
278  G4double x1= x/(1 + x);
279  G4double xlog = G4Log(1.0 + x);
280  xSection = xlog - x1 - fb*(x + x1 - 2*xlog);
281  }
282 
283  if(xSection < 0.0) {
284  ++nwarnings;
285  if(nwarnings < nwarnlimit) {
286  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom"
287  << " scattering on e- <0"
288  << G4endl;
289  G4cout << "cross= " << xSection
290  << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
291  << " Z= " << targetZ << " "
292  << particle->GetParticleName() << G4endl;
293  G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
294  << " x= " << x << G4endl;
295  }
296  xSection = 0.0;
297  }
298  }
299  /*
300  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom: \n"
301  << " Z= " << targetZ
302  << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
303  << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
304  << " zmaxN= " << (1.0 - cosThetaMax)/screenZ
305  << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
306  */
307  // scattering off nucleus
308  if(cosTMax < 1.0) {
309  G4double x = (1.0 - cosTMax)/screenZ;
310  G4double y;
311  if(x < numlimit) {
312  G4double x2 = 0.5*x*x;
313  y = x2*((1.0 - 1.3333333*x + 3*x2) - fb*x*(0.6666667 - x));
314  } else {
315  G4double x1= x/(1 + x);
316  G4double xlog = G4Log(1.0 + x);
317  y = xlog - x1 - fb*(x + x1 - 2*xlog);
318  }
319 
320  if(y < 0.0) {
321  ++nwarnings;
322  if(nwarnings < nwarnlimit) {
323  G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom"
324  << " scattering on nucleus <0"
325  << G4endl;
326  G4cout << "y= " << y
327  << " e(MeV)= " << tkin << " Z= " << targetZ << " "
328  << particle->GetParticleName() << G4endl;
329  G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
330  << " x= " << x <<G4endl;
331  }
332  y = 0.0;
333  }
334  xSection += y*targetZ;
335  }
336  xSection *= kinFactor;
337 
338  /*
339  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
340  << " screenZ= " << screenZ << " formF= " << formfactA
341  << " for " << particle->GetParticleName()
342  << " m= " << mass << " 1/v= " << sqrt(invbeta2)
343  << " p= " << sqrt(mom2)
344  << " x= " << x << G4endl;
345  */
346  return xSection;
347 }
348 
349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
350 
351 G4ThreeVector&
352 G4WentzelOKandVIxSection::SampleSingleScattering(G4double cosTMin,
353  G4double cosTMax,
354  G4double elecRatio)
355 {
356  temp.set(0.0,0.0,1.0);
357  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
358 
359  G4double formf = formfactA;
360  G4double cost1 = cosTMin;
361  G4double cost2 = cosTMax;
362  if(elecRatio > 0.0) {
363  if(rndmEngineMod->flat() <= elecRatio) {
364  formf = 0.0;
365  cost1 = std::max(cost1,cosTetMaxElec);
366  cost2 = std::max(cost2,cosTetMaxElec);
367  }
368  }
369  if(cost1 > cost2) {
370 
371  G4double w1 = 1. - cost1 + screenZ;
372  G4double w2 = 1. - cost2 + screenZ;
373  G4double z1 = w1*w2/(w1 + rndmEngineMod->flat()*(w2 - w1)) - screenZ;
374 
375  G4double fm = 1.0;
376  if(fNucFormfactor == fExponentialNF) {
377  fm += formf*z1;
378  fm = 1.0/(fm*fm);
379  } else if(fNucFormfactor == fGaussianNF) {
380  fm = G4Exp(-2*formf*z1);
381  } else if(fNucFormfactor == fFlatNF) {
382  static const G4double ccoef = 0.00508/MeV;
383  G4double x = std::sqrt(2.*mom2*z1)*ccoef*2.;
384  fm = FlatFormfactor(x);
385  fm *= FlatFormfactor(x*0.6
386  *fG4pow->A13(fNistManager->GetAtomicMassAmu(targetZ)));
387  }
388  G4double grej;
389  if(fMottXSection) {
390  fMottXSection->SetupKinematic(tkin, targetZ);
391  grej = fMottXSection->RatioMottRutherfordCosT(std::sqrt(z1))*fm*fm;
392  } else {
393  grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2. - z1))
394  *fm*fm/(1.0 + z1*factD);
395  }
396  // G4cout << "SampleSingleScattering: E= " << tkin << " z1= "
397  // << z1 << " grej= "<< grej << " mottFact= "<< fMottFactor<< G4endl;
398  if(fMottFactor*rndmEngineMod->flat() <= grej ) {
399  // exclude "false" scattering due to formfactor and spin effect
400  G4double cost = 1.0 - z1;
401  if(cost > 1.0) { cost = 1.0; }
402  else if(cost < -1.0) { cost =-1.0; }
403  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
404  //G4cout << "sint= " << sint << G4endl;
405  G4double phi = twopi*rndmEngineMod->flat();
406  temp.set(sint*cos(phi),sint*sin(phi),cost);
407  }
408  }
409  return temp;
410 }
411 
412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
413 
414 void
415 G4WentzelOKandVIxSection::ComputeMaxElectronScattering(G4double cutEnergy)
416 {
417  if(mass > MeV) {
418  G4double ratio = electron_mass_c2/mass;
419  G4double tau = tkin/mass;
420  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
421  (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
422  cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
423  } else {
424 
425  G4double tmax = (particle == theElectron) ? 0.5*tkin : tkin;
426  G4double t = std::min(cutEnergy, tmax);
427  G4double mom21 = t*(t + 2.0*electron_mass_c2);
428  G4double t1 = tkin - t;
429  //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
430  //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
431  if(t1 > 0.0) {
432  G4double mom22 = t1*(t1 + 2.0*mass);
433  G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
434  if(ctm < 1.0) { cosTetMaxElec = ctm; }
435  if(particle == theElectron && cosTetMaxElec < 0.0) {
436  cosTetMaxElec = 0.0;
437  }
438  }
439  }
440 }
441 
442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
443 
444 G4double
445 G4WentzelOKandVIxSection::ComputeSecondTransportMoment(G4double /*CosThetaMax*/)
446 {
447  return 0.0;
448 }
449 
450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......