Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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
31#include "G4SystemOfUnits.hh"
34
35//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36
37using namespace std;
38
39//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40
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 }
73}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
78{}
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81
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
117
118 if (isInitialised) { return; }
120 isInitialised = true;
121
122}
123
124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
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
171void 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 {
193
194 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
196 level,
197 theIncomingTrack);
198 }
199}
200
201//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
202
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
246G4int 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
276G4double 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
@ eExcitedMolecule
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
G4double PartialCrossSection(G4double energy, G4int level)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4DNAEmfietzoglouExcitationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAEmfietzoglouExcitationModel")
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &= *(new G4DataVector()))
static G4DNAMolecularMaterial * Instance()
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:215
size_t GetIndex() const
Definition: G4Material.hh:261
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:576
const G4Track * GetCurrentTrack() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:585
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:109
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:529
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:522
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:592
void ProposeLocalEnergyDeposit(G4double anEnergyPart)