source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNAScreenedRutherfordElasticModel.cc@ 1316

Last change on this file since 1316 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 15.2 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4DNAScreenedRutherfordElasticModel.cc,v 1.10 2010/01/07 18:10:50 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29
30#include "G4DNAScreenedRutherfordElasticModel.hh"
31
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
33
34using namespace std;
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37
38G4DNAScreenedRutherfordElasticModel::G4DNAScreenedRutherfordElasticModel
39(const G4ParticleDefinition*, const G4String& nam)
40:G4VEmModel(nam),isInitialised(false)
41{
42
43 killBelowEnergy = 8.23*eV; // Minimum e- energy for energy loss by excitation
44 lowEnergyLimit = 0 * eV;
45 lowEnergyLimitOfModel = 7 * eV; // The model lower energy is 7 eV
46 intermediateEnergyLimit = 200 * eV; // Switch between two final state models
47 highEnergyLimit = 10 * MeV;
48 SetLowEnergyLimit(lowEnergyLimit);
49 SetHighEnergyLimit(highEnergyLimit);
50
51 verboseLevel= 0;
52 // Verbosity scale:
53 // 0 = nothing
54 // 1 = warning for energy non-conservation
55 // 2 = details of energy budget
56 // 3 = calculation of cross sections, file openings, sampling of atoms
57 // 4 = entering in methods
58
59 if( verboseLevel>0 )
60 {
61 G4cout << "Screened Rutherford Elastic model is constructed " << G4endl
62 << "Energy range: "
63 << lowEnergyLimit / eV << " eV - "
64 << highEnergyLimit / MeV << " MeV"
65 << G4endl;
66 }
67
68}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72G4DNAScreenedRutherfordElasticModel::~G4DNAScreenedRutherfordElasticModel()
73{}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77void G4DNAScreenedRutherfordElasticModel::Initialise(const G4ParticleDefinition* /*particle*/,
78 const G4DataVector& /*cuts*/)
79{
80
81 if (verboseLevel > 3)
82 G4cout << "Calling G4DNAScreenedRutherfordElasticModel::Initialise()" << G4endl;
83
84 // Energy limits
85
86 if (LowEnergyLimit() < lowEnergyLimit)
87 {
88 G4cout << "G4DNAScreenedRutherfordElasticModel: low energy limit increased from " <<
89 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
90 SetLowEnergyLimit(lowEnergyLimit);
91 }
92
93 if (HighEnergyLimit() > highEnergyLimit)
94 {
95 G4cout << "G4DNAScreenedRutherfordElasticModel: high energy limit decreased from " <<
96 HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
97 SetHighEnergyLimit(highEnergyLimit);
98 }
99
100 // Constants for final stae by Brenner & Zaider
101
102 betaCoeff.push_back(7.51525);
103 betaCoeff.push_back(-0.41912);
104 betaCoeff.push_back(7.2017E-3);
105 betaCoeff.push_back(-4.646E-5);
106 betaCoeff.push_back(1.02897E-7);
107
108 deltaCoeff.push_back(2.9612);
109 deltaCoeff.push_back(-0.26376);
110 deltaCoeff.push_back(4.307E-3);
111 deltaCoeff.push_back(-2.6895E-5);
112 deltaCoeff.push_back(5.83505E-8);
113
114 gamma035_10Coeff.push_back(-1.7013);
115 gamma035_10Coeff.push_back(-1.48284);
116 gamma035_10Coeff.push_back(0.6331);
117 gamma035_10Coeff.push_back(-0.10911);
118 gamma035_10Coeff.push_back(8.358E-3);
119 gamma035_10Coeff.push_back(-2.388E-4);
120
121 gamma10_100Coeff.push_back(-3.32517);
122 gamma10_100Coeff.push_back(0.10996);
123 gamma10_100Coeff.push_back(-4.5255E-3);
124 gamma10_100Coeff.push_back(5.8372E-5);
125 gamma10_100Coeff.push_back(-2.4659E-7);
126
127 gamma100_200Coeff.push_back(2.4775E-2);
128 gamma100_200Coeff.push_back(-2.96264E-5);
129 gamma100_200Coeff.push_back(-1.20655E-7);
130
131 //
132
133 if( verboseLevel>0 )
134 {
135 G4cout << "Screened Rutherford elastic model is initialized " << G4endl
136 << "Energy range: "
137 << LowEnergyLimit() / eV << " eV - "
138 << HighEnergyLimit() / MeV << " MeV"
139 << G4endl;
140 }
141
142 if(!isInitialised)
143 {
144 isInitialised = true;
145
146 if(pParticleChange)
147 fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
148 else
149 fParticleChangeForGamma = new G4ParticleChangeForGamma();
150 }
151
152 // InitialiseElementSelectors(particle,cuts);
153
154}
155
156//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157
158G4double G4DNAScreenedRutherfordElasticModel::CrossSectionPerVolume(const G4Material* material,
159 const G4ParticleDefinition*,
160 G4double ekin,
161 G4double,
162 G4double)
163{
164 if (verboseLevel > 3)
165 G4cout << "Calling CrossSectionPerVolume() of G4DNAScreenedRutherfordElasticModel" << G4endl;
166
167 // Calculate total cross section for model
168
169 G4double sigma=0;
170
171 if (material->GetName() == "G4_WATER")
172 {
173
174 if (ekin < highEnergyLimit)
175 {
176
177 //SI : XS must not be zero otherwise sampling of secondaries method ignored
178 if (ekin < lowEnergyLimitOfModel) ekin = lowEnergyLimitOfModel;
179 //
180
181 G4double z = 10.;
182 G4double n = ScreeningFactor(ekin,z);
183 G4double crossSection = RutherfordCrossSection(ekin, z);
184 sigma = pi * crossSection / (n * (n + 1.));
185 }
186
187 if (verboseLevel > 3)
188 {
189 G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
190 G4cout << " - Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
191 G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
192 }
193
194 }
195
196 return sigma*material->GetAtomicNumDensityVector()[1];
197}
198
199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
200
201G4double G4DNAScreenedRutherfordElasticModel::RutherfordCrossSection(G4double k, G4double z)
202{
203 //
204 // e^4 / K + m_e c^2 \^2
205 // sigma_Ruth(K) = Z (Z+1) -------------------- | --------------------- |
206 // (4 pi epsilon_0)^2 \ K * (K + 2 m_e c^2) /
207 //
208 // Where K is the electron non-relativistic kinetic energy
209 //
210 // NIM 155, pp. 145-156, 1978
211
212 G4double length =(e_squared * (k + electron_mass_c2)) / (4 * pi *epsilon0 * k * ( k + 2 * electron_mass_c2));
213 G4double cross = z * ( z + 1) * length * length;
214
215 return cross;
216}
217
218//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219
220G4double G4DNAScreenedRutherfordElasticModel::ScreeningFactor(G4double k, G4double z)
221{
222 //
223 // alpha_1 + beta_1 ln(K/eV) constK Z^(2/3)
224 // n(T) = -------------------------- -----------------
225 // K/(m_e c^2) 2 + K/(m_e c^2)
226 //
227 // Where K is the electron non-relativistic kinetic energy
228 //
229 // n(T) > 0 for T < ~ 400 MeV
230 //
231 // NIM 155, pp. 145-156, 1978
232 // Formulae (2) and (5)
233
234 const G4double alpha_1(1.64);
235 const G4double beta_1(-0.0825);
236 const G4double constK(1.7E-5);
237
238 G4double numerator = (alpha_1 + beta_1 * std::log(k/eV)) * constK * std::pow(z, 2./3.);
239
240 k /= electron_mass_c2;
241
242 G4double denominator = k * (2 + k);
243
244 G4double value = 0.;
245 if (denominator > 0.) value = numerator / denominator;
246
247 return value;
248}
249
250//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251
252void G4DNAScreenedRutherfordElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
253 const G4MaterialCutsCouple* /*couple*/,
254 const G4DynamicParticle* aDynamicElectron,
255 G4double,
256 G4double)
257{
258
259 if (verboseLevel > 3)
260 G4cout << "Calling SampleSecondaries() of G4DNAScreenedRutherfordElasticModel" << G4endl;
261
262 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
263
264 if (electronEnergy0 < killBelowEnergy)
265 {
266 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
267 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
268 return ;
269 }
270
271 G4double cosTheta = 0.;
272
273 if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
274 {
275 if (electronEnergy0<intermediateEnergyLimit)
276 {
277 if (verboseLevel > 3) G4cout << "---> Using Brenner & Zaider model" << G4endl;
278 cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
279 }
280
281 if (electronEnergy0>=intermediateEnergyLimit)
282 {
283 if (verboseLevel > 3) G4cout << "---> Using Screened Rutherford model" << G4endl;
284 G4double z = 10.;
285 cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
286 }
287
288 G4double phi = 2. * pi * G4UniformRand();
289
290 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
291 G4ThreeVector xVers = zVers.orthogonal();
292 G4ThreeVector yVers = zVers.cross(xVers);
293
294 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
295 G4double yDir = xDir;
296 xDir *= std::cos(phi);
297 yDir *= std::sin(phi);
298
299 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
300
301 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
302
303 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
304 }
305
306}
307
308//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
309
310G4double G4DNAScreenedRutherfordElasticModel::BrennerZaiderRandomizeCosTheta(G4double k)
311{
312 // d sigma_el 1 beta(K)
313 // ------------ (K) ~ --------------------------------- + ---------------------------------
314 // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2
315 //
316 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
317 //
318 // Phys. Med. Biol. 29 N.4 (1983) 443-447
319
320 // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
321
322 k /= eV;
323
324 G4double beta = std::exp(CalculatePolynomial(k,betaCoeff));
325 G4double delta = std::exp(CalculatePolynomial(k,deltaCoeff));
326 G4double gamma;
327
328 if (k > 100.)
329 {
330 gamma = CalculatePolynomial(k, gamma100_200Coeff);
331 // Only in this case it is not the exponent of the polynomial
332 }
333 else
334 {
335 if (k>10)
336 {
337 gamma = std::exp(CalculatePolynomial(k, gamma10_100Coeff));
338 }
339 else
340 {
341 gamma = std::exp(CalculatePolynomial(k, gamma035_10Coeff));
342 }
343 }
344
345 // ***** Original method
346
347 G4double oneOverMax = 1. / (1./(4.*gamma*gamma) + beta/( (2.+2.*delta)*(2.+2.*delta) ));
348
349 G4double cosTheta = 0.;
350 G4double leftDenominator = 0.;
351 G4double rightDenominator = 0.;
352 G4double fCosTheta = 0.;
353
354 do
355 {
356 cosTheta = 2. * G4UniformRand() - 1.;
357
358 leftDenominator = (1. + 2.*gamma - cosTheta);
359 rightDenominator = (1. + 2.*delta + cosTheta);
360 if ( (leftDenominator * rightDenominator) != 0. )
361 {
362 fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
363 }
364 }
365 while (fCosTheta < G4UniformRand());
366
367 return cosTheta;
368
369 // ***** Alternative method using cumulative probability
370/*
371 G4double cosTheta = -1;
372 G4double cumul = 0;
373 G4double value = 0;
374 G4double leftDenominator = 0.;
375 G4double rightDenominator = 0.;
376
377 // Number of integration steps in the -1,1 range
378 G4int iMax=200;
379
380 G4double random = G4UniformRand();
381
382 // Cumulate differential cross section
383 for (G4int i=0; i<iMax; i++)
384 {
385 cosTheta = -1 + i*2./(iMax-1);
386 leftDenominator = (1. + 2.*gamma - cosTheta);
387 rightDenominator = (1. + 2.*delta + cosTheta);
388 if ( (leftDenominator * rightDenominator) != 0. )
389 {
390 cumul = cumul + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
391 }
392 }
393
394 // Select cosTheta
395 for (G4int i=0; i<iMax; i++)
396 {
397 cosTheta = -1 + i*2./(iMax-1);
398 leftDenominator = (1. + 2.*gamma - cosTheta);
399 rightDenominator = (1. + 2.*delta + cosTheta);
400 if (cumul !=0 && (leftDenominator * rightDenominator) != 0.)
401 value = value + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)) / cumul;
402 if (random < value) break;
403 }
404
405 return cosTheta;
406*/
407
408}
409
410//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
411
412G4double G4DNAScreenedRutherfordElasticModel::CalculatePolynomial(G4double k, std::vector<G4double>& vec)
413{
414 // Sum_{i=0}^{size-1} vector_i k^i
415 //
416 // Phys. Med. Biol. 29 N.4 (1983) 443-447
417
418 G4double result = 0.;
419 size_t size = vec.size();
420
421 while (size>0)
422 {
423 size--;
424
425 result *= k;
426 result += vec[size];
427 }
428
429 return result;
430}
431
432//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
433
434G4double G4DNAScreenedRutherfordElasticModel::ScreenedRutherfordRandomizeCosTheta(G4double k, G4double z)
435{
436
437 // d sigma_el sigma_Ruth(K)
438 // ------------ (K) ~ -----------------------------
439 // d Omega (1 + 2 n(K) - cos(theta))^2
440 //
441 // We extract cos(theta) distributed as (1 + 2 n(K) - cos(theta))^-2
442 //
443 // 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)
444 //
445 // Phys. Med. Biol. 45 (2000) 3171-3194
446
447 // ***** Original method
448
449 G4double n = ScreeningFactor(k, z);
450
451 G4double oneOverMax = (4.*n*n);
452
453 G4double cosTheta = 0.;
454 G4double fCosTheta;
455
456 do
457 {
458 cosTheta = 2. * G4UniformRand() - 1.;
459 fCosTheta = (1 + 2.*n - cosTheta);
460 if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta);
461 }
462 while (fCosTheta < G4UniformRand());
463
464 return cosTheta;
465
466 // ***** Alternative method using cumulative probability
467/*
468 G4double cosTheta = -1;
469 G4double cumul = 0;
470 G4double value = 0;
471 G4double n = ScreeningFactor(k, z);
472 G4double fCosTheta;
473
474 // Number of integration steps in the -1,1 range
475 G4int iMax=200;
476
477 G4double random = G4UniformRand();
478
479 // Cumulate differential cross section
480 for (G4int i=0; i<iMax; i++)
481 {
482 cosTheta = -1 + i*2./(iMax-1);
483 fCosTheta = (1 + 2.*n - cosTheta);
484 if (fCosTheta !=0.) cumul = cumul + 1./(fCosTheta*fCosTheta);
485 }
486
487 // Select cosTheta
488 for (G4int i=0; i<iMax; i++)
489 {
490 cosTheta = -1 + i*2./(iMax-1);
491 fCosTheta = (1 + 2.*n - cosTheta);
492 if (cumul !=0.) value = value + (1./(fCosTheta*fCosTheta)) / cumul;
493 if (random < value) break;
494 }
495 return cosTheta;
496*/
497}
498
Note: See TracBrowser for help on using the repository browser.