| 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 // Based on the work described in
27 // Rad Res 163, 98-111 (2005) 27 // Rad Res 163, 98-111 (2005)
28 // D. Emfietzoglou, H. Nikjoo 28 // D. Emfietzoglou, H. Nikjoo
29 // 29 //
30 // Authors of the class (2014): 30 // Authors of the class (2014):
31 // I. Kyriakou (kyriak@cc.uoi.gr) 31 // I. Kyriakou (kyriak@cc.uoi.gr)
32 // D. Emfietzoglou (demfietz@cc.uoi.gr) 32 // D. Emfietzoglou (demfietz@cc.uoi.gr)
33 // S. Incerti (incerti@cenbg.in2p3.fr) 33 // S. Incerti (incerti@cenbg.in2p3.fr)
34 // 34 //
35 35
36 #include "G4DNAEmfietzoglouIonisationModel.hh" 36 #include "G4DNAEmfietzoglouIonisationModel.hh"
37 #include "G4PhysicalConstants.hh" 37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh" 38 #include "G4SystemOfUnits.hh"
39 #include "G4UAtomicDeexcitation.hh" 39 #include "G4UAtomicDeexcitation.hh"
40 #include "G4LossTableManager.hh" 40 #include "G4LossTableManager.hh"
41 #include "G4DNAChemistryManager.hh" 41 #include "G4DNAChemistryManager.hh"
42 #include "G4DNAMolecularMaterial.hh" 42 #include "G4DNAMolecularMaterial.hh"
43 #include "G4DNABornAngle.hh" 43 #include "G4DNABornAngle.hh"
44 #include "G4DeltaAngle.hh" 44 #include "G4DeltaAngle.hh"
45 45
46 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 47
48 using namespace std; 48 using namespace std;
49 49
50 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 51
52 G4DNAEmfietzoglouIonisationModel::G4DNAEmfietz 52 G4DNAEmfietzoglouIonisationModel::G4DNAEmfietzoglouIonisationModel(const G4ParticleDefinition*,
53 53 const G4String& nam) :
54 G4VEmModel(nam) 54 G4VEmModel(nam)
55 { 55 {
56 verboseLevel = 0; 56 verboseLevel = 0;
57 // Verbosity scale: 57 // Verbosity scale:
58 // 0 = nothing 58 // 0 = nothing
59 // 1 = warning for energy non-conservation 59 // 1 = warning for energy non-conservation
60 // 2 = details of energy budget 60 // 2 = details of energy budget
61 // 3 = calculation of cross sections, file o 61 // 3 = calculation of cross sections, file openings, sampling of atoms
62 // 4 = entering in methods 62 // 4 = entering in methods
63 63
64 if(verboseLevel > 0) 64 if(verboseLevel > 0)
65 { 65 {
66 G4cout << "Emfietzoglou ionisation model i 66 G4cout << "Emfietzoglou ionisation model is constructed " << G4endl;
67 } 67 }
68 68
69 // Mark this model as "applicable" for atomi 69 // Mark this model as "applicable" for atomic deexcitation
70 SetDeexcitationFlag(true); 70 SetDeexcitationFlag(true);
71 fAtomDeexcitation = nullptr; 71 fAtomDeexcitation = nullptr;
72 fParticleChangeForGamma = nullptr; 72 fParticleChangeForGamma = nullptr;
73 fpMolWaterDensity = nullptr; 73 fpMolWaterDensity = nullptr;
74 74
75 // Define default angular generator 75 // Define default angular generator
76 SetAngularDistribution(new G4DNABornAngle()) 76 SetAngularDistribution(new G4DNABornAngle());
77 77
78 SetLowEnergyLimit(10. * eV); 78 SetLowEnergyLimit(10. * eV);
79 SetHighEnergyLimit(10. * keV); 79 SetHighEnergyLimit(10. * keV);
80 80
81 // Selection of computation method 81 // Selection of computation method
82 82
83 fasterCode = false; 83 fasterCode = false;
84 84
85 // Selection of stationary mode 85 // Selection of stationary mode
86 86
87 statCode = false; 87 statCode = false;
88 } 88 }
89 89
90 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 91
92 G4DNAEmfietzoglouIonisationModel::~G4DNAEmfiet 92 G4DNAEmfietzoglouIonisationModel::~G4DNAEmfietzoglouIonisationModel()
93 { 93 {
94 // Cross section 94 // Cross section
95 95
96 std::map<G4String, G4DNACrossSectionDataSet* 96 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
97 for(pos = tableData.begin(); pos != tableDat 97 for(pos = tableData.begin(); pos != tableData.end(); ++pos)
98 { 98 {
99 G4DNACrossSectionDataSet* table = pos->sec 99 G4DNACrossSectionDataSet* table = pos->second;
100 delete table; 100 delete table;
101 } 101 }
102 102
103 // Final state 103 // Final state
104 104
105 eVecm.clear(); 105 eVecm.clear();
106 106
107 } 107 }
108 108
109 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 110
111 void G4DNAEmfietzoglouIonisationModel::Initial 111 void G4DNAEmfietzoglouIonisationModel::Initialise(const G4ParticleDefinition* particle,
112 112 const G4DataVector& /*cuts*/)
113 { 113 {
114 114
115 if(verboseLevel > 3) 115 if(verboseLevel > 3)
116 { 116 {
117 G4cout << "Calling G4DNAEmfietzoglouIonisa 117 G4cout << "Calling G4DNAEmfietzoglouIonisationModel::Initialise()" << G4endl;
118 } 118 }
119 119
120 // Energy limits 120 // Energy limits
121 121
122 G4String fileElectron("dna/sigma_ionisation_ 122 G4String fileElectron("dna/sigma_ionisation_e_emfietzoglou");
123 123
124 G4ParticleDefinition* electronDef = G4Electr 124 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
125 125
126 G4String electron; 126 G4String electron;
127 127
128 G4double scaleFactor = (1.e-22 / 3.343) * m* 128 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
129 129
130 const char *path = G4FindDataDir("G4LEDATA") 130 const char *path = G4FindDataDir("G4LEDATA");
131 131
132 // *** ELECTRON 132 // *** ELECTRON
133 133
134 electron = electronDef->GetParticleName(); 134 electron = electronDef->GetParticleName();
135 135
136 tableFile[electron] = fileElectron; 136 tableFile[electron] = fileElectron;
137 137
138 // Cross section 138 // Cross section
139 139
140 auto tableE = new G4DNACrossSectionDataSet( 140 auto tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
141 tableE->LoadData(fileElectron); 141 tableE->LoadData(fileElectron);
142 142
143 tableData[electron] = tableE; 143 tableData[electron] = tableE;
144 144
145 // Final state 145 // Final state
146 146
147 std::ostringstream eFullFileName; 147 std::ostringstream eFullFileName;
148 148
149 if (fasterCode) eFullFileName << path << "/d 149 if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_emfietzoglou.dat";
150 if (!fasterCode) eFullFileName << path << "/ 150 if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_emfietzoglou.dat";
151 151
152 std::ifstream eDiffCrossSection(eFullFileNam 152 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
153 153
154 if (!eDiffCrossSection) 154 if (!eDiffCrossSection)
155 { 155 {
156 if (fasterCode) G4Exception("G4DNAEmfietzo 156 if (fasterCode) G4Exception("G4DNAEmfietzoglouIonisationModel::Initialise","em0003",
157 FatalException,"Missing data file:/dna 157 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_emfietzoglou.dat");
158 158
159 if (!fasterCode) G4Exception("G4DNAEmfietz 159 if (!fasterCode) G4Exception("G4DNAEmfietzoglouIonisationModel::Initialise","em0003",
160 FatalException,"Missing data file:/dna 160 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_emfietzoglou.dat");
161 } 161 }
162 162
163 // 163 //
164 164
165 // Clear the arrays for re-initialization ca 165 // Clear the arrays for re-initialization case (MT mode)
166 // March 25th, 2014 - Vaclav Stepan, Sebasti 166 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
167 167
168 eTdummyVec.clear(); 168 eTdummyVec.clear();
169 169
170 eVecm.clear(); 170 eVecm.clear();
171 171
172 eProbaShellMap->clear(); 172 eProbaShellMap->clear();
173 173
174 eDiffCrossSectionData->clear(); 174 eDiffCrossSectionData->clear();
175 175
176 eNrjTransfData->clear(); 176 eNrjTransfData->clear();
177 177
178 // 178 //
179 179
180 eTdummyVec.push_back(0.); 180 eTdummyVec.push_back(0.);
181 while(!eDiffCrossSection.eof()) 181 while(!eDiffCrossSection.eof())
182 { 182 {
183 G4double tDummy; 183 G4double tDummy;
184 G4double eDummy; 184 G4double eDummy;
185 eDiffCrossSection>>tDummy>>eDummy; 185 eDiffCrossSection>>tDummy>>eDummy;
186 if (tDummy != eTdummyVec.back()) eTdummyVe 186 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
187 for (G4int j=0; j<5; j++) 187 for (G4int j=0; j<5; j++)
188 { 188 {
189 eDiffCrossSection>>eDiffCrossSectionData 189 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
190 190
191 if (fasterCode) 191 if (fasterCode)
192 { 192 {
193 eNrjTransfData[j][tDummy][eDiffCrossSe 193 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
194 eProbaShellMap[j][tDummy].push_back(eD 194 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
195 } 195 }
196 196
197 // SI - only if eof is not reached 197 // SI - only if eof is not reached
198 if (!eDiffCrossSection.eof() && !fasterC 198 if (!eDiffCrossSection.eof() && !fasterCode)
199 { 199 {
200 eDiffCrossSectionData[j][tDummy][eDumm 200 eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
201 } 201 }
202 202
203 if (!fasterCode) eVecm[tDummy].push_back 203 if (!fasterCode) eVecm[tDummy].push_back(eDummy);
204 204
205 } 205 }
206 } 206 }
207 207
208 // 208 //
209 209
210 if( verboseLevel>0 ) 210 if( verboseLevel>0 )
211 { 211 {
212 G4cout << "Emfietzoglou ionisation model i 212 G4cout << "Emfietzoglou ionisation model is initialized " << G4endl
213 << "Energy range: " 213 << "Energy range: "
214 << LowEnergyLimit() / eV << " eV - " 214 << LowEnergyLimit() / eV << " eV - "
215 << HighEnergyLimit() / keV << " keV for " 215 << HighEnergyLimit() / keV << " keV for "
216 << particle->GetParticleName() 216 << particle->GetParticleName()
217 << G4endl; 217 << G4endl;
218 } 218 }
219 219
220 // Initialize water density pointer 220 // Initialize water density pointer
221 221
222 fpMolWaterDensity = 222 fpMolWaterDensity =
223 G4DNAMolecularMaterial::Instance()-> 223 G4DNAMolecularMaterial::Instance()->
224 GetNumMolPerVolTableFor(G4Material::Ge 224 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
225 225
226 // AD 226 // AD
227 227
228 fAtomDeexcitation = G4LossTableManager::Inst 228 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
229 229
230 if (isInitialised) 230 if (isInitialised)
231 { return;} 231 { return;}
232 fParticleChangeForGamma = GetParticleChangeF 232 fParticleChangeForGamma = GetParticleChangeForGamma();
233 isInitialised = true; 233 isInitialised = true;
234 } 234 }
235 235
236 //....oooOO0OOooo........oooOO0OOooo........oo 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
237 237
238 G4double G4DNAEmfietzoglouIonisationModel:: 238 G4double G4DNAEmfietzoglouIonisationModel::
239 CrossSectionPerVolume(const G4Material* materi 239 CrossSectionPerVolume(const G4Material* material,
240 const G4ParticleDefiniti 240 const G4ParticleDefinition* particleDefinition,
241 G4double ekin, 241 G4double ekin,
242 G4double, 242 G4double,
243 G4double) 243 G4double)
244 { 244 {
245 if(verboseLevel > 3) 245 if(verboseLevel > 3)
246 { 246 {
247 G4cout 247 G4cout
248 << "Calling CrossSectionPerVolume() of 248 << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouIonisationModel"
249 << G4endl; 249 << G4endl;
250 } 250 }
251 251
252 if (particleDefinition != G4Electron::Electr 252 if (particleDefinition != G4Electron::ElectronDefinition()) return 0; // necessary ??
253 253
254 // Calculate total cross section for model 254 // Calculate total cross section for model
255 255
256 G4double sigma=0; 256 G4double sigma=0;
257 257
258 G4double waterDensity = (*fpMolWaterDensity) 258 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
259 259
260 const G4String& particleName = particleDefin 260 const G4String& particleName = particleDefinition->GetParticleName();
261 261
262 if (ekin >= LowEnergyLimit() && ekin <= High 262 if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit())
263 { 263 {
264 std::map< G4String,G4DNACrossSectionDataSe 264 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
265 pos = tableData.find(particleName); 265 pos = tableData.find(particleName);
266 266
267 if (pos != tableData.end()) 267 if (pos != tableData.end())
268 { 268 {
269 G4DNACrossSectionDataSet* table = pos->s 269 G4DNACrossSectionDataSet* table = pos->second;
270 if (table != nullptr) 270 if (table != nullptr)
271 { 271 {
272 sigma = table->FindValue(ekin); 272 sigma = table->FindValue(ekin);
273 } 273 }
274 } 274 }
275 else 275 else
276 { 276 {
277 G4Exception("G4DNAEmfietzoglouIonisation 277 G4Exception("G4DNAEmfietzoglouIonisationModel::CrossSectionPerVolume","em0002",
278 FatalException,"Model not applicable 278 FatalException,"Model not applicable to particle type.");
279 } 279 }
280 } 280 }
281 281
282 if (verboseLevel > 2) 282 if (verboseLevel > 2)
283 { 283 {
284 G4cout << "_______________________________ 284 G4cout << "__________________________________" << G4endl;
285 G4cout << "G4DNAEmfietzoglouIonisationMode 285 G4cout << "G4DNAEmfietzoglouIonisationModel - XS INFO START" << G4endl;
286 G4cout << "Kinetic energy(eV)=" << ekin/eV 286 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
287 G4cout << "Cross section per water molecul 287 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
288 G4cout << "Cross section per water molecul 288 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
289 G4cout << "G4DNAEmfietzoglouIonisationMode 289 G4cout << "G4DNAEmfietzoglouIonisationModel - XS INFO END" << G4endl;
290 } 290 }
291 291
292 return sigma*waterDensity; 292 return sigma*waterDensity;
293 } 293 }
294 294
295 //....oooOO0OOooo........oooOO0OOooo........oo 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296 296
297 void G4DNAEmfietzoglouIonisationModel:: 297 void G4DNAEmfietzoglouIonisationModel::
298 SampleSecondaries(std::vector<G4DynamicParticl 298 SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
299 const G4MaterialCutsCouple* 299 const G4MaterialCutsCouple* couple,
300 const G4DynamicParticle* par 300 const G4DynamicParticle* particle,
301 G4double, 301 G4double,
302 G4double) 302 G4double)
303 { 303 {
304 304
305 if(verboseLevel > 3) 305 if(verboseLevel > 3)
306 { 306 {
307 G4cout << "Calling SampleSecondaries() of 307 G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouIonisationModel"
308 << G4endl; 308 << G4endl;
309 } 309 }
310 310
311 G4double k = particle->GetKineticEnergy(); 311 G4double k = particle->GetKineticEnergy();
312 312
313 const G4String& particleName = particle->Get 313 const G4String& particleName = particle->GetDefinition()->GetParticleName();
314 314
315 if (k >= LowEnergyLimit() && k <= HighEnergy 315 if (k >= LowEnergyLimit() && k <= HighEnergyLimit())
316 { 316 {
317 G4ParticleMomentum primaryDirection = part 317 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
318 G4double particleMass = particle->GetDefin 318 G4double particleMass = particle->GetDefinition()->GetPDGMass();
319 G4double totalEnergy = k + particleMass; 319 G4double totalEnergy = k + particleMass;
320 G4double pSquare = k * (totalEnergy + part 320 G4double pSquare = k * (totalEnergy + particleMass);
321 G4double totalMomentum = std::sqrt(pSquare 321 G4double totalMomentum = std::sqrt(pSquare);
322 322
323 G4int ionizationShell = 0; 323 G4int ionizationShell = 0;
324 324
325 ionizationShell = RandomSelect(k,particleN 325 ionizationShell = RandomSelect(k,particleName);
326 326
327 G4double bindingEnergy = 0; 327 G4double bindingEnergy = 0;
328 bindingEnergy = waterStructure.IonisationE 328 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
329 329
330 // SI : additional protection if tcs inter 330 // SI : additional protection if tcs interpolation method is modified
331 if (k<bindingEnergy) return; 331 if (k<bindingEnergy) return;
332 // 332 //
333 333
334 G4double secondaryKinetic=-1000*eV; 334 G4double secondaryKinetic=-1000*eV;
335 335
336 if (!fasterCode) secondaryKinetic = Random 336 if (!fasterCode) secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
337 337
338 if (fasterCode) 338 if (fasterCode)
339 secondaryKinetic = RandomizeEjectedElectro 339 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
340 340
341 // SI - For atom. deexc. tagging - 23/05/2 341 // SI - For atom. deexc. tagging - 23/05/2017
342 342
343 G4int Z = 8; 343 G4int Z = 8;
344 344
345 G4ThreeVector deltaDirection = 345 G4ThreeVector deltaDirection =
346 GetAngularDistribution()->SampleDirectionF 346 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
347 Z, ionizationShell, 347 Z, ionizationShell,
348 couple->GetMaterial()); 348 couple->GetMaterial());
349 349
350 if (secondaryKinetic>0) 350 if (secondaryKinetic>0)
351 { 351 {
352 auto dp = new G4DynamicParticle (G4Elec 352 auto dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
353 fvect->push_back(dp); 353 fvect->push_back(dp);
354 } 354 }
355 355
356 G4double deltaTotalMomentum = std::sqrt(se 356 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
357 357
358 G4double finalPx = totalMomentum*primaryDi 358 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
359 G4double finalPy = totalMomentum*primaryDi 359 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
360 G4double finalPz = totalMomentum*primaryDi 360 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
361 G4double finalMomentum = std::sqrt(finalPx 361 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
362 finalPx /= finalMomentum; 362 finalPx /= finalMomentum;
363 finalPy /= finalMomentum; 363 finalPy /= finalMomentum;
364 finalPz /= finalMomentum; 364 finalPz /= finalMomentum;
365 365
366 G4ThreeVector direction; 366 G4ThreeVector direction;
367 direction.set(finalPx,finalPy,finalPz); 367 direction.set(finalPx,finalPy,finalPz);
368 368
369 fParticleChangeForGamma->ProposeMomentumDi 369 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
370 370
371 // AM: sample deexcitation 371 // AM: sample deexcitation
372 // here we assume that H_{2}O electronic l 372 // here we assume that H_{2}O electronic levels are the same as Oxygen.
373 // this can be considered true with a roug 373 // this can be considered true with a rough 10% error in energy on K-shell,
374 374
375 size_t secNumberInit = 0;// need to know a 375 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
376 size_t secNumberFinal = 0;// So I'll make 376 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
377 377
378 G4double scatteredEnergy = k-bindingEnergy 378 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
379 379
380 // SI: only atomic deexcitation from K she 380 // SI: only atomic deexcitation from K shell is considered
381 if((fAtomDeexcitation != nullptr) && ioniz 381 if((fAtomDeexcitation != nullptr) && ionizationShell == 4)
382 { 382 {
383 const G4AtomicShell* shell 383 const G4AtomicShell* shell
384 = fAtomDeexcitation->GetAtomicShell(Z, 384 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
385 secNumberInit = fvect->size(); 385 secNumberInit = fvect->size();
386 fAtomDeexcitation->GenerateParticles(fve 386 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
387 secNumberFinal = fvect->size(); 387 secNumberFinal = fvect->size();
388 388
389 if(secNumberFinal > secNumberInit) { 389 if(secNumberFinal > secNumberInit) {
390 for (size_t i=secNumberInit; i<secNumberFina 390 for (size_t i=secNumberInit; i<secNumberFinal; ++i) {
391 //Check if there is enough residual 391 //Check if there is enough residual energy
392 if (bindingEnergy >= ((*fvect)[i])-> 392 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
393 { 393 {
394 //Ok, this is a valid secondary: 394 //Ok, this is a valid secondary: keep it
395 bindingEnergy -= ((*fvect)[i])->GetKine 395 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
396 } 396 }
397 else 397 else
398 { 398 {
399 //Invalid secondary: not enough energy 399 //Invalid secondary: not enough energy to create it!
400 //Keep its energy in the local deposit 400 //Keep its energy in the local deposit
401 delete (*fvect)[i]; 401 delete (*fvect)[i];
402 (*fvect)[i]=nullptr; 402 (*fvect)[i]=nullptr;
403 } 403 }
404 } 404 }
405 } 405 }
406 406
407 } 407 }
408 408
409 //This should never happen 409 //This should never happen
410 if(bindingEnergy < 0.0) 410 if(bindingEnergy < 0.0)
411 G4Exception("G4DNAEmfietzoglouIonisatioMo 411 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
412 "em2050",FatalException,"Nega 412 "em2050",FatalException,"Negative local energy deposit");
413 413
414 //bindingEnergy has been decreased 414 //bindingEnergy has been decreased
415 //by the amount of energy taken away by de 415 //by the amount of energy taken away by deexc. products
416 if (!statCode) 416 if (!statCode)
417 { 417 {
418 fParticleChangeForGamma->SetProposedKine 418 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
419 fParticleChangeForGamma->ProposeLocalEne 419 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
420 } 420 }
421 else 421 else
422 { 422 {
423 fParticleChangeForGamma->SetProposedKine 423 fParticleChangeForGamma->SetProposedKineticEnergy(k);
424 fParticleChangeForGamma->ProposeLocalEne 424 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
425 } 425 }
426 // TEST ////////////////////////// 426 // TEST //////////////////////////
427 // if (secondaryKinetic<0) abort(); 427 // if (secondaryKinetic<0) abort();
428 // if (scatteredEnergy<0) abort(); 428 // if (scatteredEnergy<0) abort();
429 // if (k-scatteredEnergy-secondaryKinetic- 429 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
430 // if (k-scatteredEnergy<0) abort(); 430 // if (k-scatteredEnergy<0) abort();
431 ///////////////////////////////// 431 /////////////////////////////////
432 432
433 const G4Track * theIncomingTrack = fPartic 433 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
434 G4DNAChemistryManager::Instance()->CreateW 434 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
435 ionizationShell, 435 ionizationShell,
436 theIncomingTrack); 436 theIncomingTrack);
437 } 437 }
438 438
439 } 439 }
440 440
441 //....oooOO0OOooo........oooOO0OOooo........oo 441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
442 442
443 G4double 443 G4double
444 G4DNAEmfietzoglouIonisationModel:: 444 G4DNAEmfietzoglouIonisationModel::
445 RandomizeEjectedElectronEnergy(G4ParticleDefin 445 RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
446 G4double k, 446 G4double k,
447 G4int shell) 447 G4int shell)
448 { 448 {
449 // G4cout << "*** SLOW computation for " 449 // G4cout << "*** SLOW computation for "
450 // << " " << particleDefinition->GetP 450 // << " " << particleDefinition->GetParticleName() << G4endl;
451 451
452 if(particleDefinition == G4Electron::Electro 452 if(particleDefinition == G4Electron::ElectronDefinition())
453 { 453 {
454 G4double maximumEnergyTransfer = 0.; 454 G4double maximumEnergyTransfer = 0.;
455 if((k + waterStructure.IonisationEnergy(sh 455 if((k + waterStructure.IonisationEnergy(shell)) / 2. > k)
456 maximumEnergyTransfer = k; 456 maximumEnergyTransfer = k;
457 else 457 else
458 maximumEnergyTransfer = (k + waterStruct 458 maximumEnergyTransfer = (k + waterStructure.IonisationEnergy(shell))/ 2.;
459 459
460 // SI : original method 460 // SI : original method
461 /* 461 /*
462 G4double crossSectionMaximum = 0.; 462 G4double crossSectionMaximum = 0.;
463 for(G4double value=waterStructure.Ionisat 463 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
464 { 464 {
465 G4double differentialCrossSection = Diffe 465 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
466 if(differentialCrossSection >= crossSecti 466 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
467 } 467 }
468 */ 468 */
469 469
470 // SI : alternative method 470 // SI : alternative method
471 G4double crossSectionMaximum = 0.; 471 G4double crossSectionMaximum = 0.;
472 472
473 G4double minEnergy = waterStructure.Ionisa 473 G4double minEnergy = waterStructure.IonisationEnergy(shell);
474 G4double maxEnergy = maximumEnergyTransfer 474 G4double maxEnergy = maximumEnergyTransfer;
475 G4int nEnergySteps = 50; 475 G4int nEnergySteps = 50;
476 476
477 G4double value(minEnergy); 477 G4double value(minEnergy);
478 G4double stpEnergy(std::pow(maxEnergy / va 478 G4double stpEnergy(std::pow(maxEnergy / value,
479 1. / static_ca 479 1. / static_cast<G4double>(nEnergySteps - 1)));
480 G4int step(nEnergySteps); 480 G4int step(nEnergySteps);
481 while(step > 0) 481 while(step > 0)
482 { 482 {
483 step--; 483 step--;
484 G4double differentialCrossSection = 484 G4double differentialCrossSection =
485 DifferentialCrossSection(particleDef 485 DifferentialCrossSection(particleDefinition,
486 k / eV, 486 k / eV,
487 value / eV, 487 value / eV,
488 shell); 488 shell);
489 if(differentialCrossSection >= crossSect 489 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum =
490 differentialCrossSection; 490 differentialCrossSection;
491 value *= stpEnergy; 491 value *= stpEnergy;
492 } 492 }
493 // 493 //
494 494
495 G4double secondaryElectronKineticEnergy = 495 G4double secondaryElectronKineticEnergy = 0.;
496 do 496 do
497 { 497 {
498 secondaryElectronKineticEnergy = G4Unifo 498 secondaryElectronKineticEnergy = G4UniformRand()* (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
499 }while(G4UniformRand()*crossSectionMaximum 499 }while(G4UniformRand()*crossSectionMaximum >
500 DifferentialCrossSection(particleDefin 500 DifferentialCrossSection(particleDefinition, k/eV,
501 (secondaryElectronKineticEnergy+wa 501 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
502 502
503 return secondaryElectronKineticEnergy; 503 return secondaryElectronKineticEnergy;
504 504
505 } 505 }
506 506
507 return 0; 507 return 0;
508 } 508 }
509 509
510 //....oooOO0OOooo........oooOO0OOooo........oo 510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
511 511
512 // The following section is not used anymore b 512 // The following section is not used anymore but is kept for memory
513 // GetAngularDistribution()->SampleDirectionFo 513 // GetAngularDistribution()->SampleDirectionForShell is used instead
514 514
515 /* 515 /*
516 void G4DNAEmfietzoglouIonisationModel::Random 516 void G4DNAEmfietzoglouIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
517 G4double k, 517 G4double k,
518 G4double secKinetic, 518 G4double secKinetic,
519 G4double & cosTheta, 519 G4double & cosTheta,
520 G4double & phi ) 520 G4double & phi )
521 { 521 {
522 if (particleDefinition == G4Electron::Electro 522 if (particleDefinition == G4Electron::ElectronDefinition())
523 { 523 {
524 phi = twopi * G4UniformRand(); 524 phi = twopi * G4UniformRand();
525 if (secKinetic < 50.*eV) cosTheta = (2.*G4Uni 525 if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
526 else if (secKinetic <= 200.*eV) 526 else if (secKinetic <= 200.*eV)
527 { 527 {
528 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4 528 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
529 else cosTheta = G4UniformRand()*(std::sqrt(2. 529 else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
530 } 530 }
531 else 531 else
532 { 532 {
533 G4double sin2O = (1.-secKinetic/k) / (1.+secK 533 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
534 cosTheta = std::sqrt(1.-sin2O); 534 cosTheta = std::sqrt(1.-sin2O);
535 } 535 }
536 } 536 }
537 537
538 else if (particleDefinition == G4Proton::Prot 538 else if (particleDefinition == G4Proton::ProtonDefinition())
539 { 539 {
540 G4double maxSecKinetic = 4.* (electron_mass_c 540 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
541 phi = twopi * G4UniformRand(); 541 phi = twopi * G4UniformRand();
542 542
543 // cosTheta = std::sqrt(secKinetic / maxSecKi 543 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
544 544
545 // Restriction below 100 eV from Emfietzoglou 545 // Restriction below 100 eV from Emfietzoglou (2000)
546 546
547 if (secKinetic>100*eV) cosTheta = std::sqrt(s 547 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
548 else cosTheta = (2.*G4UniformRand())-1.; 548 else cosTheta = (2.*G4UniformRand())-1.;
549 549
550 } 550 }
551 } 551 }
552 */ 552 */
553 553
554 //....oooOO0OOooo........oooOO0OOooo........oo 554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
555 G4double G4DNAEmfietzoglouIonisationModel::Dif 555 G4double G4DNAEmfietzoglouIonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
556 556 G4double k,
557 557 G4double energyTransfer,
558 558 G4int ionizationLevelIndex)
559 { 559 {
560 G4double sigma = 0.; 560 G4double sigma = 0.;
561 561
562 if(energyTransfer >= waterStructure.Ionisati 562 if(energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
563 { 563 {
564 G4double valueT1 = 0; 564 G4double valueT1 = 0;
565 G4double valueT2 = 0; 565 G4double valueT2 = 0;
566 G4double valueE21 = 0; 566 G4double valueE21 = 0;
567 G4double valueE22 = 0; 567 G4double valueE22 = 0;
568 G4double valueE12 = 0; 568 G4double valueE12 = 0;
569 G4double valueE11 = 0; 569 G4double valueE11 = 0;
570 570
571 G4double xs11 = 0; 571 G4double xs11 = 0;
572 G4double xs12 = 0; 572 G4double xs12 = 0;
573 G4double xs21 = 0; 573 G4double xs21 = 0;
574 G4double xs22 = 0; 574 G4double xs22 = 0;
575 575
576 if(particleDefinition == G4Electron::Elect 576 if(particleDefinition == G4Electron::ElectronDefinition())
577 { 577 {
578 // Protection against out of boundary ac 578 // Protection against out of boundary access
579 if (k==eTdummyVec.back()) k=k*(1.-1e-12) 579 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
580 // 580 //
581 581
582 // k should be in eV and energy transfer 582 // k should be in eV and energy transfer eV also
583 583
584 auto t2 = std::upper_bound(eTdummyVec.be 584 auto t2 = std::upper_bound(eTdummyVec.begin(),
585 585 eTdummyVec.end(),
586 586 k);
587 587
588 auto t1 = t2 - 1; 588 auto t1 = t2 - 1;
589 589
590 // SI : the following condition avoids s 590 // SI : the following condition avoids situations where energyTransfer >last vector element
591 // added strict limitations (09/08/2017) 591 // added strict limitations (09/08/2017)
592 if(energyTransfer < eVecm[(*t1)].back() 592 if(energyTransfer < eVecm[(*t1)].back() &&
593 energyTransfer < eVecm[(*t2)].back()) 593 energyTransfer < eVecm[(*t2)].back())
594 { 594 {
595 auto e12 = 595 auto e12 =
596 std::upper_bound(eVecm[(*t1)].begi 596 std::upper_bound(eVecm[(*t1)].begin(),
597 eVecm[(*t1)].end( 597 eVecm[(*t1)].end(),
598 energyTransfer); 598 energyTransfer);
599 auto e11 = e12 - 1; 599 auto e11 = e12 - 1;
600 600
601 auto e22 = 601 auto e22 =
602 std::upper_bound(eVecm[(*t2)].begi 602 std::upper_bound(eVecm[(*t2)].begin(),
603 eVecm[(*t2)].end( 603 eVecm[(*t2)].end(),
604 energyTransfer); 604 energyTransfer);
605 auto e21 = e22 - 1; 605 auto e21 = e22 - 1;
606 606
607 valueT1 = *t1; 607 valueT1 = *t1;
608 valueT2 = *t2; 608 valueT2 = *t2;
609 valueE21 = *e21; 609 valueE21 = *e21;
610 valueE22 = *e22; 610 valueE22 = *e22;
611 valueE12 = *e12; 611 valueE12 = *e12;
612 valueE11 = *e11; 612 valueE11 = *e11;
613 613
614 xs11 = eDiffCrossSectionData[ionizatio 614 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
615 xs12 = eDiffCrossSectionData[ionizatio 615 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
616 xs21 = eDiffCrossSectionData[ionizatio 616 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
617 xs22 = eDiffCrossSectionData[ionizatio 617 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
618 618
619 //G4cout << "-------------------" << G 619 //G4cout << "-------------------" << G4endl;
620 //G4cout << "ionizationLevelIndex=" << 620 //G4cout << "ionizationLevelIndex=" << ionizationLevelIndex << G4endl;
621 //G4cout << "valueT1/eV=" << valueT1 < 621 //G4cout << "valueT1/eV=" << valueT1 << " valueT2/eV=" << valueT2 << G4endl;
622 //G4cout << "valueE11/eV=" << valueE11 622 //G4cout << "valueE11/eV=" << valueE11 << " valueE12/eV=" << valueE12
623 // << " valueE21/eV=" << valueE2 623 // << " valueE21/eV=" << valueE21 << " valueE22/eV=" << valueE22 << G4endl;
624 //G4cout << "xs11=" << xs11 / ((1.e-22 624 //G4cout << "xs11=" << xs11 / ((1.e-22 / 3.343) * m*m) << G4endl;
625 //G4cout << "xs12=" << xs12 / ((1.e-22 625 //G4cout << "xs12=" << xs12 / ((1.e-22 / 3.343) * m*m) << G4endl;
626 //G4cout << "xs21=" << xs21 / ((1.e-22 626 //G4cout << "xs21=" << xs21 / ((1.e-22 / 3.343) * m*m) << G4endl;
627 //G4cout << "xs22=" << xs22 / ((1.e-22 627 //G4cout << "xs22=" << xs22 / ((1.e-22 / 3.343) * m*m) << G4endl;
628 //G4cout << "###################" << G 628 //G4cout << "###################" << G4endl;
629 629
630 } 630 }
631 631
632 } 632 }
633 633
634 G4double xsProduct = xs11 * xs12 * xs21 * 634 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
635 if(xsProduct != 0.) 635 if(xsProduct != 0.)
636 { 636 {
637 sigma = QuadInterpolator(valueE11, 637 sigma = QuadInterpolator(valueE11,
638 valueE12, 638 valueE12,
639 valueE21, 639 valueE21,
640 valueE22, 640 valueE22,
641 xs11, 641 xs11,
642 xs12, 642 xs12,
643 xs21, 643 xs21,
644 xs22, 644 xs22,
645 valueT1, 645 valueT1,
646 valueT2, 646 valueT2,
647 k, 647 k,
648 energyTransfer) 648 energyTransfer);
649 } 649 }
650 650
651 } 651 }
652 652
653 return sigma; 653 return sigma;
654 } 654 }
655 655
656 //....oooOO0OOooo........oooOO0OOooo........oo 656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
657 657
658 G4double G4DNAEmfietzoglouIonisationModel::Int 658 G4double G4DNAEmfietzoglouIonisationModel::Interpolate(G4double e1,
659 659 G4double e2,
660 660 G4double e,
661 661 G4double xs1,
662 662 G4double xs2)
663 { 663 {
664 664
665 G4double value = 0.; 665 G4double value = 0.;
666 666
667 // Log-log interpolation by default 667 // Log-log interpolation by default
668 668
669 if(e1 != 0 && e2 != 0 && (std::log10(e2) - s 669 if(e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
670 && !fasterCode) 670 && !fasterCode)
671 { 671 {
672 G4double a = (std::log10(xs2) - std::log10 672 G4double a = (std::log10(xs2) - std::log10(xs1))
673 / (std::log10(e2) - std::log10(e1)); 673 / (std::log10(e2) - std::log10(e1));
674 G4double b = std::log10(xs2) - a * std::lo 674 G4double b = std::log10(xs2) - a * std::log10(e2);
675 G4double sigma = a * std::log10(e) + b; 675 G4double sigma = a * std::log10(e) + b;
676 value = (std::pow(10., sigma)); 676 value = (std::pow(10., sigma));
677 } 677 }
678 678
679 // Switch to lin-lin interpolation 679 // Switch to lin-lin interpolation
680 /* 680 /*
681 if ((e2-e1)!=0) 681 if ((e2-e1)!=0)
682 { 682 {
683 G4double d1 = xs1; 683 G4double d1 = xs1;
684 G4double d2 = xs2; 684 G4double d2 = xs2;
685 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1) 685 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
686 } 686 }
687 */ 687 */
688 688
689 // Switch to log-lin interpolation for faste 689 // Switch to log-lin interpolation for faster code
690 if((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && 690 if((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
691 { 691 {
692 G4double d1 = std::log10(xs1); 692 G4double d1 = std::log10(xs1);
693 G4double d2 = std::log10(xs2); 693 G4double d2 = std::log10(xs2);
694 value = std::pow(10., (d1 + (d2 - d1) * (e 694 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
695 } 695 }
696 696
697 // Switch to lin-lin interpolation for faste 697 // Switch to lin-lin interpolation for faster code
698 // in case one of xs1 or xs2 (=cum proba) va 698 // in case one of xs1 or xs2 (=cum proba) value is zero
699 699
700 if((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) 700 if((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
701 { 701 {
702 G4double d1 = xs1; 702 G4double d1 = xs1;
703 G4double d2 = xs2; 703 G4double d2 = xs2;
704 value = (d1 + (d2 - d1) * (e - e1) / (e2 - 704 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
705 } 705 }
706 706
707 /* 707 /*
708 G4cout 708 G4cout
709 << e1 << " " 709 << e1 << " "
710 << e2 << " " 710 << e2 << " "
711 << e << " " 711 << e << " "
712 << xs1 << " " 712 << xs1 << " "
713 << xs2 << " " 713 << xs2 << " "
714 << value 714 << value
715 << G4endl; 715 << G4endl;
716 */ 716 */
717 717
718 return value; 718 return value;
719 } 719 }
720 720
721 //....oooOO0OOooo........oooOO0OOooo........oo 721 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
722 722
723 G4double G4DNAEmfietzoglouIonisationModel::Qua 723 G4double G4DNAEmfietzoglouIonisationModel::QuadInterpolator(G4double e11,
724 724 G4double e12,
725 725 G4double e21,
726 726 G4double e22,
727 727 G4double xs11,
728 728 G4double xs12,
729 729 G4double xs21,
730 730 G4double xs22,
731 731 G4double t1,
732 732 G4double t2,
733 733 G4double t,
734 734 G4double e)
735 { 735 {
736 G4double interpolatedvalue1 = Interpolate(e1 736 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
737 G4double interpolatedvalue2 = Interpolate(e2 737 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
738 G4double value = Interpolate(t1, 738 G4double value = Interpolate(t1,
739 t2, 739 t2,
740 t, 740 t,
741 interpolatedval 741 interpolatedvalue1,
742 interpolatedval 742 interpolatedvalue2);
743 743
744 return value; 744 return value;
745 } 745 }
746 746
747 //....oooOO0OOooo........oooOO0OOooo........oo 747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
748 748
749 G4int G4DNAEmfietzoglouIonisationModel::Random 749 G4int G4DNAEmfietzoglouIonisationModel::RandomSelect(G4double k,
750 750 const G4String& particle)
751 { 751 {
752 G4int level = 0; 752 G4int level = 0;
753 753
754 auto pos = tableData.find(particle); 754 auto pos = tableData.find(particle);
755 755
756 if(pos != tableData.cend()) 756 if(pos != tableData.cend())
757 { 757 {
758 G4DNACrossSectionDataSet* table = pos->sec 758 G4DNACrossSectionDataSet* table = pos->second;
759 759
760 if(table != nullptr) 760 if(table != nullptr)
761 { 761 {
762 auto valuesBuffer = new G4double[table- 762 auto valuesBuffer = new G4double[table->NumberOfComponents()];
763 const auto n = (G4int)table->NumberOfCo 763 const auto n = (G4int)table->NumberOfComponents();
764 G4int i(n); 764 G4int i(n);
765 G4double value = 0.; 765 G4double value = 0.;
766 766
767 while(i > 0) 767 while(i > 0)
768 { 768 {
769 i--; 769 i--;
770 valuesBuffer[i] = table->GetComponent( 770 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
771 value += valuesBuffer[i]; 771 value += valuesBuffer[i];
772 } 772 }
773 773
774 value *= G4UniformRand(); 774 value *= G4UniformRand();
775 775
776 i = n; 776 i = n;
777 777
778 while(i > 0) 778 while(i > 0)
779 { 779 {
780 i--; 780 i--;
781 781
782 if(valuesBuffer[i] > value) 782 if(valuesBuffer[i] > value)
783 { 783 {
784 delete[] valuesBuffer; 784 delete[] valuesBuffer;
785 return i; 785 return i;
786 } 786 }
787 value -= valuesBuffer[i]; 787 value -= valuesBuffer[i];
788 } 788 }
789 789
790 delete[] valuesBuffer; 790 delete[] valuesBuffer;
791 791
792 } 792 }
793 } 793 }
794 else 794 else
795 { 795 {
796 G4Exception("G4DNAEmfietzoglouIonisationMo 796 G4Exception("G4DNAEmfietzoglouIonisationModel::RandomSelect",
797 "em0002", 797 "em0002",
798 FatalException, 798 FatalException,
799 "Model not applicable to parti 799 "Model not applicable to particle type.");
800 } 800 }
801 801
802 return level; 802 return level;
803 } 803 }
804 804
805 //....oooOO0OOooo........oooOO0OOooo........oo 805 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
806 806
807 G4double G4DNAEmfietzoglouIonisationModel::Ran 807 G4double G4DNAEmfietzoglouIonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
808 808 G4double k,
809 809 G4int shell)
810 { 810 {
811 //G4cout << "*** FAST computation for " << " 811 //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
812 812
813 G4double secondaryElectronKineticEnergy = 0. 813 G4double secondaryElectronKineticEnergy = 0.;
814 814
815 secondaryElectronKineticEnergy = RandomTrans 815 secondaryElectronKineticEnergy = RandomTransferedEnergy(particleDefinition,
816 816 k / eV,
817 817 shell)
818 * eV 818 * eV
819 - waterStru 819 - waterStructure.IonisationEnergy(shell);
820 820
821 //G4cout << RandomTransferedEnergy(particleD 821 //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
822 if(secondaryElectronKineticEnergy < 0.) retu 822 if(secondaryElectronKineticEnergy < 0.) return 0.;
823 823
824 return secondaryElectronKineticEnergy; 824 return secondaryElectronKineticEnergy;
825 } 825 }
826 826
827 //....oooOO0OOooo........oooOO0OOooo........oo 827 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
828 828
829 G4double G4DNAEmfietzoglouIonisationModel::Ran 829 G4double G4DNAEmfietzoglouIonisationModel::RandomTransferedEnergy(G4ParticleDefinition* particleDefinition,
830 830 G4double k,
831 831 G4int ionizationLevelIndex)
832 { 832 {
833 833
834 G4double random = G4UniformRand(); 834 G4double random = G4UniformRand();
835 835
836 G4double nrj = 0.; 836 G4double nrj = 0.;
837 837
838 G4double valueK1 = 0; 838 G4double valueK1 = 0;
839 G4double valueK2 = 0; 839 G4double valueK2 = 0;
840 G4double valuePROB21 = 0; 840 G4double valuePROB21 = 0;
841 G4double valuePROB22 = 0; 841 G4double valuePROB22 = 0;
842 G4double valuePROB12 = 0; 842 G4double valuePROB12 = 0;
843 G4double valuePROB11 = 0; 843 G4double valuePROB11 = 0;
844 844
845 G4double nrjTransf11 = 0; 845 G4double nrjTransf11 = 0;
846 G4double nrjTransf12 = 0; 846 G4double nrjTransf12 = 0;
847 G4double nrjTransf21 = 0; 847 G4double nrjTransf21 = 0;
848 G4double nrjTransf22 = 0; 848 G4double nrjTransf22 = 0;
849 849
850 if (particleDefinition == G4Electron::Electr 850 if (particleDefinition == G4Electron::ElectronDefinition())
851 { 851 {
852 // Protection against out of boundary acce 852 // Protection against out of boundary access
853 if (k==eTdummyVec.back()) k=k*(1.-1e-12); 853 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
854 // 854 //
855 855
856 // k should be in eV 856 // k should be in eV
857 auto k2 = std::upper_bound(eTdummyVec.begi 857 auto k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
858 858
859 auto k1 = k2-1; 859 auto k1 = k2-1;
860 860
861 /* 861 /*
862 G4cout << "----> k=" << k 862 G4cout << "----> k=" << k
863 << " " << *k1 863 << " " << *k1
864 << " " << *k2 864 << " " << *k2
865 << " " << random 865 << " " << random
866 << " " << ionizationLevelIndex 866 << " " << ionizationLevelIndex
867 << " " << eProbaShellMap[ionizationLevelI 867 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
868 << " " << eProbaShellMap[ionizationLevelI 868 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
869 << G4endl; 869 << G4endl;
870 */ 870 */
871 871
872 // SI : the following condition avoids sit 872 // SI : the following condition avoids situations where random >last vector element
873 if ( random <= eProbaShellMap[ionizationLe 873 if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
874 && random <= eProbaShellMap[ionization 874 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
875 875
876 { 876 {
877 auto prob12 = std::upper_bound(eProbaShe 877 auto prob12 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
878 eProbaShellMap[ionizationLevelIndex] 878 eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
879 879
880 auto prob11 = prob12-1; 880 auto prob11 = prob12-1;
881 881
882 auto prob22 = std::upper_bound(eProbaShe 882 auto prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
883 eProbaShellMap[ionizationLevelIndex] 883 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
884 884
885 auto prob21 = prob22-1; 885 auto prob21 = prob22-1;
886 886
887 valueK1 =*k1; 887 valueK1 =*k1;
888 valueK2 =*k2; 888 valueK2 =*k2;
889 valuePROB21 =*prob21; 889 valuePROB21 =*prob21;
890 valuePROB22 =*prob22; 890 valuePROB22 =*prob22;
891 valuePROB12 =*prob12; 891 valuePROB12 =*prob12;
892 valuePROB11 =*prob11; 892 valuePROB11 =*prob11;
893 893
894 /* 894 /*
895 G4cout << " " << random << " " < 895 G4cout << " " << random << " " << valuePROB11 << " "
896 << valuePROB12 << " " << valuePROB21 << 896 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
897 */ 897 */
898 898
899 nrjTransf11 = eNrjTransfData[ionizationL 899 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
900 nrjTransf12 = eNrjTransfData[ionizationL 900 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
901 nrjTransf21 = eNrjTransfData[ionizationL 901 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
902 nrjTransf22 = eNrjTransfData[ionizationL 902 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
903 903
904 /* 904 /*
905 G4cout << " " << ionizationLevel 905 G4cout << " " << ionizationLevelIndex << " "
906 << random << " " <<valueK1 << " " << va 906 << random << " " <<valueK1 << " " << valueK2 << G4endl;
907 907
908 G4cout << " " << random << " " < 908 G4cout << " " << random << " " << nrjTransf11 << " "
909 << nrjTransf12 << " " << nrjTransf21 << 909 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
910 */ 910 */
911 911
912 } 912 }
913 913
914 // Avoids cases where cum xs is zero for k 914 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
915 915
916 if ( random > eProbaShellMap[ionizationLev 916 if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
917 917
918 { 918 {
919 auto prob22 = std::upper_bound(eProbaShe 919 auto prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
920 eProbaShellMap[ionizationLevelIndex] 920 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
921 921
922 auto prob21 = prob22-1; 922 auto prob21 = prob22-1;
923 923
924 valueK1 =*k1; 924 valueK1 =*k1;
925 valueK2 =*k2; 925 valueK2 =*k2;
926 valuePROB21 =*prob21; 926 valuePROB21 =*prob21;
927 valuePROB22 =*prob22; 927 valuePROB22 =*prob22;
928 928
929 //G4cout << " " << random << " " 929 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
930 930
931 nrjTransf21 = eNrjTransfData[ionizationL 931 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
932 nrjTransf22 = eNrjTransfData[ionizationL 932 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
933 933
934 G4double interpolatedvalue2 = Interpolat 934 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
935 935
936 // zeros are explicitly set 936 // zeros are explicitly set
937 937
938 G4double value = Interpolate(valueK1, va 938 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
939 939
940 /* 940 /*
941 G4cout << " " << ionizationLevel 941 G4cout << " " << ionizationLevelIndex << " "
942 << random << " " <<valueK1 << " " << va 942 << random << " " <<valueK1 << " " << valueK2 << G4endl;
943 943
944 G4cout << " " << random << " " < 944 G4cout << " " << random << " " << nrjTransf11 << " "
945 << nrjTransf12 << " " << nrjTransf21 << 945 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
946 946
947 G4cout << "ici" << " " << value << G4en 947 G4cout << "ici" << " " << value << G4endl;
948 */ 948 */
949 949
950 return value; 950 return value;
951 } 951 }
952 952
953 } 953 }
954 954
955 // End electron 955 // End electron
956 956
957 G4double nrjTransfProduct = nrjTransf11 * nr 957 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
958 958
959 //G4cout << "nrjTransfProduct=" << nrjTransf 959 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
960 960
961 if (nrjTransfProduct != 0.) 961 if (nrjTransfProduct != 0.)
962 { 962 {
963 nrj = QuadInterpolator( valuePROB11, value 963 nrj = QuadInterpolator( valuePROB11, valuePROB12,
964 valuePROB21, valuePROB22, 964 valuePROB21, valuePROB22,
965 nrjTransf11, nrjTransf12, 965 nrjTransf11, nrjTransf12,
966 nrjTransf21, nrjTransf22, 966 nrjTransf21, nrjTransf22,
967 valueK1, valueK2, 967 valueK1, valueK2,
968 k, random); 968 k, random);
969 } 969 }
970 970
971 //G4cout << nrj << endl; 971 //G4cout << nrj << endl;
972 972
973 return nrj; 973 return nrj;
974 } 974 }
975 975