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

Last change on this file since 1346 was 1340, checked in by garnier, 15 years ago

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