| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4ScreeningMottCrossSection.cc 26 // G4ScreeningMottCrossSection.cc
27 //-------------------------------------------- 27 //-------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class header file 29 // GEANT4 Class header file
30 // 30 //
31 // File name: G4ScreeningMottCrossSection 31 // File name: G4ScreeningMottCrossSection
32 // 32 //
33 // Author: Cristina Consolandi 33 // Author: Cristina Consolandi
34 // 34 //
35 // Creation date: 20.10.2011 << 35 // Creation date: 20.10.2011
36 // 36 //
37 // Modifications: 37 // Modifications:
38 // 27-05-2012 Added Analytic Fitting to the Mo 38 // 27-05-2012 Added Analytic Fitting to the Mott Cross Section by means of G4MottCoefficients class.
39 // 39 //
40 // 40 //
41 // Class Description: 41 // Class Description:
42 // Computation of electron Coulomb Scattering 42 // Computation of electron Coulomb Scattering Cross Section.
43 // Suitable for high energy electrons and lig << 43 // Suitable for high energy electrons and light target materials.
44 // 44 //
45 // Reference: 45 // Reference:
46 // M.J. Boschini et al. 46 // M.J. Boschini et al.
47 // "Non Ionizing Energy Loss induced by El 47 // "Non Ionizing Energy Loss induced by Electrons in the Space Environment"
48 // Proc. of the 13th International Confer << 48 // Proc. of the 13th International Conference on Particle Physics and Advanced Technology
49 // (13th ICPPAT, Como 3-7/10/2011), World 49 // (13th ICPPAT, Como 3-7/10/2011), World Scientific (Singapore).
50 // Available at: http://arxiv.org/abs/1111.40 50 // Available at: http://arxiv.org/abs/1111.4042v4
51 // 51 //
52 // 1) Mott Differential Cross Section App << 52 // 1) Mott Differential Cross Section Approximation:
53 // For Target material up to Z=92 (U): 53 // For Target material up to Z=92 (U):
54 // As described in http://arxiv.org/ab << 54 // As described in http://arxiv.org/abs/1111.4042v4
55 // par. 2.1 , eq. (16)-(17) 55 // par. 2.1 , eq. (16)-(17)
56 // Else (Z>92): 56 // Else (Z>92):
57 // W. A. McKinley and H. Fashbach, Phys. R 57 // W. A. McKinley and H. Fashbach, Phys. Rev. 74, (1948) 1759.
58 // 2) Screening coefficient: << 58 // 2) Screening coefficient:
59 // vomn G. Moliere, Z. Naturforsh A2 (194 59 // vomn G. Moliere, Z. Naturforsh A2 (1947), 133-145; A3 (1948), 78.
60 // 3) Nuclear Form Factor: << 60 // 3) Nuclear Form Factor:
61 // A.V. Butkevich et al. Nucl. Instr. Met << 61 // A.V. Butkevich et al. Nucl. Instr. and Meth. in Phys. Res. A 488 (2002), 282-294.
62 // 62 //
63 // ------------------------------------------- << 63 // -------------------------------------------------------------------------------------
64 // 64 //
65 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 66
67 #include "G4ScreeningMottCrossSection.hh" 67 #include "G4ScreeningMottCrossSection.hh"
68 #include "G4PhysicalConstants.hh" 68 #include "G4PhysicalConstants.hh"
69 #include "G4SystemOfUnits.hh" 69 #include "G4SystemOfUnits.hh"
>> 70 #include "G4MottCoefficients.hh"
70 #include "Randomize.hh" 71 #include "Randomize.hh"
71 #include "G4Proton.hh" 72 #include "G4Proton.hh"
72 #include "G4LossTableManager.hh" 73 #include "G4LossTableManager.hh"
73 #include "G4NucleiProperties.hh" 74 #include "G4NucleiProperties.hh"
74 #include "G4Element.hh" 75 #include "G4Element.hh"
75 #include "G4UnitsTable.hh" 76 #include "G4UnitsTable.hh"
76 #include "G4NistManager.hh" <<
77 #include "G4ThreeVector.hh" <<
78 #include "G4Pow.hh" <<
79 #include "G4MottData.hh" <<
80 77
81 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 79
83 static G4double invlog10 = 1./std::log(10.); <<
84 80
85 static const G4double angle[DIMMOTT]={ << 81 using namespace std;
86 1e-07,1.02329e-07,1.04713e-07,1.07152e-07,1.09 <<
87 <<
88 //using namespace std; <<
89 82
90 G4ScreeningMottCrossSection::G4ScreeningMottCr 83 G4ScreeningMottCrossSection::G4ScreeningMottCrossSection():
91 cosThetaMin(1.0), 84 cosThetaMin(1.0),
92 cosThetaMax(-1.0), 85 cosThetaMax(-1.0),
93 alpha(fine_structure_const), 86 alpha(fine_structure_const),
94 htc2(hbarc_squared), 87 htc2(hbarc_squared),
95 e2(electron_mass_c2*classic_electr_radius) << 88 e2(electron_mass_c2*classic_electr_radius)
96 { 89 {
97 fTotalCross=0; <<
98 90
99 fNistManager = G4NistManager::Instance(); << 91 TotalCross=0;
100 fG4pow = G4Pow::GetInstance(); <<
101 particle=nullptr; <<
102 92
103 spin = mass = mu_rel = 0.0; << 93 fNistManager = G4NistManager::Instance();
104 tkinLab = momLab2 = invbetaLab2 = 0.0; << 94 particle=0;
105 tkin = mom2 = invbeta2 = beta = gamma = 0.0; <<
106 95
107 targetMass = As = etag = ecut = 0.0; << 96 spin = mass = mu_rel=0;
>> 97 tkinLab = momLab2 = invbetaLab2=0;
>> 98 tkin = mom2 = invbeta2=beta=gamma=0;
108 99
109 targetZ = targetA = 0; << 100 Trec=targetZ = targetMass = As =0;
>> 101 etag = ecut = 0.0;
110 102
111 cosTetMinNuc = cosTetMaxNuc = 0.0; << 103 targetA = 0;
112 } << 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 113
114 //....Ooooo0ooooo........oooOO0OOooo........oo <<
115 114
116 G4ScreeningMottCrossSection::~G4ScreeningMottC <<
117 115
>> 116
>> 117 }
>> 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 119
>> 120 G4ScreeningMottCrossSection::~G4ScreeningMottCrossSection()
>> 121 {}
118 //....oooOO0OOooo........oooOO0OOooo........oo 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 123
120 void G4ScreeningMottCrossSection::Initialise(c 124 void G4ScreeningMottCrossSection::Initialise(const G4ParticleDefinition* p,
121 G << 125 G4double CosThetaLim)
122 { 126 {
123 SetupParticle(p); <<
124 tkin = mom2 = 0.0; <<
125 ecut = etag = DBL_MAX; <<
126 particle = p; <<
127 cosThetaMin = cosThetaLim; <<
128 } <<
129 127
>> 128 SetupParticle(p);
>> 129 tkin = targetZ = mom2 = DBL_MIN;
>> 130 ecut = etag = DBL_MAX;
>> 131 particle = p;
>> 132 cosThetaMin = CosThetaLim;
>> 133
>> 134 }
130 //....oooOO0OOooo........oooOO0OOooo........oo 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;
131 152
132 void G4ScreeningMottCrossSection::SetupKinemat <<
133 { <<
134 //...Target <<
135 const G4int iz = std::min(92, Z); <<
136 const G4int ia = G4lrint(fNistManager->GetAt <<
137 <<
138 targetZ = iz; <<
139 targetA = ia; <<
140 targetMass = G4NucleiProperties::GetNuclearM <<
141 <<
142 //G4cout<<"......... targetA "<< targetA <<G <<
143 //G4cout<<"......... targetMass "<< targetMa <<
144 <<
145 // incident particle lab <<
146 tkinLab = ekin; <<
147 momLab2 = tkinLab*(tkinLab + 2.0*mass); <<
148 invbetaLab2 = 1.0 + mass*mass/momLab2; <<
149 <<
150 const G4double etot = tkinLab + mass; <<
151 const G4double ptot = std::sqrt(momLab2); <<
152 const G4double m12 = mass*mass; <<
153 <<
154 // relativistic reduced mass from publucatio <<
155 // A.P. Martynenko, R.N. Faustov, Teoret. ma <<
156 <<
157 //incident particle & target nucleus <<
158 const G4double Ecm = std::sqrt(m12 + targetM <<
159 mu_rel = mass*targetMass/Ecm; <<
160 const G4double momCM= ptot*targetMass/Ecm; <<
161 // relative system <<
162 mom2 = momCM*momCM; <<
163 const G4double x = mu_rel*mu_rel/mom2; <<
164 invbeta2 = 1.0 + x; <<
165 tkin = momCM*std::sqrt(invbeta2) - mu_rel;// <<
166 const G4double beta2 = 1./invbeta2; <<
167 beta = std::sqrt(beta2) ; <<
168 const G4double gamma2= invbeta2/x; <<
169 gamma = std::sqrt(gamma2); <<
170 <<
171 //Thomas-Fermi screening length <<
172 const G4double alpha2 = alpha*alpha; <<
173 const G4double aU = 0.88534*CLHEP::Bohr_radi <<
174 const G4double twoR2 = aU*aU; <<
175 As = 0.25*htc2*(1.13 + 3.76*targetZ*targetZ* <<
176 <<
177 //Integration Angles definition <<
178 cosTetMinNuc = cosThetaMin; <<
179 cosTetMaxNuc = cosThetaMax; <<
180 } 153 }
181 154
182 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183 156
184 G4double G4ScreeningMottCrossSection::FormFact << 157 G4double G4ScreeningMottCrossSection::GetScreeningAngle(){
185 { <<
186 G4double M=targetMass; <<
187 G4double E=tkinLab; <<
188 G4double Etot=E+mass; <<
189 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+ <<
190 G4double T=Tmax*sin2t2; <<
191 G4double q2=T*(T+2.*M); <<
192 q2/=htc2;//1/cm2 <<
193 G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targ <<
194 G4double xN= (RN*RN*q2); <<
195 G4double den=(1.+xN/12.); <<
196 G4double FN=1./(den*den); <<
197 G4double form2=(FN*FN); <<
198 158
199 return form2; << 159 SetScreeningCoefficient();
200 } << 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;
201 170
>> 171
>> 172 return screenangle;
>> 173
>> 174 }
202 //....oooOO0OOooo........oooOO0OOooo........oo 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203 176
204 G4double G4ScreeningMottCrossSection::FormFact << 177 void G4ScreeningMottCrossSection::SetupKinematic(G4double ekin, G4double Z )
205 { 178 {
206 G4double M=targetMass; <<
207 G4double E=tkinLab; <<
208 G4double Etot=E+mass; <<
209 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+ <<
210 G4double T=Tmax*sin2t2; <<
211 G4double q2=T*(T+2.*M); <<
212 q2/=htc2;//1/cm2 <<
213 G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targ <<
214 G4double xN= (RN*RN*q2); <<
215 179
216 G4double expo=(-xN/6.); << 180 //...Target
217 G4double FN=G4Exp(expo); << 181 G4int iz = G4int(Z);
>> 182 G4double A = fNistManager->GetAtomicMassAmu(iz);
>> 183 G4int ia = G4int(A);
>> 184 G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);
218 185
219 G4double form2=(FN*FN); << 186 targetZ = Z;
>> 187 targetA = fNistManager->GetAtomicMassAmu(iz);
>> 188 targetMass= mass2;
220 189
221 return form2; << 190 mottcoeff->SetMottCoeff(targetZ, coeffb);
222 } << 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);
223 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 }
224 //....oooOO0OOooo........oooOO0OOooo........oo 235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
225 236
226 G4double G4ScreeningMottCrossSection::FormFact << 237 G4double G4ScreeningMottCrossSection::FormFactor2ExpHof(G4double angle){
227 { <<
228 G4double M=targetMass; <<
229 G4double E=tkinLab; <<
230 G4double Etot=E+mass; <<
231 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+ <<
232 G4double T=Tmax*sin2t2; <<
233 G4double q2=T*(T+2.*M); <<
234 q2=q2/(htc2*0.01);//1/cm2 <<
235 <<
236 G4double q=std::sqrt(q2); <<
237 G4double R0=1.2E-13*fG4pow->Z13(targetA); <<
238 G4double R1=2.0E-13; <<
239 <<
240 G4double x0=q*R0; <<
241 G4double F0=(3./fG4pow->powN(x0,3))*(std::si <<
242 238
243 G4double x1=q*R1; <<
244 G4double F1=(3./fG4pow->powN(x1,3))*(std::si <<
245 239
246 G4double F=F0*F1; << 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);
247 252
248 G4double form2=(F*F); <<
249 253
250 return form2; 254 return form2;
>> 255
>> 256 //cout<<"..................... form2 "<< form2<<endl;
251 } 257 }
252 258
253 //....oooOO0OOooo........oooOO0OOooo........oo 259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254 260
255 G4double G4ScreeningMottCrossSection::McFcorre << 261 G4double G4ScreeningMottCrossSection::McFcorrection(G4double angle ){
256 { <<
257 const G4double sint = std::sqrt(sin2t2); <<
258 return 1.-beta*beta*sin2t2 + targetZ*alpha*b <<
259 } <<
260 262
261 //....oooOO0OOooo........oooOO0OOooo........oo <<
262 263
263 G4double G4ScreeningMottCrossSection::RatioMot << 264 G4double beta2=1./invbeta2;
264 { << 265 G4double sintmezzi=std::sin(angle/2.);
265 return RatioMottRutherfordCosT(std::sqrt(1. << 266 G4double sin2tmezzi = sintmezzi*sintmezzi;
266 } << 267 G4double R=1.-beta2*sin2tmezzi + targetZ*alpha*beta*pi*sintmezzi*(1.-sintmezzi);
267 268
268 //....oooOO0OOooo........oooOO0OOooo........oo <<
269 269
270 G4double G4ScreeningMottCrossSection::RatioMot << 270 return R;
271 { <<
272 G4double R(0.); <<
273 const G4double shift = 0.7181228; <<
274 G4double beta0 = beta - shift; <<
275 <<
276 G4double b0 = 1.0; <<
277 G4double b[6]; <<
278 for(G4int i=0; i<6; ++i) { <<
279 b[i] = b0; <<
280 b0 *= beta0; <<
281 } <<
282 <<
283 b0 = 1.0; <<
284 for(G4int j=0; j<=4; ++j) { <<
285 G4double a = 0.0; <<
286 for(G4int k=0; k<=5; ++k){ <<
287 a += fMottCoef[targetZ][j][k]*b[k]; <<
288 } <<
289 R += a*b0; <<
290 b0 *= fcost; <<
291 } <<
292 return R; <<
293 } 271 }
294 <<
295 //....oooOO0OOooo........oooOO0OOooo........oo 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 296
297 G4double G4ScreeningMottCrossSection::GetTrans <<
298 { <<
299 G4double x = G4Log(tkinLab)*invlog10; <<
300 G4double res(0.), x0(1.0); <<
301 for(G4int i=0; i<11; ++i) { <<
302 res += fPRM[targetZ][i]*x0; <<
303 x0 *= x; <<
304 } <<
305 //G4cout << "Z= " << targetZ << " x0= " << x <<
306 return res; <<
307 } <<
308 297
>> 298 return R;
>> 299
>> 300 }
309 //....oooOO0OOooo........oooOO0OOooo........oo 301 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
310 302
311 G4double << 303 G4double G4ScreeningMottCrossSection::NuclearCrossSection()
312 G4ScreeningMottCrossSection::DifferentialXSect <<
313 { 304 {
314 const G4double ang = angle[idx]; << 305 if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
315 const G4double z1 = 1.0 - std::cos(ang); <<
316 G4double step; <<
317 if(0 == idx) { <<
318 step = 0.5*(angle[1] + angle[0]); <<
319 } else if(DIMMOTT == idx + 1) { <<
320 step = CLHEP::pi - 0.5*(angle[DIMMOTT-2] + <<
321 } else { <<
322 step = 0.5*(angle[idx+1] - angle[idx-1]); <<
323 } <<
324 <<
325 G4double F2 = 1.; <<
326 switch (form) { <<
327 case 1: <<
328 F2 = FormFactor2ExpHof(z1*0.5); <<
329 break; <<
330 case 2: <<
331 F2 = FormFactor2Gauss(z1*0.5); <<
332 break; <<
333 case 3: <<
334 F2 = FormFactor2UniformHelm(z1*0.5); <<
335 break; <<
336 } <<
337 <<
338 const G4double R = RatioMottRutherfordCosT(s <<
339 <<
340 G4double den = 2.*As + z1; <<
341 G4double func = 1./(den*den); <<
342 <<
343 G4double fatt = (G4double)targetZ/(mu_rel*ga <<
344 G4double sigma= e2*e2*fatt*fatt*func; <<
345 G4double pi2sintet = CLHEP::twopi*std::sqrt( <<
346 306
347 G4double Xsec = std::max(pi2sintet*F2*R*sigm << 307 TotalCross=0;
348 return Xsec; << 308
349 } << 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);
350 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 }
351 //....oooOO0OOooo........oooOO0OOooo........oo 392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 393 G4double G4ScreeningMottCrossSection::AngleDistribution(G4double anglein){
352 394
353 G4double <<
354 G4ScreeningMottCrossSection::NuclearCrossSecti <<
355 { <<
356 fTotalCross=0.; <<
357 if(cosTetMaxNuc >= cosTetMinNuc) return fTot <<
358 if(0 == cross.size()) { cross.resize(DIMMOTT <<
359 <<
360 //G4cout<<"MODEL: "<<fast<<G4endl; <<
361 //************ PRECISE COMPUTATION <<
362 if(fast == 0){ <<
363 <<
364 for(G4int i=0; i<DIMMOTT; ++i){ <<
365 G4double y = DifferentialXSection(i, for <<
366 fTotalCross += y; <<
367 cross[i] = fTotalCross; <<
368 if(fTotalCross*1.e-9 > y) { <<
369 for(G4int j=i+1; j<DIMMOTT; ++j) { <<
370 cross[j] = fTotalCross; <<
371 } <<
372 break; <<
373 } <<
374 } <<
375 //**************** FAST COMPUTATION <<
376 } else if(fast == 1) { <<
377 395
378 G4double p0 = electron_mass_c2*classic_ele << 396 G4double total=TotalCross ;
379 G4double coeff = twopi*p0*p0; << 397 G4double fatt= e2*targetZ/(mu_rel*gamma*beta*beta);
>> 398 G4double fatt2=fatt*fatt;
>> 399 total/=fatt2;
380 400
381 G4double fac = coeff*targetZ*targetZ*invbe << 401 G4double R=0;
>> 402 if (coeffb[0][0]!=0){
>> 403 // cout<<" Mott....targetZ "<< targetZ<<endl;
>> 404 R=RatioMottRutherford(anglein);
>> 405 }
382 406
383 G4double x = 1.0 - cosTetMinNuc; << 407 else if (coeffb[0][0]==0){
384 G4double x1 = x + 2*As; << 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;
385 417
386 // scattering with nucleus <<
387 fTotalCross = fac*(cosTetMinNuc - cosTetMa <<
388 (x1*(1.0 - cosTetMaxNuc + 2*As)); <<
389 } <<
390 /* <<
391 G4cout << "Energy(MeV): " << tkinLab/CLHEP:: <<
392 << " Total Cross(mb): " << fTotalCros <<
393 << " form= " << form << " fast= " << <<
394 */ <<
395 return fTotalCross; <<
396 } 418 }
397 419
398 //....oooOO0OOooo........oooOO0OOooo........oo 420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
399 421
400 G4double << 422 G4double G4ScreeningMottCrossSection::GetScatteringAngle()
401 G4ScreeningMottCrossSection::GetScatteringAngl <<
402 { 423 {
403 G4double scattangle = 0.0; << 424
404 G4double r = G4UniformRand(); << 425
405 //************ PRECISE COMPUTATION << 426 //cout<<" G4ScreeningMottCrossSection::SampleCosineTheta ............."<<endl;
406 if(fast == 0){ << 427
407 //G4cout << "r= " << r << " tot= " << fTot << 428 if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
408 r *= fTotalCross; << 429
409 for(G4int i=0; i<DIMMOTT; ++i){ << 430 G4double anglemin=std::acos(cosTetMinNuc);
410 //G4cout << i << ". r= " << r << " xs= " << 431 G4double anglemax= std::acos(cosTetMaxNuc);
411 if(r <= cross[i]) { << 432
412 scattangle = ComputeAngle(i, r); << 433 const G4double limit = 1.e-9;
413 break; << 434 if(anglemin < limit) {
414 } << 435 anglemin = GetScreeningAngle()/10.;
415 } << 436 if(anglemin < limit) { anglemin = limit; }
416 << 437 }
417 //**************** FAST COMPUTATION << 438
418 } else if(fast == 1) { << 439 // cout<<"................ tkinLab "<< G4BestUnit(tkinLab,"Energy") << " anglemin= "<<anglemin<<endl;
419 << 440 //cout<<"anglemax= "<<anglemax<<endl;
420 G4double limit = GetTransitionRandom(); << 441 G4double r =G4UniformRand();
421 if(limit > 0.0) { << 442
422 G4double Sz = 2*As; << 443 G4double loganglemin=log10(anglemin);
423 G4double x = Sz-((Sz*(2+Sz))/(Sz+2*limit << 444 G4double loganglemax=log10(anglemax);
424 //G4cout << "limit= " << limit << " Sz= << 445 G4double logdangle=0.01;
425 G4double angle_limit = (std::abs(x) < 1. << 446
426 //G4cout<<"ANGLE LIMIT: "<<angle_limit<< << 447 G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);
427 << 448
428 if(r > limit && angle_limit != 0.0) { << 449 std::vector<G4double> angle;
429 x = Sz-((Sz*(2+Sz))/(Sz+2*r))+1.0; << 450 std::vector<G4double> tet;
430 scattangle = (x >= 1.0) ? 0.0 : ((x <= -1.0) << 451 std::vector<G4double> dangle;
431 //G4cout<<"FAST: scattangle= "<< scattangle << 452
432 } << 453 for(G4int i=0; i<=bins; i++ ){
433 } else { << 454 tet.push_back(0);
434 //G4cout<<"SLOW: total= "<<fTotalCross<< << 455 dangle.push_back(0);
435 r *= fTotalCross; << 456 angle.push_back(pow(10.,loganglemin+logdangle*i));
436 G4double tot = 0.0; << 457 }
437 for(G4int i=0; i<DIMMOTT; ++i) { << 458
438 G4double xs = DifferentialXSection(i, form); << 459
439 tot += xs; << 460 G4int dim = tet.size();
440 cross[i] = tot; << 461 G4double scattangle=0;
441 if(r <= tot) { << 462 G4double y=0;
442 scattangle = ComputeAngle(i, r); << 463 G4double dy=0;
443 break; << 464 G4double area=0;
444 } << 465
445 } << 466 for(G4int i=0; i<dim;i++){
446 } << 467
447 } << 468 if(i!=dim-1){
448 //****************************************** << 469 dangle[i]=(angle[i+1]-angle[i]);
449 //G4cout<<"Energy(MeV): "<<tkinLab/MeV<<" SC << 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
450 return scattangle; 489 return scattangle;
>> 490
>> 491
451 } 492 }
452 493
453 G4double G4ScreeningMottCrossSection::ComputeA << 494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
454 { << 495
455 G4double x1, x2, y; << 496 G4ThreeVector G4ScreeningMottCrossSection::GetNewDirection(){
456 if(0 == i) { << 497
457 x1 = 0.0; << 498 G4ThreeVector dir(0.0,0.0,1.0);
458 x2 = 0.5*(angle[0] + angle[1]); << 499
459 y = cross[0]; << 500
460 } else if(i == DIMMOTT-1) { << 501 G4double z1=GetScatteringAngle();
461 x1 = 0.5*(angle[i] + angle[i-1]); << 502
462 x2 = CLHEP::pi; << 503 G4double cost = cos(z1);
463 y = cross[i] - cross[i-1]; << 504 G4double sint = sin(z1);
464 r -= cross[i-1]; << 505 G4double phi = twopi* G4UniformRand();
465 } else { << 506 G4double dirx = sint*cos(phi);
466 x1 = 0.5*(angle[i] + angle[i-1]); << 507 G4double diry = sint*sin(phi);
467 x2 = 0.5*(angle[i] + angle[i+1]); << 508 G4double dirz = cost;
468 y = cross[i] - cross[i-1]; << 509
469 r -= cross[i-1]; << 510
470 } << 511 //.......set Trc
471 //G4cout << i << ". r= " << r << " y= " << y << 512 G4double etot=tkinLab+mass;
472 // << " x1= " << " x2= " << x2 << G4endl; << 513 G4double mass2=targetMass;
473 return x1 + (x2 - x1)*r/y; << 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
474 } 526 }
>> 527
>> 528
475 529