Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4WentzelVIRelXSection.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: G4WentzelVIRelXSection.cc 96934 2016-05-18 09:10:41Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4WentzelVIRelXSection
34 //
35 // Author: V.Ivanchenko
36 //
37 // Creation date: 08.06.2012 from G4WentzelOKandVIxSection
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 "G4EmParameters.hh"
57 #include "G4Log.hh"
58 #include "G4Exp.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 G4double G4WentzelVIRelXSection::ScreenRSquare[] = {0.0};
63 G4double G4WentzelVIRelXSection::FormFactor[] = {0.0};
64 
65 using namespace std;
66 
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 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115 
117 {}
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120 
122  G4double cosThetaLim)
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 }
137 
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139 
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 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155 
156 G4double
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 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 
183 G4double
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 }
277 
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 
282  G4double cosTMax,
283  G4double elecRatio)
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 }
325 
326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
327 
328 void
329 G4WentzelVIRelXSection::ComputeMaxElectronScattering(G4double cutEnergy)
330 {
331  if(mass > MeV) {
332  G4double ratio = electron_mass_c2/mass;
333  G4double tau = tkin/mass;
334  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
335  (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
336  //tmax = std::min(tmax, targetZ*targetZ*10*eV);
337  cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
338  } else {
339 
340  G4double tmax = tkin;
341  if(particle == theElectron) { tmax *= 0.5; }
342  //tmax = std::min(tmax, targetZ*targetZ*10*eV);
343  G4double t = std::min(cutEnergy, tmax);
344  G4double mom21 = t*(t + 2.0*electron_mass_c2);
345  G4double t1 = tkin - t;
346  //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
347  //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
348  if(t1 > 0.0) {
349  G4double mom22 = t1*(t1 + 2.0*mass);
350  G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
351  if(ctm < 1.0) { cosTetMaxElec = ctm; }
352  if(particle == theElectron && cosTetMaxElec < 0.0) {
353  cosTetMaxElec = 0.0;
354  }
355  }
356  }
357 }
358 
359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void set(double x, double y, double z)
const G4double a0
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4WentzelVIRelXSection(G4bool combined=true)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static constexpr double hbarc
const char * p
Definition: xmltok.h:285
G4double GetA27(G4int Z) const
tuple x
Definition: test.py:50
virtual double flat()=0
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
const G4String & GetParticleName() const
static constexpr double twopi
Definition: G4SIunits.hh:76
void SetupParticle(const G4ParticleDefinition *)
G4GLOB_DLL std::ostream G4cout
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
G4double SetupTarget(G4int Z, G4double cut)
bool G4bool
Definition: G4Types.hh:79
static constexpr double eplus
Definition: G4SIunits.hh:199
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
G4double G4Log(G4double x)
Definition: G4Log.hh:230
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
tuple t1
Definition: plottest35.py:33
static G4EmParameters * Instance()
static constexpr double fermi
Definition: SystemOfUnits.h:83
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetPDGSpin() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
void Initialise(const G4ParticleDefinition *, G4double CosThetaLim)
double G4double
Definition: G4Types.hh:76
G4ThreeVector & SampleSingleScattering(G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio)
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
static constexpr double pi
Definition: SystemOfUnits.h:54
G4double FactorForAngleLimit() const