| 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 // G4ScreeningMottCrossSection.cc
27 //--------------------------------------------
28 //
29 // GEANT4 Class header file
30 //
31 // File name: G4ScreeningMottCrossSection
32 //
33 // Author: Cristina Consolandi
34 //
35 // Creation date: 20.10.2011
36 //
37 // Modifications:
38 // 27-05-2012 Added Analytic Fitting to the Mo
39 //
40 //
41 // Class Description:
42 // Computation of electron Coulomb Scattering
43 // Suitable for high energy electrons and lig
44 //
45 // Reference:
46 // M.J. Boschini et al.
47 // "Non Ionizing Energy Loss induced by El
48 // Proc. of the 13th International Confer
49 // (13th ICPPAT, Como 3-7/10/2011), World
50 // Available at: http://arxiv.org/abs/1111.40
51 //
52 // 1) Mott Differential Cross Section App
53 // For Target material up to Z=92 (U):
54 // As described in http://arxiv.org/ab
55 // par. 2.1 , eq. (16)-(17)
56 // Else (Z>92):
57 // W. A. McKinley and H. Fashbach, Phys. R
58 // 2) Screening coefficient:
59 // vomn G. Moliere, Z. Naturforsh A2 (194
60 // 3) Nuclear Form Factor:
61 // A.V. Butkevich et al. Nucl. Instr. Met
62 //
63 // -------------------------------------------
64 //
65 //....oooOO0OOooo........oooOO0OOooo........oo
66
67 #include "G4ScreeningMottCrossSection.hh"
68 #include "G4PhysicalConstants.hh"
69 #include "G4SystemOfUnits.hh"
70 #include "Randomize.hh"
71 #include "G4Proton.hh"
72 #include "G4LossTableManager.hh"
73 #include "G4NucleiProperties.hh"
74 #include "G4Element.hh"
75 #include "G4UnitsTable.hh"
76 #include "G4NistManager.hh"
77 #include "G4ThreeVector.hh"
78 #include "G4Pow.hh"
79 #include "G4MottData.hh"
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 static G4double invlog10 = 1./std::log(10.);
84
85 static const G4double angle[DIMMOTT]={
86 1e-07,1.02329e-07,1.04713e-07,1.07152e-07,1.09
87
88 //using namespace std;
89
90 G4ScreeningMottCrossSection::G4ScreeningMottCr
91 cosThetaMin(1.0),
92 cosThetaMax(-1.0),
93 alpha(fine_structure_const),
94 htc2(hbarc_squared),
95 e2(electron_mass_c2*classic_electr_radius)
96 {
97 fTotalCross=0;
98
99 fNistManager = G4NistManager::Instance();
100 fG4pow = G4Pow::GetInstance();
101 particle=nullptr;
102
103 spin = mass = mu_rel = 0.0;
104 tkinLab = momLab2 = invbetaLab2 = 0.0;
105 tkin = mom2 = invbeta2 = beta = gamma = 0.0;
106
107 targetMass = As = etag = ecut = 0.0;
108
109 targetZ = targetA = 0;
110
111 cosTetMinNuc = cosTetMaxNuc = 0.0;
112 }
113
114 //....Ooooo0ooooo........oooOO0OOooo........oo
115
116 G4ScreeningMottCrossSection::~G4ScreeningMottC
117
118 //....oooOO0OOooo........oooOO0OOooo........oo
119
120 void G4ScreeningMottCrossSection::Initialise(c
121 G
122 {
123 SetupParticle(p);
124 tkin = mom2 = 0.0;
125 ecut = etag = DBL_MAX;
126 particle = p;
127 cosThetaMin = cosThetaLim;
128 }
129
130 //....oooOO0OOooo........oooOO0OOooo........oo
131
132 void G4ScreeningMottCrossSection::SetupKinemat
133 {
134 //...Target
135 const G4int iz = std::min(92, Z);
136 const G4int ia = G4lrint(fNistManager->GetAt
137
138 targetZ = iz;
139 targetA = ia;
140 targetMass = G4NucleiProperties::GetNuclearM
141
142 //G4cout<<"......... targetA "<< targetA <<G
143 //G4cout<<"......... targetMass "<< targetMa
144
145 // incident particle lab
146 tkinLab = ekin;
147 momLab2 = tkinLab*(tkinLab + 2.0*mass);
148 invbetaLab2 = 1.0 + mass*mass/momLab2;
149
150 const G4double etot = tkinLab + mass;
151 const G4double ptot = std::sqrt(momLab2);
152 const G4double m12 = mass*mass;
153
154 // relativistic reduced mass from publucatio
155 // A.P. Martynenko, R.N. Faustov, Teoret. ma
156
157 //incident particle & target nucleus
158 const G4double Ecm = std::sqrt(m12 + targetM
159 mu_rel = mass*targetMass/Ecm;
160 const G4double momCM= ptot*targetMass/Ecm;
161 // relative system
162 mom2 = momCM*momCM;
163 const G4double x = mu_rel*mu_rel/mom2;
164 invbeta2 = 1.0 + x;
165 tkin = momCM*std::sqrt(invbeta2) - mu_rel;//
166 const G4double beta2 = 1./invbeta2;
167 beta = std::sqrt(beta2) ;
168 const G4double gamma2= invbeta2/x;
169 gamma = std::sqrt(gamma2);
170
171 //Thomas-Fermi screening length
172 const G4double alpha2 = alpha*alpha;
173 const G4double aU = 0.88534*CLHEP::Bohr_radi
174 const G4double twoR2 = aU*aU;
175 As = 0.25*htc2*(1.13 + 3.76*targetZ*targetZ*
176
177 //Integration Angles definition
178 cosTetMinNuc = cosThetaMin;
179 cosTetMaxNuc = cosThetaMax;
180 }
181
182 //....oooOO0OOooo........oooOO0OOooo........oo
183
184 G4double G4ScreeningMottCrossSection::FormFact
185 {
186 G4double M=targetMass;
187 G4double E=tkinLab;
188 G4double Etot=E+mass;
189 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+
190 G4double T=Tmax*sin2t2;
191 G4double q2=T*(T+2.*M);
192 q2/=htc2;//1/cm2
193 G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targ
194 G4double xN= (RN*RN*q2);
195 G4double den=(1.+xN/12.);
196 G4double FN=1./(den*den);
197 G4double form2=(FN*FN);
198
199 return form2;
200 }
201
202 //....oooOO0OOooo........oooOO0OOooo........oo
203
204 G4double G4ScreeningMottCrossSection::FormFact
205 {
206 G4double M=targetMass;
207 G4double E=tkinLab;
208 G4double Etot=E+mass;
209 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+
210 G4double T=Tmax*sin2t2;
211 G4double q2=T*(T+2.*M);
212 q2/=htc2;//1/cm2
213 G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targ
214 G4double xN= (RN*RN*q2);
215
216 G4double expo=(-xN/6.);
217 G4double FN=G4Exp(expo);
218
219 G4double form2=(FN*FN);
220
221 return form2;
222 }
223
224 //....oooOO0OOooo........oooOO0OOooo........oo
225
226 G4double G4ScreeningMottCrossSection::FormFact
227 {
228 G4double M=targetMass;
229 G4double E=tkinLab;
230 G4double Etot=E+mass;
231 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+
232 G4double T=Tmax*sin2t2;
233 G4double q2=T*(T+2.*M);
234 q2=q2/(htc2*0.01);//1/cm2
235
236 G4double q=std::sqrt(q2);
237 G4double R0=1.2E-13*fG4pow->Z13(targetA);
238 G4double R1=2.0E-13;
239
240 G4double x0=q*R0;
241 G4double F0=(3./fG4pow->powN(x0,3))*(std::si
242
243 G4double x1=q*R1;
244 G4double F1=(3./fG4pow->powN(x1,3))*(std::si
245
246 G4double F=F0*F1;
247
248 G4double form2=(F*F);
249
250 return form2;
251 }
252
253 //....oooOO0OOooo........oooOO0OOooo........oo
254
255 G4double G4ScreeningMottCrossSection::McFcorre
256 {
257 const G4double sint = std::sqrt(sin2t2);
258 return 1.-beta*beta*sin2t2 + targetZ*alpha*b
259 }
260
261 //....oooOO0OOooo........oooOO0OOooo........oo
262
263 G4double G4ScreeningMottCrossSection::RatioMot
264 {
265 return RatioMottRutherfordCosT(std::sqrt(1.
266 }
267
268 //....oooOO0OOooo........oooOO0OOooo........oo
269
270 G4double G4ScreeningMottCrossSection::RatioMot
271 {
272 G4double R(0.);
273 const G4double shift = 0.7181228;
274 G4double beta0 = beta - shift;
275
276 G4double b0 = 1.0;
277 G4double b[6];
278 for(G4int i=0; i<6; ++i) {
279 b[i] = b0;
280 b0 *= beta0;
281 }
282
283 b0 = 1.0;
284 for(G4int j=0; j<=4; ++j) {
285 G4double a = 0.0;
286 for(G4int k=0; k<=5; ++k){
287 a += fMottCoef[targetZ][j][k]*b[k];
288 }
289 R += a*b0;
290 b0 *= fcost;
291 }
292 return R;
293 }
294
295 //....oooOO0OOooo........oooOO0OOooo........oo
296
297 G4double G4ScreeningMottCrossSection::GetTrans
298 {
299 G4double x = G4Log(tkinLab)*invlog10;
300 G4double res(0.), x0(1.0);
301 for(G4int i=0; i<11; ++i) {
302 res += fPRM[targetZ][i]*x0;
303 x0 *= x;
304 }
305 //G4cout << "Z= " << targetZ << " x0= " << x
306 return res;
307 }
308
309 //....oooOO0OOooo........oooOO0OOooo........oo
310
311 G4double
312 G4ScreeningMottCrossSection::DifferentialXSect
313 {
314 const G4double ang = angle[idx];
315 const G4double z1 = 1.0 - std::cos(ang);
316 G4double step;
317 if(0 == idx) {
318 step = 0.5*(angle[1] + angle[0]);
319 } else if(DIMMOTT == idx + 1) {
320 step = CLHEP::pi - 0.5*(angle[DIMMOTT-2] +
321 } else {
322 step = 0.5*(angle[idx+1] - angle[idx-1]);
323 }
324
325 G4double F2 = 1.;
326 switch (form) {
327 case 1:
328 F2 = FormFactor2ExpHof(z1*0.5);
329 break;
330 case 2:
331 F2 = FormFactor2Gauss(z1*0.5);
332 break;
333 case 3:
334 F2 = FormFactor2UniformHelm(z1*0.5);
335 break;
336 }
337
338 const G4double R = RatioMottRutherfordCosT(s
339
340 G4double den = 2.*As + z1;
341 G4double func = 1./(den*den);
342
343 G4double fatt = (G4double)targetZ/(mu_rel*ga
344 G4double sigma= e2*e2*fatt*fatt*func;
345 G4double pi2sintet = CLHEP::twopi*std::sqrt(
346
347 G4double Xsec = std::max(pi2sintet*F2*R*sigm
348 return Xsec;
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oo
352
353 G4double
354 G4ScreeningMottCrossSection::NuclearCrossSecti
355 {
356 fTotalCross=0.;
357 if(cosTetMaxNuc >= cosTetMinNuc) return fTot
358 if(0 == cross.size()) { cross.resize(DIMMOTT
359
360 //G4cout<<"MODEL: "<<fast<<G4endl;
361 //************ PRECISE COMPUTATION
362 if(fast == 0){
363
364 for(G4int i=0; i<DIMMOTT; ++i){
365 G4double y = DifferentialXSection(i, for
366 fTotalCross += y;
367 cross[i] = fTotalCross;
368 if(fTotalCross*1.e-9 > y) {
369 for(G4int j=i+1; j<DIMMOTT; ++j) {
370 cross[j] = fTotalCross;
371 }
372 break;
373 }
374 }
375 //**************** FAST COMPUTATION
376 } else if(fast == 1) {
377
378 G4double p0 = electron_mass_c2*classic_ele
379 G4double coeff = twopi*p0*p0;
380
381 G4double fac = coeff*targetZ*targetZ*invbe
382
383 G4double x = 1.0 - cosTetMinNuc;
384 G4double x1 = x + 2*As;
385
386 // scattering with nucleus
387 fTotalCross = fac*(cosTetMinNuc - cosTetMa
388 (x1*(1.0 - cosTetMaxNuc + 2*As));
389 }
390 /*
391 G4cout << "Energy(MeV): " << tkinLab/CLHEP::
392 << " Total Cross(mb): " << fTotalCros
393 << " form= " << form << " fast= " <<
394 */
395 return fTotalCross;
396 }
397
398 //....oooOO0OOooo........oooOO0OOooo........oo
399
400 G4double
401 G4ScreeningMottCrossSection::GetScatteringAngl
402 {
403 G4double scattangle = 0.0;
404 G4double r = G4UniformRand();
405 //************ PRECISE COMPUTATION
406 if(fast == 0){
407 //G4cout << "r= " << r << " tot= " << fTot
408 r *= fTotalCross;
409 for(G4int i=0; i<DIMMOTT; ++i){
410 //G4cout << i << ". r= " << r << " xs= "
411 if(r <= cross[i]) {
412 scattangle = ComputeAngle(i, r);
413 break;
414 }
415 }
416
417 //**************** FAST COMPUTATION
418 } else if(fast == 1) {
419
420 G4double limit = GetTransitionRandom();
421 if(limit > 0.0) {
422 G4double Sz = 2*As;
423 G4double x = Sz-((Sz*(2+Sz))/(Sz+2*limit
424 //G4cout << "limit= " << limit << " Sz=
425 G4double angle_limit = (std::abs(x) < 1.
426 //G4cout<<"ANGLE LIMIT: "<<angle_limit<<
427
428 if(r > limit && angle_limit != 0.0) {
429 x = Sz-((Sz*(2+Sz))/(Sz+2*r))+1.0;
430 scattangle = (x >= 1.0) ? 0.0 : ((x <= -1.0)
431 //G4cout<<"FAST: scattangle= "<< scattangle
432 }
433 } else {
434 //G4cout<<"SLOW: total= "<<fTotalCross<<
435 r *= fTotalCross;
436 G4double tot = 0.0;
437 for(G4int i=0; i<DIMMOTT; ++i) {
438 G4double xs = DifferentialXSection(i, form);
439 tot += xs;
440 cross[i] = tot;
441 if(r <= tot) {
442 scattangle = ComputeAngle(i, r);
443 break;
444 }
445 }
446 }
447 }
448 //******************************************
449 //G4cout<<"Energy(MeV): "<<tkinLab/MeV<<" SC
450 return scattangle;
451 }
452
453 G4double G4ScreeningMottCrossSection::ComputeA
454 {
455 G4double x1, x2, y;
456 if(0 == i) {
457 x1 = 0.0;
458 x2 = 0.5*(angle[0] + angle[1]);
459 y = cross[0];
460 } else if(i == DIMMOTT-1) {
461 x1 = 0.5*(angle[i] + angle[i-1]);
462 x2 = CLHEP::pi;
463 y = cross[i] - cross[i-1];
464 r -= cross[i-1];
465 } else {
466 x1 = 0.5*(angle[i] + angle[i-1]);
467 x2 = 0.5*(angle[i] + angle[i+1]);
468 y = cross[i] - cross[i-1];
469 r -= cross[i-1];
470 }
471 //G4cout << i << ". r= " << r << " y= " << y
472 // << " x1= " << " x2= " << x2 << G4endl;
473 return x1 + (x2 - x1)*r/y;
474 }
475