source: trunk/source/processes/electromagnetic/standard/src/G4WentzelOKandVIxSection.cc@ 1330

Last change on this file since 1330 was 1316, checked in by garnier, 15 years ago

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