| 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 // 26 //
27 27
28 #include "G4DNAScreenedRutherfordElasticModel. 28 #include "G4DNAScreenedRutherfordElasticModel.hh"
29 #include "G4PhysicalConstants.hh" 29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 #include "G4DNAMolecularMaterial.hh" 31 #include "G4DNAMolecularMaterial.hh"
32 #include "G4Exp.hh" 32 #include "G4Exp.hh"
33 #include "G4Log.hh" 33 #include "G4Log.hh"
34 34
35 //....oooOO0OOooo........oooOO0OOooo........oo 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36 36
37 using namespace std; 37 using namespace std;
38 38
39 //....oooOO0OOooo........oooOO0OOooo........oo 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 40
41 G4DNAScreenedRutherfordElasticModel:: 41 G4DNAScreenedRutherfordElasticModel::
42 G4DNAScreenedRutherfordElasticModel(const G4Pa 42 G4DNAScreenedRutherfordElasticModel(const G4ParticleDefinition*,
43 const G4St 43 const G4String& nam) :
44 G4VEmModel(nam) << 44 G4VEmModel(nam), isInitialised(false)
45 { 45 {
46 fpWaterDensity = nullptr; << 46 fpWaterDensity = 0;
47 47
48 lowEnergyLimit = 0 * eV; 48 lowEnergyLimit = 0 * eV;
49 intermediateEnergyLimit = 200 * eV; // Switc 49 intermediateEnergyLimit = 200 * eV; // Switch between two final state models
50 highEnergyLimit = 1. * MeV; 50 highEnergyLimit = 1. * MeV;
51 51
52 SetLowEnergyLimit(lowEnergyLimit); 52 SetLowEnergyLimit(lowEnergyLimit);
53 SetHighEnergyLimit(highEnergyLimit); 53 SetHighEnergyLimit(highEnergyLimit);
54 54
55 verboseLevel = 0; 55 verboseLevel = 0;
56 // Verbosity scale: 56 // Verbosity scale:
57 // 0 = nothing 57 // 0 = nothing
58 // 1 = warning for energy non-conservation 58 // 1 = warning for energy non-conservation
59 // 2 = details of energy budget 59 // 2 = details of energy budget
60 // 3 = calculation of cross sections, file o 60 // 3 = calculation of cross sections, file openings, sampling of atoms
61 // 4 = entering in methods 61 // 4 = entering in methods
62 62
63 #ifdef SR_VERBOSE 63 #ifdef SR_VERBOSE
64 if (verboseLevel > 0) 64 if (verboseLevel > 0)
65 { 65 {
66 G4cout << "Screened Rutherford Elastic mod 66 G4cout << "Screened Rutherford Elastic model is constructed "
67 << G4endl 67 << G4endl
68 << "Energy range: " 68 << "Energy range: "
69 << lowEnergyLimit / eV << " eV - " 69 << lowEnergyLimit / eV << " eV - "
70 << highEnergyLimit / MeV << " MeV" 70 << highEnergyLimit / MeV << " MeV"
71 << G4endl; 71 << G4endl;
72 } 72 }
73 #endif 73 #endif
74 fParticleChangeForGamma = nullptr; << 74 fParticleChangeForGamma = 0;
75 75
76 // Selection of computation method 76 // Selection of computation method
77 // We do not recommend "true" usage with the 77 // We do not recommend "true" usage with the current cumul. proba. settings
78 fasterCode = false; 78 fasterCode = false;
79 } 79 }
80 80
81 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82 82
83 G4DNAScreenedRutherfordElasticModel::~G4DNAScr 83 G4DNAScreenedRutherfordElasticModel::~G4DNAScreenedRutherfordElasticModel()
84 = default; << 84 {
>> 85 }
85 86
86 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87 88
88 void G4DNAScreenedRutherfordElasticModel:: 89 void G4DNAScreenedRutherfordElasticModel::
89 Initialise(const G4ParticleDefinition* particl 90 Initialise(const G4ParticleDefinition* particle,
90 const G4DataVector& /*cuts*/) 91 const G4DataVector& /*cuts*/)
91 { 92 {
92 #ifdef SR_VERBOSE 93 #ifdef SR_VERBOSE
93 if (verboseLevel > 3) 94 if (verboseLevel > 3)
94 { 95 {
95 G4cout << "Calling G4DNAScreenedRutherford 96 G4cout << "Calling G4DNAScreenedRutherfordElasticModel::Initialise()"
96 << G4endl; 97 << G4endl;
97 } 98 }
98 #endif 99 #endif
99 100
100 if(particle->GetParticleName() != "e-") 101 if(particle->GetParticleName() != "e-")
101 { 102 {
102 G4Exception ("*** WARNING: the G4DNAScreen 103 G4Exception ("*** WARNING: the G4DNAScreenedRutherfordElasticModel is not "
103 "intented to be used with ano 104 "intented to be used with another particle than the electron",
104 "",FatalException,"") ; 105 "",FatalException,"") ;
105 } 106 }
106 107
107 // Energy limits 108 // Energy limits
108 if (LowEnergyLimit() < 9*eV) 109 if (LowEnergyLimit() < 9*eV)
109 { 110 {
110 G4Exception("*** WARNING: the G4DNAScreene 111 G4Exception("*** WARNING: the G4DNAScreenedRutherfordElasticModel class is "
111 "not validated below 9 eV", 112 "not validated below 9 eV",
112 "",JustWarning,"") ; 113 "",JustWarning,"") ;
113 } 114 }
114 115
115 if (HighEnergyLimit() > 1*MeV) 116 if (HighEnergyLimit() > 1*MeV)
116 { 117 {
117 G4Exception("*** WARNING: the G4DNAScreene 118 G4Exception("*** WARNING: the G4DNAScreenedRutherfordElasticModel class is "
118 "not validated above 1 MeV", 119 "not validated above 1 MeV",
119 "",JustWarning,"") ; 120 "",JustWarning,"") ;
120 } 121 }
121 122
122 // 123 //
123 #ifdef SR_VERBOSE 124 #ifdef SR_VERBOSE
124 if( verboseLevel>0 ) 125 if( verboseLevel>0 )
125 { 126 {
126 G4cout << "Screened Rutherford elastic mod 127 G4cout << "Screened Rutherford elastic model is initialized " << G4endl
127 << "Energy range: " 128 << "Energy range: "
128 << LowEnergyLimit() / eV << " eV - 129 << LowEnergyLimit() / eV << " eV - "
129 << HighEnergyLimit() / MeV << " MeV 130 << HighEnergyLimit() / MeV << " MeV"
130 << G4endl; 131 << G4endl;
131 } 132 }
132 #endif 133 #endif
133 134
134 if (isInitialised) { return; } // return her 135 if (isInitialised) { return; } // return here, prevent reinit consts + pointer
135 136
136 // Initialize water density pointer 137 // Initialize water density pointer
137 fpWaterDensity = G4DNAMolecularMaterial::Ins 138 fpWaterDensity = G4DNAMolecularMaterial::Instance()->
138 GetNumMolPerVolTableFor(G4Ma 139 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
139 140
140 fParticleChangeForGamma = GetParticleChangeF 141 fParticleChangeForGamma = GetParticleChangeForGamma();
141 isInitialised = true; 142 isInitialised = true;
142 143
143 // Constants for final state by Brenner & Za 144 // Constants for final state by Brenner & Zaider
144 // note: if called after if(isInitialised) n 145 // note: if called after if(isInitialised) no need for clear and resetting
145 // the values at every call 146 // the values at every call
146 147
147 betaCoeff= 148 betaCoeff=
148 { 149 {
149 7.51525, 150 7.51525,
150 -0.41912, 151 -0.41912,
151 7.2017E-3, 152 7.2017E-3,
152 -4.646E-5, 153 -4.646E-5,
153 1.02897E-7}; 154 1.02897E-7};
154 155
155 deltaCoeff= 156 deltaCoeff=
156 { 157 {
157 2.9612, 158 2.9612,
158 -0.26376, 159 -0.26376,
159 4.307E-3, 160 4.307E-3,
160 -2.6895E-5, 161 -2.6895E-5,
161 5.83505E-8}; 162 5.83505E-8};
162 163
163 gamma035_10Coeff = 164 gamma035_10Coeff =
164 { 165 {
165 -1.7013, 166 -1.7013,
166 -1.48284, 167 -1.48284,
167 0.6331, 168 0.6331,
168 -0.10911, 169 -0.10911,
169 8.358E-3, 170 8.358E-3,
170 -2.388E-4}; 171 -2.388E-4};
171 172
172 gamma10_100Coeff = 173 gamma10_100Coeff =
173 { 174 {
174 -3.32517, 175 -3.32517,
175 0.10996, 176 0.10996,
176 -4.5255E-3, 177 -4.5255E-3,
177 5.8372E-5, 178 5.8372E-5,
178 -2.4659E-7}; 179 -2.4659E-7};
179 180
180 gamma100_200Coeff = 181 gamma100_200Coeff =
181 { 182 {
182 2.4775E-2, 183 2.4775E-2,
183 -2.96264E-5, 184 -2.96264E-5,
184 -1.20655E-7}; 185 -1.20655E-7};
185 } 186 }
186 187
187 //....oooOO0OOooo........oooOO0OOooo........oo 188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188 189
189 G4double G4DNAScreenedRutherfordElasticModel:: 190 G4double G4DNAScreenedRutherfordElasticModel::
190 CrossSectionPerVolume(const G4Material* materi 191 CrossSectionPerVolume(const G4Material* material,
191 #ifdef SR_VERBOSE 192 #ifdef SR_VERBOSE
192 const G4ParticleDefiniti 193 const G4ParticleDefinition* particleDefinition,
193 #else 194 #else
194 const G4ParticleDefiniti 195 const G4ParticleDefinition*,
195 #endif 196 #endif
196 G4double ekin, 197 G4double ekin,
197 G4double, 198 G4double,
198 G4double) 199 G4double)
199 { 200 {
200 #ifdef SR_VERBOSE 201 #ifdef SR_VERBOSE
201 if (verboseLevel > 3) 202 if (verboseLevel > 3)
202 { 203 {
203 G4cout << "Calling CrossSectionPerVolume() 204 G4cout << "Calling CrossSectionPerVolume() of "
204 "G4DNAScreenedRutherfordElasticMod 205 "G4DNAScreenedRutherfordElasticModel"
205 << G4endl; 206 << G4endl;
206 } 207 }
207 #endif 208 #endif
208 209
209 // Calculate total cross section for model 210 // Calculate total cross section for model
210 211
211 G4double sigma=0.; 212 G4double sigma=0.;
212 G4double waterDensity = (*fpWaterDensity)[ma 213 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
213 214
214 if(ekin <= HighEnergyLimit() && ekin >= LowE 215 if(ekin <= HighEnergyLimit() && ekin >= LowEnergyLimit())
215 { 216 {
216 G4double z = 10.; 217 G4double z = 10.;
217 G4double n = ScreeningFactor(ekin,z); 218 G4double n = ScreeningFactor(ekin,z);
218 G4double crossSection = RutherfordCrossSec 219 G4double crossSection = RutherfordCrossSection(ekin, z);
219 sigma = pi * crossSection / (n * (n + 1.)) 220 sigma = pi * crossSection / (n * (n + 1.));
220 } 221 }
221 222
222 #ifdef SR_VERBOSE 223 #ifdef SR_VERBOSE
223 if (verboseLevel > 2) 224 if (verboseLevel > 2)
224 { 225 {
225 G4cout << "_______________________________ 226 G4cout << "__________________________________" << G4endl;
226 G4cout << "=== G4DNAScreenedRutherfordElas 227 G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO START"
227 << G4endl; 228 << G4endl;
228 G4cout << "=== Kinetic energy(eV)=" << eki 229 G4cout << "=== Kinetic energy(eV)=" << ekin/eV
229 << " particle : " << particleDefini 230 << " particle : " << particleDefinition->GetParticleName()
230 << G4endl; 231 << G4endl;
231 G4cout << "=== Cross section per water mol 232 G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm
232 << G4endl; 233 << G4endl;
233 G4cout << "=== Cross section per water mol 234 G4cout << "=== Cross section per water molecule (cm^-1)="
234 << sigma*waterDensity/(1./cm) << G4 235 << sigma*waterDensity/(1./cm) << G4endl;
235 G4cout << "=== G4DNAScreenedRutherfordElas 236 G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO END"
236 << G4endl; 237 << G4endl;
237 } 238 }
238 #endif 239 #endif
239 240
240 return sigma*waterDensity; 241 return sigma*waterDensity;
241 } 242 }
242 243
243 //....oooOO0OOooo........oooOO0OOooo........oo 244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
244 245
245 G4double G4DNAScreenedRutherfordElasticModel:: 246 G4double G4DNAScreenedRutherfordElasticModel::RutherfordCrossSection(G4double k,
246 247 G4double z)
247 { 248 {
248 // 249 //
249 // e^4 250 // e^4 / K + m_e c^2 \^2
250 // sigma_Ruth(K) = Z (Z+1) ----------------- 251 // sigma_Ruth(K) = Z (Z+1) -------------------- | --------------------- |
251 // (4 pi epsilon_0) 252 // (4 pi epsilon_0)^2 \ K * (K + 2 m_e c^2) /
252 // 253 //
253 // Where K is the electron non-relativistic 254 // Where K is the electron non-relativistic kinetic energy
254 // 255 //
255 // NIM 155, pp. 145-156, 1978 256 // NIM 155, pp. 145-156, 1978
256 257
257 G4double length = (e_squared * (k + electron 258 G4double length = (e_squared * (k + electron_mass_c2))
258 / (4 * pi * epsilon0 * k * (k + 2 * elec 259 / (4 * pi * epsilon0 * k * (k + 2 * electron_mass_c2));
259 G4double cross = z * (z + 1) * length * leng 260 G4double cross = z * (z + 1) * length * length;
260 261
261 return cross; 262 return cross;
262 } 263 }
263 264
264 //....oooOO0OOooo........oooOO0OOooo........oo 265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
265 266
266 G4double G4DNAScreenedRutherfordElasticModel:: 267 G4double G4DNAScreenedRutherfordElasticModel::ScreeningFactor(G4double k,
267 268 G4double z)
268 { 269 {
269 // 270 //
270 // alpha_1 + beta_1 ln(K/eV) const 271 // alpha_1 + beta_1 ln(K/eV) constK Z^(2/3)
271 // n(T) = -------------------------- ------- 272 // n(T) = -------------------------- -----------------
272 // K/(m_e c^2) 2 + K 273 // K/(m_e c^2) 2 + K/(m_e c^2)
273 // 274 //
274 // Where K is the electron non-relativistic 275 // Where K is the electron non-relativistic kinetic energy
275 // 276 //
276 // n(T) > 0 for T < ~ 400 MeV 277 // n(T) > 0 for T < ~ 400 MeV
277 // 278 //
278 // NIM 155, pp. 145-156, 1978 279 // NIM 155, pp. 145-156, 1978
279 // Formulae (2) and (5) 280 // Formulae (2) and (5)
280 281
281 const G4double alpha_1(1.64); 282 const G4double alpha_1(1.64);
282 const G4double beta_1(-0.0825); 283 const G4double beta_1(-0.0825);
283 const G4double constK(1.7E-5); 284 const G4double constK(1.7E-5);
284 285
285 G4double numerator = (alpha_1 + beta_1 * G4L 286 G4double numerator = (alpha_1 + beta_1 * G4Log(k / eV)) * constK
286 * std::pow(z, 2. / 3.); 287 * std::pow(z, 2. / 3.);
287 288
288 k /= electron_mass_c2; 289 k /= electron_mass_c2;
289 290
290 G4double denominator = k * (2 + k); 291 G4double denominator = k * (2 + k);
291 292
292 G4double value = 0.; 293 G4double value = 0.;
293 if (denominator > 0.) value = numerator / de 294 if (denominator > 0.) value = numerator / denominator;
294 295
295 return value; 296 return value;
296 } 297 }
297 298
298 //....oooOO0OOooo........oooOO0OOooo........oo 299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299 300
300 void G4DNAScreenedRutherfordElasticModel:: 301 void G4DNAScreenedRutherfordElasticModel::
301 SampleSecondaries(std::vector<G4DynamicParticl 302 SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
302 const G4MaterialCutsCouple* 303 const G4MaterialCutsCouple* /*couple*/,
303 const G4DynamicParticle* aDy 304 const G4DynamicParticle* aDynamicElectron,
304 G4double, 305 G4double,
305 G4double) 306 G4double)
306 { 307 {
307 #ifdef SR_VERBOSE 308 #ifdef SR_VERBOSE
308 if (verboseLevel > 3) 309 if (verboseLevel > 3)
309 { 310 {
310 G4cout << "Calling SampleSecondaries() of 311 G4cout << "Calling SampleSecondaries() of "
311 "G4DNAScreenedRutherfordElasticM 312 "G4DNAScreenedRutherfordElasticModel"
312 << G4endl; 313 << G4endl;
313 } 314 }
314 #endif 315 #endif
315 316
316 G4double electronEnergy0 = aDynamicElectron- 317 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
317 G4double cosTheta = 0.; 318 G4double cosTheta = 0.;
318 319
319 if (electronEnergy0<intermediateEnergyLimit) 320 if (electronEnergy0<intermediateEnergyLimit)
320 { 321 {
321 #ifdef SR_VERBOSE 322 #ifdef SR_VERBOSE
322 if (verboseLevel > 3) 323 if (verboseLevel > 3)
323 {G4cout << "---> Using Brenner & Zaider 324 {G4cout << "---> Using Brenner & Zaider model" << G4endl;}
324 #endif 325 #endif
325 cosTheta = BrennerZaiderRandomizeCosTheta( 326 cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
326 } 327 }
327 328
328 if (electronEnergy0>=intermediateEnergyLimit 329 if (electronEnergy0>=intermediateEnergyLimit)
329 { 330 {
330 #ifdef SR_VERBOSE 331 #ifdef SR_VERBOSE
331 if (verboseLevel > 3) 332 if (verboseLevel > 3)
332 {G4cout << "---> Using Screened Rutherfo 333 {G4cout << "---> Using Screened Rutherford model" << G4endl;}
333 #endif 334 #endif
334 G4double z = 10.; 335 G4double z = 10.;
335 cosTheta = ScreenedRutherfordRandomizeCosT 336 cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
336 } 337 }
337 338
338 G4double phi = 2. * pi * G4UniformRand(); 339 G4double phi = 2. * pi * G4UniformRand();
339 340
340 G4ThreeVector zVers = aDynamicElectron->GetM 341 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
341 G4ThreeVector xVers = zVers.orthogonal(); 342 G4ThreeVector xVers = zVers.orthogonal();
342 G4ThreeVector yVers = zVers.cross(xVers); 343 G4ThreeVector yVers = zVers.cross(xVers);
343 344
344 G4double xDir = std::sqrt(1. - cosTheta*cosT 345 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
345 G4double yDir = xDir; 346 G4double yDir = xDir;
346 xDir *= std::cos(phi); 347 xDir *= std::cos(phi);
347 yDir *= std::sin(phi); 348 yDir *= std::sin(phi);
348 349
349 G4ThreeVector zPrimeVers((xDir*xVers + yDir* 350 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
350 351
351 fParticleChangeForGamma->ProposeMomentumDire 352 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
352 353
353 fParticleChangeForGamma->SetProposedKineticE 354 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
354 } 355 }
355 356
356 //....oooOO0OOooo........oooOO0OOooo........oo 357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
357 358
358 G4double G4DNAScreenedRutherfordElasticModel:: 359 G4double G4DNAScreenedRutherfordElasticModel::
359 BrennerZaiderRandomizeCosTheta(G4double k) 360 BrennerZaiderRandomizeCosTheta(G4double k)
360 { 361 {
361 // d sigma_el 1 362 // d sigma_el 1 beta(K)
362 // ------------ (K) ~ ---------------------- 363 // ------------ (K) ~ --------------------------------- + ---------------------------------
363 // d Omega (1 + 2 gamma(K) - cos 364 // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2
364 // 365 //
365 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/( 366 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
366 // 367 //
367 // Phys. Med. Biol. 29 N.4 (1983) 443-447 368 // Phys. Med. Biol. 29 N.4 (1983) 443-447
368 369
369 // gamma(K), beta(K) and delta(K) are polyno 370 // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
370 371
371 k /= eV; 372 k /= eV;
372 373
373 G4double beta = G4Exp(CalculatePolynomial(k, 374 G4double beta = G4Exp(CalculatePolynomial(k, betaCoeff));
374 G4double delta = G4Exp(CalculatePolynomial(k 375 G4double delta = G4Exp(CalculatePolynomial(k, deltaCoeff));
375 G4double gamma; 376 G4double gamma;
376 377
377 if (k > 100.) 378 if (k > 100.)
378 { 379 {
379 gamma = CalculatePolynomial(k, gamma100_20 380 gamma = CalculatePolynomial(k, gamma100_200Coeff);
380 // Only in this case it is not the exponen 381 // Only in this case it is not the exponent of the polynomial
381 } 382 }
382 else 383 else
383 { 384 {
384 if (k > 10) 385 if (k > 10)
385 { 386 {
386 gamma = G4Exp(CalculatePolynomial(k, gam 387 gamma = G4Exp(CalculatePolynomial(k, gamma10_100Coeff));
387 } 388 }
388 else 389 else
389 { 390 {
390 gamma = G4Exp(CalculatePolynomial(k, gam 391 gamma = G4Exp(CalculatePolynomial(k, gamma035_10Coeff));
391 } 392 }
392 } 393 }
393 394
394 // ***** Original method 395 // ***** Original method
395 396
396 if (!fasterCode) 397 if (!fasterCode)
397 { 398 {
398 399
399 G4double oneOverMax = 1. 400 G4double oneOverMax = 1.
400 / (1. / (4. * gamma * gamma) + beta 401 / (1. / (4. * gamma * gamma) + beta
401 / ((2. + 2. * delta) * (2. + 2. * de 402 / ((2. + 2. * delta) * (2. + 2. * delta)));
402 403
403 G4double cosTheta = 0.; 404 G4double cosTheta = 0.;
404 G4double leftDenominator = 0.; 405 G4double leftDenominator = 0.;
405 G4double rightDenominator = 0.; 406 G4double rightDenominator = 0.;
406 G4double fCosTheta = 0.; 407 G4double fCosTheta = 0.;
407 408
408 do 409 do
409 { 410 {
410 cosTheta = 2. * G4UniformRand()- 1.; 411 cosTheta = 2. * G4UniformRand()- 1.;
411 412
412 leftDenominator = (1. + 2.*gamma - cosThe 413 leftDenominator = (1. + 2.*gamma - cosTheta);
413 rightDenominator = (1. + 2.*delta + cosTh 414 rightDenominator = (1. + 2.*delta + cosTheta);
414 if ( (leftDenominator * rightDenominator) 415 if ( (leftDenominator * rightDenominator) != 0. )
415 { 416 {
416 fCosTheta = oneOverMax * (1./(leftDenom 417 fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator)
417 + beta/(right 418 + beta/(rightDenominator*rightDenominator));
418 } 419 }
419 } 420 }
420 while (fCosTheta < G4UniformRand()); 421 while (fCosTheta < G4UniformRand());
421 422
422 return cosTheta; 423 return cosTheta;
423 } 424 }
424 425
425 // ***** Alternative method using cumulative 426 // ***** Alternative method using cumulative probability
426 427
427 if (fasterCode) 428 if (fasterCode)
428 { 429 {
429 430
430 // 431 //
431 // modified by Shogo OKADA @ KEK, JP, 2016. 432 // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
432 // 433 //
433 // An integral of differential cross-sectio 434 // An integral of differential cross-section formula shown above this member function
434 // (integral variable: cos(theta), integral 435 // (integral variable: cos(theta), integral interval: [-1, x]) is as follows:
435 // 436 //
436 // 1.0 + x beta * ( 437 // 1.0 + x beta * (1 + x)
437 // I = --------------------- + ------------ 438 // I = --------------------- + ---------------------- (1)
438 // (a - x) * (a + 1.0) (b + x) * 439 // (a - x) * (a + 1.0) (b + x) * (b - 1.0)
439 // 440 //
440 // where a = 1.0 + 2.0 * gamma(K), b = 1.0 441 // where a = 1.0 + 2.0 * gamma(K), b = 1.0 + 2.0 * delta(K)
441 // 442 //
442 // Then, a cumulative probability (cp) is a 443 // Then, a cumulative probability (cp) is as follows:
443 // 444 //
444 // cp 1.0 + x beta * 445 // cp 1.0 + x beta * (1 + x)
445 // ---- = --------------------- + --------- 446 // ---- = --------------------- + ---------------------- (2)
446 // S (a - x) * (a + 1.0) (b + x) 447 // S (a - x) * (a + 1.0) (b + x) * (b - 1.0)
447 // 448 //
448 // where 1/S is the integral of differnetic 449 // where 1/S is the integral of differnetical cross-section (1) on interval [-1, 1]
449 // 450 //
450 // 1 2.0 2. 451 // 1 2.0 2.0 * beta
451 // --- = ----------------------- + ------- 452 // --- = ----------------------- + ----------------------- (3)
452 // S (a - 1.0) * (a + 1.0) (b + 1 453 // S (a - 1.0) * (a + 1.0) (b + 1.0) * (b - 1.0)
453 // 454 //
454 // x is calculated from the quadratic equat 455 // x is calculated from the quadratic equation derived from (2) and (3):
455 // 456 //
456 // A * x^2 + B * x + C = 0 457 // A * x^2 + B * x + C = 0
457 // 458 //
458 // where A, B, anc C are coefficients of th 459 // where A, B, anc C are coefficients of the equation:
459 // A = S * {(b - 1.0) - beta * (a + 1.0)} 460 // A = S * {(b - 1.0) - beta * (a + 1.0)} + cp * (a + 1.0) * (b - 1.0),
460 // B = S * {(b - 1.0) * (b + 1.0) + beta * 461 // B = S * {(b - 1.0) * (b + 1.0) + beta * (a - 1.0) * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * (a - b)
461 // C = S * {b * (b - 1.0) + beta * a * (a 462 // C = S * {b * (b - 1.0) + beta * a * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * ab
462 // 463 //
463 464
464 // sampling cumulative probability 465 // sampling cumulative probability
465 G4double cp = G4UniformRand(); 466 G4double cp = G4UniformRand();
466 467
467 G4double a = 1.0 + 2.0 * gamma; 468 G4double a = 1.0 + 2.0 * gamma;
468 G4double b = 1.0 + 2.0 * delta; 469 G4double b = 1.0 + 2.0 * delta;
469 G4double a1 = a - 1.0; 470 G4double a1 = a - 1.0;
470 G4double a2 = a + 1.0; 471 G4double a2 = a + 1.0;
471 G4double b1 = b - 1.0; 472 G4double b1 = b - 1.0;
472 G4double b2 = b + 1.0; 473 G4double b2 = b + 1.0;
473 G4double c1 = a - b; 474 G4double c1 = a - b;
474 G4double c2 = a * b; 475 G4double c2 = a * b;
475 476
476 G4double S = 2.0 / (a1 * a2) + 2.0 * beta / 477 G4double S = 2.0 / (a1 * a2) + 2.0 * beta / (b1 * b2); S = 1.0 / S;
477 478
478 // coefficients for the quadratic equation 479 // coefficients for the quadratic equation
479 G4double A = S * (b1 - beta * a2) + cp * a2 480 G4double A = S * (b1 - beta * a2) + cp * a2 * b1;
480 G4double B = S * (b1 * b2 + beta * a1 * a2) 481 G4double B = S * (b1 * b2 + beta * a1 * a2) - cp * a2 * b1 * c1;
481 G4double C = S * (b * b1 + beta * a * a2) - 482 G4double C = S * (b * b1 + beta * a * a2) - cp * a2 * b1 * c2;
482 483
483 // calculate cos(theta) 484 // calculate cos(theta)
484 return (-1.0 * B + std::sqrt(B * B - 4.0 * 485 return (-1.0 * B + std::sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
485 486
486 /* 487 /*
487 G4double cosTheta = -1; 488 G4double cosTheta = -1;
488 G4double cumul = 0; 489 G4double cumul = 0;
489 G4double value = 0; 490 G4double value = 0;
490 G4double leftDenominator = 0.; 491 G4double leftDenominator = 0.;
491 G4double rightDenominator = 0.; 492 G4double rightDenominator = 0.;
492 493
493 // Number of integration steps in the -1,1 494 // Number of integration steps in the -1,1 range
494 G4int iMax=200; 495 G4int iMax=200;
495 496
496 G4double random = G4UniformRand(); 497 G4double random = G4UniformRand();
497 498
498 // Cumulate differential cross section 499 // Cumulate differential cross section
499 for (G4int i=0; i<iMax; i++) 500 for (G4int i=0; i<iMax; i++)
500 { 501 {
501 cosTheta = -1 + i*2./(iMax-1); 502 cosTheta = -1 + i*2./(iMax-1);
502 leftDenominator = (1. + 2.*gamma - cosThet 503 leftDenominator = (1. + 2.*gamma - cosTheta);
503 rightDenominator = (1. + 2.*delta + cosThe 504 rightDenominator = (1. + 2.*delta + cosTheta);
504 if ( (leftDenominator * rightDenominator) 505 if ( (leftDenominator * rightDenominator) != 0. )
505 { 506 {
506 cumul = cumul + (1./(leftDenominator*lef 507 cumul = cumul + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
507 } 508 }
508 } 509 }
509 510
510 // Select cosTheta 511 // Select cosTheta
511 for (G4int i=0; i<iMax; i++) 512 for (G4int i=0; i<iMax; i++)
512 { 513 {
513 cosTheta = -1 + i*2./(iMax-1); 514 cosTheta = -1 + i*2./(iMax-1);
514 leftDenominator = (1. + 2.*gamma - cosThe 515 leftDenominator = (1. + 2.*gamma - cosTheta);
515 rightDenominator = (1. + 2.*delta + cosTh 516 rightDenominator = (1. + 2.*delta + cosTheta);
516 if (cumul !=0 && (leftDenominator * right 517 if (cumul !=0 && (leftDenominator * rightDenominator) != 0.)
517 value = value + (1./(leftDenominator*lef 518 value = value + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)) / cumul;
518 if (random < value) break; 519 if (random < value) break;
519 } 520 }
520 521
521 return cosTheta; 522 return cosTheta;
522 */ 523 */
523 } 524 }
524 525
525 return 0.; 526 return 0.;
526 } 527 }
527 528
528 //....oooOO0OOooo........oooOO0OOooo........oo 529 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
529 530
530 G4double G4DNAScreenedRutherfordElasticModel:: 531 G4double G4DNAScreenedRutherfordElasticModel::
531 CalculatePolynomial(G4double k, 532 CalculatePolynomial(G4double k,
532 std::vector<G4double>& vec 533 std::vector<G4double>& vec)
533 { 534 {
534 // Sum_{i=0}^{size-1} vector_i k^i 535 // Sum_{i=0}^{size-1} vector_i k^i
535 // 536 //
536 // Phys. Med. Biol. 29 N.4 (1983) 443-447 537 // Phys. Med. Biol. 29 N.4 (1983) 443-447
537 538
538 G4double result = 0.; 539 G4double result = 0.;
539 size_t size = vec.size(); 540 size_t size = vec.size();
540 541
541 while (size > 0) 542 while (size > 0)
542 { 543 {
543 size--; 544 size--;
544 545
545 result *= k; 546 result *= k;
546 result += vec[size]; 547 result += vec[size];
547 } 548 }
548 549
549 return result; 550 return result;
550 } 551 }
551 552
552 //....oooOO0OOooo........oooOO0OOooo........oo 553 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553 554
554 G4double G4DNAScreenedRutherfordElasticModel:: 555 G4double G4DNAScreenedRutherfordElasticModel::
555 ScreenedRutherfordRandomizeCosTheta(G4double k 556 ScreenedRutherfordRandomizeCosTheta(G4double k,
556 G4double z 557 G4double z)
557 { 558 {
558 559
559 // d sigma_el sigma_Ruth(K) 560 // d sigma_el sigma_Ruth(K)
560 // ------------ (K) ~ ---------------------- 561 // ------------ (K) ~ -----------------------------
561 // d Omega (1 + 2 n(K) - cos(the 562 // d Omega (1 + 2 n(K) - cos(theta))^2
562 // 563 //
563 // We extract cos(theta) distributed as (1 + 564 // We extract cos(theta) distributed as (1 + 2 n(K) - cos(theta))^-2
564 // 565 //
565 // Maximum is for theta=0: 1/(4 n(K)^2) (Whe 566 // Maximum is for theta=0: 1/(4 n(K)^2) (When n(K) is positive, that is always satisfied within the validity of the process)
566 // 567 //
567 // Phys. Med. Biol. 45 (2000) 3171-3194 568 // Phys. Med. Biol. 45 (2000) 3171-3194
568 569
569 // ***** Original method 570 // ***** Original method
570 571
571 if (!fasterCode) 572 if (!fasterCode)
572 { 573 {
573 574
574 G4double n = ScreeningFactor(k, z); 575 G4double n = ScreeningFactor(k, z);
575 576
576 G4double oneOverMax = (4. * n * n); 577 G4double oneOverMax = (4. * n * n);
577 578
578 G4double cosTheta = 0.; 579 G4double cosTheta = 0.;
579 G4double fCosTheta; 580 G4double fCosTheta;
580 581
581 do 582 do
582 { 583 {
583 cosTheta = 2. * G4UniformRand()- 1.; 584 cosTheta = 2. * G4UniformRand()- 1.;
584 fCosTheta = (1 + 2.*n - cosTheta); 585 fCosTheta = (1 + 2.*n - cosTheta);
585 if (fCosTheta !=0.) fCosTheta = oneOverMa 586 if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta);
586 } 587 }
587 while (fCosTheta < G4UniformRand()); 588 while (fCosTheta < G4UniformRand());
588 589
589 return cosTheta; 590 return cosTheta;
590 } 591 }
591 592
592 // ***** Alternative method using cumulative 593 // ***** Alternative method using cumulative probability
>> 594 else
>> 595 {
593 596
594 <<
595 // 597 //
596 // modified by Shogo OKADA @ KEK, JP, 2016. 598 // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
597 // 599 //
598 // The cumulative probability (cp) is calcu 600 // The cumulative probability (cp) is calculated by integrating
599 // the differential cross-section fomula wi 601 // the differential cross-section fomula with cos(theta):
600 // 602 //
601 // n(K) * (1.0 + cos(theta)) 603 // n(K) * (1.0 + cos(theta))
602 // cp = --------------------------------- 604 // cp = ---------------------------------
603 // 1.0 + 2.0 * n(K) - cos(theta) 605 // 1.0 + 2.0 * n(K) - cos(theta)
604 // 606 //
605 // Then, cos(theta) is as follows: 607 // Then, cos(theta) is as follows:
606 // 608 //
607 // cp * (1.0 + 2.0 * n(K)) - 609 // cp * (1.0 + 2.0 * n(K)) - n(K)
608 // cos(theta) = --------------------------- 610 // cos(theta) = --------------------------------
609 // n(k) + cp 611 // n(k) + cp
610 // 612 //
611 // where, K is kinetic energy, n(K) is scre 613 // where, K is kinetic energy, n(K) is screeing factor, and cp is cumulative probability
612 // 614 //
613 615
614 G4double n = ScreeningFactor(k, z); 616 G4double n = ScreeningFactor(k, z);
615 G4double cp = G4UniformRand(); 617 G4double cp = G4UniformRand();
616 G4double numerator = cp * (1.0 + 2.0 * n) - 618 G4double numerator = cp * (1.0 + 2.0 * n) - n;
617 G4double denominator = n + cp; 619 G4double denominator = n + cp;
618 return numerator / denominator; 620 return numerator / denominator;
619 621
620 /* 622 /*
621 G4double cosTheta = -1; 623 G4double cosTheta = -1;
622 G4double cumul = 0; 624 G4double cumul = 0;
623 G4double value = 0; 625 G4double value = 0;
624 G4double n = ScreeningFactor(k, z); 626 G4double n = ScreeningFactor(k, z);
625 G4double fCosTheta; 627 G4double fCosTheta;
626 628
627 // Number of integration steps in the -1,1 629 // Number of integration steps in the -1,1 range
628 G4int iMax=200; 630 G4int iMax=200;
629 631
630 G4double random = G4UniformRand(); 632 G4double random = G4UniformRand();
631 633
632 // Cumulate differential cross section 634 // Cumulate differential cross section
633 for (G4int i=0; i<iMax; i++) 635 for (G4int i=0; i<iMax; i++)
634 { 636 {
635 cosTheta = -1 + i*2./(iMax-1); 637 cosTheta = -1 + i*2./(iMax-1);
636 fCosTheta = (1 + 2.*n - cosTheta); 638 fCosTheta = (1 + 2.*n - cosTheta);
637 if (fCosTheta !=0.) cumul = cumul + 1./(f 639 if (fCosTheta !=0.) cumul = cumul + 1./(fCosTheta*fCosTheta);
638 } 640 }
639 641
640 // Select cosTheta 642 // Select cosTheta
641 for (G4int i=0; i<iMax; i++) 643 for (G4int i=0; i<iMax; i++)
642 { 644 {
643 cosTheta = -1 + i*2./(iMax-1); 645 cosTheta = -1 + i*2./(iMax-1);
644 fCosTheta = (1 + 2.*n - cosTheta); 646 fCosTheta = (1 + 2.*n - cosTheta);
645 if (cumul !=0.) value = value + (1./(fCos 647 if (cumul !=0.) value = value + (1./(fCosTheta*fCosTheta)) / cumul;
646 if (random < value) break; 648 if (random < value) break;
647 } 649 }
648 return cosTheta; 650 return cosTheta;
649 */ 651 */
650 << 652 }
651 653
652 //return 0.; 654 //return 0.;
653 } 655 }
654 656
655 657
656 658