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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 16.0 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.9 2009/08/13 11:32:47 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-03-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 // Test if water material
155
156 flagMaterialIsWater= false;
157 densityWater = 0;
158
159 const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
160
161 if(theCoupleTable)
162 {
163 G4int numOfCouples = theCoupleTable->GetTableSize();
164
165 if(numOfCouples>0)
166 {
167 for (G4int i=0; i<numOfCouples; i++)
168 {
169 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
170 const G4Material* material = couple->GetMaterial();
171
172 if (material->GetName() == "G4_WATER")
173 {
174 G4double density = material->GetAtomicNumDensityVector()[1];
175 flagMaterialIsWater = true;
176 densityWater = density;
177
178 if (verboseLevel > 3)
179 G4cout << "****** Water material is found with density(cm^-3)=" << density/(cm*cm*cm) << G4endl;
180 }
181
182 }
183
184 } // if(numOfCouples>0)
185
186 } // if (theCoupleTable)
187
188}
189
190//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
191
192G4double G4DNAScreenedRutherfordElasticModel::CrossSectionPerVolume(const G4Material*,
193 const G4ParticleDefinition*,
194 G4double ekin,
195 G4double,
196 G4double)
197{
198 if (verboseLevel > 3)
199 G4cout << "Calling CrossSectionPerVolume() of G4DNAScreenedRutherfordElasticModel" << G4endl;
200
201 // Calculate total cross section for model
202
203 G4double sigma=0;
204
205 if (flagMaterialIsWater)
206 {
207
208 if (ekin < highEnergyLimit)
209 {
210
211 //SI : XS must not be zero otherwise sampling of secondaries method ignored
212 if (ekin < lowEnergyLimitOfModel) ekin = lowEnergyLimitOfModel;
213 //
214
215 G4double z = 10.;
216 G4double n = ScreeningFactor(ekin,z);
217 G4double crossSection = RutherfordCrossSection(ekin, z);
218 sigma = pi * crossSection / (n * (n + 1.));
219 }
220
221 if (verboseLevel > 3)
222 {
223 G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
224 G4cout << " - Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
225 G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*densityWater/(1./cm) << G4endl;
226 }
227
228 } // if (flagMaterialIsWater)
229
230 return sigma*densityWater;
231}
232
233//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
234
235G4double G4DNAScreenedRutherfordElasticModel::RutherfordCrossSection(G4double k, G4double z)
236{
237 //
238 // e^4 / K + m_e c^2 \^2
239 // sigma_Ruth(K) = Z (Z+1) -------------------- | --------------------- |
240 // (4 pi epsilon_0)^2 \ K * (K + 2 m_e c^2) /
241 //
242 // Where K is the electron non-relativistic kinetic energy
243 //
244 // NIM 155, pp. 145-156, 1978
245
246 G4double length =(e_squared * (k + electron_mass_c2)) / (4 * pi *epsilon0 * k * ( k + 2 * electron_mass_c2));
247 G4double cross = z * ( z + 1) * length * length;
248
249 return cross;
250}
251
252//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
253
254G4double G4DNAScreenedRutherfordElasticModel::ScreeningFactor(G4double k, G4double z)
255{
256 //
257 // alpha_1 + beta_1 ln(K/eV) constK Z^(2/3)
258 // n(T) = -------------------------- -----------------
259 // K/(m_e c^2) 2 + K/(m_e c^2)
260 //
261 // Where K is the electron non-relativistic kinetic energy
262 //
263 // n(T) > 0 for T < ~ 400 MeV
264 //
265 // NIM 155, pp. 145-156, 1978
266 // Formulae (2) and (5)
267
268 const G4double alpha_1(1.64);
269 const G4double beta_1(-0.0825);
270 const G4double constK(1.7E-5);
271
272 G4double numerator = (alpha_1 + beta_1 * std::log(k/eV)) * constK * std::pow(z, 2./3.);
273
274 k /= electron_mass_c2;
275
276 G4double denominator = k * (2 + k);
277
278 G4double value = 0.;
279 if (denominator > 0.) value = numerator / denominator;
280
281 return value;
282}
283
284//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
285
286void G4DNAScreenedRutherfordElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
287 const G4MaterialCutsCouple* /*couple*/,
288 const G4DynamicParticle* aDynamicElectron,
289 G4double,
290 G4double)
291{
292
293 if (verboseLevel > 3)
294 G4cout << "Calling SampleSecondaries() of G4DNAScreenedRutherfordElasticModel" << G4endl;
295
296 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
297
298 if (electronEnergy0 < killBelowEnergy)
299 {
300 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
301 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
302 return ;
303 }
304
305 G4double cosTheta = 0.;
306
307 if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
308 {
309 if (electronEnergy0<intermediateEnergyLimit)
310 {
311 if (verboseLevel > 3) G4cout << "---> Using Brenner & Zaider model" << G4endl;
312 cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
313 }
314
315 if (electronEnergy0>=intermediateEnergyLimit)
316 {
317 if (verboseLevel > 3) G4cout << "---> Using Screened Rutherford model" << G4endl;
318 G4double z = 10.;
319 cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
320 }
321
322 G4double phi = 2. * pi * G4UniformRand();
323
324 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
325 G4ThreeVector xVers = zVers.orthogonal();
326 G4ThreeVector yVers = zVers.cross(xVers);
327
328 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
329 G4double yDir = xDir;
330 xDir *= std::cos(phi);
331 yDir *= std::sin(phi);
332
333 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
334
335 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
336
337 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
338 }
339
340}
341
342//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
343
344G4double G4DNAScreenedRutherfordElasticModel::BrennerZaiderRandomizeCosTheta(G4double k)
345{
346 // d sigma_el 1 beta(K)
347 // ------------ (K) ~ --------------------------------- + ---------------------------------
348 // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2
349 //
350 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
351 //
352 // Phys. Med. Biol. 29 N.4 (1983) 443-447
353
354 // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
355
356 k /= eV;
357
358 G4double beta = std::exp(CalculatePolynomial(k,betaCoeff));
359 G4double delta = std::exp(CalculatePolynomial(k,deltaCoeff));
360 G4double gamma;
361
362 if (k > 100.)
363 {
364 gamma = CalculatePolynomial(k, gamma100_200Coeff);
365 // Only in this case it is not the exponent of the polynomial
366 }
367 else
368 {
369 if (k>10)
370 {
371 gamma = std::exp(CalculatePolynomial(k, gamma10_100Coeff));
372 }
373 else
374 {
375 gamma = std::exp(CalculatePolynomial(k, gamma035_10Coeff));
376 }
377 }
378
379 // ***** Original method
380
381 G4double oneOverMax = 1. / (1./(4.*gamma*gamma) + beta/( (2.+2.*delta)*(2.+2.*delta) ));
382
383 G4double cosTheta = 0.;
384 G4double leftDenominator = 0.;
385 G4double rightDenominator = 0.;
386 G4double fCosTheta = 0.;
387
388 do
389 {
390 cosTheta = 2. * G4UniformRand() - 1.;
391
392 leftDenominator = (1. + 2.*gamma - cosTheta);
393 rightDenominator = (1. + 2.*delta + cosTheta);
394 if ( (leftDenominator * rightDenominator) != 0. )
395 {
396 fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
397 }
398 }
399 while (fCosTheta < G4UniformRand());
400
401 return cosTheta;
402
403 // ***** Alternative method using cumulative probability
404/*
405 G4double cosTheta = -1;
406 G4double cumul = 0;
407 G4double value = 0;
408 G4double leftDenominator = 0.;
409 G4double rightDenominator = 0.;
410
411 // Number of integration steps in the -1,1 range
412 G4int iMax=200;
413
414 G4double random = G4UniformRand();
415
416 // Cumulate differential cross section
417 for (G4int i=0; i<iMax; i++)
418 {
419 cosTheta = -1 + i*2./(iMax-1);
420 leftDenominator = (1. + 2.*gamma - cosTheta);
421 rightDenominator = (1. + 2.*delta + cosTheta);
422 if ( (leftDenominator * rightDenominator) != 0. )
423 {
424 cumul = cumul + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
425 }
426 }
427
428 // Select cosTheta
429 for (G4int i=0; i<iMax; i++)
430 {
431 cosTheta = -1 + i*2./(iMax-1);
432 leftDenominator = (1. + 2.*gamma - cosTheta);
433 rightDenominator = (1. + 2.*delta + cosTheta);
434 if (cumul !=0 && (leftDenominator * rightDenominator) != 0.)
435 value = value + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)) / cumul;
436 if (random < value) break;
437 }
438
439 return cosTheta;
440*/
441
442}
443
444//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
445
446G4double G4DNAScreenedRutherfordElasticModel::CalculatePolynomial(G4double k, std::vector<G4double>& vec)
447{
448 // Sum_{i=0}^{size-1} vector_i k^i
449 //
450 // Phys. Med. Biol. 29 N.4 (1983) 443-447
451
452 G4double result = 0.;
453 size_t size = vec.size();
454
455 while (size>0)
456 {
457 size--;
458
459 result *= k;
460 result += vec[size];
461 }
462
463 return result;
464}
465
466//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
467
468G4double G4DNAScreenedRutherfordElasticModel::ScreenedRutherfordRandomizeCosTheta(G4double k, G4double z)
469{
470
471 // d sigma_el sigma_Ruth(K)
472 // ------------ (K) ~ -----------------------------
473 // d Omega (1 + 2 n(K) - cos(theta))^2
474 //
475 // We extract cos(theta) distributed as (1 + 2 n(K) - cos(theta))^-2
476 //
477 // 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)
478 //
479 // Phys. Med. Biol. 45 (2000) 3171-3194
480
481 // ***** Original method
482
483 G4double n = ScreeningFactor(k, z);
484
485 G4double oneOverMax = (4.*n*n);
486
487 G4double cosTheta = 0.;
488 G4double fCosTheta;
489
490 do
491 {
492 cosTheta = 2. * G4UniformRand() - 1.;
493 fCosTheta = (1 + 2.*n - cosTheta);
494 if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta);
495 }
496 while (fCosTheta < G4UniformRand());
497
498 return cosTheta;
499
500 // ***** Alternative method using cumulative probability
501/*
502 G4double cosTheta = -1;
503 G4double cumul = 0;
504 G4double value = 0;
505 G4double n = ScreeningFactor(k, z);
506 G4double fCosTheta;
507
508 // Number of integration steps in the -1,1 range
509 G4int iMax=200;
510
511 G4double random = G4UniformRand();
512
513 // Cumulate differential cross section
514 for (G4int i=0; i<iMax; i++)
515 {
516 cosTheta = -1 + i*2./(iMax-1);
517 fCosTheta = (1 + 2.*n - cosTheta);
518 if (fCosTheta !=0.) cumul = cumul + 1./(fCosTheta*fCosTheta);
519 }
520
521 // Select cosTheta
522 for (G4int i=0; i<iMax; i++)
523 {
524 cosTheta = -1 + i*2./(iMax-1);
525 fCosTheta = (1 + 2.*n - cosTheta);
526 if (cumul !=0.) value = value + (1./(fCosTheta*fCosTheta)) / cumul;
527 if (random < value) break;
528 }
529 return cosTheta;
530*/
531}
532
Note: See TracBrowser for help on using the repository browser.