Geant4  9.6.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$
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 
49 #include "G4WentzelVIRelXSection.hh"
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 "G4LossTableManager.hh"
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 
60 G4double G4WentzelVIRelXSection::ScreenRSquare[] = {0.0};
61 G4double G4WentzelVIRelXSection::FormFactor[] = {0.0};
62 
63 using namespace std;
64 
65 G4WentzelVIRelXSection::G4WentzelVIRelXSection() :
66  numlimit(0.1),
67  nwarnings(0),
68  nwarnlimit(50),
69  alpha2(fine_structure_const*fine_structure_const)
70 {
71  fNistManager = G4NistManager::Instance();
72  fG4pow = G4Pow::GetInstance();
73  theElectron = G4Electron::Electron();
74  thePositron = G4Positron::Positron();
75  theProton = G4Proton::Proton();
76  lowEnergyLimit = 1.0*eV;
77  G4double p0 = electron_mass_c2*classic_electr_radius;
78  coeff = twopi*p0*p0;
79  particle = 0;
80 
81  // Thomas-Fermi screening radii
82  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
83 
84  if(0.0 == ScreenRSquare[0]) {
85  G4double a0 = electron_mass_c2/0.88534;
86  G4double constn = 6.937e-6/(MeV*MeV);
87 
88  ScreenRSquare[0] = alpha2*a0*a0;
89  for(G4int j=1; j<100; ++j) {
90  G4double x = a0*fG4pow->Z13(j);
91  //ScreenRSquare[j] = 0.5*(1 + exp(-j*j*0.001))*alpha2*x*x;
92  ScreenRSquare[j] = 0.5*alpha2*x*x;
93  x = fNistManager->GetA27(j);
94  FormFactor[j] = constn*x*x;
95  }
96  }
97  currentMaterial = 0;
98  elecXSRatio = factB = factD = formfactA = screenZ = 0.0;
99  cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;
100 
101  factB1= 0.5*CLHEP::pi*fine_structure_const;
102 
103  Initialise(theElectron, 1.0);
104  targetMass = proton_mass_c2;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
109 G4WentzelVIRelXSection::~G4WentzelVIRelXSection()
110 {}
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113 
114 void G4WentzelVIRelXSection::Initialise(const G4ParticleDefinition* p,
115  G4double CosThetaLim)
116 {
117  SetupParticle(p);
118  tkin = mom2 = momCM2 = 0.0;
119  ecut = etag = DBL_MAX;
120  targetZ = 0;
121  cosThetaMax = CosThetaLim;
122  G4double a =
123  G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
124  factorA2 = 0.5*a*a;
125  currentMaterial = 0;
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 
130 void G4WentzelVIRelXSection::SetupParticle(const G4ParticleDefinition* p)
131 {
132  particle = p;
133  mass = particle->GetPDGMass();
134  spin = particle->GetPDGSpin();
135  if(0.0 != spin) { spin = 0.5; }
136  G4double q = std::fabs(particle->GetPDGCharge()/eplus);
137  chargeSquare = q*q;
138  charge3 = chargeSquare*q;
139  tkin = 0.0;
140  currentMaterial = 0;
141  targetZ = 0;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 
146 G4double
147 G4WentzelVIRelXSection::SetupTarget(G4int Z, G4double cut)
148 {
149  G4double cosTetMaxNuc2 = cosTetMaxNuc;
150  if(Z != targetZ || tkin != etag) {
151  etag = tkin;
152  targetZ = Z;
153  if(targetZ > 99) { targetZ = 99; }
154  SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2);
155  //G4double tmass2 = targetMass*targetMass;
156  //G4double etot = tkin + mass;
157  //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass;
158  //momCM2 = mom2*tmass2/invmass2;
159  //gam0pcmp = (etot + targetMass)*targetMass/invmass2;
160  //pcmp2 = tmass2/invmass2;
161 
162  kinFactor = coeff*targetZ*chargeSquare*invbeta2/mom2;
163 
164  screenZ = ScreenRSquare[targetZ]/mom2;
165  if(Z > 1) {
166  screenZ *= std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare);
167  /*
168  if(mass > MeV) {
169  screenZ *= std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare);
170  } else {
171  G4double tau = tkin/mass;
172  // screenZ *= std::min(Z*invbeta2,
173  screenZ *= std::min(Z*1.13,
174  (1.13 +3.76*Z*Z*invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(Z)))));
175  }
176  */
177  }
178  if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
179  cosTetMaxNuc2 = 0.0;
180  }
181  formfactA = FormFactor[targetZ]*mom2;
182 
183  cosTetMaxElec = 1.0;
184  ComputeMaxElectronScattering(cut);
185  }
186  return cosTetMaxNuc2;
187 }
188 
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
190 
191 G4double
192 G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax)
193 {
194  G4double xsec = 0.0;
195  if(cosTMax >= 1.0) { return xsec; }
196 
197  G4double xSection = 0.0;
198  G4double x = 0;
199  G4double y = 0;
200  G4double x1= 0;
201  G4double x2= 0;
202  G4double xlog = 0.0;
203 
204  G4double costm = std::max(cosTMax,cosTetMaxElec);
205  G4double fb = screenZ*factB;
206 
207  // scattering off electrons
208  if(costm < 1.0) {
209  x = (1.0 - costm)/screenZ;
210  if(x < numlimit) {
211  x2 = 0.5*x*x;
212  y = x2*(1.0 - 1.3333333*x + 3*x2);
213  if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
214  } else {
215  x1= x/(1 + x);
216  xlog = log(1.0 + x);
217  y = xlog - x1;
218  if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
219  }
220 
221  if(y < 0.0) {
222  ++nwarnings;
223  if(nwarnings < nwarnlimit) {
224  G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
225  << G4endl;
226  G4cout << "y= " << y
227  << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
228  << " Z= " << targetZ << " "
229  << particle->GetParticleName() << G4endl;
230  G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
231  << " x= " << x << G4endl;
232  }
233  y = 0.0;
234  }
235  xSection = y;
236  }
237  /*
238  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
239  << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
240  << " cut(MeV)= " << ecut/MeV
241  << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
242  << " zmaxN= " << (1.0 - cosThetaMax)/screenZ
243  << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
244  */
245  // scattering off nucleus
246  if(cosTMax < 1.0) {
247  x = (1.0 - cosTMax)/screenZ;
248  if(x < numlimit) {
249  x2 = 0.5*x*x;
250  y = x2*(1.0 - 1.3333333*x + 3*x2);
251  if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
252  } else {
253  x1= x/(1 + x);
254  xlog = log(1.0 + x);
255  y = xlog - x1;
256  if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
257  }
258 
259  if(y < 0.0) {
260  ++nwarnings;
261  if(nwarnings < nwarnlimit) {
262  G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
263  << G4endl;
264  G4cout << "y= " << y
265  << " e(MeV)= " << tkin << " Z= " << targetZ << " "
266  << particle->GetParticleName() << G4endl;
267  G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
268  << " x= " << " x1= " << x1 <<G4endl;
269  }
270  y = 0.0;
271  }
272  xSection += y*targetZ;
273  }
274  xSection *= kinFactor;
275  /*
276  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
277  << " screenZ= " << screenZ << " formF= " << formfactA
278  << " for " << particle->GetParticleName()
279  << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
280  << " x= " << x
281  << G4endl;
282  */
283  return xSection;
284 }
285 
286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
287 
288 G4ThreeVector
289 G4WentzelVIRelXSection::SampleSingleScattering(G4double cosTMin,
290  G4double cosTMax,
291  G4double elecRatio)
292 {
293  G4ThreeVector v(0.0,0.0,1.0);
294 
295  G4double formf = formfactA;
296  G4double cost1 = cosTMin;
297  G4double cost2 = cosTMax;
298  if(elecRatio > 0.0) {
299  if(G4UniformRand() <= elecRatio) {
300  formf = 0.0;
301  cost1 = std::max(cost1,cosTetMaxElec);
302  cost2 = std::max(cost2,cosTetMaxElec);
303  }
304  }
305  if(cost1 < cost2) { return v; }
306 
307  G4double w1 = 1. - cost1 + screenZ;
308  G4double w2 = 1. - cost2 + screenZ;
309  G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;
310 
311  //G4double fm = 1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass);
312  G4double fm = 1.0 + formf*z1;
313  //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm );
314  G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm );
315  // "false" scattering
316  if( G4UniformRand() > grej ) { return v; }
317  // }
318  G4double cost = 1.0 - z1;
319 
320  if(cost > 1.0) { cost = 1.0; }
321  else if(cost < -1.0) { cost =-1.0; }
322  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
323  //G4cout << "sint= " << sint << G4endl;
324  G4double phi = twopi*G4UniformRand();
325  G4double vx1 = sint*cos(phi);
326  G4double vy1 = sint*sin(phi);
327 
328  // only direction is changed
329  v.set(vx1,vy1,cost);
330  return v;
331 }
332 
333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
334 
335 void
336 G4WentzelVIRelXSection::ComputeMaxElectronScattering(G4double cutEnergy)
337 {
338  if(mass > MeV) {
339  G4double ratio = electron_mass_c2/mass;
340  G4double tau = tkin/mass;
341  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
342  (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
343  //tmax = std::min(tmax, targetZ*targetZ*10*eV);
344  cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
345  } else {
346 
347  G4double tmax = tkin;
348  if(particle == theElectron) { tmax *= 0.5; }
349  //tmax = std::min(tmax, targetZ*targetZ*10*eV);
350  G4double t = std::min(cutEnergy, tmax);
351  G4double mom21 = t*(t + 2.0*electron_mass_c2);
352  G4double t1 = tkin - t;
353  //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
354  //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
355  if(t1 > 0.0) {
356  G4double mom22 = t1*(t1 + 2.0*mass);
357  G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
358  if(ctm < 1.0) { cosTetMaxElec = ctm; }
359  if(particle == theElectron && cosTetMaxElec < 0.0) { cosTetMaxElec = 0.0; }
360  }
361  }
362 }
363 
364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:33:38 for Geant4 by   doxygen 1.8.4