| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Based on the work described in << 26 // $Id$
27 // Rad Res 163, 98-111 (2005) <<
28 // D. Emfietzoglou, H. Nikjoo <<
29 // <<
30 // Authors of the class (2014): <<
31 // I. Kyriakou (kyriak@cc.uoi.gr) <<
32 // D. Emfietzoglou (demfietz@cc.uoi.gr) <<
33 // S. Incerti (incerti@cenbg.in2p3.fr) <<
34 // 27 //
35 28
36 #include "G4DNAEmfietzoglouExcitationModel.hh" 29 #include "G4DNAEmfietzoglouExcitationModel.hh"
>> 30 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh"
38 #include "G4DNAChemistryManager.hh" 32 #include "G4DNAChemistryManager.hh"
39 #include "G4DNAMolecularMaterial.hh" 33 #include "G4DNAMolecularMaterial.hh"
40 34
41 //....oooOO0OOooo........oooOO0OOooo........oo 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 36
43 using namespace std; 37 using namespace std;
44 38
45 //....oooOO0OOooo........oooOO0OOooo........oo 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 40
47 G4DNAEmfietzoglouExcitationModel::G4DNAEmfietz 41 G4DNAEmfietzoglouExcitationModel::G4DNAEmfietzoglouExcitationModel(const G4ParticleDefinition*,
48 << 42 const G4String& nam)
49 :G4VEmModel(nam) << 43 :G4VEmModel(nam),isInitialised(false)
50 { 44 {
51 fpMolWaterDensity = nullptr; << 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();
52 54
53 verboseLevel= 0; 55 verboseLevel= 0;
54 // Verbosity scale: 56 // Verbosity scale:
55 // 0 = nothing 57 // 0 = nothing
56 // 1 = warning for energy non-conservation 58 // 1 = warning for energy non-conservation
57 // 2 = details of energy budget 59 // 2 = details of energy budget
58 // 3 = calculation of cross sections, file 60 // 3 = calculation of cross sections, file openings, sampling of atoms
59 // 4 = entering in methods 61 // 4 = entering in methods
60 62
61 if( verboseLevel>0 ) << 63 if (verboseLevel > 3)
62 { << 64 if( verboseLevel>0 )
63 G4cout << "Emfietzoglou excitation model << 65 {
64 } << 66 G4cout << "Emfietzoglou Excitation model is constructed " << G4endl
65 fParticleChangeForGamma = nullptr; << 67 << "Energy range: "
66 << 68 << lowEnergyLimit / eV << " eV - "
67 SetLowEnergyLimit(8.*eV); << 69 << highEnergyLimit / MeV << " MeV"
68 SetHighEnergyLimit(10.*keV); << 70 << G4endl;
69 << 71 }
70 // Selection of stationary mode << 72 fParticleChangeForGamma = 0;
71 statCode = false; <<
72 } 73 }
73 74
74 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 76
76 G4DNAEmfietzoglouExcitationModel::~G4DNAEmfiet 77 G4DNAEmfietzoglouExcitationModel::~G4DNAEmfietzoglouExcitationModel()
77 { << 78 {}
78 // Cross section <<
79 <<
80 std::map< G4String,G4DNACrossSectionDataSe <<
81 for (pos = tableData.begin(); pos != table <<
82 { <<
83 G4DNACrossSectionDataSet* table = pos- <<
84 delete table; <<
85 } <<
86 <<
87 } <<
88 79
89 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 81
91 void G4DNAEmfietzoglouExcitationModel::Initial << 82 void G4DNAEmfietzoglouExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
92 cons << 83 const G4DataVector& /*cuts*/)
93 { 84 {
94 85
95 if (verboseLevel > 3) 86 if (verboseLevel > 3)
96 G4cout << "Calling G4DNAEmfietzoglouEx 87 G4cout << "Calling G4DNAEmfietzoglouExcitationModel::Initialise()" << G4endl;
97 88
98 G4String fileElectron("dna/sigma_excitatio << 89 // Energy limits
99 <<
100 G4ParticleDefinition* electronDef = G4Elec <<
101 <<
102 G4String electron; <<
103 90
104 G4double scaleFactor = (1.e-22 / 3.343) * << 91 if (LowEnergyLimit() < lowEnergyLimit)
105 << 92 {
106 // *** ELECTRON << 93 G4cout << "G4DNAEmfietzoglouExcitationModel: low energy limit increased from " <<
107 << 94 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
108 electron = electronDef->GetParticleName(); << 95 SetLowEnergyLimit(lowEnergyLimit);
109 << 96 }
110 tableFile[electron] = fileElectron; <<
111 <<
112 // Cross section <<
113 <<
114 auto tableE = new G4DNACrossSectionDataSe <<
115 tableE->LoadData(fileElectron); <<
116 97
117 tableData[electron] = tableE; << 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 }
118 104
119 // 105 //
120 <<
121 if( verboseLevel>0 ) 106 if( verboseLevel>0 )
122 { 107 {
123 G4cout << "Emfietzoglou excitation model << 108 G4cout << "Emfietzoglou Excitation model is initialized " << G4endl
124 << "Energy range: " << 109 << "Energy range: "
125 << LowEnergyLimit() / eV << " eV << 110 << LowEnergyLimit() / eV << " eV - "
126 << HighEnergyLimit() / keV << " k << 111 << HighEnergyLimit() / MeV << " MeV"
127 << particle->GetParticleName() << 112 << G4endl;
128 << G4endl; <<
129 } 113 }
130 114
131 // Initialize water density pointer 115 // Initialize water density pointer
132 fpMolWaterDensity = G4DNAMolecularMaterial << 116 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
133 117
134 if (isInitialised) return; << 118 if (isInitialised) { return; }
135 fParticleChangeForGamma = GetParticleChang 119 fParticleChangeForGamma = GetParticleChangeForGamma();
136 isInitialised = true; 120 isInitialised = true;
>> 121
137 } 122 }
138 123
139 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140 125
141 G4double G4DNAEmfietzoglouExcitationModel::Cro 126 G4double G4DNAEmfietzoglouExcitationModel::CrossSectionPerVolume(const G4Material* material,
142 << 127 const G4ParticleDefinition* particleDefinition,
143 << 128 G4double ekin,
144 << 129 G4double,
145 << 130 G4double)
146 { 131 {
147 if (verboseLevel > 3) 132 if (verboseLevel > 3)
148 G4cout << "Calling CrossSectionPerVolu 133 G4cout << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouExcitationModel" << G4endl;
149 134
150 if (particleDefinition != G4Electron::Elec <<
151 <<
152 // Calculate total cross section for model 135 // Calculate total cross section for model
153 136
154 G4double sigma=0; 137 G4double sigma=0;
155 138
156 G4double waterDensity = (*fpMolWaterDensit << 139 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
157 140
158 const G4String& particleName = particleDef << 141 if(waterDensity!= 0.0)
159 << 142 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
160 if (ekin >= LowEnergyLimit() && ekin <= Hi <<
161 { 143 {
162 std::map< G4String,G4DNACrossSectionData <<
163 pos = tableData.find(particleName); <<
164 144
165 if (pos != tableData.end()) << 145 if (particleDefinition == G4Electron::ElectronDefinition())
166 { << 146 {
167 G4DNACrossSectionDataSet* table = pos- << 147 if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
168 if (table != nullptr) sigma = table->F << 148 {
169 } << 149 sigma = Sum(ekin);
170 else << 150 }
171 { << 151 }
172 G4Exception("G4DNAEmfietzoglouExcitati << 152
173 FatalException,"Mo << 153 if (verboseLevel > 2)
174 } << 154 {
175 } << 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 }
176 163
177 if (verboseLevel > 2) <<
178 { <<
179 G4cout << "_____________________________ <<
180 G4cout << "G4DNAEmfietzoglouExcitationMo <<
181 G4cout << "Kinetic energy(eV)=" << ekin/ <<
182 G4cout << "Cross section per water molec <<
183 G4cout << "Cross section per water molec <<
184 //G4cout << " Cross section per water <<
185 ///sigma*material->GetAtomicNumDensityVe <<
186 G4cout << "G4DNAEmfietzoglouExcitationMo <<
187 } 164 }
188 165
189 return sigma*waterDensity; << 166 return sigma*material->GetAtomicNumDensityVector()[1];
190 } 167 }
191 168
192 //....oooOO0OOooo........oooOO0OOooo........oo 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
193 170
194 void G4DNAEmfietzoglouExcitationModel::SampleS 171 void G4DNAEmfietzoglouExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
195 << 172 const G4MaterialCutsCouple* /*couple*/,
196 << 173 const G4DynamicParticle* aDynamicElectron,
197 << 174 G4double,
198 << 175 G4double)
199 { 176 {
200 177
201 if (verboseLevel > 3) 178 if (verboseLevel > 3)
202 G4cout << "Calling SampleSecondaries() 179 G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouExcitationModel" << G4endl;
203 180
204 G4double k = aDynamicParticle->GetKineticE << 181 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
205 182
206 const G4String& particleName = aDynamicPar << 183 G4int level = RandomSelect(electronEnergy0);
207 184
208 G4int level = RandomSelect(k,particleName) << 185 G4double excitationEnergy = waterExcitation.ExcitationEnergy(level);
209 G4double excitationEnergy = waterStructure << 186 G4double newEnergy = electronEnergy0 - excitationEnergy;
210 G4double newEnergy = k - excitationEnergy; <<
211 187
212 if (newEnergy > 0) << 188 if (electronEnergy0 < highEnergyLimit)
213 { 189 {
214 fParticleChangeForGamma->ProposeMoment << 190 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
215 << 191 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
216 if (!statCode) fParticleChangeForGamma << 192 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
217 else fParticleChangeForGamma->SetPropo << 193
218 << 194 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
219 fParticleChangeForGamma->ProposeLocalE << 195 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
>> 196 level,
>> 197 theIncomingTrack);
220 } 198 }
221 <<
222 const G4Track * theIncomingTrack = fPartic <<
223 G4DNAChemistryManager::Instance()->CreateW <<
224 <<
225 <<
226 } 199 }
227 200
228 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
229 202
230 G4int G4DNAEmfietzoglouExcitationModel::Random << 203 G4double G4DNAEmfietzoglouExcitationModel::PartialCrossSection(G4double t, G4int level)
231 { 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;
232 220
233 G4int level = 0; << 221 const G4double density = 3.34192e+19 * mm3;
234 222
235 std::map< G4String,G4DNACrossSectionDataSe << 223 const G4double aj[]={0.0205, 0.0209, 0.0130, 0.0026, 0.0025};
236 pos = tableData.find(particle); << 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;
237 227
238 if (pos != tableData.end()) << 228 G4double sigma = 0.;
239 { <<
240 G4DNACrossSectionDataSet* table = pos- <<
241 229
242 if (table != nullptr) << 230 G4double exc = waterExcitation.ExcitationEnergy(level);
243 { <<
244 auto valuesBuffer = new G4double[ <<
245 const auto n = (G4int)table->Numb <<
246 G4int i(n); <<
247 G4double value = 0.; <<
248 <<
249 //Check reading of initial xs file <<
250 //G4cout << table->GetComponent(0)-> <<
251 //G4cout << table->GetComponent(1) <<
252 //G4cout << table->GetComponent(2) <<
253 //G4cout << table->GetComponent(3) <<
254 //G4cout << table->GetComponent(4) <<
255 //G4cout << table->GetComponent(5) <<
256 //G4cout << table->GetComponent(6) <<
257 //abort(); <<
258 231
259 while (i>0) << 232 if (t >= exc)
260 { << 233 {
261 i--; << 234 G4double excitationSigma = ( aj[level] / (2.*pi*Bohr_radius))
262 valuesBuffer[i] = table->GetCo << 235 * (exc / t)
263 value += valuesBuffer[i]; << 236 * std::log(cj[level]*(t/r))
264 } << 237 * std::pow((1.- (exc/t)), pj[level]);
>> 238 sigma = excitationSigma / density;
>> 239 }
265 240
266 value *= G4UniformRand(); << 241 return sigma;
>> 242 }
267 243
268 i = n; << 244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
269 245
270 while (i > 0) << 246 G4int G4DNAEmfietzoglouExcitationModel::RandomSelect(G4double k)
271 { << 247 {
272 i--; << 248 G4int i = nLevels;
>> 249 G4double value = 0.;
>> 250 std::deque<double> values;
273 251
274 if (valuesBuffer[i] > value) << 252 while (i > 0)
275 { << 253 {
276 delete[] valuesBuffer; << 254 i--;
277 return i; << 255 G4double partial = PartialCrossSection(k,i);
278 } << 256 values.push_front(partial);
279 value -= valuesBuffer[i]; << 257 value += partial;
280 } << 258 }
281 259
282 delete[] valuesBuffer; << 260 value *= G4UniformRand();
>> 261
>> 262 i = nLevels;
283 263
284 } << 264 while (i > 0)
>> 265 {
>> 266 i--;
>> 267 if (values[i] > value) return i;
>> 268 value -= values[i];
285 } 269 }
286 else << 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++)
287 { 281 {
288 G4Exception("G4DNAEmfietzoglouExcitati << 282 totalCrossSection += PartialCrossSection(k,i);
289 FatalException,"Model not <<
290 } 283 }
291 return level; << 284 return totalCrossSection;
292 } 285 }
>> 286
293 287