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

Last change on this file since 1005 was 968, checked in by garnier, 17 years ago

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