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

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
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-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.