| 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 // $Id: G4DNABornIonisationModel1.cc 87631 2014-12-14 12:42:05Z matkara $
26 // 27 //
27 28
28 #include "G4DNABornIonisationModel1.hh" 29 #include "G4DNABornIonisationModel1.hh"
29 #include "G4PhysicalConstants.hh" 30 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh"
31 #include "G4UAtomicDeexcitation.hh" 32 #include "G4UAtomicDeexcitation.hh"
32 #include "G4LossTableManager.hh" 33 #include "G4LossTableManager.hh"
33 #include "G4DNAChemistryManager.hh" 34 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMolecularMaterial.hh" 35 #include "G4DNAMolecularMaterial.hh"
35 #include "G4DNABornAngle.hh" 36 #include "G4DNABornAngle.hh"
36 #include "G4DeltaAngle.hh" 37 #include "G4DeltaAngle.hh"
37 #include "G4Exp.hh" 38 #include "G4Exp.hh"
38 39
39 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 41
41 using namespace std; 42 using namespace std;
42 43
43 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 45
45 G4DNABornIonisationModel1::G4DNABornIonisation 46 G4DNABornIonisationModel1::G4DNABornIonisationModel1(const G4ParticleDefinition*,
46 47 const G4String& nam) :
47 G4VEmModel(nam) << 48 G4VEmModel(nam), isInitialised(false)
48 { 49 {
49 verboseLevel = 0; 50 verboseLevel = 0;
50 // Verbosity scale: 51 // Verbosity scale:
51 // 0 = nothing 52 // 0 = nothing
52 // 1 = warning for energy non-conservation 53 // 1 = warning for energy non-conservation
53 // 2 = details of energy budget 54 // 2 = details of energy budget
54 // 3 = calculation of cross sections, file o 55 // 3 = calculation of cross sections, file openings, sampling of atoms
55 // 4 = entering in methods 56 // 4 = entering in methods
56 57
57 if (verboseLevel > 0) 58 if (verboseLevel > 0)
58 { 59 {
59 G4cout << "Born ionisation model is constr 60 G4cout << "Born ionisation model is constructed " << G4endl;
60 } 61 }
61 62
62 // Mark this model as "applicable" for atomi 63 // Mark this model as "applicable" for atomic deexcitation
63 SetDeexcitationFlag(true); 64 SetDeexcitationFlag(true);
64 fAtomDeexcitation = nullptr; << 65 fAtomDeexcitation = 0;
65 fParticleChangeForGamma = nullptr; << 66 fParticleChangeForGamma = 0;
66 fpMolWaterDensity = nullptr; << 67 fpMolWaterDensity = 0;
67 68
68 // Define default angular generator 69 // Define default angular generator
69 SetAngularDistribution(new G4DNABornAngle()) 70 SetAngularDistribution(new G4DNABornAngle());
70 71
71 // Selection of computation method 72 // Selection of computation method
72 73
73 fasterCode = false; 74 fasterCode = false;
74 75
75 // Selection of stationary mode 76 // Selection of stationary mode
76 77
77 statCode = false; 78 statCode = false;
78 79
79 // Selection of SP scaling 80 // Selection of SP scaling
80 81
81 spScaling = true; 82 spScaling = true;
82 } 83 }
83 84
84 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 86
86 G4DNABornIonisationModel1::~G4DNABornIonisatio 87 G4DNABornIonisationModel1::~G4DNABornIonisationModel1()
87 { 88 {
88 // Cross section 89 // Cross section
89 90
90 std::map<G4String, G4DNACrossSectionDataSet* 91 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
91 for (pos = tableData.begin(); pos != tableDa 92 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
92 { 93 {
93 G4DNACrossSectionDataSet* table = pos->sec 94 G4DNACrossSectionDataSet* table = pos->second;
94 delete table; 95 delete table;
95 } 96 }
96 97
97 // Final state 98 // Final state
98 99
99 eVecm.clear(); 100 eVecm.clear();
100 pVecm.clear(); 101 pVecm.clear();
101 } 102 }
102 103
103 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104 105
105 void G4DNABornIonisationModel1::Initialise(con 106 void G4DNABornIonisationModel1::Initialise(const G4ParticleDefinition* particle,
106 con 107 const G4DataVector& /*cuts*/)
107 { 108 {
108 109
109 if (verboseLevel > 3) 110 if (verboseLevel > 3)
110 { 111 {
111 G4cout << "Calling G4DNABornIonisationMode 112 G4cout << "Calling G4DNABornIonisationModel1::Initialise()" << G4endl;
112 } 113 }
113 114
114 // Energy limits 115 // Energy limits
115 116
116 G4String fileElectron("dna/sigma_ionisation_ 117 G4String fileElectron("dna/sigma_ionisation_e_born");
117 G4String fileProton("dna/sigma_ionisation_p_ 118 G4String fileProton("dna/sigma_ionisation_p_born");
118 119
119 G4ParticleDefinition* electronDef = G4Electr 120 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
120 G4ParticleDefinition* protonDef = G4Proton:: 121 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
121 122
122 G4String electron; 123 G4String electron;
123 G4String proton; 124 G4String proton;
124 125
125 G4double scaleFactor = (1.e-22 / 3.343) * m* 126 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
126 127
127 const char *path = G4FindDataDir("G4LEDATA") << 128 char *path = getenv("G4LEDATA");
128 129
129 // *** ELECTRON 130 // *** ELECTRON
130 131
131 electron = electronDef->GetParticleName(); 132 electron = electronDef->GetParticleName();
132 133
133 tableFile[electron] = fileElectron; 134 tableFile[electron] = fileElectron;
134 135
135 lowEnergyLimit[electron] = 11. * eV; 136 lowEnergyLimit[electron] = 11. * eV;
136 highEnergyLimit[electron] = 1. * MeV; 137 highEnergyLimit[electron] = 1. * MeV;
137 138
138 // Cross section 139 // Cross section
139 140
140 auto tableE = new G4DNACrossSectionDataSet( << 141 G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
141 tableE->LoadData(fileElectron); 142 tableE->LoadData(fileElectron);
142 143
143 tableData[electron] = tableE; 144 tableData[electron] = tableE;
144 145
145 // Final state 146 // Final state
146 147
147 std::ostringstream eFullFileName; 148 std::ostringstream eFullFileName;
148 149
149 if (fasterCode) eFullFileName << path << "/d 150 if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
150 if (!fasterCode) eFullFileName << path << "/ 151 if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
151 152
152 std::ifstream eDiffCrossSection(eFullFileNam 153 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
153 154
154 if (!eDiffCrossSection) 155 if (!eDiffCrossSection)
155 { 156 {
156 if (fasterCode) G4Exception("G4DNABornIoni 157 if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
157 FatalException,"Missing data file:/dna 158 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat");
158 159
159 if (!fasterCode) G4Exception("G4DNABornIon 160 if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
160 FatalException,"Missing data file:/dna 161 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
161 } 162 }
162 163
163 // Clear the arrays for re-initialization ca 164 // Clear the arrays for re-initialization case (MT mode)
164 // March 25th, 2014 - Vaclav Stepan, Sebasti 165 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
165 166
166 eTdummyVec.clear(); 167 eTdummyVec.clear();
167 pTdummyVec.clear(); 168 pTdummyVec.clear();
168 169
169 eVecm.clear(); 170 eVecm.clear();
170 pVecm.clear(); 171 pVecm.clear();
171 172
172 for (G4int j=0; j<5; j++) 173 for (G4int j=0; j<5; j++)
173 { 174 {
174 eProbaShellMap[j].clear(); 175 eProbaShellMap[j].clear();
175 pProbaShellMap[j].clear(); 176 pProbaShellMap[j].clear();
176 177
177 eDiffCrossSectionData[j].clear(); 178 eDiffCrossSectionData[j].clear();
178 pDiffCrossSectionData[j].clear(); 179 pDiffCrossSectionData[j].clear();
179 180
180 eNrjTransfData[j].clear(); 181 eNrjTransfData[j].clear();
181 pNrjTransfData[j].clear(); 182 pNrjTransfData[j].clear();
182 } 183 }
183 184
184 // 185 //
185 186
186 eTdummyVec.push_back(0.); 187 eTdummyVec.push_back(0.);
187 while(!eDiffCrossSection.eof()) 188 while(!eDiffCrossSection.eof())
188 { 189 {
189 G4double tDummy; 190 G4double tDummy;
190 G4double eDummy; 191 G4double eDummy;
191 eDiffCrossSection>>tDummy>>eDummy; 192 eDiffCrossSection>>tDummy>>eDummy;
192 if (tDummy != eTdummyVec.back()) eTdummyVe 193 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
193 194
194 G4double tmp; 195 G4double tmp;
195 for (G4int j=0; j<5; j++) 196 for (G4int j=0; j<5; j++)
196 { 197 {
197 eDiffCrossSection>> tmp; 198 eDiffCrossSection>> tmp;
198 199
199 eDiffCrossSectionData[j][tDummy][eDummy] 200 eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
200 201
201 if (fasterCode) 202 if (fasterCode)
202 { 203 {
203 eNrjTransfData[j][tDummy][eDiffCrossSe 204 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
204 eProbaShellMap[j][tDummy].push_back(eD 205 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
205 } 206 }
206 207
207 // SI - only if eof is not reached 208 // SI - only if eof is not reached
208 if (!eDiffCrossSection.eof() && !fasterC 209 if (!eDiffCrossSection.eof() && !fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
209 210
210 if (!fasterCode) eVecm[tDummy].push_back 211 if (!fasterCode) eVecm[tDummy].push_back(eDummy);
211 212
212 } 213 }
213 } 214 }
214 215
215 // *** PROTON 216 // *** PROTON
216 217
217 proton = protonDef->GetParticleName(); 218 proton = protonDef->GetParticleName();
218 219
219 tableFile[proton] = fileProton; 220 tableFile[proton] = fileProton;
220 221
221 lowEnergyLimit[proton] = 500. * keV; 222 lowEnergyLimit[proton] = 500. * keV;
222 highEnergyLimit[proton] = 100. * MeV; 223 highEnergyLimit[proton] = 100. * MeV;
223 224
224 // Cross section 225 // Cross section
225 226
226 auto tableP = new G4DNACrossSectionDataSet( << 227 G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
227 tableP->LoadData(fileProton); 228 tableP->LoadData(fileProton);
228 229
229 tableData[proton] = tableP; 230 tableData[proton] = tableP;
230 231
231 // Final state 232 // Final state
232 233
233 std::ostringstream pFullFileName; 234 std::ostringstream pFullFileName;
234 235
235 if (fasterCode) pFullFileName << path << "/d 236 if (fasterCode) pFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
236 237
237 if (!fasterCode) pFullFileName << path << "/ 238 if (!fasterCode) pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
238 239
239 std::ifstream pDiffCrossSection(pFullFileNam 240 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
240 241
241 if (!pDiffCrossSection) 242 if (!pDiffCrossSection)
242 { 243 {
243 if (fasterCode) G4Exception("G4DNABornIoni 244 if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
244 FatalException,"Missing data file:/dna 245 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat");
245 246
246 if (!fasterCode) G4Exception("G4DNABornIon 247 if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
247 FatalException,"Missing data file:/dna 248 FatalException,"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
248 } 249 }
249 250
250 pTdummyVec.push_back(0.); 251 pTdummyVec.push_back(0.);
251 while(!pDiffCrossSection.eof()) 252 while(!pDiffCrossSection.eof())
252 { 253 {
253 G4double tDummy; 254 G4double tDummy;
254 G4double eDummy; 255 G4double eDummy;
255 pDiffCrossSection>>tDummy>>eDummy; 256 pDiffCrossSection>>tDummy>>eDummy;
256 if (tDummy != pTdummyVec.back()) pTdummyVe 257 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
257 for (G4int j=0; j<5; j++) 258 for (G4int j=0; j<5; j++)
258 { 259 {
259 pDiffCrossSection>>pDiffCrossSectionData 260 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
260 261
261 if (fasterCode) 262 if (fasterCode)
262 { 263 {
263 pNrjTransfData[j][tDummy][pDiffCrossSe 264 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
264 pProbaShellMap[j][tDummy].push_back(pD 265 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
265 } 266 }
266 267
267 // SI - only if eof is not reached ! 268 // SI - only if eof is not reached !
268 if (!pDiffCrossSection.eof() && !fasterC 269 if (!pDiffCrossSection.eof() && !fasterCode) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
269 270
270 if (!fasterCode) pVecm[tDummy].push_back 271 if (!fasterCode) pVecm[tDummy].push_back(eDummy);
271 } 272 }
272 } 273 }
273 274
274 // 275 //
275 276
276 if (particle==electronDef) 277 if (particle==electronDef)
277 { 278 {
278 SetLowEnergyLimit(lowEnergyLimit[electron] 279 SetLowEnergyLimit(lowEnergyLimit[electron]);
279 SetHighEnergyLimit(highEnergyLimit[electro 280 SetHighEnergyLimit(highEnergyLimit[electron]);
280 } 281 }
281 282
282 if (particle==protonDef) 283 if (particle==protonDef)
283 { 284 {
284 SetLowEnergyLimit(lowEnergyLimit[proton]); 285 SetLowEnergyLimit(lowEnergyLimit[proton]);
285 SetHighEnergyLimit(highEnergyLimit[proton] 286 SetHighEnergyLimit(highEnergyLimit[proton]);
286 } 287 }
287 288
288 if( verboseLevel>0 ) 289 if( verboseLevel>0 )
289 { 290 {
290 G4cout << "Born ionisation model is initia 291 G4cout << "Born ionisation model is initialized " << G4endl
291 << "Energy range: " 292 << "Energy range: "
292 << LowEnergyLimit() / eV << " eV - " 293 << LowEnergyLimit() / eV << " eV - "
293 << HighEnergyLimit() / keV << " keV for " 294 << HighEnergyLimit() / keV << " keV for "
294 << particle->GetParticleName() 295 << particle->GetParticleName()
295 << G4endl; 296 << G4endl;
296 } 297 }
297 298
298 // Initialize water density pointer 299 // Initialize water density pointer
299 300
300 fpMolWaterDensity = G4DNAMolecularMaterial:: 301 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
301 GetNumMolPerVolTableFor(G4Material::GetMater 302 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
302 303
303 // AD 304 // AD
304 305
305 fAtomDeexcitation = G4LossTableManager::Inst 306 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
306 307
307 // 308 //
308 309
309 if (isInitialised) 310 if (isInitialised)
310 { return;} 311 { return;}
311 fParticleChangeForGamma = GetParticleChangeF 312 fParticleChangeForGamma = GetParticleChangeForGamma();
312 isInitialised = true; 313 isInitialised = true;
313 } 314 }
314 315
315 //....oooOO0OOooo........oooOO0OOooo........oo 316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
316 317
317 G4double G4DNABornIonisationModel1::CrossSecti 318 G4double G4DNABornIonisationModel1::CrossSectionPerVolume(const G4Material* material,
318 319 const G4ParticleDefinition* particleDefinition,
319 320 G4double ekin,
320 321 G4double,
321 322 G4double)
322 { 323 {
323 if (verboseLevel > 3) 324 if (verboseLevel > 3)
324 { 325 {
325 G4cout << "Calling CrossSectionPerVolume() 326 G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel1"
326 << G4endl; 327 << G4endl;
327 } 328 }
328 329
329 if ( 330 if (
330 particleDefinition != G4Proton::ProtonDe 331 particleDefinition != G4Proton::ProtonDefinition()
331 && 332 &&
332 particleDefinition != G4Electron::Electr 333 particleDefinition != G4Electron::ElectronDefinition()
333 ) 334 )
334 335
335 return 0; 336 return 0;
336 337
337 // Calculate total cross section for model 338 // Calculate total cross section for model
338 339
339 G4double lowLim = 0; 340 G4double lowLim = 0;
340 G4double highLim = 0; 341 G4double highLim = 0;
341 G4double sigma=0; 342 G4double sigma=0;
342 343
343 G4double waterDensity = (*fpMolWaterDensity) 344 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
344 345
345 const G4String& particleName = particleDefin << 346 if(waterDensity!= 0.0)
346 << 347 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
347 std::map< G4String,G4double,std::less<G4Stri <<
348 pos1 = lowEnergyLimit.find(particleName); <<
349 if (pos1 != lowEnergyLimit.end()) <<
350 { 348 {
351 lowLim = pos1->second; << 349 const G4String& particleName = particleDefinition->GetParticleName();
352 } <<
353 350
354 std::map< G4String,G4double,std::less<G4Stri << 351 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
355 pos2 = highEnergyLimit.find(particleName); << 352 pos1 = lowEnergyLimit.find(particleName);
356 if (pos2 != highEnergyLimit.end()) << 353 if (pos1 != lowEnergyLimit.end())
357 { << 354 {
358 highLim = pos2->second; << 355 lowLim = pos1->second;
359 } << 356 }
360 357
361 if (ekin >= lowLim && ekin <= highLim) << 358 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
362 { << 359 pos2 = highEnergyLimit.find(particleName);
363 std::map< G4String,G4DNACrossSectionDataSe << 360 if (pos2 != highEnergyLimit.end())
364 pos = tableData.find(particleName); << 361 {
>> 362 highLim = pos2->second;
>> 363 }
365 364
366 if (pos != tableData.end()) << 365 if (ekin >= lowLim && ekin < highLim)
367 { 366 {
368 G4DNACrossSectionDataSet* table = pos->s << 367 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
369 if (table != nullptr) << 368 pos = tableData.find(particleName);
>> 369
>> 370 if (pos != tableData.end())
370 { 371 {
371 sigma = table->FindValue(ekin); << 372 G4DNACrossSectionDataSet* table = pos->second;
>> 373 if (table != 0)
>> 374 {
>> 375 sigma = table->FindValue(ekin);
372 376
373 // ICRU49 electronic SP scaling - ZF, << 377 // ICRU49 electronic SP scaling - ZF, SI
374 378
375 if (particleDefinition == G4Proton::Pr << 379 if (particleDefinition == G4Proton::ProtonDefinition() && ekin < 70*MeV && spScaling)
376 { << 380 {
377 G4double A = 1.39241700556072800000E- << 381 G4double A = 1.39241700556072800000E-009 ;
378 G4double B = -8.52610412942622630000E << 382 G4double B = -8.52610412942622630000E-002 ;
379 sigma = sigma * G4Exp(A*(ekin/eV)+B); << 383 sigma = sigma * G4Exp(A*(ekin/eV)+B);
380 } << 384 }
381 // << 385 //
382 386
>> 387 }
>> 388 }
>> 389 else
>> 390 {
>> 391 G4Exception("G4DNABornIonisationModel1::CrossSectionPerVolume","em0002",
>> 392 FatalException,"Model not applicable to particle type.");
383 } 393 }
384 } 394 }
385 else << 395
>> 396 if (verboseLevel > 2)
386 { 397 {
387 G4Exception("G4DNABornIonisationModel1:: << 398 G4cout << "__________________________________" << G4endl;
388 FatalException,"Model not applicable << 399 G4cout << "G4DNABornIonisationModel1 - XS INFO START" << G4endl;
>> 400 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
>> 401 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
>> 402 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
>> 403 G4cout << "G4DNABornIonisationModel1 - XS INFO END" << G4endl;
389 } 404 }
390 } << 405 } // if (waterMaterial)
391 <<
392 if (verboseLevel > 2) <<
393 { <<
394 G4cout << "_______________________________ <<
395 G4cout << "G4DNABornIonisationModel1 - XS <<
396 G4cout << "Kinetic energy(eV)=" << ekin/eV <<
397 G4cout << "Cross section per water molecul <<
398 G4cout << "Cross section per water molecul <<
399 G4cout << "G4DNABornIonisationModel1 - XS <<
400 } <<
401 406
402 return sigma*waterDensity; 407 return sigma*waterDensity;
403 } 408 }
404 409
405 //....oooOO0OOooo........oooOO0OOooo........oo 410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406 411
407 void G4DNABornIonisationModel1::SampleSecondar 412 void G4DNABornIonisationModel1::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
408 413 const G4MaterialCutsCouple* couple,
409 414 const G4DynamicParticle* particle,
410 415 G4double,
411 416 G4double)
412 { 417 {
413 418
414 if (verboseLevel > 3) 419 if (verboseLevel > 3)
415 { 420 {
416 G4cout << "Calling SampleSecondaries() of 421 G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel1"
417 << G4endl; 422 << G4endl;
418 } 423 }
419 424
420 G4double lowLim = 0; 425 G4double lowLim = 0;
421 G4double highLim = 0; 426 G4double highLim = 0;
422 427
423 G4double k = particle->GetKineticEnergy(); 428 G4double k = particle->GetKineticEnergy();
424 429
425 const G4String& particleName = particle->Get 430 const G4String& particleName = particle->GetDefinition()->GetParticleName();
426 431
427 std::map< G4String,G4double,std::less<G4Stri 432 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
428 pos1 = lowEnergyLimit.find(particleName); 433 pos1 = lowEnergyLimit.find(particleName);
429 434
430 if (pos1 != lowEnergyLimit.end()) 435 if (pos1 != lowEnergyLimit.end())
431 { 436 {
432 lowLim = pos1->second; 437 lowLim = pos1->second;
433 } 438 }
434 439
435 std::map< G4String,G4double,std::less<G4Stri 440 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
436 pos2 = highEnergyLimit.find(particleName); 441 pos2 = highEnergyLimit.find(particleName);
437 442
438 if (pos2 != highEnergyLimit.end()) 443 if (pos2 != highEnergyLimit.end())
439 { 444 {
440 highLim = pos2->second; 445 highLim = pos2->second;
441 } 446 }
442 447
443 if (k >= lowLim && k <= highLim) << 448 if (k >= lowLim && k < highLim)
444 { 449 {
445 G4ParticleMomentum primaryDirection = part 450 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
446 G4double particleMass = particle->GetDefin 451 G4double particleMass = particle->GetDefinition()->GetPDGMass();
447 G4double totalEnergy = k + particleMass; 452 G4double totalEnergy = k + particleMass;
448 G4double pSquare = k * (totalEnergy + part 453 G4double pSquare = k * (totalEnergy + particleMass);
449 G4double totalMomentum = std::sqrt(pSquare 454 G4double totalMomentum = std::sqrt(pSquare);
450 455
451 G4int ionizationShell = 0; 456 G4int ionizationShell = 0;
452 457
453 if (!fasterCode) ionizationShell = RandomS 458 if (!fasterCode) ionizationShell = RandomSelect(k,particleName);
454 459
455 // SI: The following protection is necessa 460 // SI: The following protection is necessary to avoid infinite loops :
456 // sigmadiff_ionisation_e_born.dat has no 461 // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
457 // sigmadiff_cumulated_ionisation_e_born. 462 // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
458 // this is due to the fact that the max a 463 // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
459 // strictly above this value have non zer 464 // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
460 465
461 if (fasterCode) 466 if (fasterCode)
462 do 467 do
463 { 468 {
464 ionizationShell = RandomSelect(k,particl 469 ionizationShell = RandomSelect(k,particleName);
465 } while (k<19*eV && ionizationShell==2 && << 470 }while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());
466 471
467 G4double bindingEnergy = 0; 472 G4double bindingEnergy = 0;
468 bindingEnergy = waterStructure.IonisationE 473 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
469 474
470 // SI: additional protection if tcs interp 475 // SI: additional protection if tcs interpolation method is modified
471 if (k<bindingEnergy) return; 476 if (k<bindingEnergy) return;
472 // 477 //
473 478
474 G4double secondaryKinetic=-1000*eV; 479 G4double secondaryKinetic=-1000*eV;
475 480
476 if (!fasterCode) << 481 if (fasterCode == false)
477 { 482 {
478 secondaryKinetic = RandomizeEjectedElect 483 secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
479 } 484 }
480 else 485 else
481 { 486 {
482 secondaryKinetic = RandomizeEjectedElect 487 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
483 } 488 }
484 // 489 //
485 490
486 G4int Z = 8; 491 G4int Z = 8;
487 492
488 G4ThreeVector deltaDirection = 493 G4ThreeVector deltaDirection =
489 GetAngularDistribution()->SampleDirectionF 494 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
490 Z, ionizationShell, 495 Z, ionizationShell,
491 couple->GetMaterial()); 496 couple->GetMaterial());
492 497
493 if (secondaryKinetic>0) 498 if (secondaryKinetic>0)
494 { 499 {
495 auto dp = new G4DynamicParticle (G4Elec << 500 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
496 fvect->push_back(dp); 501 fvect->push_back(dp);
497 } 502 }
498 503
499 if (particle->GetDefinition() == G4Electro 504 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
500 { 505 {
501 G4double deltaTotalMomentum = std::sqrt( 506 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
502 507
503 G4double finalPx = totalMomentum*primary 508 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
504 G4double finalPy = totalMomentum*primary 509 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
505 G4double finalPz = totalMomentum*primary 510 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
506 G4double finalMomentum = std::sqrt(final 511 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
507 finalPx /= finalMomentum; 512 finalPx /= finalMomentum;
508 finalPy /= finalMomentum; 513 finalPy /= finalMomentum;
509 finalPz /= finalMomentum; 514 finalPz /= finalMomentum;
510 515
511 G4ThreeVector direction; 516 G4ThreeVector direction;
512 direction.set(finalPx,finalPy,finalPz); 517 direction.set(finalPx,finalPy,finalPz);
513 518
514 fParticleChangeForGamma->ProposeMomentum 519 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
515 } 520 }
516 521
517 else fParticleChangeForGamma->ProposeMomen 522 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
518 523
519 // AM: sample deexcitation 524 // AM: sample deexcitation
520 // here we assume that H_{2}O electronic l 525 // here we assume that H_{2}O electronic levels are the same as Oxygen.
521 // this can be considered true with a roug 526 // this can be considered true with a rough 10% error in energy on K-shell,
522 527
523 std::size_t secNumberInit = 0;// need to k << 528 G4int secNumberInit = 0;// need to know at a certain point the energy of secondaries
524 std::size_t secNumberFinal = 0;// So I'll << 529 G4int secNumberFinal = 0;// So I'll make the diference and then sum the energies
525 530
526 G4double scatteredEnergy = k-bindingEnergy << 531 if(fAtomDeexcitation)
527 <<
528 // SI: only atomic deexcitation from K she <<
529 if((fAtomDeexcitation != nullptr) && ioniz <<
530 { 532 {
531 const G4AtomicShell* shell = << 533 G4AtomicShellEnumerator as = fKShell;
532 fAtomDeexcitation->GetAtomicShell(Z, G <<
533 secNumberInit = fvect->size(); <<
534 fAtomDeexcitation->GenerateParticles(fve <<
535 secNumberFinal = fvect->size(); <<
536 <<
537 //TEST <<
538 //G4cout << "ionizationShell=" << ioniza <<
539 //G4cout << "bindingEnergy=" << bindingE <<
540 534
541 if(secNumberFinal > secNumberInit) << 535 if (ionizationShell <5 && ionizationShell >1)
542 { 536 {
543 for (std::size_t i=secNumberInit; i<secNumbe << 537 as = G4AtomicShellEnumerator(4-ionizationShell);
544 { << 538 }
545 //Check if there is enough residual << 539 else if (ionizationShell <2)
546 if (bindingEnergy >= ((*fvect)[i])-> << 540 {
547 { << 541 as = G4AtomicShellEnumerator(3);
548 //Ok, this is a valid secondary: <<
549 bindingEnergy -= ((*fvect)[i])->GetKine <<
550 //G4cout << "--deex nrj=" << ((*fvect)[ <<
551 //<< G4endl; <<
552 } <<
553 else <<
554 { <<
555 //Invalid secondary: not enough energy <<
556 //Keep its energy in the local deposit <<
557 delete (*fvect)[i]; <<
558 (*fvect)[i]=nullptr; <<
559 } <<
560 } <<
561 } 542 }
562 543
563 //TEST << 544 // FOR DEBUG ONLY
564 //G4cout << "k=" << k/eV<< G4endl; << 545 // if (ionizationShell == 4) {
565 //G4cout << "secondaryKinetic=" << seconda << 546 //
566 //G4cout << "scatteredEnergy=" << scattere << 547 // G4cout << "Z: " << Z << " as: " << as
567 //G4cout << "deposited energy=" << binding << 548 // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
568 // << 549 // G4cout << "Press <Enter> key to continue..." << G4endl;
>> 550 // G4cin.ignore();
>> 551 // }
569 552
>> 553 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
>> 554 secNumberInit = fvect->size();
>> 555 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
>> 556 secNumberFinal = fvect->size();
>> 557 }
>> 558
>> 559 // Note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
>> 560
>> 561 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
>> 562 G4double deexSecEnergy = 0;
>> 563 for (G4int j=secNumberInit; j < secNumberFinal; j++)
>> 564 {
>> 565 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
570 } 566 }
571 567
572 //This should never happen <<
573 if(bindingEnergy < 0.0) <<
574 G4Exception("G4DNABornIonisatioModel1::Sa <<
575 "em2050",FatalException,"Nega <<
576 <<
577 //bindingEnergy has been decreased <<
578 //by the amount of energy taken away by de <<
579 if (!statCode) 568 if (!statCode)
580 { 569 {
581 fParticleChangeForGamma->SetProposedKine 570 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
582 fParticleChangeForGamma->ProposeLocalEne << 571 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
583 } 572 }
584 else 573 else
585 { 574 {
586 fParticleChangeForGamma->SetProposedKine 575 fParticleChangeForGamma->SetProposedKineticEnergy(k);
587 fParticleChangeForGamma->ProposeLocalEne 576 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
588 } 577 }
589 578
590 // TEST ////////////////////////// 579 // TEST //////////////////////////
591 // if (secondaryKinetic<0) abort(); 580 // if (secondaryKinetic<0) abort();
592 // if (scatteredEnergy<0) abort(); 581 // if (scatteredEnergy<0) abort();
593 // if (k-scatteredEnergy-secondaryKinetic- 582 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
594 // if (k-scatteredEnergy<0) abort(); 583 // if (k-scatteredEnergy<0) abort();
595 ///////////////////////////////// 584 /////////////////////////////////
596 585
597 const G4Track * theIncomingTrack = fPartic 586 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
598 G4DNAChemistryManager::Instance()->CreateW 587 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
599 ionizationShell, 588 ionizationShell,
600 theIncomingTrack); 589 theIncomingTrack);
601 } 590 }
602 } 591 }
603 592
604 //....oooOO0OOooo........oooOO0OOooo........oo 593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
605 594
606 G4double G4DNABornIonisationModel1::RandomizeE 595 G4double G4DNABornIonisationModel1::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
607 596 G4double k,
608 597 G4int shell)
609 { 598 {
610 // G4cout << "*** SLOW computation for " << 599 // G4cout << "*** SLOW computation for " << " " << particleDefinition->GetParticleName() << G4endl;
611 600
612 if (particleDefinition == G4Electron::Electr 601 if (particleDefinition == G4Electron::ElectronDefinition())
613 { 602 {
614 G4double maximumEnergyTransfer = 0.; 603 G4double maximumEnergyTransfer = 0.;
615 if ((k + waterStructure.IonisationEnergy(s 604 if ((k + waterStructure.IonisationEnergy(shell)) / 2. > k)
616 maximumEnergyTransfer = k; 605 maximumEnergyTransfer = k;
617 else 606 else
618 maximumEnergyTransfer = (k + waterStruct 607 maximumEnergyTransfer = (k + waterStructure.IonisationEnergy(shell)) / 2.;
619 608
620 // SI : original method 609 // SI : original method
621 /* 610 /*
622 G4double crossSectionMaximum = 0.; << 611 G4double crossSectionMaximum = 0.;
623 for(G4double value=waterStructure.Ionisati << 612 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
624 { << 613 {
625 G4double differentialCrossSection = Diffe 614 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
626 if(differentialCrossSection >= crossSecti 615 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
627 } << 616 }
628 */ << 617 */
629 618
630 // SI : alternative method 619 // SI : alternative method
631 G4double crossSectionMaximum = 0.; 620 G4double crossSectionMaximum = 0.;
632 621
633 G4double minEnergy = waterStructure.Ionisa 622 G4double minEnergy = waterStructure.IonisationEnergy(shell);
634 G4double maxEnergy = maximumEnergyTransfer 623 G4double maxEnergy = maximumEnergyTransfer;
635 G4int nEnergySteps = 50; 624 G4int nEnergySteps = 50;
636 625
637 G4double value(minEnergy); 626 G4double value(minEnergy);
638 G4double stpEnergy(std::pow(maxEnergy / va 627 G4double stpEnergy(std::pow(maxEnergy / value,
639 1. / static_ca 628 1. / static_cast<G4double>(nEnergySteps - 1)));
640 G4int step(nEnergySteps); 629 G4int step(nEnergySteps);
641 while (step > 0) 630 while (step > 0)
642 { 631 {
643 step--; 632 step--;
644 G4double differentialCrossSection = 633 G4double differentialCrossSection =
645 DifferentialCrossSection(particleDef 634 DifferentialCrossSection(particleDefinition,
646 k / eV, 635 k / eV,
647 value / eV, 636 value / eV,
648 shell); 637 shell);
649 if (differentialCrossSection >= crossSec 638 if (differentialCrossSection >= crossSectionMaximum)
650 crossSectionMaximum = differentialCros 639 crossSectionMaximum = differentialCrossSection;
651 value *= stpEnergy; 640 value *= stpEnergy;
652 } 641 }
653 // 642 //
654 643
655 G4double secondaryElectronKineticEnergy = 644 G4double secondaryElectronKineticEnergy = 0.;
656 do 645 do
657 { 646 {
658 secondaryElectronKineticEnergy = G4Unifo 647 secondaryElectronKineticEnergy = G4UniformRand()* (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
659 } while(G4UniformRand()*crossSectionMaximu << 648 }while(G4UniformRand()*crossSectionMaximum >
660 DifferentialCrossSection(particleDefin 649 DifferentialCrossSection(particleDefinition, k/eV,
661 (secondaryElectronKineticEnergy+wa 650 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
662 651
663 return secondaryElectronKineticEnergy; 652 return secondaryElectronKineticEnergy;
664 653
665 } 654 }
666 655
667 if (particleDefinition == G4Proton::ProtonDe << 656 else if (particleDefinition == G4Proton::ProtonDefinition())
668 { 657 {
669 G4double maximumKineticEnergyTransfer = 4. 658 G4double maximumKineticEnergyTransfer = 4.
670 * (electron_mass_c2 / proton_mass_c2) 659 * (electron_mass_c2 / proton_mass_c2) * k;
671 660
672 G4double crossSectionMaximum = 0.; 661 G4double crossSectionMaximum = 0.;
673 for (G4double value = waterStructure.Ionis 662 for (G4double value = waterStructure.IonisationEnergy(shell);
674 value <= 4. * waterStructure.Ionisatio 663 value <= 4. * waterStructure.IonisationEnergy(shell); value += 0.1 * eV)
675 { 664 {
676 G4double differentialCrossSection = 665 G4double differentialCrossSection =
677 DifferentialCrossSection(particleDef 666 DifferentialCrossSection(particleDefinition,
678 k / eV, 667 k / eV,
679 value / eV, 668 value / eV,
680 shell); 669 shell);
681 if (differentialCrossSection >= crossSec 670 if (differentialCrossSection >= crossSectionMaximum)
682 crossSectionMaximum = differentialCros 671 crossSectionMaximum = differentialCrossSection;
683 } 672 }
684 673
685 G4double secondaryElectronKineticEnergy = 674 G4double secondaryElectronKineticEnergy = 0.;
686 do 675 do
687 { 676 {
688 secondaryElectronKineticEnergy = G4Unifo 677 secondaryElectronKineticEnergy = G4UniformRand()* maximumKineticEnergyTransfer;
689 } while(G4UniformRand()*crossSectionMaximu << 678 }while(G4UniformRand()*crossSectionMaximum >=
690 DifferentialCrossSection(particleDefin 679 DifferentialCrossSection(particleDefinition, k/eV,
691 (secondaryElectronKineticEnergy+wa 680 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
692 681
693 return secondaryElectronKineticEnergy; 682 return secondaryElectronKineticEnergy;
694 } 683 }
695 684
696 return 0; 685 return 0;
697 } 686 }
698 687
699 //....oooOO0OOooo........oooOO0OOooo........oo 688 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
700 689
701 // The following section is not used anymore b 690 // The following section is not used anymore but is kept for memory
702 // GetAngularDistribution()->SampleDirectionFo 691 // GetAngularDistribution()->SampleDirectionForShell is used instead
703 692
704 /* 693 /*
705 void G4DNABornIonisationModel1::RandomizeEjec 694 void G4DNABornIonisationModel1::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
706 G4double k, 695 G4double k,
707 G4double secKinetic, 696 G4double secKinetic,
708 G4double & cosTheta, 697 G4double & cosTheta,
709 G4double & phi ) 698 G4double & phi )
710 { 699 {
711 if (particleDefinition == G4Electron::Electro 700 if (particleDefinition == G4Electron::ElectronDefinition())
712 { 701 {
713 phi = twopi * G4UniformRand(); 702 phi = twopi * G4UniformRand();
714 if (secKinetic < 50.*eV) cosTheta = (2.*G4Uni 703 if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
715 else if (secKinetic <= 200.*eV) 704 else if (secKinetic <= 200.*eV)
716 { 705 {
717 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4 706 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
718 else cosTheta = G4UniformRand()*(std::sqrt(2. 707 else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
719 } 708 }
720 else 709 else
721 { 710 {
722 G4double sin2O = (1.-secKinetic/k) / (1.+secK 711 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
723 cosTheta = std::sqrt(1.-sin2O); 712 cosTheta = std::sqrt(1.-sin2O);
724 } 713 }
725 } 714 }
726 715
727 else if (particleDefinition == G4Proton::Prot 716 else if (particleDefinition == G4Proton::ProtonDefinition())
728 { 717 {
729 G4double maxSecKinetic = 4.* (electron_mass_c 718 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
730 phi = twopi * G4UniformRand(); 719 phi = twopi * G4UniformRand();
731 720
732 // cosTheta = std::sqrt(secKinetic / maxSecKi 721 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
733 722
734 // Restriction below 100 eV from Emfietzoglou 723 // Restriction below 100 eV from Emfietzoglou (2000)
735 724
736 if (secKinetic>100*eV) cosTheta = std::sqrt(s 725 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
737 else cosTheta = (2.*G4UniformRand())-1.; 726 else cosTheta = (2.*G4UniformRand())-1.;
738 727
739 } 728 }
740 } 729 }
741 */ 730 */
742 731
743 //....oooOO0OOooo........oooOO0OOooo........oo 732 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
744 G4double G4DNABornIonisationModel1::Differenti 733 G4double G4DNABornIonisationModel1::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
745 734 G4double k,
746 735 G4double energyTransfer,
747 736 G4int ionizationLevelIndex)
748 { 737 {
749 G4double sigma = 0.; 738 G4double sigma = 0.;
750 739
751 if (energyTransfer >= waterStructure.Ionisat 740 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
752 { 741 {
753 G4double valueT1 = 0; 742 G4double valueT1 = 0;
754 G4double valueT2 = 0; 743 G4double valueT2 = 0;
755 G4double valueE21 = 0; 744 G4double valueE21 = 0;
756 G4double valueE22 = 0; 745 G4double valueE22 = 0;
757 G4double valueE12 = 0; 746 G4double valueE12 = 0;
758 G4double valueE11 = 0; 747 G4double valueE11 = 0;
759 748
760 G4double xs11 = 0; 749 G4double xs11 = 0;
761 G4double xs12 = 0; 750 G4double xs12 = 0;
762 G4double xs21 = 0; 751 G4double xs21 = 0;
763 G4double xs22 = 0; 752 G4double xs22 = 0;
764 753
765 if (particleDefinition == G4Electron::Elec 754 if (particleDefinition == G4Electron::ElectronDefinition())
766 { 755 {
767 <<
768 // Protection against out of boundary ac <<
769 if (k==eTdummyVec.back()) k=k*(1.-1e-12) <<
770 // <<
771 <<
772 // k should be in eV and energy transfer 756 // k should be in eV and energy transfer eV also
773 757
774 auto t2 = std::upper_bound(eTdummyVec.be << 758 std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
775 759 eTdummyVec.end(),
776 760 k);
777 761
778 auto t1 = t2 - 1; << 762 std::vector<G4double>::iterator t1 = t2 - 1;
779 763
780 // SI : the following condition avoids s 764 // SI : the following condition avoids situations where energyTransfer >last vector element
781 if (energyTransfer <= eVecm[(*t1)].back( 765 if (energyTransfer <= eVecm[(*t1)].back()
782 && energyTransfer <= eVecm[(*t2)].ba 766 && energyTransfer <= eVecm[(*t2)].back())
783 { 767 {
784 auto e12 = << 768 std::vector<G4double>::iterator e12 =
785 std::upper_bound(eVecm[(*t1)].begi 769 std::upper_bound(eVecm[(*t1)].begin(),
786 eVecm[(*t1)].end( 770 eVecm[(*t1)].end(),
787 energyTransfer); 771 energyTransfer);
788 auto e11 = e12 - 1; << 772 std::vector<G4double>::iterator e11 = e12 - 1;
789 773
790 auto e22 = << 774 std::vector<G4double>::iterator e22 =
791 std::upper_bound(eVecm[(*t2)].begi 775 std::upper_bound(eVecm[(*t2)].begin(),
792 eVecm[(*t2)].end( 776 eVecm[(*t2)].end(),
793 energyTransfer); 777 energyTransfer);
794 auto e21 = e22 - 1; << 778 std::vector<G4double>::iterator e21 = e22 - 1;
795 779
796 valueT1 = *t1; 780 valueT1 = *t1;
797 valueT2 = *t2; 781 valueT2 = *t2;
798 valueE21 = *e21; 782 valueE21 = *e21;
799 valueE22 = *e22; 783 valueE22 = *e22;
800 valueE12 = *e12; 784 valueE12 = *e12;
801 valueE11 = *e11; 785 valueE11 = *e11;
802 786
803 xs11 = eDiffCrossSectionData[ionizatio 787 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
804 xs12 = eDiffCrossSectionData[ionizatio 788 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
805 xs21 = eDiffCrossSectionData[ionizatio 789 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
806 xs22 = eDiffCrossSectionData[ionizatio 790 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
807 791
808 } 792 }
809 793
810 } 794 }
811 795
812 if (particleDefinition == G4Proton::Proton 796 if (particleDefinition == G4Proton::ProtonDefinition())
813 { 797 {
814 // Protection against out of boundary ac <<
815 if (k==pTdummyVec.back()) k=k*(1.-1e-12) <<
816 // <<
817 <<
818 // k should be in eV and energy transfer 798 // k should be in eV and energy transfer eV also
819 auto t2 = std::upper_bound(pTdummyVec.be << 799 std::vector<G4double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),
820 800 pTdummyVec.end(),
821 801 k);
822 auto t1 = t2 - 1; << 802 std::vector<G4double>::iterator t1 = t2 - 1;
823 803
824 auto e12 = std::upper_bound(pVecm[(*t1)] << 804 std::vector<G4double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),
825 805 pVecm[(*t1)].end(),
826 806 energyTransfer);
827 auto e11 = e12 - 1; << 807 std::vector<G4double>::iterator e11 = e12 - 1;
828 808
829 auto e22 = std::upper_bound(pVecm[(*t2)] << 809 std::vector<G4double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),
830 810 pVecm[(*t2)].end(),
831 811 energyTransfer);
832 auto e21 = e22 - 1; << 812 std::vector<G4double>::iterator e21 = e22 - 1;
833 813
834 valueT1 = *t1; 814 valueT1 = *t1;
835 valueT2 = *t2; 815 valueT2 = *t2;
836 valueE21 = *e21; 816 valueE21 = *e21;
837 valueE22 = *e22; 817 valueE22 = *e22;
838 valueE12 = *e12; 818 valueE12 = *e12;
839 valueE11 = *e11; 819 valueE11 = *e11;
840 820
841 xs11 = pDiffCrossSectionData[ionizationL 821 xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
842 xs12 = pDiffCrossSectionData[ionizationL 822 xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
843 xs21 = pDiffCrossSectionData[ionizationL 823 xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
844 xs22 = pDiffCrossSectionData[ionizationL 824 xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
845 825
846 } 826 }
847 827
848 G4double xsProduct = xs11 * xs12 * xs21 * 828 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
849 if (xsProduct != 0.) 829 if (xsProduct != 0.)
850 { 830 {
851 sigma = QuadInterpolator(valueE11, 831 sigma = QuadInterpolator(valueE11,
852 valueE12, 832 valueE12,
853 valueE21, 833 valueE21,
854 valueE22, 834 valueE22,
855 xs11, 835 xs11,
856 xs12, 836 xs12,
857 xs21, 837 xs21,
858 xs22, 838 xs22,
859 valueT1, 839 valueT1,
860 valueT2, 840 valueT2,
861 k, 841 k,
862 energyTransfer) 842 energyTransfer);
863 } 843 }
864 } 844 }
865 845
866 return sigma; 846 return sigma;
867 } 847 }
868 848
869 //....oooOO0OOooo........oooOO0OOooo........oo 849 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
870 850
871 G4double G4DNABornIonisationModel1::Interpolat 851 G4double G4DNABornIonisationModel1::Interpolate(G4double e1,
872 852 G4double e2,
873 853 G4double e,
874 854 G4double xs1,
875 855 G4double xs2)
876 { 856 {
877 G4double value = 0.; 857 G4double value = 0.;
878 858
879 // Log-log interpolation by default 859 // Log-log interpolation by default
880 860
881 if (e1 != 0 && e2 != 0 && (std::log10(e2) - 861 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
882 && !fasterCode) 862 && !fasterCode)
883 { 863 {
884 G4double a = (std::log10(xs2) - std::log10 864 G4double a = (std::log10(xs2) - std::log10(xs1))
885 / (std::log10(e2) - std::log10(e1)); 865 / (std::log10(e2) - std::log10(e1));
886 G4double b = std::log10(xs2) - a * std::lo 866 G4double b = std::log10(xs2) - a * std::log10(e2);
887 G4double sigma = a * std::log10(e) + b; 867 G4double sigma = a * std::log10(e) + b;
888 value = (std::pow(10., sigma)); 868 value = (std::pow(10., sigma));
889 } 869 }
890 870
891 // Switch to lin-lin interpolation 871 // Switch to lin-lin interpolation
892 /* 872 /*
893 if ((e2-e1)!=0) 873 if ((e2-e1)!=0)
894 { 874 {
895 G4double d1 = xs1; 875 G4double d1 = xs1;
896 G4double d2 = xs2; 876 G4double d2 = xs2;
897 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1) 877 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
898 } 878 }
899 */ 879 */
900 880
901 // Switch to log-lin interpolation for faste 881 // Switch to log-lin interpolation for faster code
902 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 & 882 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
903 { 883 {
904 G4double d1 = std::log10(xs1); 884 G4double d1 = std::log10(xs1);
905 G4double d2 = std::log10(xs2); 885 G4double d2 = std::log10(xs2);
906 value = std::pow(10., (d1 + (d2 - d1) * (e 886 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
907 } 887 }
908 888
909 // Switch to lin-lin interpolation for faste 889 // Switch to lin-lin interpolation for faster code
910 // in case one of xs1 or xs2 (=cum proba) va 890 // in case one of xs1 or xs2 (=cum proba) value is zero
911 891
912 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) 892 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
913 { 893 {
914 G4double d1 = xs1; 894 G4double d1 = xs1;
915 G4double d2 = xs2; 895 G4double d2 = xs2;
916 value = (d1 + (d2 - d1) * (e - e1) / (e2 - 896 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
917 } 897 }
918 898
919 /* 899 /*
920 G4cout 900 G4cout
921 << e1 << " " 901 << e1 << " "
922 << e2 << " " 902 << e2 << " "
923 << e << " " 903 << e << " "
924 << xs1 << " " 904 << xs1 << " "
925 << xs2 << " " 905 << xs2 << " "
926 << value 906 << value
927 << G4endl; 907 << G4endl;
928 */ 908 */
929 909
930 return value; 910 return value;
931 } 911 }
932 912
933 //....oooOO0OOooo........oooOO0OOooo........oo 913 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
934 914
935 G4double G4DNABornIonisationModel1::QuadInterp 915 G4double G4DNABornIonisationModel1::QuadInterpolator(G4double e11,
936 916 G4double e12,
937 917 G4double e21,
938 918 G4double e22,
939 919 G4double xs11,
940 920 G4double xs12,
941 921 G4double xs21,
942 922 G4double xs22,
943 923 G4double t1,
944 924 G4double t2,
945 925 G4double t,
946 926 G4double e)
947 { 927 {
948 G4double interpolatedvalue1 = Interpolate(e1 928 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
949 G4double interpolatedvalue2 = Interpolate(e2 929 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
950 G4double value = Interpolate(t1, 930 G4double value = Interpolate(t1,
951 t2, 931 t2,
952 t, 932 t,
953 interpolatedval 933 interpolatedvalue1,
954 interpolatedval 934 interpolatedvalue2);
955 935
956 return value; 936 return value;
957 } 937 }
958 938
959 G4double G4DNABornIonisationModel1::GetPartial 939 G4double G4DNABornIonisationModel1::GetPartialCrossSection(const G4Material* /*material*/,
960 940 G4int level,
961 941 const G4ParticleDefinition* particle,
962 942 G4double kineticEnergy)
963 { 943 {
964 std::map<G4String, G4DNACrossSectionDataSet* 944 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
965 pos = tableData.find(particle->GetParticleNa 945 pos = tableData.find(particle->GetParticleName());
966 946
967 if (pos != tableData.end()) 947 if (pos != tableData.end())
968 { 948 {
969 G4DNACrossSectionDataSet* table = pos->sec 949 G4DNACrossSectionDataSet* table = pos->second;
970 return table->GetComponent(level)->FindVal 950 return table->GetComponent(level)->FindValue(kineticEnergy);
971 } 951 }
972 952
973 return 0; 953 return 0;
974 } 954 }
975 955
976 //....oooOO0OOooo........oooOO0OOooo........oo 956 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
977 957
978 G4int G4DNABornIonisationModel1::RandomSelect( 958 G4int G4DNABornIonisationModel1::RandomSelect(G4double k,
979 959 const G4String& particle)
980 { 960 {
981 G4int level = 0; 961 G4int level = 0;
982 962
983 std::map<G4String, G4DNACrossSectionDataSet* 963 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
984 pos = tableData.find(particle); 964 pos = tableData.find(particle);
985 965
986 if (pos != tableData.end()) 966 if (pos != tableData.end())
987 { 967 {
988 G4DNACrossSectionDataSet* table = pos->sec 968 G4DNACrossSectionDataSet* table = pos->second;
989 969
990 if (table != nullptr) << 970 if (table != 0)
991 { 971 {
992 auto valuesBuffer = new G4double[table- << 972 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
993 const auto n = (G4int)table->NumberOfCo << 973 const size_t n(table->NumberOfComponents());
994 G4int i(n); << 974 size_t i(n);
995 G4double value = 0.; 975 G4double value = 0.;
996 976
997 while (i > 0) 977 while (i > 0)
998 { 978 {
999 i--; 979 i--;
1000 valuesBuffer[i] = table->GetComponent 980 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1001 value += valuesBuffer[i]; 981 value += valuesBuffer[i];
1002 } 982 }
1003 983
1004 value *= G4UniformRand(); 984 value *= G4UniformRand();
1005 985
1006 i = n; 986 i = n;
1007 987
1008 while (i > 0) 988 while (i > 0)
1009 { 989 {
1010 i--; 990 i--;
1011 991
1012 if (valuesBuffer[i] > value) 992 if (valuesBuffer[i] > value)
1013 { 993 {
1014 delete[] valuesBuffer; 994 delete[] valuesBuffer;
1015 return i; 995 return i;
1016 } 996 }
1017 value -= valuesBuffer[i]; 997 value -= valuesBuffer[i];
1018 } 998 }
1019 999
1020 << 1000 if (valuesBuffer)
1021 delete[] valuesBuffer; 1001 delete[] valuesBuffer;
1022 1002
1023 } 1003 }
1024 } else 1004 } else
1025 { 1005 {
1026 G4Exception("G4DNABornIonisationModel1::R 1006 G4Exception("G4DNABornIonisationModel1::RandomSelect",
1027 "em0002", 1007 "em0002",
1028 FatalException, 1008 FatalException,
1029 "Model not applicable to part 1009 "Model not applicable to particle type.");
1030 } 1010 }
1031 1011
1032 return level; 1012 return level;
1033 } 1013 }
1034 1014
1035 //....oooOO0OOooo........oooOO0OOooo........o 1015 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1036 1016
1037 G4double G4DNABornIonisationModel1::Randomize 1017 G4double G4DNABornIonisationModel1::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
1038 1018 G4double k,
1039 1019 G4int shell)
1040 { 1020 {
1041 //G4cout << "*** FAST computation for " << 1021 //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
1042 1022
1043 G4double secondaryElectronKineticEnergy = 0 1023 G4double secondaryElectronKineticEnergy = 0.;
1044 1024
1045 G4double random = G4UniformRand(); 1025 G4double random = G4UniformRand();
1046 1026
1047 secondaryElectronKineticEnergy = Transfered 1027 secondaryElectronKineticEnergy = TransferedEnergy(particleDefinition,
1048 1028 k / eV,
1049 1029 shell,
1050 1030 random) * eV
1051 - waterStructure.IonisationEnergy(shell 1031 - waterStructure.IonisationEnergy(shell);
1052 1032
1053 //G4cout << RandomTransferedEnergy(particle 1033 //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
1054 <<
1055 if (secondaryElectronKineticEnergy < 0.) 1034 if (secondaryElectronKineticEnergy < 0.)
1056 return 0.; 1035 return 0.;
1057 // 1036 //
1058 1037
1059 return secondaryElectronKineticEnergy; 1038 return secondaryElectronKineticEnergy;
1060 } 1039 }
1061 1040
1062 //....oooOO0OOooo........oooOO0OOooo........o 1041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1063 1042
1064 G4double G4DNABornIonisationModel1::Transfere 1043 G4double G4DNABornIonisationModel1::TransferedEnergy(G4ParticleDefinition* particleDefinition,
1065 1044 G4double k,
1066 1045 G4int ionizationLevelIndex,
1067 1046 G4double random)
1068 { 1047 {
1069 G4double nrj = 0.; 1048 G4double nrj = 0.;
1070 1049
1071 G4double valueK1 = 0; 1050 G4double valueK1 = 0;
1072 G4double valueK2 = 0; 1051 G4double valueK2 = 0;
1073 G4double valuePROB21 = 0; 1052 G4double valuePROB21 = 0;
1074 G4double valuePROB22 = 0; 1053 G4double valuePROB22 = 0;
1075 G4double valuePROB12 = 0; 1054 G4double valuePROB12 = 0;
1076 G4double valuePROB11 = 0; 1055 G4double valuePROB11 = 0;
1077 1056
1078 G4double nrjTransf11 = 0; 1057 G4double nrjTransf11 = 0;
1079 G4double nrjTransf12 = 0; 1058 G4double nrjTransf12 = 0;
1080 G4double nrjTransf21 = 0; 1059 G4double nrjTransf21 = 0;
1081 G4double nrjTransf22 = 0; 1060 G4double nrjTransf22 = 0;
1082 1061
1083 if (particleDefinition == G4Electron::Elect 1062 if (particleDefinition == G4Electron::ElectronDefinition())
1084 { 1063 {
1085 // Protection against out of boundary acc <<
1086 if (k==eTdummyVec.back()) k=k*(1.-1e-12); <<
1087 // <<
1088 <<
1089 // k should be in eV 1064 // k should be in eV
1090 auto k2 = std::upper_bound(eTdummyVec.beg << 1065 std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),
1091 1066 eTdummyVec.end(),
1092 1067 k);
1093 auto k1 = k2 - 1; << 1068 std::vector<G4double>::iterator k1 = k2 - 1;
1094 1069
1095 /* 1070 /*
1096 G4cout << "----> k=" << k 1071 G4cout << "----> k=" << k
1097 << " " << *k1 1072 << " " << *k1
1098 << " " << *k2 1073 << " " << *k2
1099 << " " << random 1074 << " " << random
1100 << " " << ionizationLevelIndex 1075 << " " << ionizationLevelIndex
1101 << " " << eProbaShellMap[ionizationLevel 1076 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1102 << " " << eProbaShellMap[ionizationLevel 1077 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
1103 << G4endl; 1078 << G4endl;
1104 */ 1079 */
1105 1080
1106 // SI : the following condition avoids si 1081 // SI : the following condition avoids situations where random >last vector element
1107 if (random <= eProbaShellMap[ionizationLe 1082 if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1108 && random <= eProbaShellMap[ionizatio 1083 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
1109 { 1084 {
1110 auto prob12 = << 1085 std::vector<G4double>::iterator prob12 =
1111 std::upper_bound(eProbaShellMap[ion 1086 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1112 eProbaShellMap[ion 1087 eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1113 random); 1088 random);
1114 1089
1115 auto prob11 = prob12 - 1; << 1090 std::vector<G4double>::iterator prob11 = prob12 - 1;
1116 1091
1117 auto prob22 = << 1092 std::vector<G4double>::iterator prob22 =
1118 std::upper_bound(eProbaShellMap[ion 1093 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1119 eProbaShellMap[ion 1094 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1120 random); 1095 random);
1121 1096
1122 auto prob21 = prob22 - 1; << 1097 std::vector<G4double>::iterator prob21 = prob22 - 1;
1123 1098
1124 valueK1 = *k1; 1099 valueK1 = *k1;
1125 valueK2 = *k2; 1100 valueK2 = *k2;
1126 valuePROB21 = *prob21; 1101 valuePROB21 = *prob21;
1127 valuePROB22 = *prob22; 1102 valuePROB22 = *prob22;
1128 valuePROB12 = *prob12; 1103 valuePROB12 = *prob12;
1129 valuePROB11 = *prob11; 1104 valuePROB11 = *prob11;
1130 1105
1131 /* 1106 /*
1132 G4cout << " " << random << " " 1107 G4cout << " " << random << " " << valuePROB11 << " "
1133 << valuePROB12 << " " << valuePROB21 < 1108 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1134 */ 1109 */
1135 1110
1136 nrjTransf11 = eNrjTransfData[ionization 1111 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1137 nrjTransf12 = eNrjTransfData[ionization 1112 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1138 nrjTransf21 = eNrjTransfData[ionization 1113 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1139 nrjTransf22 = eNrjTransfData[ionization 1114 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1140 1115
1141 /* 1116 /*
1142 G4cout << " " << ionizationLeve 1117 G4cout << " " << ionizationLevelIndex << " "
1143 << random << " " <<valueK1 << " " << v 1118 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1144 1119
1145 G4cout << " " << random << " " 1120 G4cout << " " << random << " " << nrjTransf11 << " "
1146 << nrjTransf12 << " " << nrjTransf21 < 1121 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1147 */ 1122 */
1148 1123
1149 } 1124 }
1150 <<
1151 // Avoids cases where cum xs is zero for 1125 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1152 if (random > eProbaShellMap[ionizationLev 1126 if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
1153 { 1127 {
1154 auto prob22 = << 1128 std::vector<G4double>::iterator prob22 =
1155 std::upper_bound(eProbaShellMap[ion 1129 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1156 eProbaShellMap[ion 1130 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1157 random); 1131 random);
1158 1132
1159 auto prob21 = prob22 - 1; << 1133 std::vector<G4double>::iterator prob21 = prob22 - 1;
1160 1134
1161 valueK1 = *k1; 1135 valueK1 = *k1;
1162 valueK2 = *k2; 1136 valueK2 = *k2;
1163 valuePROB21 = *prob21; 1137 valuePROB21 = *prob21;
1164 valuePROB22 = *prob22; 1138 valuePROB22 = *prob22;
1165 1139
1166 //G4cout << " " << random << " " 1140 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1167 1141
1168 nrjTransf21 = eNrjTransfData[ionization 1142 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1169 nrjTransf22 = eNrjTransfData[ionization 1143 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1170 1144
1171 G4double interpolatedvalue2 = Interpola 1145 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1172 1146 valuePROB22,
1173 1147 random,
1174 1148 nrjTransf21,
1175 1149 nrjTransf22);
1176 1150
1177 // zeros are explicitly set << 1151 // zeros are explicitely set
1178 1152
1179 G4double value = Interpolate(valueK1, v 1153 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1180 1154
1181 /* 1155 /*
1182 G4cout << " " << ionizationLeve 1156 G4cout << " " << ionizationLevelIndex << " "
1183 << random << " " <<valueK1 << " " << v 1157 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1184 1158
1185 G4cout << " " << random << " " 1159 G4cout << " " << random << " " << nrjTransf11 << " "
1186 << nrjTransf12 << " " << nrjTransf21 < 1160 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1187 1161
1188 G4cout << "ici" << " " << value << G4e 1162 G4cout << "ici" << " " << value << G4endl;
1189 */ 1163 */
1190 1164
1191 return value; 1165 return value;
1192 } 1166 }
1193 } 1167 }
1194 // 1168 //
1195 else if (particleDefinition == G4Proton::Pr 1169 else if (particleDefinition == G4Proton::ProtonDefinition())
1196 { 1170 {
1197 // Protection against out of boundary acc <<
1198 if (k==pTdummyVec.back()) k=k*(1.-1e-12); <<
1199 // <<
1200 <<
1201 // k should be in eV 1171 // k should be in eV
1202 1172
1203 auto k2 = std::upper_bound(pTdummyVec.beg << 1173 std::vector<G4double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),
1204 1174 pTdummyVec.end(),
1205 1175 k);
1206 1176
1207 auto k1 = k2 - 1; << 1177 std::vector<G4double>::iterator k1 = k2 - 1;
1208 1178
1209 /* 1179 /*
1210 G4cout << "----> k=" << k 1180 G4cout << "----> k=" << k
1211 << " " << *k1 1181 << " " << *k1
1212 << " " << *k2 1182 << " " << *k2
1213 << " " << random 1183 << " " << random
1214 << " " << ionizationLevelIndex 1184 << " " << ionizationLevelIndex
1215 << " " << pProbaShellMap[ionizationLevel 1185 << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1216 << " " << pProbaShellMap[ionizationLevel 1186 << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1217 << G4endl; 1187 << G4endl;
1218 */ 1188 */
1219 1189
1220 // SI : the following condition avoids si 1190 // SI : the following condition avoids situations where random > last vector element,
1221 // for eg. when the last element is 1191 // for eg. when the last element is zero
1222 if (random <= pProbaShellMap[ionizationLe 1192 if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1223 && random <= pProbaShellMap[ionizatio 1193 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1224 { 1194 {
1225 auto prob12 = << 1195 std::vector<G4double>::iterator prob12 =
1226 std::upper_bound(pProbaShellMap[ion 1196 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1227 pProbaShellMap[ion 1197 pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1228 random); 1198 random);
1229 1199
1230 auto prob11 = prob12 - 1; << 1200 std::vector<G4double>::iterator prob11 = prob12 - 1;
1231 1201
1232 auto prob22 = << 1202 std::vector<G4double>::iterator prob22 =
1233 std::upper_bound(pProbaShellMap[ion 1203 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1234 pProbaShellMap[ion 1204 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1235 random); 1205 random);
1236 1206
1237 auto prob21 = prob22 - 1; << 1207 std::vector<G4double>::iterator prob21 = prob22 - 1;
1238 1208
1239 valueK1 = *k1; 1209 valueK1 = *k1;
1240 valueK2 = *k2; 1210 valueK2 = *k2;
1241 valuePROB21 = *prob21; 1211 valuePROB21 = *prob21;
1242 valuePROB22 = *prob22; 1212 valuePROB22 = *prob22;
1243 valuePROB12 = *prob12; 1213 valuePROB12 = *prob12;
1244 valuePROB11 = *prob11; 1214 valuePROB11 = *prob11;
1245 1215
1246 /* 1216 /*
1247 G4cout << " " << random << " " 1217 G4cout << " " << random << " " << valuePROB11 << " "
1248 << valuePROB12 << " " << valuePROB21 < 1218 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1249 */ 1219 */
1250 1220
1251 nrjTransf11 = pNrjTransfData[ionization 1221 nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1252 nrjTransf12 = pNrjTransfData[ionization 1222 nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1253 nrjTransf21 = pNrjTransfData[ionization 1223 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1254 nrjTransf22 = pNrjTransfData[ionization 1224 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1255 1225
1256 /* 1226 /*
1257 G4cout << " " << ionizationLeve 1227 G4cout << " " << ionizationLevelIndex << " "
1258 << random << " " <<valueK1 << " " << v 1228 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1259 1229
1260 G4cout << " " << random << " " 1230 G4cout << " " << random << " " << nrjTransf11 << " "
1261 << nrjTransf12 << " " << nrjTransf21 < 1231 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1262 */ 1232 */
1263 } 1233 }
1264 1234
1265 // Avoids cases where cum xs is zero for 1235 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1266 1236
1267 if (random > pProbaShellMap[ionizationLev 1237 if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1268 { 1238 {
1269 auto prob22 = << 1239 std::vector<G4double>::iterator prob22 =
1270 std::upper_bound(pProbaShellMap[ion 1240 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1271 pProbaShellMap[ion 1241 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1272 random); 1242 random);
1273 1243
1274 auto prob21 = prob22 - 1; << 1244 std::vector<G4double>::iterator prob21 = prob22 - 1;
1275 1245
1276 valueK1 = *k1; 1246 valueK1 = *k1;
1277 valueK2 = *k2; 1247 valueK2 = *k2;
1278 valuePROB21 = *prob21; 1248 valuePROB21 = *prob21;
1279 valuePROB22 = *prob22; 1249 valuePROB22 = *prob22;
1280 1250
1281 //G4cout << " " << random << " " 1251 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1282 1252
1283 nrjTransf21 = pNrjTransfData[ionization 1253 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1284 nrjTransf22 = pNrjTransfData[ionization 1254 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1285 1255
1286 G4double interpolatedvalue2 = Interpola 1256 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1287 1257 valuePROB22,
1288 1258 random,
1289 1259 nrjTransf21,
1290 1260 nrjTransf22);
1291 1261
1292 // zeros are explicitly set << 1262 // zeros are explicitely set
1293 1263
1294 G4double value = Interpolate(valueK1, v 1264 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1295 1265
1296 /* 1266 /*
1297 G4cout << " " << ionizationLeve 1267 G4cout << " " << ionizationLevelIndex << " "
1298 << random << " " <<valueK1 << " " << v 1268 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1299 1269
1300 G4cout << " " << random << " " 1270 G4cout << " " << random << " " << nrjTransf11 << " "
1301 << nrjTransf12 << " " << nrjTransf21 < 1271 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1302 1272
1303 G4cout << "ici" << " " << value << G4e 1273 G4cout << "ici" << " " << value << G4endl;
1304 */ 1274 */
1305 1275
1306 return value; 1276 return value;
1307 } 1277 }
1308 } 1278 }
1309 // End electron and proton cases 1279 // End electron and proton cases
1310 1280
1311 G4double nrjTransfProduct = nrjTransf11 * n 1281 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1312 * nrjTransf22; 1282 * nrjTransf22;
1313 //G4cout << "nrjTransfProduct=" << nrjTrans 1283 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1314 1284
1315 if (nrjTransfProduct != 0.) 1285 if (nrjTransfProduct != 0.)
1316 { 1286 {
1317 nrj = QuadInterpolator(valuePROB11, 1287 nrj = QuadInterpolator(valuePROB11,
1318 valuePROB12, 1288 valuePROB12,
1319 valuePROB21, 1289 valuePROB21,
1320 valuePROB22, 1290 valuePROB22,
1321 nrjTransf11, 1291 nrjTransf11,
1322 nrjTransf12, 1292 nrjTransf12,
1323 nrjTransf21, 1293 nrjTransf21,
1324 nrjTransf22, 1294 nrjTransf22,
1325 valueK1, 1295 valueK1,
1326 valueK2, 1296 valueK2,
1327 k, 1297 k,
1328 random); 1298 random);
1329 } 1299 }
1330 //G4cout << nrj << endl; 1300 //G4cout << nrj << endl;
1331 1301
1332 return nrj; 1302 return nrj;
1333 } 1303 }
1334 1304