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

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

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