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