Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAEmfietzoglouExcitationModel.cc
Go to the documentation of this file.
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$
27 //
28 
29 #include "G4DNAEmfietzoglouExcitationModel.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4DNAChemistryManager.hh"
33 #include "G4DNAMolecularMaterial.hh"
34 
35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36 
37 using namespace std;
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 
41 G4DNAEmfietzoglouExcitationModel::G4DNAEmfietzoglouExcitationModel(const G4ParticleDefinition*,
42  const G4String& nam)
43  :G4VEmModel(nam),isInitialised(false)
44 {
45  // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
46  fpWaterDensity = 0;
47 
48  lowEnergyLimit = 8.23 * eV;
49  highEnergyLimit = 10 * MeV;
50  SetLowEnergyLimit(lowEnergyLimit);
51  SetHighEnergyLimit(highEnergyLimit);
52 
53  nLevels = waterExcitation.NumberOfLevels();
54 
55  verboseLevel= 0;
56  // Verbosity scale:
57  // 0 = nothing
58  // 1 = warning for energy non-conservation
59  // 2 = details of energy budget
60  // 3 = calculation of cross sections, file openings, sampling of atoms
61  // 4 = entering in methods
62 
63  if (verboseLevel > 3)
64  if( verboseLevel>0 )
65  {
66  G4cout << "Emfietzoglou Excitation model is constructed " << G4endl
67  << "Energy range: "
68  << lowEnergyLimit / eV << " eV - "
69  << highEnergyLimit / MeV << " MeV"
70  << G4endl;
71  }
72  fParticleChangeForGamma = 0;
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
77 G4DNAEmfietzoglouExcitationModel::~G4DNAEmfietzoglouExcitationModel()
78 {}
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
82 void G4DNAEmfietzoglouExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
83  const G4DataVector& /*cuts*/)
84 {
85 
86  if (verboseLevel > 3)
87  G4cout << "Calling G4DNAEmfietzoglouExcitationModel::Initialise()" << G4endl;
88 
89  // Energy limits
90 
91  if (LowEnergyLimit() < lowEnergyLimit)
92  {
93  G4cout << "G4DNAEmfietzoglouExcitationModel: low energy limit increased from " <<
94  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
95  SetLowEnergyLimit(lowEnergyLimit);
96  }
97 
98  if (HighEnergyLimit() > highEnergyLimit)
99  {
100  G4cout << "G4DNAEmfietzoglouExcitationModel: high energy limit decreased from " <<
101  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
102  SetHighEnergyLimit(highEnergyLimit);
103  }
104 
105  //
106  if( verboseLevel>0 )
107  {
108  G4cout << "Emfietzoglou Excitation model is initialized " << G4endl
109  << "Energy range: "
110  << LowEnergyLimit() / eV << " eV - "
111  << HighEnergyLimit() / MeV << " MeV"
112  << G4endl;
113  }
114 
115  // Initialize water density pointer
116  fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
117 
118  if (isInitialised) { return; }
119  fParticleChangeForGamma = GetParticleChangeForGamma();
120  isInitialised = true;
121 
122 }
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125 
126 G4double G4DNAEmfietzoglouExcitationModel::CrossSectionPerVolume(const G4Material* material,
127  const G4ParticleDefinition* particleDefinition,
128  G4double ekin,
129  G4double,
130  G4double)
131 {
132  if (verboseLevel > 3)
133  G4cout << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouExcitationModel" << G4endl;
134 
135  // Calculate total cross section for model
136 
137  G4double sigma=0;
138 
139  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
140 
141  if(waterDensity!= 0.0)
142  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
143  {
144 
145  if (particleDefinition == G4Electron::ElectronDefinition())
146  {
147  if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
148  {
149  sigma = Sum(ekin);
150  }
151  }
152 
153  if (verboseLevel > 2)
154  {
155  G4cout << "__________________________________" << G4endl;
156  G4cout << "°°° G4DNAEmfietzoglouExcitationModel - XS INFO START" << G4endl;
157  G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
158  G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
159  G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
160  // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
161  G4cout << "°°° G4DNAEmfietzoglouExcitationModel - XS INFO END" << G4endl;
162  }
163 
164  }
165 
166  return sigma*material->GetAtomicNumDensityVector()[1];
167 }
168 
169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
170 
171 void G4DNAEmfietzoglouExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
172  const G4MaterialCutsCouple* /*couple*/,
173  const G4DynamicParticle* aDynamicElectron,
174  G4double,
175  G4double)
176 {
177 
178  if (verboseLevel > 3)
179  G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouExcitationModel" << G4endl;
180 
181  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
182 
183  G4int level = RandomSelect(electronEnergy0);
184 
185  G4double excitationEnergy = waterExcitation.ExcitationEnergy(level);
186  G4double newEnergy = electronEnergy0 - excitationEnergy;
187 
188  if (electronEnergy0 < highEnergyLimit)
189  {
190  fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
191  fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
192  fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
193 
194  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
195  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
196  level,
197  theIncomingTrack);
198  }
199 }
200 
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
202 
203 G4double G4DNAEmfietzoglouExcitationModel::PartialCrossSection(G4double t, G4int level)
204 {
205  // Aj T
206  // Sigma(T) = ------------- (Bj / T) ln(Cj ---) [1 - Bj / T]^Pj
207  // 2 pi alpha0 R
208  //
209  // Sigma is the macroscopic cross section = N sigma, where N = number of target particles per unit volume
210  // and sigma is the microscopic cross section
211  // T is the incoming electron kinetic energy
212  // alpha0 is the Bohr Radius (Bohr_radius)
213  // Aj, Bj, Cj & Pj are parameters that can be found in Emfietzoglou's papers
214  //
215  // From Phys. Med. Biol. 48 (2003) 2355-2371, D.Emfietzoglou,
216  // Monte Carlo Simulation of the energy loss of low energy electrons in liquid Water
217  //
218  // Scaling for macroscopic cross section: number of water moleculs per unit volume
219  // const G4double sigma0 = (10. / 3.343e22) * cm2;
220 
221  const G4double density = 3.34192e+19 * mm3;
222 
223  const G4double aj[]={0.0205, 0.0209, 0.0130, 0.0026, 0.0025};
224  const G4double cj[]={4.9801, 3.3850, 2.8095, 1.9242, 3.4624};
225  const G4double pj[]={0.4757, 0.3483, 0.4443, 0.3429, 0.4379};
226  const G4double r = 13.6 * eV;
227 
228  G4double sigma = 0.;
229 
230  G4double exc = waterExcitation.ExcitationEnergy(level);
231 
232  if (t >= exc)
233  {
234  G4double excitationSigma = ( aj[level] / (2.*pi*Bohr_radius))
235  * (exc / t)
236  * std::log(cj[level]*(t/r))
237  * std::pow((1.- (exc/t)), pj[level]);
238  sigma = excitationSigma / density;
239  }
240 
241  return sigma;
242 }
243 
244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
245 
246 G4int G4DNAEmfietzoglouExcitationModel::RandomSelect(G4double k)
247 {
248  G4int i = nLevels;
249  G4double value = 0.;
250  std::deque<double> values;
251 
252  while (i > 0)
253  {
254  i--;
255  G4double partial = PartialCrossSection(k,i);
256  values.push_front(partial);
257  value += partial;
258  }
259 
260  value *= G4UniformRand();
261 
262  i = nLevels;
263 
264  while (i > 0)
265  {
266  i--;
267  if (values[i] > value) return i;
268  value -= values[i];
269  }
270 
271  return 0;
272 }
273 
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
275 
276 G4double G4DNAEmfietzoglouExcitationModel::Sum(G4double k)
277 {
278  G4double totalCrossSection = 0.;
279 
280  for (G4int i=0; i<nLevels; i++)
281  {
282  totalCrossSection += PartialCrossSection(k,i);
283  }
284  return totalCrossSection;
285 }
286 
Generated on Sat May 25 2013 14:33:30 for Geant4 by   doxygen 1.8.4