Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4ScreeningMottCrossSection.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// 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 Mott Cross Section by means of G4MottCoefficients class.
39//
40//
41// Class Description:
42// Computation of electron Coulomb Scattering Cross Section.
43// Suitable for high energy electrons and light target materials.
44//
45// Reference:
46// M.J. Boschini et al.
47// "Non Ionizing Energy Loss induced by Electrons in the Space Environment"
48// Proc. of the 13th International Conference on Particle Physics and Advanced Technology
49// (13th ICPPAT, Como 3-7/10/2011), World Scientific (Singapore).
50// Available at: http://arxiv.org/abs/1111.4042v4
51//
52// 1) Mott Differential Cross Section Approximation:
53// For Target material up to Z=92 (U):
54// As described in http://arxiv.org/abs/1111.4042v4
55// par. 2.1 , eq. (16)-(17)
56// Else (Z>92):
57// W. A. McKinley and H. Fashbach, Phys. Rev. 74, (1948) 1759.
58// 2) Screening coefficient:
59// vomn G. Moliere, Z. Naturforsh A2 (1947), 133-145; A3 (1948), 78.
60// 3) Nuclear Form Factor:
61// A.V. Butkevich et al. Nucl. Instr. and Meth. in Phys. Res. A 488 (2002), 282-294.
62//
63// -------------------------------------------------------------------------------------
64//
65//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66
67#include "G4ScreeningMottCrossSection.hh"
68#include "G4PhysicalConstants.hh"
69#include "G4SystemOfUnits.hh"
70#include "G4MottCoefficients.hh"
71#include "Randomize.hh"
72#include "G4Proton.hh"
73#include "G4LossTableManager.hh"
74#include "G4NucleiProperties.hh"
75#include "G4Element.hh"
76#include "G4UnitsTable.hh"
77
78//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79
80
81using namespace std;
82
83G4ScreeningMottCrossSection::G4ScreeningMottCrossSection():
84 cosThetaMin(1.0),
85 cosThetaMax(-1.0),
86 alpha(fine_structure_const),
87 htc2(hbarc_squared),
88 e2(electron_mass_c2*classic_electr_radius)
89{
90
91 TotalCross=0;
92
93 fNistManager = G4NistManager::Instance();
94 particle=0;
95
96 spin = mass = mu_rel=0;
97 tkinLab = momLab2 = invbetaLab2=0;
98 tkin = mom2 = invbeta2=beta=gamma=0;
99
100 Trec=targetZ = targetMass = As =0;
101 etag = ecut = 0.0;
102
103 targetA = 0;
104
105 cosTetMinNuc=0;
106 cosTetMaxNuc=0;
107
108 for(G4int i=0 ; i<5; i++){
109 for(G4int j=0; j< 6; j++){
110 coeffb[i][j]=0;
111 }
112 }
113
114
115
116
117}
118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119
120G4ScreeningMottCrossSection::~G4ScreeningMottCrossSection()
121{}
122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123
124void G4ScreeningMottCrossSection::Initialise(const G4ParticleDefinition* p,
125 G4double CosThetaLim)
126{
127
128 SetupParticle(p);
129 tkin = targetZ = mom2 = DBL_MIN;
130 ecut = etag = DBL_MAX;
131 particle = p;
132 cosThetaMin = CosThetaLim;
133
134}
135//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
136void G4ScreeningMottCrossSection::SetScreeningCoefficient(){
137
138 G4double alpha2=alpha*alpha;
139 //Bohr radius
140 G4double a0= Bohr_radius ;//0.529e-8*cm;
141 //Thomas-Fermi screening length
142 G4double aU=0.88534*a0/pow(targetZ,1./3.);
143 G4double twoR2=aU*aU;
144
145 G4double factor= 1.13 + 3.76*targetZ*targetZ*invbeta2*alpha2;
146 As=0.25*(htc2)/(twoR2*mom2)*factor;
147
148
149
150
151//cout<<"0k .........................As "<<As<<endl;
152
153}
154
155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156
157G4double G4ScreeningMottCrossSection::GetScreeningAngle(){
158
159 SetScreeningCoefficient();
160
161 //cout<<" As "<<As<<endl;
162 if(As < 0.0) { As = 0.0; }
163 else if(As > 1.0) { As = 1.0; }
164
165 G4double screenangle=2.*asin(sqrt(As));
166
167 // cout<<" screenangle "<< screenangle <<endl;
168
169 if(screenangle>=pi) screenangle=pi;
170
171
172return screenangle;
173
174}
175//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
176
177void G4ScreeningMottCrossSection::SetupKinematic(G4double ekin, G4double Z )
178{
179
180 //...Target
181 G4int iz = G4int(Z);
182 G4double A = fNistManager->GetAtomicMassAmu(iz);
183 G4int ia = G4int(A);
184 G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);
185
186 targetZ = Z;
187 targetA = fNistManager->GetAtomicMassAmu(iz);
188 targetMass= mass2;
189
190 mottcoeff->SetMottCoeff(targetZ, coeffb);
191
192 //cout<<"......... targetA "<< targetA <<endl;
193 //cout<<"......... targetMass "<< targetMass/MeV <<endl;
194
195 // incident particle lab
196 tkinLab = ekin;
197 momLab2 = tkinLab*(tkinLab + 2.0*mass);
198 invbetaLab2 = 1.0 + mass*mass/momLab2;
199
200 G4double etot = tkinLab + mass;
201 G4double ptot = sqrt(momLab2);
202 G4double m12 = mass*mass;
203
204
205 // relativistic reduced mass from publucation
206 // A.P. Martynenko, R.N. Faustov, Teoret. mat. Fiz. 64 (1985) 179
207
208 //incident particle & target nucleus
209 G4double Ecm=sqrt(m12 + mass2*mass2 + 2.0*etot*mass2);
210 mu_rel=mass*mass2/Ecm;
211 G4double momCM= ptot*mass2/Ecm;
212 // relative system
213 mom2 = momCM*momCM;
214 invbeta2 = 1.0 + mu_rel*mu_rel/mom2;
215 tkin = momCM*sqrt(invbeta2) - mu_rel;//Ekin of mu_rel
216 G4double beta2=1./invbeta2;
217 beta=std::sqrt(beta2) ;
218 G4double gamma2= 1./(1.-beta2);
219 gamma=std::sqrt(gamma2);
220
221
222
223 //.........................................................
224
225
226 G4double screenangle=GetScreeningAngle()/10.;
227 //cout<<" screenangle [rad] "<<screenangle/rad <<endl;
228
229 cosTetMinNuc =min( cosThetaMin ,cos(screenangle));
230 cosTetMaxNuc =cosThetaMax;
231
232//cout<<"ok..............mu_rel "<<mu_rel<<endl;
233
234}
235//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236
237G4double G4ScreeningMottCrossSection::FormFactor2ExpHof(G4double angle){
238
239
240 G4double M=targetMass;
241 G4double E=tkinLab;
242 G4double Etot=E+mass;
243 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
244 G4double T=Tmax*pow(sin(angle/2.),2.);
245 G4double q2=T*(T+2.*M);
246 q2/=htc2;//1/cm2
247 G4double RN=1.27e-13*pow(targetA,0.27)*cm;
248 G4double xN= (RN*RN*q2);
249 G4double den=(1.+xN/12.);
250 G4double FN=1./(den*den);
251 G4double form2=(FN*FN);
252
253
254 return form2;
255
256//cout<<"..................... form2 "<< form2<<endl;
257}
258
259//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260
261G4double G4ScreeningMottCrossSection::McFcorrection(G4double angle ){
262
263
264 G4double beta2=1./invbeta2;
265 G4double sintmezzi=std::sin(angle/2.);
266 G4double sin2tmezzi = sintmezzi*sintmezzi;
267 G4double R=1.-beta2*sin2tmezzi + targetZ*alpha*beta*pi*sintmezzi*(1.-sintmezzi);
268
269
270return R;
271}
272//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
273G4double G4ScreeningMottCrossSection::RatioMottRutherford(G4double angle){
274
275
276 G4double R=0;
277 G4double fcost=std::sqrt((1. -cos(angle)));
278 G4double a[5];
279 G4double shift=0.7181228;
280 G4double beta0= beta -shift;
281
282 for(G4int j=0 ;j<=4;j++){
283 a[j]=0;
284 }
285
286 for(G4int j=0 ;j<=4;j++){
287 for(G4int k=0;k<=5;k++ ){
288 a[j]+=coeffb[j][k]*pow(beta0,k);
289 }
290 }
291
292 for(G4int j=0 ;j<=4 ;j++){
293 R+=a[j]* pow(fcost,j);
294 }
295
296
297
298return R;
299
300}
301//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
302
303G4double G4ScreeningMottCrossSection::NuclearCrossSection()
304{
305 if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
306
307 TotalCross=0;
308
309 G4double anglemin =std::acos(cosTetMinNuc);
310 G4double anglemax =std::acos(cosTetMaxNuc);
311
312 const G4double limit = 1.e-9;
313 if(anglemin < limit) {
314 anglemin = GetScreeningAngle()/10.;
315 if(anglemin < limit) { anglemin = limit; }
316 }
317
318 //cout<<" anglemin "<< anglemin <<endl;
319
320 G4double loganglemin=log10(anglemin);
321 G4double loganglemax=log10(anglemax);
322 G4double logdangle=0.01;
323
324 G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);
325
326
327 vector<G4double> angle;
328 vector<G4double> tet;
329 vector<G4double> dangle;
330 vector<G4double> cross;
331
332 for(G4int i=0; i<=bins; i++ ){
333 tet.push_back(0);
334 dangle.push_back(0);
335 angle.push_back(pow(10.,loganglemin+logdangle*i));
336 cross.push_back(0);
337 }
338
339
340 G4int dim = tet.size();
341 //cout<<"dim--- "<<dim<<endl;
342
343
344 for(G4int i=0; i<dim;i++){
345
346 if(i!=dim-1){
347 dangle[i]=(angle[i+1]-angle[i]);
348 tet[i]=(angle[i+1]+angle[i])/2.;
349 }else if(i==dim-1){
350 break;
351 }
352
353
354 G4double R=0;
355 G4double F2=FormFactor2ExpHof(tet[i]);
356
357 if (coeffb[0][0]!=0){
358 //cout<<" Mott....targetZ "<< targetZ<<endl;
359 R=RatioMottRutherford(tet[i]);
360 }
361
362 else if (coeffb[0][0]==0){
363 // cout<<" McF.... targetZ "<< targetZ<<endl;
364 R=McFcorrection(tet[i]);
365 }
366
367
368//cout<<"----------------- R "<<R<<" F2 "<<F2<<endl;
369
370
371// cout<<"angle "<<tet[i] << " F2 "<<F2<<endl;
372
373 G4double e4=e2*e2;
374 G4double den=2.*As+2.*pow(sin(tet[i]/2.),2.);
375 G4double func=1./(den*den);
376
377 G4double fatt= targetZ/(mu_rel*gamma*beta*beta);
378 G4double sigma=e4*fatt*fatt*func;
379 cross[i]=F2*R*sigma;
380 G4double pi2sintet=2.*pi*sin(tet[i]);
381
382 TotalCross+=pi2sintet*cross[i]*dangle[i];
383 }//end integral
384
385
386//cout<< "ok ......... TotalCross "<<TotalCross<<endl;
387
388
389return TotalCross;
390
391}
392//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
393G4double G4ScreeningMottCrossSection::AngleDistribution(G4double anglein){
394
395
396 G4double total=TotalCross ;
397 G4double fatt= e2*targetZ/(mu_rel*gamma*beta*beta);
398 G4double fatt2=fatt*fatt;
399 total/=fatt2;
400
401 G4double R=0;
402 if (coeffb[0][0]!=0){
403 // cout<<" Mott....targetZ "<< targetZ<<endl;
404 R=RatioMottRutherford(anglein);
405 }
406
407 else if (coeffb[0][0]==0){
408 // cout<<" McF.... targetZ "<< targetZ<<endl;
409 R=McFcorrection(anglein);
410 }
411
412
413 G4double y=2.*pi*sin(anglein)*R*FormFactor2ExpHof(anglein)/
414 ((2*As+2.*pow(sin(anglein/2.),2.))*(2.*As+2.*pow(sin(anglein/2.),2.) ));
415
416return y/total;
417
418}
419
420//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
421
422G4double G4ScreeningMottCrossSection::GetScatteringAngle()
423{
424
425
426//cout<<" G4ScreeningMottCrossSection::SampleCosineTheta ............."<<endl;
427
428 if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
429
430 G4double anglemin=std::acos(cosTetMinNuc);
431 G4double anglemax= std::acos(cosTetMaxNuc);
432
433 const G4double limit = 1.e-9;
434 if(anglemin < limit) {
435 anglemin = GetScreeningAngle()/10.;
436 if(anglemin < limit) { anglemin = limit; }
437 }
438
439// cout<<"................ tkinLab "<< G4BestUnit(tkinLab,"Energy") << " anglemin= "<<anglemin<<endl;
440 //cout<<"anglemax= "<<anglemax<<endl;
441 G4double r =G4UniformRand();
442
443 G4double loganglemin=log10(anglemin);
444 G4double loganglemax=log10(anglemax);
445 G4double logdangle=0.01;
446
447 G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);
448
449 std::vector<G4double> angle;
450 std::vector<G4double> tet;
451 std::vector<G4double> dangle;
452
453 for(G4int i=0; i<=bins; i++ ){
454 tet.push_back(0);
455 dangle.push_back(0);
456 angle.push_back(pow(10.,loganglemin+logdangle*i));
457 }
458
459
460 G4int dim = tet.size();
461 G4double scattangle=0;
462 G4double y=0;
463 G4double dy=0;
464 G4double area=0;
465
466 for(G4int i=0; i<dim;i++){
467
468 if(i!=dim-1){
469 dangle[i]=(angle[i+1]-angle[i]);
470 tet[i]=(angle[i+1]+angle[i])/2.;
471 }else if(i==dim-1){
472 break;
473 }
474
475 y+=AngleDistribution(tet[i])*dangle[i];
476 dy= y-area ;
477 area=y;
478
479 if(r >=y-dy && r<=y+dy ){
480 scattangle= angle[i] +G4UniformRand()*dangle[i];
481 //cout<<"y "<<y <<" r "<< r << " y/r "<<y/r << endl;
482 break;
483
484 }
485
486 }
487
488
489 return scattangle;
490
491
492}
493
494//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
495
496G4ThreeVector G4ScreeningMottCrossSection::GetNewDirection(){
497
498G4ThreeVector dir(0.0,0.0,1.0);
499
500
501 G4double z1=GetScatteringAngle();
502
503 G4double cost = cos(z1);
504 G4double sint = sin(z1);
505 G4double phi = twopi* G4UniformRand();
506 G4double dirx = sint*cos(phi);
507 G4double diry = sint*sin(phi);
508 G4double dirz = cost;
509
510
511 //.......set Trc
512 G4double etot=tkinLab+mass;
513 G4double mass2=targetMass;
514 Trec=(1.0 - cost)* mass2*(etot*etot - mass*mass )/
515 (mass*mass + mass2*mass2+ 2.*mass2*etot);
516
517
518
519 dir.set(dirx,diry,dirz);
520
521
522
523
524return dir;
525
526}
527
528
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4MottCoefficients::SetMottCoeff
void SetMottCoeff(G4double targetZ, G4double coeff[5][6])
Definition: G4MottCoefficients.cc:81
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:343
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:53
G4ScreeningMottCrossSection::SetupKinematic
void SetupKinematic(G4double kinEnergy, G4double Z)
Definition: G4ScreeningMottCrossSection.cc:177
G4ScreeningMottCrossSection::Initialise
void Initialise(const G4ParticleDefinition *, G4double cosThetaLim)
Definition: G4ScreeningMottCrossSection.cc:124
G4ScreeningMottCrossSection::SetupParticle
void SetupParticle(const G4ParticleDefinition *)
Definition: G4ScreeningMottCrossSection.hh:184
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83