| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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........oo
44 //....oooOO0OOooo........oooOO0OOooo........oo
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 #include "G4AutoLock.hh"
58
59 //....oooOO0OOooo........oooOO0OOooo........oo
60
61 G4double G4WentzelOKandVIxSection::ScreenRSqua
62 G4double G4WentzelOKandVIxSection::ScreenRSqua
63 G4double G4WentzelOKandVIxSection::FormFactor[
64
65 namespace
66 {
67 G4Mutex theWOKVIMutex = G4MUTEX_INITIALIZER;
68 }
69
70 const G4double alpha2 = CLHEP::fine_structure_
71 const G4double factB1= 0.5*CLHEP::pi*CLHEP::fi
72 const G4double numlimit = 0.1;
73 const G4int nwarnlimit = 50;
74
75 using namespace std;
76
77 G4WentzelOKandVIxSection::G4WentzelOKandVIxSec
78 temp(0.,0.,0.),
79 isCombined(comb)
80 {
81 fNistManager = G4NistManager::Instance();
82 fG4pow = G4Pow::GetInstance();
83
84 theElectron = G4Electron::Electron();
85 thePositron = G4Positron::Positron();
86 theProton = G4Proton::Proton();
87
88 G4double p0 = CLHEP::electron_mass_c2*CLHEP:
89 coeff = CLHEP::twopi*p0*p0;
90 targetMass = CLHEP::proton_mass_c2;
91 }
92
93 //....oooOO0OOooo........oooOO0OOooo........oo
94
95 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSe
96 {
97 delete fMottXSection;
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oo
101
102 void G4WentzelOKandVIxSection::Initialise(cons
103 G4do
104 {
105 SetupParticle(p);
106 tkin = mom2 = momCM2 = 0.0;
107 ecut = etag = DBL_MAX;
108 targetZ = 0;
109
110 // cosThetaMax is below 1.0 only when MSC is
111 if(isCombined) { cosThetaMax = cosThetaLim;
112 G4EmParameters* param = G4EmParameters::Inst
113 G4double a = param->FactorForAngleLimit()*CL
114 factorA2 = 0.5*a*a;
115 currentMaterial = nullptr;
116
117 fNucFormfactor = param->NuclearFormfactorTyp
118 if(0.0 == ScreenRSquare[0]) { InitialiseA();
119
120 // Mott corrections always added
121 if((p == theElectron || p == thePositron) &&
122 fMottXSection = new G4ScreeningMottCrossSe
123 fMottXSection->Initialise(p, 1.0);
124 }
125 /*
126 G4cout << "G4WentzelOKandVIxSection::Initial
127 << p->GetParticleName() << " cosThetaMax= "
128 << " " << ScreenRSquare[0] << " coeff= " <
129 */
130 }
131
132 //....oooOO0OOooo........oooOO0OOooo........oo
133
134 void G4WentzelOKandVIxSection::InitialiseA()
135 {
136 // Thomas-Fermi screening radii
137 // Formfactors from A.V. Butkevich et al., N
138 if(0.0 != ScreenRSquare[0]) { return; }
139 G4AutoLock l(&theWOKVIMutex);
140 if(0.0 == ScreenRSquare[0]) {
141 const G4double invmev2 = 1./(CLHEP::MeV*CL
142 G4double a0 = CLHEP::electron_mass_c2/0.88
143 G4double constn = 6.937e-6*invmev2;
144 G4double fct = G4EmParameters::Instance()-
145
146 G4double afact = 0.5*fct*alpha2*a0*a0;
147 ScreenRSquare[0] = afact;
148 ScreenRSquare[1] = afact;
149 ScreenRSquareElec[1] = afact;
150 FormFactor[1] = 3.097e-6*invmev2;
151
152 for(G4int j=2; j<100; ++j) {
153 G4double x = fG4pow->Z13(j);
154 ScreenRSquare[j] = afact*(1 + G4Exp(-j*j
155 ScreenRSquareElec[j] = afact*x*x;
156 x = fNistManager->GetA27(j);
157 FormFactor[j] = constn*x*x;
158 }
159 }
160 l.unlock();
161 }
162
163 //....oooOO0OOooo........oooOO0OOooo........oo
164
165 void G4WentzelOKandVIxSection::SetupParticle(c
166 {
167 particle = p;
168 mass = particle->GetPDGMass();
169 spin = particle->GetPDGSpin();
170 if(0.0 != spin) { spin = 0.5; }
171 G4double q = std::abs(particle->GetPDGCharge
172 chargeSquare = q*q;
173 charge3 = chargeSquare*q;
174 tkin = 0.0;
175 currentMaterial = nullptr;
176 targetZ = 0;
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oo
180
181 G4double
182 G4WentzelOKandVIxSection::SetupKinematic(G4dou
183 {
184 if(ekin != tkin || mat != currentMaterial) {
185 currentMaterial = mat;
186 tkin = ekin;
187 mom2 = tkin*(tkin + 2.0*mass);
188 invbeta2 = 1.0 + mass*mass/mom2;
189 factB = spin/invbeta2;
190 cosTetMaxNuc = isCombined ?
191 std::max(cosThetaMax, 1.-factorA2*mat->G
192 : cosThetaMax;
193 }
194 return cosTetMaxNuc;
195 }
196
197 //....oooOO0OOooo........oooOO0OOooo........oo
198
199 G4double
200 G4WentzelOKandVIxSection::SetupTarget(G4int Z,
201 {
202 G4double cosTetMaxNuc2 = cosTetMaxNuc;
203 if(Z != targetZ || tkin != etag) {
204 etag = tkin;
205 targetZ = std::min(Z, 99);
206 G4double massT = (1 == Z) ? CLHEP::proton_
207 fNistManager->GetAtomicMassAmu(Z)*CLHEP:
208 SetTargetMass(massT);
209
210 kinFactor = coeff*Z*chargeSquare*invbeta2/
211 if(particle == theElectron && fMottXSectio
212 fMottFactor = (1.0 + 2.0e-4*Z*Z);
213 }
214
215 if(1 == Z) {
216 screenZ = ScreenRSquare[targetZ]/mom2;
217 } else if(mass > MeV) {
218 screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z
219 ScreenRSquare[targetZ]/mom2;
220 } else {
221 G4double tau = tkin/mass;
222 screenZ = std::min(Z*1.13,(1.13 +3.76*Z*
223 *invbeta2*alpha2*std::sqrt(tau/(tau
224 ScreenRSquareElec[targetZ]/mom2;
225 }
226 if(targetZ == 1 && particle == theProton &
227 cosTetMaxNuc2 = 0.0;
228 }
229 formfactA = FormFactor[targetZ]*mom2;
230
231 cosTetMaxElec = 1.0;
232 ComputeMaxElectronScattering(cut);
233 }
234 //G4cout << "SetupTarget: Z= " << targetZ <
235 // << " fMottFactor= " << fMottFactor << "
236 return cosTetMaxNuc2;
237 }
238
239 //....oooOO0OOooo........oooOO0OOooo........oo
240
241 G4double
242 G4WentzelOKandVIxSection::ComputeTransportCros
243 {
244 G4double xSection = 0.0;
245 if(cosTMax >= 1.0) { return xSection; }
246
247 G4double costm = std::max(cosTMax,cosTetMaxE
248 G4double fb = screenZ*factB;
249
250 // scattering off electrons
251 if(costm < 1.0) {
252 G4double x = (1.0 - costm)/screenZ;
253 if(x < numlimit) {
254 G4double x2 = 0.5*x*x;
255 xSection = x2*((1.0 - 1.3333333*x + 3*x2
256 } else {
257 G4double x1= x/(1 + x);
258 G4double xlog = G4Log(1.0 + x);
259 xSection = xlog - x1 - fb*(x + x1 - 2*xl
260 }
261
262 if(xSection < 0.0) {
263 ++nwarnings;
264 if(nwarnings < nwarnlimit) {
265 G4cout << "G4WentzelOKandVIxSection::C
266 << " scattering on e- <0"
267 << G4endl;
268 G4cout << "cross= " << xSection
269 << " e(MeV)= " << tkin << " p(M
270 << " Z= " << targetZ << " "
271 << particle->GetParticleName()
272 G4cout << " 1-costm= " << 1.0-costm <<
273 << " x= " << x << G4endl;
274 }
275 xSection = 0.0;
276 }
277 }
278 /*
279 G4cout << "G4WentzelOKandVIxSection::Com
280 << " Z= " << targetZ
281 << " e(MeV)= " << tkin/MeV << " XSel= "
282 << " zmaxE= " << (1.0 - cosTetMaxElec)/s
283 << " zmaxN= " << (1.0 - cosThetaMax)/scr
284 << " 1-costm= " << 1.0 - cosThetaMax <<
285 */
286 // scattering off nucleus
287 if(cosTMax < 1.0) {
288 G4double x = (1.0 - cosTMax)/screenZ;
289 G4double y;
290 if(x < numlimit) {
291 G4double x2 = 0.5*x*x;
292 y = x2*((1.0 - 1.3333333*x + 3*x2) - fb*
293 } else {
294 G4double x1= x/(1 + x);
295 G4double xlog = G4Log(1.0 + x);
296 y = xlog - x1 - fb*(x + x1 - 2*xlog);
297 }
298
299 if(y < 0.0) {
300 ++nwarnings;
301 if(nwarnings < nwarnlimit) {
302 G4cout << "G4WentzelOKandVIxSection::C
303 << " scattering on nucleus <0"
304 << G4endl;
305 G4cout << "y= " << y
306 << " e(MeV)= " << tkin << " Z=
307 << particle->GetParticleName()
308 G4cout << " formfactA= " << formfactA
309 << " x= " << x <<G4endl;
310 }
311 y = 0.0;
312 }
313 xSection += y*targetZ;
314 }
315 xSection *= kinFactor;
316
317 /*
318 G4cout << "Z= " << targetZ << " XStot= " <<
319 << " screenZ= " << screenZ << " formF
320 << " for " << particle->GetParticleNa
321 << " m= " << mass << " 1/v= " << sqrt(invbe
322 << " p= " << sqrt(mom2)
323 << " x= " << x << G4endl;
324 */
325 return xSection;
326 }
327
328 //....oooOO0OOooo........oooOO0OOooo........oo
329
330 G4ThreeVector&
331 G4WentzelOKandVIxSection::SampleSingleScatteri
332
333
334 {
335 temp.set(0.0,0.0,1.0);
336 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra
337
338 G4double formf = formfactA;
339 G4double cost1 = cosTMin;
340 G4double cost2 = cosTMax;
341 if(elecRatio > 0.0) {
342 if(rndmEngineMod->flat() <= elecRatio) {
343 formf = 0.0;
344 cost1 = std::max(cost1,cosTetMaxElec);
345 cost2 = std::max(cost2,cosTetMaxElec);
346 }
347 }
348 if(cost1 > cost2) {
349 G4double w1 = 1. - cost1;
350 G4double w2 = 1. - cost2;
351 G4double w3 = rndmEngineMod->flat()*(w2 -
352 G4double z1 = ((w2 - w3)*screenZ + w1*w2)/
353 G4double fm = 1.0;
354
355 if(fNucFormfactor == fExponentialNF) {
356 fm += formf*z1;
357 fm = 1.0/(fm*fm);
358 } else if(fNucFormfactor == fGaussianNF) {
359 fm = G4Exp(-2*formf*z1);
360 } else if(fNucFormfactor == fFlatNF) {
361 static const G4double ccoef = 0.00508/CL
362 G4double x = std::sqrt(2.*mom2*z1)*ccoef
363 fm = FlatFormfactor(x);
364 fm *= FlatFormfactor(x*0.6*fG4pow->A13(f
365 }
366 // G4cout << " fm=" << fm << " " << fMott
367 G4double grej;
368 if(nullptr != fMottXSection) {
369 fMottXSection->SetupKinematic(tkin, targ
370 grej = fMottXSection->RatioMottRutherfor
371 } else {
372 grej = (1. - z1*factB + factB1*targetZ*s
373 *fm/(1.0 + z1*factD);
374 }
375 if(fMottFactor*rndmEngineMod->flat() <= gr
376 // exclude "false" scattering due to for
377 G4double cost = 1.0 - z1;
378 if(cost > 1.0) { cost = 1.0; }
379 else if(cost < -1.0) { cost =-1.0; }
380 G4double sint = sqrt((1.0 - cost)*(1.0 +
381 //G4cout << "sint= " << sint << G4endl;
382 G4double phi = twopi*rndmEngineMod->fla
383 temp.set(sint*cos(phi),sint*sin(phi),cos
384 }
385 }
386 return temp;
387 }
388
389 //....oooOO0OOooo........oooOO0OOooo........oo
390
391 void
392 G4WentzelOKandVIxSection::ComputeMaxElectronSc
393 {
394 if(mass > MeV) {
395 G4double ratio = electron_mass_c2/mass;
396 G4double tau = tkin/mass;
397 G4double tmax = 2.0*electron_mass_c2*tau*(
398 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*rat
399 cosTetMaxElec = 1.0 - std::min(cutEnergy,
400 } else {
401
402 G4double tmax = (particle == theElectron)
403 G4double t = std::min(cutEnergy, tmax);
404 G4double mom21 = t*(t + 2.0*electron_mass_
405 G4double t1 = tkin - t;
406 //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax
407 //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4end
408 if(t1 > 0.0) {
409 G4double mom22 = t1*(t1 + 2.0*mass);
410 G4double ctm = (mom2 + mom22 - mom21)*0.
411 if(ctm < 1.0) { cosTetMaxElec = ctm; }
412 if(particle == theElectron && cosTetMaxE
413 cosTetMaxElec = 0.0;
414 }
415 }
416 }
417 }
418
419 //....oooOO0OOooo........oooOO0OOooo........oo
420
421 G4double
422 G4WentzelOKandVIxSection::ComputeSecondTranspo
423 {
424 return 0.0;
425 }
426
427 //....oooOO0OOooo........oooOO0OOooo........oo
428