source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNAEmfietzoglouExcitationModel.cc@ 1346

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

tag geant4.9.4 beta 1 + modifs locales
File size: 9.0 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: G4DNAEmfietzoglouExcitationModel.cc,v 1.10 2010/06/08 21:50:00 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-01 $
28//
29
30#include "G4DNAEmfietzoglouExcitationModel.hh"
31
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
33
34using namespace std;
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37
38G4DNAEmfietzoglouExcitationModel::G4DNAEmfietzoglouExcitationModel(const G4ParticleDefinition*,
39 const G4String& nam)
40:G4VEmModel(nam),isInitialised(false)
41{
42
43 lowEnergyLimit = 8.23 * eV;
44 highEnergyLimit = 10 * MeV;
45 SetLowEnergyLimit(lowEnergyLimit);
46 SetHighEnergyLimit(highEnergyLimit);
47
48 nLevels = waterExcitation.NumberOfLevels();
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 > 3)
59
60 if( verboseLevel>0 )
61 {
62 G4cout << "Emfietzoglou Excitation model is constructed " << G4endl
63 << "Energy range: "
64 << lowEnergyLimit / eV << " eV - "
65 << highEnergyLimit / MeV << " MeV"
66 << G4endl;
67 }
68}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72G4DNAEmfietzoglouExcitationModel::~G4DNAEmfietzoglouExcitationModel()
73{}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77void G4DNAEmfietzoglouExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
78 const G4DataVector& /*cuts*/)
79{
80
81 if (verboseLevel > 3)
82 G4cout << "Calling G4DNAEmfietzoglouExcitationModel::Initialise()" << G4endl;
83
84 // Energy limits
85
86 if (LowEnergyLimit() < lowEnergyLimit)
87 {
88 G4cout << "G4DNAEmfietzoglouExcitationModel: 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 << "G4DNAEmfietzoglouExcitationModel: high energy limit decreased from " <<
96 HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
97 SetHighEnergyLimit(highEnergyLimit);
98 }
99
100 //
101 if( verboseLevel>0 )
102 {
103 G4cout << "Emfietzoglou Excitation model is initialized " << G4endl
104 << "Energy range: "
105 << LowEnergyLimit() / eV << " eV - "
106 << HighEnergyLimit() / MeV << " MeV"
107 << G4endl;
108 }
109
110 if(!isInitialised)
111 {
112 isInitialised = true;
113
114 if(pParticleChange)
115 fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
116 else
117 fParticleChangeForGamma = new G4ParticleChangeForGamma();
118 }
119
120 // InitialiseElementSelectors(particle,cuts);
121
122}
123
124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
126G4double 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 if (material->GetName() == "G4_WATER")
140 {
141
142 if (particleDefinition == G4Electron::ElectronDefinition())
143 {
144 if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
145 {
146 sigma = Sum(ekin);
147 }
148 }
149
150 if (verboseLevel > 3)
151 {
152 G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
153 G4cout << " - Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
154 G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
155 }
156
157 }
158
159 return sigma*material->GetAtomicNumDensityVector()[1];
160}
161
162//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
163
164void G4DNAEmfietzoglouExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
165 const G4MaterialCutsCouple* /*couple*/,
166 const G4DynamicParticle* aDynamicElectron,
167 G4double,
168 G4double)
169{
170
171 if (verboseLevel > 3)
172 G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouExcitationModel" << G4endl;
173
174 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
175
176 G4int level = RandomSelect(electronEnergy0);
177
178 G4double excitationEnergy = waterExcitation.ExcitationEnergy(level);
179 G4double newEnergy = electronEnergy0 - excitationEnergy;
180
181 if (electronEnergy0 < highEnergyLimit)
182 {
183 if (newEnergy >= lowEnergyLimit)
184 {
185 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
186 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
187 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
188 }
189
190 else
191 {
192 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
193 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
194 }
195 }
196}
197
198//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
199
200G4double G4DNAEmfietzoglouExcitationModel::PartialCrossSection(G4double t, G4int level)
201{
202 // Aj T
203 // Sigma(T) = ------------- (Bj / T) ln(Cj ---) [1 - Bj / T]^Pj
204 // 2 pi alpha0 R
205 //
206 // Sigma is the macroscopic cross section = N sigma, where N = number of target particles per unit volume
207 // and sigma is the microscopic cross section
208 // T is the incoming electron kinetic energy
209 // alpha0 is the Bohr Radius (Bohr_radius)
210 // Aj, Bj, Cj & Pj are parameters that can be found in Emfietzoglou's papers
211 //
212 // From Phys. Med. Biol. 48 (2003) 2355-2371, D.Emfietzoglou,
213 // Monte Carlo Simulation of the energy loss of low energy electrons in liquid Water
214 //
215 // Scaling for macroscopic cross section: number of water moleculs per unit volume
216 // const G4double sigma0 = (10. / 3.343e22) * cm2;
217
218 const G4double density = 3.34192e+19 * mm3;
219
220 const G4double aj[]={0.0205, 0.0209, 0.0130, 0.0026, 0.0025};
221 const G4double cj[]={4.9801, 3.3850, 2.8095, 1.9242, 3.4624};
222 const G4double pj[]={0.4757, 0.3483, 0.4443, 0.3429, 0.4379};
223 const G4double r = 13.6 * eV;
224
225 G4double sigma = 0.;
226
227 G4double exc = waterExcitation.ExcitationEnergy(level);
228
229 if (t >= exc)
230 {
231 G4double excitationSigma = ( aj[level] / (2.*pi*Bohr_radius))
232 * (exc / t)
233 * std::log(cj[level]*(t/r))
234 * std::pow((1.- (exc/t)), pj[level]);
235 sigma = excitationSigma / density;
236 }
237
238 return sigma;
239}
240
241//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242
243G4int G4DNAEmfietzoglouExcitationModel::RandomSelect(G4double k)
244{
245 G4int i = nLevels;
246 G4double value = 0.;
247 std::deque<double> values;
248
249 while (i > 0)
250 {
251 i--;
252 G4double partial = PartialCrossSection(k,i);
253 values.push_front(partial);
254 value += partial;
255 }
256
257 value *= G4UniformRand();
258
259 i = nLevels;
260
261 while (i > 0)
262 {
263 i--;
264 if (values[i] > value) return i;
265 value -= values[i];
266 }
267
268 return 0;
269}
270
271//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
272
273G4double G4DNAEmfietzoglouExcitationModel::Sum(G4double k)
274{
275 G4double totalCrossSection = 0.;
276
277 for (G4int i=0; i<nLevels; i++)
278 {
279 totalCrossSection += PartialCrossSection(k,i);
280 }
281 return totalCrossSection;
282}
283
Note: See TracBrowser for help on using the repository browser.