Geant4  10.01.p02
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 
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 
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  TotalCross=0;
91 
93  particle=0;
94 
95  spin = mass = mu_rel=0;
97  tkin = mom2 = invbeta2=beta=gamma=0;
98 
99  Trec=targetZ = targetMass = As =0;
100  etag = ecut = 0.0;
101 
102  targetA = 0;
103 
104  cosTetMinNuc=0;
105  cosTetMaxNuc=0;
106 
107  for(G4int i=0 ; i<5; i++){
108  for(G4int j=0; j< 6; j++){
109  coeffb[i][j]=0;
110  }
111  }
112 
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 
119 {
120  delete mottcoeff;
121 }
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 
125  G4double CosThetaLim)
126 {
127  SetupParticle(p);
128  tkin = targetZ = mom2 = DBL_MIN;
129  ecut = etag = DBL_MAX;
130  particle = p;
131  cosThetaMin = CosThetaLim;
132 
133 }
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
136 {
137  G4double alpha2=alpha*alpha;
138  //Bohr radius
139  G4double a0= Bohr_radius ;//0.529e-8*cm;
140  //Thomas-Fermi screening length
141  G4double aU=0.88534*a0/pow(targetZ,1./3.);
142  G4double twoR2=aU*aU;
143 
144  G4double factor= 1.13 + 3.76*targetZ*targetZ*invbeta2*alpha2;
145  As=0.25*(htc2)/(twoR2*mom2)*factor;
146  //cout<<"0k .........................As "<<As<<endl;
147 }
148 
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 
152 {
154 
155  //cout<<" As "<<As<<endl;
156  if(As < 0.0) { As = 0.0; }
157  else if(As > 1.0) { As = 1.0; }
158 
159  G4double screenangle=2.*asin(sqrt(As));
160  // cout<<" screenangle "<< screenangle <<endl;
161  if(screenangle>=pi) screenangle=pi;
162 
163  return screenangle;
164 }
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166 
168 {
169  //...Target
170  G4int iz = G4int(Z);
172  G4int ia = G4int(A);
174 
175  targetZ = Z;
177  targetMass= mass2;
178 
180 
181  //cout<<"......... targetA "<< targetA <<endl;
182  //cout<<"......... targetMass "<< targetMass/MeV <<endl;
183 
184  // incident particle lab
185  tkinLab = ekin;
186  momLab2 = tkinLab*(tkinLab + 2.0*mass);
187  invbetaLab2 = 1.0 + mass*mass/momLab2;
188 
189  G4double etot = tkinLab + mass;
190  G4double ptot = sqrt(momLab2);
191  G4double m12 = mass*mass;
192 
193  // relativistic reduced mass from publucation
194  // A.P. Martynenko, R.N. Faustov, Teoret. mat. Fiz. 64 (1985) 179
195 
196  //incident particle & target nucleus
197  G4double Ecm=sqrt(m12 + mass2*mass2 + 2.0*etot*mass2);
198  mu_rel=mass*mass2/Ecm;
199  G4double momCM= ptot*mass2/Ecm;
200  // relative system
201  mom2 = momCM*momCM;
202  invbeta2 = 1.0 + mu_rel*mu_rel/mom2;
203  tkin = momCM*sqrt(invbeta2) - mu_rel;//Ekin of mu_rel
204  G4double beta2=1./invbeta2;
205  beta=std::sqrt(beta2) ;
206  G4double gamma2= 1./(1.-beta2);
207  gamma=std::sqrt(gamma2);
208 
209  //.........................................................
210 
211  G4double screenangle=GetScreeningAngle()/10.;
212  //cout<<" screenangle [rad] "<<screenangle/rad <<endl;
213 
214  cosTetMinNuc =min( cosThetaMin ,cos(screenangle));
216 
217  //cout<<"ok..............mu_rel "<<mu_rel<<endl;
218 }
219 
220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
221 
223 {
224  G4double M=targetMass;
225  G4double E=tkinLab;
226  G4double Etot=E+mass;
227  G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
228  G4double T=Tmax*pow(sin(angle/2.),2.);
229  G4double q2=T*(T+2.*M);
230  q2/=htc2;//1/cm2
231  G4double RN=1.27e-13*pow(targetA,0.27)*cm;
232  G4double xN= (RN*RN*q2);
233  G4double den=(1.+xN/12.);
234  G4double FN=1./(den*den);
235  G4double form2=(FN*FN);
236 
237  return form2;
238 
239  //cout<<"..................... form2 "<< form2<<endl;
240 }
241 
242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
243 
245 {
246  G4double beta2=1./invbeta2;
247  G4double sintmezzi=std::sin(angle/2.);
248  G4double sin2tmezzi = sintmezzi*sintmezzi;
249  G4double R=1.-beta2*sin2tmezzi + targetZ*alpha*beta*pi*sintmezzi*(1.-sintmezzi);
250  return R;
251 }
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254 {
255  G4double R=0;
256  G4double fcost=std::sqrt((1. -cos(angle)));
257  G4double a[5];
258  G4double shift=0.7181228;
259  G4double beta0= beta -shift;
260 
261  for(G4int j=0 ;j<=4;j++){
262  a[j]=0;
263  }
264 
265  for(G4int j=0 ;j<=4;j++){
266  for(G4int k=0;k<=5;k++ ){
267  a[j]+=coeffb[j][k]*pow(beta0,k);
268  }
269  }
270 
271  for(G4int j=0 ;j<=4 ;j++){
272  R+=a[j]* pow(fcost,j);
273  }
274  return R;
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
278 
280 {
281  if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
282 
283  TotalCross=0;
284 
285  G4double anglemin =std::acos(cosTetMinNuc);
286  G4double anglemax =std::acos(cosTetMaxNuc);
287 
288  static const G4double limit = 1.e-9;
289  if(anglemin < limit) {
290  anglemin = GetScreeningAngle()/10.;
291  if(anglemin < limit) { anglemin = limit; }
292  }
293 
294  //cout<<" anglemin "<< anglemin <<endl;
295 
296  G4double loganglemin=log10(anglemin);
297  G4double loganglemax=log10(anglemax);
298  G4double logdangle=0.01;
299 
300  G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);
301 
302  vector<G4double> angle;
303  vector<G4double> tet;
304  vector<G4double> dangle;
305  vector<G4double> cross;
306 
307  for(G4int i=0; i<=bins; i++ ){
308  tet.push_back(0);
309  dangle.push_back(0);
310  angle.push_back(pow(10.,loganglemin+logdangle*i));
311  cross.push_back(0);
312  }
313 
314  G4int dim = tet.size();
315  //cout<<"dim--- "<<dim<<endl;
316 
317  for(G4int i=0; i<dim;i++){
318 
319  if(i!=dim-1){
320  dangle[i]=(angle[i+1]-angle[i]);
321  tet[i]=(angle[i+1]+angle[i])/2.;
322  }else if(i==dim-1){
323  break;
324  }
325 
326  G4double R=0;
327  G4double F2=FormFactor2ExpHof(tet[i]);
328 
329  if (coeffb[0][0]!=0){
330  //cout<<" Mott....targetZ "<< targetZ<<endl;
331  R=RatioMottRutherford(tet[i]);
332  } else if (coeffb[0][0]==0){
333  // cout<<" McF.... targetZ "<< targetZ<<endl;
334  R=McFcorrection(tet[i]);
335  }
336 
337  //cout<<"----------------- R "<<R<<" F2 "<<F2<<endl;
338  // cout<<"angle "<<tet[i] << " F2 "<<F2<<endl;
339 
340  G4double e4=e2*e2;
341  G4double den=2.*As+2.*pow(sin(tet[i]/2.),2.);
342  G4double func=1./(den*den);
343 
345  G4double sigma=e4*fatt*fatt*func;
346  cross[i]=F2*R*sigma;
347  G4double pi2sintet=2.*pi*sin(tet[i]);
348 
349  TotalCross+=pi2sintet*cross[i]*dangle[i];
350  }//end integral
351 
352  //cout<< "ok ......... TotalCross "<<TotalCross<<endl;
353  return TotalCross;
354 }
355 
356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
357 
359 {
362  G4double fatt2=fatt*fatt;
363  total/=fatt2;
364 
365  G4double R=0;
366  if (coeffb[0][0]!=0){
367  // cout<<" Mott....targetZ "<< targetZ<<endl;
368  R=RatioMottRutherford(anglein);
369  } else if (coeffb[0][0]==0){
370  // cout<<" McF.... targetZ "<< targetZ<<endl;
371  R=McFcorrection(anglein);
372  }
373 
374  G4double y=2.*pi*sin(anglein)*R*FormFactor2ExpHof(anglein)/
375  ((2*As+2.*pow(sin(anglein/2.),2.))*(2.*As+2.*pow(sin(anglein/2.),2.) ));
376 
377  return y/total;
378 }
379 
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381 
383 {
384  //cout<<" G4ScreeningMottCrossSection::SampleCosineTheta ............."<<endl;
385  if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
386 
387  G4double anglemin=std::acos(cosTetMinNuc);
388  G4double anglemax= std::acos(cosTetMaxNuc);
389 
390  static const G4double limit = 1.e-9;
391  if(anglemin < limit) {
392  anglemin = GetScreeningAngle()/10.;
393  if(anglemin < limit) { anglemin = limit; }
394  }
395 
396  // cout<<"................ tkinLab "<< G4BestUnit(tkinLab,"Energy") << " anglemin= "<<anglemin<<endl;
397  //cout<<"anglemax= "<<anglemax<<endl;
398  G4double r =G4UniformRand();
399 
400  G4double loganglemin=log10(anglemin);
401  G4double loganglemax=log10(anglemax);
402  G4double logdangle=0.01;
403 
404  G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);
405 
406  std::vector<G4double> angle;
407  std::vector<G4double> tet;
408  std::vector<G4double> dangle;
409 
410  for(G4int i=0; i<=bins; i++ ){
411  tet.push_back(0);
412  dangle.push_back(0);
413  angle.push_back(pow(10.,loganglemin+logdangle*i));
414  }
415 
416  G4int dim = tet.size();
417  G4double scattangle=0;
418  G4double y=0;
419  G4double dy=0;
420  G4double area=0;
421 
422  for(G4int i=0; i<dim;i++){
423 
424  if(i!=dim-1){
425  dangle[i]=(angle[i+1]-angle[i]);
426  tet[i]=(angle[i+1]+angle[i])/2.;
427  }else if(i==dim-1){
428  break;
429  }
430 
431  y+=AngleDistribution(tet[i])*dangle[i];
432  dy= y-area ;
433  area=y;
434 
435  if(r >=y-dy && r<=y+dy ){
436  scattangle= angle[i] +G4UniformRand()*dangle[i];
437  //cout<<"y "<<y <<" r "<< r << " y/r "<<y/r << endl;
438  break;
439  }
440  }
441  return scattangle;
442 }
443 
444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
445 
447 
448  G4ThreeVector dir(0.0,0.0,1.0);
449 
451 
452  G4double sint = sin(z1);
453  G4double cost = sqrt(1.0 - sint*sint);
454  G4double phi = twopi* G4UniformRand();
455  G4double dirx = sint*cos(phi);
456  G4double diry = sint*sin(phi);
457  G4double dirz = cost;
458 
459  //.......set Trc
460  G4double etot=tkinLab+mass;
461  G4double mass2=targetMass;
462  Trec=(1.0 - cost)* mass2*(etot*etot - mass*mass )/
463  (mass*mass + mass2*mass2+ 2.*mass2*etot);
464 
465  dir.set(dirx,diry,dirz);
466 
467  return dir;
468 }
469 
470 
static const double cm
Definition: G4SIunits.hh:106
const G4double a0
static G4double GetNuclearMass(const G4double A, const G4double Z)
const G4ParticleDefinition * particle
CLHEP::Hep3Vector G4ThreeVector
void Initialise(const G4ParticleDefinition *, G4double cosThetaLim)
void SetupParticle(const G4ParticleDefinition *)
const G4double pi
static const G4double e2
G4double a
Definition: TRTMaterials.hh:39
static const G4double e4
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
#define G4UniformRand()
Definition: Randomize.hh:93
G4double iz
Definition: TRTMaterials.hh:39
static const G4double A[nN]
static const G4double factor
G4double total(Particle const *const p1, Particle const *const p2)
#define DBL_MIN
Definition: templates.hh:75
G4double GetAtomicMassAmu(const G4String &symb) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetMottCoeff(G4double targetZ, G4double coeff[5][6])
double G4double
Definition: G4Types.hh:76
void SetupKinematic(G4double kinEnergy, G4double Z)
static const G4double alpha
#define DBL_MAX
Definition: templates.hh:83