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

Last change on this file since 1129 was 1058, checked in by garnier, 15 years ago

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