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