Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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