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