Geant4  10.00.p03
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 
91  TotalCross=0;
92 
94  particle=0;
95 
96  spin = mass = mu_rel=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 
115 }
116 
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 
120 {
121  delete mottcoeff;
122 }
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 
126  G4double CosThetaLim)
127 {
128 
129  SetupParticle(p);
130  tkin = targetZ = mom2 = DBL_MIN;
131  ecut = etag = DBL_MAX;
132  particle = p;
133  cosThetaMin = CosThetaLim;
134 
135 }
136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138 
139  G4double alpha2=alpha*alpha;
140  //Bohr radius
141  G4double a0= Bohr_radius ;//0.529e-8*cm;
142  //Thomas-Fermi screening length
143  G4double aU=0.88534*a0/pow(targetZ,1./3.);
144  G4double twoR2=aU*aU;
145 
146  G4double factor= 1.13 + 3.76*targetZ*targetZ*invbeta2*alpha2;
147  As=0.25*(htc2)/(twoR2*mom2)*factor;
148 
149 
150 
151 
152 //cout<<"0k .........................As "<<As<<endl;
153 
154 }
155 
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157 
159 
161 
162  //cout<<" As "<<As<<endl;
163  if(As < 0.0) { As = 0.0; }
164  else if(As > 1.0) { As = 1.0; }
165 
166  G4double screenangle=2.*asin(sqrt(As));
167 
168  // cout<<" screenangle "<< screenangle <<endl;
169 
170  if(screenangle>=pi) screenangle=pi;
171 
172 
173 return screenangle;
174 
175 }
176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
177 
179 {
180 
181  //...Target
182  G4int iz = G4int(Z);
184  G4int ia = G4int(A);
186 
187  targetZ = Z;
189  targetMass= mass2;
190 
192 
193  //cout<<"......... targetA "<< targetA <<endl;
194  //cout<<"......... targetMass "<< targetMass/MeV <<endl;
195 
196  // incident particle lab
197  tkinLab = ekin;
198  momLab2 = tkinLab*(tkinLab + 2.0*mass);
199  invbetaLab2 = 1.0 + mass*mass/momLab2;
200 
201  G4double etot = tkinLab + mass;
202  G4double ptot = sqrt(momLab2);
203  G4double m12 = mass*mass;
204 
205 
206  // relativistic reduced mass from publucation
207  // A.P. Martynenko, R.N. Faustov, Teoret. mat. Fiz. 64 (1985) 179
208 
209  //incident particle & target nucleus
210  G4double Ecm=sqrt(m12 + mass2*mass2 + 2.0*etot*mass2);
211  mu_rel=mass*mass2/Ecm;
212  G4double momCM= ptot*mass2/Ecm;
213  // relative system
214  mom2 = momCM*momCM;
215  invbeta2 = 1.0 + mu_rel*mu_rel/mom2;
216  tkin = momCM*sqrt(invbeta2) - mu_rel;//Ekin of mu_rel
217  G4double beta2=1./invbeta2;
218  beta=std::sqrt(beta2) ;
219  G4double gamma2= 1./(1.-beta2);
220  gamma=std::sqrt(gamma2);
221 
222 
223 
224  //.........................................................
225 
226 
227  G4double screenangle=GetScreeningAngle()/10.;
228  //cout<<" screenangle [rad] "<<screenangle/rad <<endl;
229 
230  cosTetMinNuc =min( cosThetaMin ,cos(screenangle));
232 
233 //cout<<"ok..............mu_rel "<<mu_rel<<endl;
234 
235 }
236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
237 
239 
240 
241  G4double M=targetMass;
242  G4double E=tkinLab;
243  G4double Etot=E+mass;
244  G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
245  G4double T=Tmax*pow(sin(angle/2.),2.);
246  G4double q2=T*(T+2.*M);
247  q2/=htc2;//1/cm2
248  G4double RN=1.27e-13*pow(targetA,0.27)*cm;
249  G4double xN= (RN*RN*q2);
250  G4double den=(1.+xN/12.);
251  G4double FN=1./(den*den);
252  G4double form2=(FN*FN);
253 
254 
255  return form2;
256 
257 //cout<<"..................... form2 "<< form2<<endl;
258 }
259 
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
261 
263 
264 
265  G4double beta2=1./invbeta2;
266  G4double sintmezzi=std::sin(angle/2.);
267  G4double sin2tmezzi = sintmezzi*sintmezzi;
268  G4double R=1.-beta2*sin2tmezzi + targetZ*alpha*beta*pi*sintmezzi*(1.-sintmezzi);
269 
270 
271 return R;
272 }
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
275 
276 
277  G4double R=0;
278  G4double fcost=std::sqrt((1. -cos(angle)));
279  G4double a[5];
280  G4double shift=0.7181228;
281  G4double beta0= beta -shift;
282 
283  for(G4int j=0 ;j<=4;j++){
284  a[j]=0;
285  }
286 
287  for(G4int j=0 ;j<=4;j++){
288  for(G4int k=0;k<=5;k++ ){
289  a[j]+=coeffb[j][k]*pow(beta0,k);
290  }
291  }
292 
293  for(G4int j=0 ;j<=4 ;j++){
294  R+=a[j]* pow(fcost,j);
295  }
296 
297 
298 
299 return R;
300 
301 }
302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303 
305 {
306  if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
307 
308  TotalCross=0;
309 
310  G4double anglemin =std::acos(cosTetMinNuc);
311  G4double anglemax =std::acos(cosTetMaxNuc);
312 
313  static const G4double limit = 1.e-9;
314  if(anglemin < limit) {
315  anglemin = GetScreeningAngle()/10.;
316  if(anglemin < limit) { anglemin = limit; }
317  }
318 
319  //cout<<" anglemin "<< anglemin <<endl;
320 
321  G4double loganglemin=log10(anglemin);
322  G4double loganglemax=log10(anglemax);
323  G4double logdangle=0.01;
324 
325  G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);
326 
327 
328  vector<G4double> angle;
329  vector<G4double> tet;
330  vector<G4double> dangle;
331  vector<G4double> cross;
332 
333  for(G4int i=0; i<=bins; i++ ){
334  tet.push_back(0);
335  dangle.push_back(0);
336  angle.push_back(pow(10.,loganglemin+logdangle*i));
337  cross.push_back(0);
338  }
339 
340 
341  G4int dim = tet.size();
342  //cout<<"dim--- "<<dim<<endl;
343 
344 
345  for(G4int i=0; i<dim;i++){
346 
347  if(i!=dim-1){
348  dangle[i]=(angle[i+1]-angle[i]);
349  tet[i]=(angle[i+1]+angle[i])/2.;
350  }else if(i==dim-1){
351  break;
352  }
353 
354 
355  G4double R=0;
356  G4double F2=FormFactor2ExpHof(tet[i]);
357 
358  if (coeffb[0][0]!=0){
359  //cout<<" Mott....targetZ "<< targetZ<<endl;
360  R=RatioMottRutherford(tet[i]);
361  }
362 
363  else if (coeffb[0][0]==0){
364  // cout<<" McF.... targetZ "<< targetZ<<endl;
365  R=McFcorrection(tet[i]);
366  }
367 
368 
369 //cout<<"----------------- R "<<R<<" F2 "<<F2<<endl;
370 
371 
372 // cout<<"angle "<<tet[i] << " F2 "<<F2<<endl;
373 
374  G4double e4=e2*e2;
375  G4double den=2.*As+2.*pow(sin(tet[i]/2.),2.);
376  G4double func=1./(den*den);
377 
379  G4double sigma=e4*fatt*fatt*func;
380  cross[i]=F2*R*sigma;
381  G4double pi2sintet=2.*pi*sin(tet[i]);
382 
383  TotalCross+=pi2sintet*cross[i]*dangle[i];
384  }//end integral
385 
386 
387 //cout<< "ok ......... TotalCross "<<TotalCross<<endl;
388 
389 
390 return TotalCross;
391 
392 }
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
395 
396 
399  G4double fatt2=fatt*fatt;
400  total/=fatt2;
401 
402  G4double R=0;
403  if (coeffb[0][0]!=0){
404  // cout<<" Mott....targetZ "<< targetZ<<endl;
405  R=RatioMottRutherford(anglein);
406  }
407 
408  else if (coeffb[0][0]==0){
409  // cout<<" McF.... targetZ "<< targetZ<<endl;
410  R=McFcorrection(anglein);
411  }
412 
413 
414  G4double y=2.*pi*sin(anglein)*R*FormFactor2ExpHof(anglein)/
415  ((2*As+2.*pow(sin(anglein/2.),2.))*(2.*As+2.*pow(sin(anglein/2.),2.) ));
416 
417 return y/total;
418 
419 }
420 
421 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
422 
424 {
425 
426 
427 //cout<<" G4ScreeningMottCrossSection::SampleCosineTheta ............."<<endl;
428 
429  if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
430 
431  G4double anglemin=std::acos(cosTetMinNuc);
432  G4double anglemax= std::acos(cosTetMaxNuc);
433 
434  static const G4double limit = 1.e-9;
435  if(anglemin < limit) {
436  anglemin = GetScreeningAngle()/10.;
437  if(anglemin < limit) { anglemin = limit; }
438  }
439 
440 // cout<<"................ tkinLab "<< G4BestUnit(tkinLab,"Energy") << " anglemin= "<<anglemin<<endl;
441  //cout<<"anglemax= "<<anglemax<<endl;
442  G4double r =G4UniformRand();
443 
444  G4double loganglemin=log10(anglemin);
445  G4double loganglemax=log10(anglemax);
446  G4double logdangle=0.01;
447 
448  G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);
449 
450  std::vector<G4double> angle;
451  std::vector<G4double> tet;
452  std::vector<G4double> dangle;
453 
454  for(G4int i=0; i<=bins; i++ ){
455  tet.push_back(0);
456  dangle.push_back(0);
457  angle.push_back(pow(10.,loganglemin+logdangle*i));
458  }
459 
460 
461  G4int dim = tet.size();
462  G4double scattangle=0;
463  G4double y=0;
464  G4double dy=0;
465  G4double area=0;
466 
467  for(G4int i=0; i<dim;i++){
468 
469  if(i!=dim-1){
470  dangle[i]=(angle[i+1]-angle[i]);
471  tet[i]=(angle[i+1]+angle[i])/2.;
472  }else if(i==dim-1){
473  break;
474  }
475 
476  y+=AngleDistribution(tet[i])*dangle[i];
477  dy= y-area ;
478  area=y;
479 
480  if(r >=y-dy && r<=y+dy ){
481  scattangle= angle[i] +G4UniformRand()*dangle[i];
482  //cout<<"y "<<y <<" r "<< r << " y/r "<<y/r << endl;
483  break;
484 
485  }
486 
487  }
488 
489 
490  return scattangle;
491 
492 
493 }
494 
495 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
496 
498 
499 G4ThreeVector dir(0.0,0.0,1.0);
500 
501 
503 
504  G4double sint = sin(z1);
505  G4double cost = sqrt(1.0 - sint*sint);
506  G4double phi = twopi* G4UniformRand();
507  G4double dirx = sint*cos(phi);
508  G4double diry = sint*sin(phi);
509  G4double dirz = cost;
510 
511 
512  //.......set Trc
513  G4double etot=tkinLab+mass;
514  G4double mass2=targetMass;
515  Trec=(1.0 - cost)* mass2*(etot*etot - mass*mass )/
516  (mass*mass + mass2*mass2+ 2.*mass2*etot);
517 
518  dir.set(dirx,diry,dirz);
519 
520 return dir;
521 
522 }
523 
524 
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:87
G4double iz
Definition: TRTMaterials.hh:39
static const G4double A[nN]
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