Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
81 using namespace std;
82 
83 G4ScreeningMottCrossSection::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 
120 G4ScreeningMottCrossSection::~G4ScreeningMottCrossSection()
121 {}
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 
124 void 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....
136 void 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 
157 G4double 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 
172 return screenangle;
173 
174 }
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
176 
177 void 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 
237 G4double 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 
261 G4double 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 
270 return R;
271 }
272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
273 G4double 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 
298 return R;
299 
300 }
301 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
302 
303 G4double 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 
389 return TotalCross;
390 
391 }
392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
393 G4double 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 
416 return y/total;
417 
418 }
419 
420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
421 
422 G4double 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 
496 G4ThreeVector G4ScreeningMottCrossSection::GetNewDirection(){
497 
498 G4ThreeVector 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 
524 return dir;
525 
526 }
527 
528 
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4