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