ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ScreeningMottCrossSection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ScreeningMottCrossSection.cc
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. Meth. A488 (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 "Randomize.hh"
71 #include "G4Proton.hh"
72 #include "G4LossTableManager.hh"
73 #include "G4NucleiProperties.hh"
74 #include "G4Element.hh"
75 #include "G4UnitsTable.hh"
76 #include "G4NistManager.hh"
77 #include "G4ThreeVector.hh"
78 #include "G4Pow.hh"
79 #include "G4MottData.hh"
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 
83 static G4double invlog10 = 1./std::log(10.);
84 
85 static const G4double angle[DIMMOTT]={
86 1e-07,1.02329e-07,1.04713e-07,1.07152e-07,1.09648e-07,1.12202e-07,1.14815e-07,1.1749e-07,1.20226e-07,1.23027e-07,1.25893e-07,1.28825e-07,1.31826e-07,1.34896e-07,1.38038e-07,1.41254e-07,1.44544e-07,1.47911e-07,1.51356e-07,1.54882e-07,1.58489e-07,1.62181e-07,1.65959e-07,1.69824e-07,1.7378e-07,1.77828e-07,1.8197e-07,1.86209e-07,1.90546e-07,1.94984e-07,1.99526e-07,2.04174e-07,2.0893e-07,2.13796e-07,2.18776e-07,2.23872e-07,2.29087e-07,2.34423e-07,2.39883e-07,2.45471e-07,2.51189e-07,2.5704e-07,2.63027e-07,2.69153e-07,2.75423e-07,2.81838e-07,2.88403e-07,2.95121e-07,3.01995e-07,3.0903e-07,3.16228e-07,3.23594e-07,3.31131e-07,3.38844e-07,3.46737e-07,3.54813e-07,3.63078e-07,3.71535e-07,3.80189e-07,3.89045e-07,3.98107e-07,4.0738e-07,4.16869e-07,4.2658e-07,4.36516e-07,4.46684e-07,4.57088e-07,4.67735e-07,4.7863e-07,4.89779e-07,5.01187e-07,5.12861e-07,5.24807e-07,5.37032e-07,5.49541e-07,5.62341e-07,5.7544e-07,5.88844e-07,6.0256e-07,6.16595e-07,6.30957e-07,6.45654e-07,6.60693e-07,6.76083e-07,6.91831e-07,7.07946e-07,7.24436e-07,7.4131e-07,7.58578e-07,7.76247e-07,7.94328e-07,8.12831e-07,8.31764e-07,8.51138e-07,8.70964e-07,8.91251e-07,9.12011e-07,9.33254e-07,9.54993e-07,9.77237e-07,1e-06,1.02329e-06,1.04713e-06,1.07152e-06,1.09648e-06,1.12202e-06,1.14815e-06,1.1749e-06,1.20226e-06,1.23027e-06,1.25893e-06,1.28825e-06,1.31826e-06,1.34896e-06,1.38038e-06,1.41254e-06,1.44544e-06,1.47911e-06,1.51356e-06,1.54882e-06,1.58489e-06,1.62181e-06,1.65959e-06,1.69824e-06,1.7378e-06,1.77828e-06,1.8197e-06,1.86209e-06,1.90546e-06,1.94984e-06,1.99526e-06,2.04174e-06,2.0893e-06,2.13796e-06,2.18776e-06,2.23872e-06,2.29087e-06,2.34423e-06,2.39883e-06,2.45471e-06,2.51189e-06,2.5704e-06,2.63027e-06,2.69153e-06,2.75423e-06,2.81838e-06,2.88403e-06,2.95121e-06,3.01995e-06,3.0903e-06,3.16228e-06,3.23594e-06,3.31131e-06,3.38844e-06,3.46737e-06,3.54813e-06,3.63078e-06,3.71535e-06,3.80189e-06,3.89045e-06,3.98107e-06,4.0738e-06,4.16869e-06,4.2658e-06,4.36516e-06,4.46684e-06,4.57088e-06,4.67735e-06,4.7863e-06,4.89779e-06,5.01187e-06,5.12861e-06,5.24807e-06,5.37032e-06,5.49541e-06,5.62341e-06,5.7544e-06,5.88844e-06,6.0256e-06,6.16595e-06,6.30957e-06,6.45654e-06,6.60693e-06,6.76083e-06,6.91831e-06,7.07946e-06,7.24436e-06,7.4131e-06,7.58578e-06,7.76247e-06,7.94328e-06,8.12831e-06,8.31764e-06,8.51138e-06,8.70964e-06,8.91251e-06,9.12011e-06,9.33254e-06,9.54993e-06,9.77237e-06,1e-05,1.02329e-05,1.04713e-05,1.07152e-05,1.09648e-05,1.12202e-05,1.14815e-05,1.1749e-05,1.20226e-05,1.23027e-05,1.25893e-05,1.28825e-05,1.31826e-05,1.34896e-05,1.38038e-05,1.41254e-05,1.44544e-05,1.47911e-05,1.51356e-05,1.54882e-05,1.58489e-05,1.62181e-05,1.65959e-05,1.69824e-05,1.7378e-05,1.77828e-05,1.8197e-05,1.86209e-05,1.90546e-05,1.94984e-05,1.99526e-05,2.04174e-05,2.0893e-05,2.13796e-05,2.18776e-05,2.23872e-05,2.29087e-05,2.34423e-05,2.39883e-05,2.45471e-05,2.51189e-05,2.5704e-05,2.63027e-05,2.69153e-05,2.75423e-05,2.81838e-05,2.88403e-05,2.95121e-05,3.01995e-05,3.0903e-05,3.16228e-05,3.23594e-05,3.31131e-05,3.38844e-05,3.46737e-05,3.54813e-05,3.63078e-05,3.71535e-05,3.80189e-05,3.89045e-05,3.98107e-05,4.0738e-05,4.16869e-05,4.2658e-05,4.36516e-05,4.46684e-05,4.57088e-05,4.67735e-05,4.7863e-05,4.89779e-05,5.01187e-05,5.12861e-05,5.24807e-05,5.37032e-05,5.49541e-05,5.62341e-05,5.7544e-05,5.88844e-05,6.0256e-05,6.16595e-05,6.30957e-05,6.45654e-05,6.60693e-05,6.76083e-05,6.91831e-05,7.07946e-05,7.24436e-05,7.4131e-05,7.58578e-05,7.76247e-05,7.94328e-05,8.12831e-05,8.31764e-05,8.51138e-05,8.70964e-05,8.91251e-05,9.12011e-05,9.33254e-05,9.54993e-05,9.77237e-05,0.0001,0.000102329,0.000104713,0.000107152,0.000109648,0.000112202,0.000114815,0.00011749,0.000120226,0.000123027,0.000125893,0.000128825,0.000131826,0.000134896,0.000138038,0.000141254,0.000144544,0.000147911,0.000151356,0.000154882,0.000158489,0.000162181,0.000165959,0.000169824,0.00017378,0.000177828,0.00018197,0.000186209,0.000190546,0.000194984,0.000199526,0.000204174,0.00020893,0.000213796,0.000218776,0.000223872,0.000229087,0.000234423,0.000239883,0.000245471,0.000251189,0.00025704,0.000263027,0.000269153,0.000275423,0.000281838,0.000288403,0.000295121,0.000301995,0.00030903,0.000316228,0.000323594,0.000331131,0.000338844,0.000346737,0.000354813,0.000363078,0.000371535,0.000380189,0.000389045,0.000398107,0.00040738,0.000416869,0.00042658,0.000436516,0.000446684,0.000457088,0.000467735,0.00047863,0.000489779,0.000501187,0.000512861,0.000524807,0.000537032,0.000549541,0.000562341,0.00057544,0.000588844,0.00060256,0.000616595,0.000630957,0.000645654,0.000660693,0.000676083,0.000691831,0.000707946,0.000724436,0.00074131,0.000758578,0.000776247,0.000794328,0.000812831,0.000831764,0.000851138,0.000870964,0.000891251,0.000912011,0.000933254,0.000954993,0.000977237,0.001,0.00102329,0.00104713,0.00107152,0.00109648,0.00112202,0.00114815,0.0011749,0.00120226,0.00123027,0.00125893,0.00128825,0.00131826,0.00134896,0.00138038,0.00141254,0.00144544,0.00147911,0.00151356,0.00154882,0.00158489,0.00162181,0.00165959,0.00169824,0.0017378,0.00177828,0.0018197,0.00186209,0.00190546,0.00194984,0.00199526,0.00204174,0.0020893,0.00213796,0.00218776,0.00223872,0.00229087,0.00234423,0.00239883,0.00245471,0.00251189,0.0025704,0.00263027,0.00269153,0.00275423,0.00281838,0.00288403,0.00295121,0.00301995,0.0030903,0.00316228,0.00323594,0.00331131,0.00338844,0.00346737,0.00354813,0.00363078,0.00371535,0.00380189,0.00389045,0.00398107,0.0040738,0.00416869,0.0042658,0.00436516,0.00446684,0.00457088,0.00467735,0.0047863,0.00489779,0.00501187,0.00512861,0.00524807,0.00537032,0.00549541,0.00562341,0.0057544,0.00588844,0.0060256,0.00616595,0.00630957,0.00645654,0.00660693,0.00676083,0.00691831,0.00707946,0.00724436,0.0074131,0.00758578,0.00776247,0.00794328,0.00812831,0.00831764,0.00851138,0.00870964,0.00891251,0.00912011,0.00933254,0.00954993,0.00977237,0.01,0.0102329,0.0104713,0.0107152,0.0109648,0.0112202,0.0114815,0.011749,0.0120226,0.0123027,0.0125893,0.0128825,0.0131826,0.0134896,0.0138038,0.0141254,0.0144544,0.0147911,0.0151356,0.0154882,0.0158489,0.0162181,0.0165959,0.0169824,0.017378,0.0177828,0.018197,0.0186209,0.0190546,0.0194984,0.0199526,0.0204174,0.020893,0.0213796,0.0218776,0.0223872,0.0229087,0.0234423,0.0239883,0.0245471,0.0251189,0.025704,0.0263027,0.0269153,0.0275423,0.0281838,0.0288403,0.0295121,0.0301995,0.030903,0.0316228,0.0323594,0.0331131,0.0338844,0.0346737,0.0354813,0.0363078,0.0371535,0.0380189,0.0389045,0.0398107,0.040738,0.0416869,0.042658,0.0436516,0.0446684,0.0457088,0.0467735,0.047863,0.0489779,0.0501187,0.0512861,0.0524807,0.0537032,0.0549541,0.0562341,0.057544,0.0588844,0.060256,0.0616595,0.0630957,0.0645654,0.0660693,0.0676083,0.0691831,0.0707946,0.0724436,0.074131,0.0758578,0.0776247,0.0794328,0.0812831,0.0831764,0.0851138,0.0870964,0.0891251,0.0912011,0.0933254,0.0954993,0.0977237,0.1,0.102329,0.104713,0.107152,0.109648,0.112202,0.114815,0.11749,0.120226,0.123027,0.125893,0.128825,0.131826,0.134896,0.138038,0.141254,0.144544,0.147911,0.151356,0.154882,0.158489,0.162181,0.165959,0.169824,0.17378,0.177828,0.18197,0.186209,0.190546,0.194984,0.199526,0.204174,0.20893,0.213796,0.218776,0.223872,0.229087,0.234423,0.239883,0.245471,0.251189,0.25704,0.263027,0.269153,0.275423,0.281838,0.288403,0.295121,0.301995,0.30903,0.316228,0.323594,0.331131,0.338844,0.346737,0.354813,0.363078,0.371535,0.380189,0.389045,0.398107,0.40738,0.416869,0.42658,0.436516,0.446684,0.457088,0.467735,0.47863,0.489779,0.501187,0.512861,0.524807,0.537032,0.549541,0.562341,0.57544,0.588844,0.60256,0.616595,0.630957,0.645654,0.660693,0.676083,0.691831,0.707946,0.724436,0.74131,0.758578,0.776247,0.794328,0.812831,0.831764,0.851138,0.870964,0.891251,0.912011,0.933254,0.954993,0.977237,1,1.02329,1.04713,1.07152,1.09648,1.12202,1.14815,1.1749,1.20226,1.23027,1.25893,1.28825,1.31826,1.34896,1.38038,1.41254,1.44544,1.47911,1.51356,1.54882,1.58489,1.62181,1.65959,1.69824,1.7378,1.77828,1.8197,1.86209,1.90546,1.94984,1.99526,2.04174,2.0893,2.13796,2.18776,2.23872,2.29087,2.34423,2.39883,2.45471,2.51189,2.5704,2.63027,2.69153,2.75423,2.81838,2.88403,2.95121,3.01995,3.0903};
87 
88 //using namespace std;
89 
90 G4ScreeningMottCrossSection::G4ScreeningMottCrossSection():
91  cosThetaMin(1.0),
92  cosThetaMax(-1.0),
93  alpha(fine_structure_const),
94  htc2(hbarc_squared),
95  e2(electron_mass_c2*classic_electr_radius)
96 {
97  fTotalCross=0;
98 
99  fNistManager = G4NistManager::Instance();
100  fG4pow = G4Pow::GetInstance();
101  particle=nullptr;
102 
103  spin = mass = mu_rel = 0.0;
104  tkinLab = momLab2 = invbetaLab2 = 0.0;
105  tkin = mom2 = invbeta2 = beta = gamma = 0.0;
106 
107  targetMass = As = etag = ecut = 0.0;
108 
109  targetZ = targetA = 0;
110 
111  cosTetMinNuc = cosTetMaxNuc = 0.0;
112 }
113 
114 //....Ooooo0ooooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115 
116 G4ScreeningMottCrossSection::~G4ScreeningMottCrossSection()
117 {}
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120 
121 void G4ScreeningMottCrossSection::Initialise(const G4ParticleDefinition* p,
122  G4double cosThetaLim)
123 {
124  SetupParticle(p);
125  tkin = mom2 = 0.0;
126  ecut = etag = DBL_MAX;
127  particle = p;
128  cosThetaMin = cosThetaLim;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132 
133 void G4ScreeningMottCrossSection::SetupKinematic(G4double ekin, G4int Z)
134 {
135  //...Target
136  const G4int iz = std::min(92, Z);
137  const G4int ia = G4lrint(fNistManager->GetAtomicMassAmu(iz));
138 
139  targetZ = iz;
140  targetA = ia;
141  targetMass = G4NucleiProperties::GetNuclearMass(ia, iz);
142 
143  //G4cout<<"......... targetA "<< targetA <<G4endl;
144  //G4cout<<"......... targetMass "<< targetMass/MeV <<G4endl;
145 
146  // incident particle lab
147  tkinLab = ekin;
148  momLab2 = tkinLab*(tkinLab + 2.0*mass);
149  invbetaLab2 = 1.0 + mass*mass/momLab2;
150 
151  const G4double etot = tkinLab + mass;
152  const G4double ptot = std::sqrt(momLab2);
153  const G4double m12 = mass*mass;
154 
155  // relativistic reduced mass from publucation
156  // A.P. Martynenko, R.N. Faustov, Teoret. mat. Fiz. 64 (1985) 179
157 
158  //incident particle & target nucleus
159  const G4double Ecm = std::sqrt(m12 + targetMass*targetMass + 2.0*etot*targetMass);
160  mu_rel = mass*targetMass/Ecm;
161  const G4double momCM= ptot*targetMass/Ecm;
162  // relative system
163  mom2 = momCM*momCM;
164  const G4double x = mu_rel*mu_rel/mom2;
165  invbeta2 = 1.0 + x;
166  tkin = momCM*std::sqrt(invbeta2) - mu_rel;//Ekin of mu_rel
167  const G4double beta2 = 1./invbeta2;
168  beta = std::sqrt(beta2) ;
169  const G4double gamma2= invbeta2/x;
170  gamma = std::sqrt(gamma2);
171 
172  //Thomas-Fermi screening length
173  const G4double alpha2 = alpha*alpha;
174  const G4double aU = 0.88534*CLHEP::Bohr_radius/fG4pow->Z13(targetZ);
175  const G4double twoR2 = aU*aU;
176  As = 0.25*htc2*(1.13 + 3.76*targetZ*targetZ*invbeta2*alpha2)/(twoR2*mom2);
177 
178  //Integration Angles definition
179  cosTetMinNuc = cosThetaMin;
180  cosTetMaxNuc = cosThetaMax;
181 }
182 
183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
184 
185 G4double G4ScreeningMottCrossSection::FormFactor2ExpHof(G4double sin2t2)
186 {
187  G4double M=targetMass;
188  G4double E=tkinLab;
189  G4double Etot=E+mass;
190  G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
191  G4double T=Tmax*sin2t2;
192  G4double q2=T*(T+2.*M);
193  q2/=htc2;//1/cm2
194  G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targetA)*0.27)*CLHEP::cm;
195  G4double xN= (RN*RN*q2);
196  G4double den=(1.+xN/12.);
197  G4double FN=1./(den*den);
198  G4double form2=(FN*FN);
199 
200  return form2;
201 }
202 
203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
204 
205 G4double G4ScreeningMottCrossSection::FormFactor2Gauss(G4double sin2t2)
206 {
207  G4double M=targetMass;
208  G4double E=tkinLab;
209  G4double Etot=E+mass;
210  G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
211  G4double T=Tmax*sin2t2;
212  G4double q2=T*(T+2.*M);
213  q2/=htc2;//1/cm2
214  G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targetA)*0.27)*CLHEP::cm;
215  G4double xN= (RN*RN*q2);
216 
217  G4double expo=(-xN/6.);
218  G4double FN=G4Exp(expo);
219 
220  G4double form2=(FN*FN);
221 
222  return form2;
223 }
224 
225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
226 
227 G4double G4ScreeningMottCrossSection::FormFactor2UniformHelm(G4double sin2t2)
228 {
229  G4double M=targetMass;
230  G4double E=tkinLab;
231  G4double Etot=E+mass;
232  G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
233  G4double T=Tmax*sin2t2;
234  G4double q2=T*(T+2.*M);
235  q2=q2/(htc2*0.01);//1/cm2
236 
237  G4double q=std::sqrt(q2);
238  G4double R0=1.2E-13*fG4pow->Z13(targetA);
239  G4double R1=2.0E-13;
240 
241  G4double x0=q*R0;
242  G4double F0=(3./fG4pow->powN(x0,3))*(std::sin(x0)-x0*std::cos(x0));
243 
244  G4double x1=q*R1;
245  G4double F1=(3./fG4pow->powN(x1,3))*(std::sin(x1)-x1*std::cos(x1));
246 
247  G4double F=F0*F1;
248 
249  G4double form2=(F*F);
250 
251  return form2;
252 }
253 
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255 
256 G4double G4ScreeningMottCrossSection::McFcorrection(G4double sin2t2)
257 {
258  const G4double sint = std::sqrt(sin2t2);
259  return 1.-beta*beta*sin2t2 + targetZ*alpha*beta*pi*sint*(1.-sint);
260 }
261 
262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263 
264 G4double G4ScreeningMottCrossSection::RatioMottRutherford(G4double angles)
265 {
266  return RatioMottRutherfordCosT(std::sqrt(1. - std::cos(angles)));
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270 
271 G4double G4ScreeningMottCrossSection::RatioMottRutherfordCosT(G4double fcost)
272 {
273  G4double R(0.);
274  const G4double shift = 0.7181228;
275  G4double beta0 = beta - shift;
276 
277  G4double b0 = 1.0;
278  G4double b[6];
279  for(G4int i=0; i<6; ++i) {
280  b[i] = b0;
281  b0 *= beta0;
282  }
283 
284  b0 = 1.0;
285  for(G4int j=0; j<=4; ++j) {
286  G4double a = 0.0;
287  for(G4int k=0; k<=5; ++k){
288  a += fMottCoef[targetZ][j][k]*b[k];
289  }
290  R += a*b0;
291  b0 *= fcost;
292  }
293  return R;
294 }
295 
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
297 
298 G4double G4ScreeningMottCrossSection::GetTransitionRandom()
299 {
300  G4double x = G4Log(tkinLab)*invlog10;
301  G4double res(0.), x0(1.0);
302  for(G4int i=0; i<11; ++i) {
303  res += fPRM[targetZ][i]*x0;
304  x0 *= x;
305  }
306  //G4cout << "Z= " << targetZ << " x0= " << x0 << " res= " << res << G4endl;
307  return res;
308 }
309 
310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
311 
312 G4double
313 G4ScreeningMottCrossSection::DifferentialXSection(G4int idx, G4int form)
314 {
315  const G4double ang = angle[idx];
316  const G4double z1 = 1.0 - std::cos(ang);
317  G4double step;
318  if(0 == idx) {
319  step = 0.5*(angle[1] + angle[0]);
320  } else if(DIMMOTT == idx + 1) {
321  step = CLHEP::pi - 0.5*(angle[DIMMOTT-2] + angle[DIMMOTT-1]);
322  } else {
323  step = 0.5*(angle[idx+1] - angle[idx-1]);
324  }
325 
326  G4double F2 = 1.;
327  switch (form) {
328  case 1:
329  F2 = FormFactor2ExpHof(z1*0.5);
330  break;
331  case 2:
332  F2 = FormFactor2Gauss(z1*0.5);
333  break;
334  case 3:
335  F2 = FormFactor2UniformHelm(z1*0.5);
336  break;
337  }
338 
339  const G4double R = RatioMottRutherfordCosT(std::sqrt(z1));
340 
341  G4double den = 2.*As + z1;
342  G4double func = 1./(den*den);
343 
344  G4double fatt = (G4double)targetZ/(mu_rel*gamma*beta*beta);
345  G4double sigma= e2*e2*fatt*fatt*func;
346  G4double pi2sintet = CLHEP::twopi*std::sqrt(z1*(2. - z1));
347 
348  G4double Xsec = std::max(pi2sintet*F2*R*sigma*step, 0.);
349  return Xsec;
350 }
351 
352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
353 
354 G4double
355 G4ScreeningMottCrossSection::NuclearCrossSection(G4int form, G4int fast)
356 {
357  fTotalCross=0.;
358  if(cosTetMaxNuc >= cosTetMinNuc) return fTotalCross;
359  if(0 == cross.size()) { cross.resize(DIMMOTT, 0.0); }
360 
361  //G4cout<<"MODEL: "<<fast<<G4endl;
362  //************ PRECISE COMPUTATION
363  if(fast == 0){
364 
365  for(G4int i=0; i<DIMMOTT; ++i){
366  G4double y = DifferentialXSection(i, form);
367  fTotalCross += y;
368  cross[i] = fTotalCross;
369  if(fTotalCross*1.e-9 > y) {
370  for(G4int j=i+1; j<DIMMOTT; ++j) {
371  cross[j] = fTotalCross;
372  }
373  break;
374  }
375  }
376  //**************** FAST COMPUTATION
377  } else if(fast == 1) {
378 
379  G4double p0 = electron_mass_c2*classic_electr_radius;
380  G4double coeff = twopi*p0*p0;
381 
382  G4double fac = coeff*targetZ*targetZ*invbeta2/mom2;
383 
384  G4double x = 1.0 - cosTetMinNuc;
385  G4double x1 = x + 2*As;
386 
387  // scattering with nucleus
388  fTotalCross = fac*(cosTetMinNuc - cosTetMaxNuc)/
389  (x1*(1.0 - cosTetMaxNuc + 2*As));
390  }
391  /*
392  G4cout << "Energy(MeV): " << tkinLab/CLHEP::MeV
393  << " Total Cross(mb): " << fTotalCross
394  << " form= " << form << " fast= " << fast << G4endl;
395  */
396  return fTotalCross;
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 
401 G4double
402 G4ScreeningMottCrossSection::GetScatteringAngle(G4int form, G4int fast)
403 {
404  G4double scattangle = 0.0;
405  G4double r = G4UniformRand();
406  //************ PRECISE COMPUTATION
407  if(fast == 0){
408  //G4cout << "r= " << r << " tot= " << fTotalCross << G4endl;
409  r *= fTotalCross;
410  for(G4int i=0; i<DIMMOTT; ++i){
411  //G4cout << i << ". r= " << r << " xs= " << cross[i] << G4endl;
412  if(r <= cross[i]) {
413  scattangle = ComputeAngle(i, r);
414  break;
415  }
416  }
417 
418  //**************** FAST COMPUTATION
419  } else if(fast == 1) {
420 
421  G4double limit = GetTransitionRandom();
422  if(limit > 0.0) {
423  G4double Sz = 2*As;
424  G4double x = Sz-((Sz*(2+Sz))/(Sz+2*limit))+1.0;
425  //G4cout << "limit= " << limit << " Sz= " << Sz << " x= " << x << G4endl;
426  G4double angle_limit = (std::abs(x) < 1.0) ? std::acos(x) : 0.0;
427  //G4cout<<"ANGLE LIMIT: "<<angle_limit<<" x= " << x << G4endl;
428 
429  if(r > limit && angle_limit != 0.0) {
430  x = Sz-((Sz*(2+Sz))/(Sz+2*r))+1.0;
431  scattangle = (x >= 1.0) ? 0.0 : ((x <= -1.0) ? pi : std::acos(x));
432  //G4cout<<"FAST: scattangle= "<< scattangle <<G4endl;
433  }
434  } else {
435  //G4cout<<"SLOW: total= "<<fTotalCross<< G4endl;
436  r *= fTotalCross;
437  G4double tot = 0.0;
438  for(G4int i=0; i<DIMMOTT; ++i) {
439  G4double xs = DifferentialXSection(i, form);
440  tot += xs;
441  cross[i] = tot;
442  if(r <= tot) {
443  scattangle = ComputeAngle(i, r);
444  break;
445  }
446  }
447  }
448  }
449  //************************************************
450  //G4cout<<"Energy(MeV): "<<tkinLab/MeV<<" SCATTANGLE: "<<scattangle<<G4endl;
451  return scattangle;
452 }
453 
454 G4double G4ScreeningMottCrossSection::ComputeAngle(G4int i, G4double& r)
455 {
456  G4double x1, x2, y;
457  if(0 == i) {
458  x1 = 0.0;
459  x2 = 0.5*(angle[0] + angle[1]);
460  y = cross[0];
461  } else if(i == DIMMOTT-1) {
462  x1 = 0.5*(angle[i] + angle[i-1]);
463  x2 = CLHEP::pi;
464  y = cross[i] - cross[i-1];
465  r -= cross[i-1];
466  } else {
467  x1 = 0.5*(angle[i] + angle[i-1]);
468  x2 = 0.5*(angle[i] + angle[i+1]);
469  y = cross[i] - cross[i-1];
470  r -= cross[i-1];
471  }
472  //G4cout << i << ". r= " << r << " y= " << y
473  // << " x1= " << " x2= " << x2 << G4endl;
474  return x1 + (x2 - x1)*r/y;
475 }