| 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 // 26 //
27 // G4MicroElecInelasticModel.cc, 2011/08/29 A. 27 // G4MicroElecInelasticModel.cc, 2011/08/29 A.Valentin, M. Raine
28 // 28 //
29 // Based on the following publications 29 // Based on the following publications
30 // 30 //
31 // - Inelastic cross-sections of low 31 // - Inelastic cross-sections of low energy electrons in silicon
32 // for the simulation of heavy ion tracks 32 // for the simulation of heavy ion tracks with theGeant4-DNA toolkit,
33 // NSS Conf. Record 2010, pp. 80-85. 33 // NSS Conf. Record 2010, pp. 80-85.
34 // - Geant4 physics processes for microdo 34 // - Geant4 physics processes for microdosimetry simulation:
35 // very low energy electromagnetic models 35 // very low energy electromagnetic models for electrons in Si,
36 // NIM B, vol. 288, pp. 66 - 73, 2012. 36 // NIM B, vol. 288, pp. 66 - 73, 2012.
37 // - Geant4 physics processes for microdo 37 // - Geant4 physics processes for microdosimetry simulation:
38 // very low energy electromagnetic models 38 // very low energy electromagnetic models for protons and
39 // heavy ions in Si, NIM B, vol. 287, pp. 39 // heavy ions in Si, NIM B, vol. 287, pp. 124 - 129, 2012.
40 // 40 //
41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 42
43 #include "G4MicroElecInelasticModel.hh" 43 #include "G4MicroElecInelasticModel.hh"
44 44
45 #include "globals.hh" 45 #include "globals.hh"
46 #include "G4PhysicalConstants.hh" 46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
48 #include "G4ios.hh" 48 #include "G4ios.hh"
49 #include "G4UnitsTable.hh" 49 #include "G4UnitsTable.hh"
50 #include "G4UAtomicDeexcitation.hh" 50 #include "G4UAtomicDeexcitation.hh"
51 #include "G4MicroElecSiStructure.hh" 51 #include "G4MicroElecSiStructure.hh"
52 #include "G4LossTableManager.hh" 52 #include "G4LossTableManager.hh"
53 #include "G4ionEffectiveCharge.hh" 53 #include "G4ionEffectiveCharge.hh"
54 #include "G4MicroElecCrossSectionDataSet.hh" 54 #include "G4MicroElecCrossSectionDataSet.hh"
55 #include "G4Electron.hh" 55 #include "G4Electron.hh"
56 #include "G4Proton.hh" 56 #include "G4Proton.hh"
57 #include "G4GenericIon.hh" 57 #include "G4GenericIon.hh"
58 #include "G4ParticleDefinition.hh" 58 #include "G4ParticleDefinition.hh"
59 #include "G4NistManager.hh" 59 #include "G4NistManager.hh"
60 #include "G4LogLogInterpolation.hh" 60 #include "G4LogLogInterpolation.hh"
61 #include "G4DeltaAngle.hh" 61 #include "G4DeltaAngle.hh"
62 62
63 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 64
65 using namespace std; 65 using namespace std;
66 66
67 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 68
69 G4MicroElecInelasticModel::G4MicroElecInelasti 69 G4MicroElecInelasticModel::G4MicroElecInelasticModel(const G4ParticleDefinition*,
70 70 const G4String& nam)
71 :G4VEmModel(nam),isInitialised(false) 71 :G4VEmModel(nam),isInitialised(false)
72 { 72 {
73 nistSi = G4NistManager::Instance()->FindOrBu 73 nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
74 74
75 verboseLevel= 0; 75 verboseLevel= 0;
76 // Verbosity scale: 76 // Verbosity scale:
77 // 0 = nothing 77 // 0 = nothing
78 // 1 = warning for energy non-conservation 78 // 1 = warning for energy non-conservation
79 // 2 = details of energy budget 79 // 2 = details of energy budget
80 // 3 = calculation of cross sections, file o 80 // 3 = calculation of cross sections, file openings, sampling of atoms
81 // 4 = entering in methods 81 // 4 = entering in methods
82 82
83 if( verboseLevel>0 ) 83 if( verboseLevel>0 )
84 { 84 {
85 G4cout << "MicroElec inelastic model is co 85 G4cout << "MicroElec inelastic model is constructed " << G4endl;
86 } 86 }
87 87
88 //Mark this model as "applicable" for atomic 88 //Mark this model as "applicable" for atomic deexcitation
89 SetDeexcitationFlag(true); 89 SetDeexcitationFlag(true);
90 fAtomDeexcitation = nullptr; 90 fAtomDeexcitation = nullptr;
91 fParticleChangeForGamma = nullptr; 91 fParticleChangeForGamma = nullptr;
92 92
93 // default generator 93 // default generator
94 SetAngularDistribution(new G4DeltaAngle()); 94 SetAngularDistribution(new G4DeltaAngle());
95 95
96 // Selection of computation method 96 // Selection of computation method
97 fasterCode = true; //false; 97 fasterCode = true; //false;
98 } 98 }
99 99
100 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 101
102 G4MicroElecInelasticModel::~G4MicroElecInelast 102 G4MicroElecInelasticModel::~G4MicroElecInelasticModel()
103 { 103 {
104 // Cross section 104 // Cross section
105 for (auto & pos : tableData) 105 for (auto & pos : tableData)
106 { 106 {
107 G4MicroElecCrossSectionDataSet* table = po 107 G4MicroElecCrossSectionDataSet* table = pos.second;
108 delete table; 108 delete table;
109 } 109 }
110 110
111 // Final state 111 // Final state
112 eVecm.clear(); 112 eVecm.clear();
113 pVecm.clear(); 113 pVecm.clear();
114 114
115 } 115 }
116 116
117 //....oooOO0OOooo........oooOO0OOooo........oo 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 118
119 void G4MicroElecInelasticModel::Initialise(con 119 void G4MicroElecInelasticModel::Initialise(const G4ParticleDefinition* particle,
120 con 120 const G4DataVector& /*cuts*/)
121 { 121 {
122 122
123 if (verboseLevel > 3) 123 if (verboseLevel > 3)
124 G4cout << "Calling G4MicroElecInelasticMod 124 G4cout << "Calling G4MicroElecInelasticModel::Initialise()" << G4endl;
125 125
126 // Energy limits 126 // Energy limits
127 127
128 G4String fileElectron("microelec/sigma_inela 128 G4String fileElectron("microelec/sigma_inelastic_e_Si");
129 G4String fileProton("microelec/sigma_inelast 129 G4String fileProton("microelec/sigma_inelastic_p_Si");
130 130
131 G4ParticleDefinition* electronDef = G4Electr 131 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
132 G4ParticleDefinition* protonDef = G4Proton:: 132 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
133 G4String electron = electronDef->GetParticle 133 G4String electron = electronDef->GetParticleName();
134 G4String proton = protonDef->GetParticleName 134 G4String proton = protonDef->GetParticleName();;
135 135
136 G4double scaleFactor = 1e-18 * cm *cm; 136 G4double scaleFactor = 1e-18 * cm *cm;
137 137
138 const char* path = G4FindDataDir("G4LEDATA") 138 const char* path = G4FindDataDir("G4LEDATA");
139 139
140 // *** ELECTRON 140 // *** ELECTRON
141 tableFile[electron] = fileElectron; 141 tableFile[electron] = fileElectron;
142 lowEnergyLimit[electron] = 16.7 * eV; 142 lowEnergyLimit[electron] = 16.7 * eV;
143 highEnergyLimit[electron] = 100.0 * MeV; 143 highEnergyLimit[electron] = 100.0 * MeV;
144 144
145 // Cross section 145 // Cross section
146 G4MicroElecCrossSectionDataSet* tableE = new 146 G4MicroElecCrossSectionDataSet* tableE = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
147 tableE->LoadData(fileElectron); 147 tableE->LoadData(fileElectron);
148 148
149 tableData[electron] = tableE; 149 tableData[electron] = tableE;
150 150
151 // Final state 151 // Final state
152 152
153 std::ostringstream eFullFileName; 153 std::ostringstream eFullFileName;
154 154
155 if (fasterCode) eFullFileName << path << "/m 155 if (fasterCode) eFullFileName << path << "/microelec/sigmadiff_cumulated_inelastic_e_Si.dat";
156 else eFullFileName << path << "/microelec/si 156 else eFullFileName << path << "/microelec/sigmadiff_inelastic_e_Si.dat";
157 157
158 std::ifstream eDiffCrossSection(eFullFileNam 158 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
159 159
160 if (!eDiffCrossSection) 160 if (!eDiffCrossSection)
161 { 161 {
162 if (fasterCode) G4Exception("G4MicroElecI 162 if (fasterCode) G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
163 FatalException,"Missing data file:/microe 163 FatalException,"Missing data file:/microelec/sigmadiff_cumulated_inelastic_e_Si.dat");
164 164
165 else G4Exception("G4MicroElecInelasticMod 165 else G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
166 FatalException,"Missing data file:/microe 166 FatalException,"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
167 } 167 }
168 168
169 // Clear the arrays for re-initialization ca 169 // Clear the arrays for re-initialization case (MT mode)
170 // Octobre 22nd, 2014 - Melanie Raine 170 // Octobre 22nd, 2014 - Melanie Raine
171 eTdummyVec.clear(); 171 eTdummyVec.clear();
172 pTdummyVec.clear(); 172 pTdummyVec.clear();
173 173
174 eVecm.clear(); 174 eVecm.clear();
175 pVecm.clear(); 175 pVecm.clear();
176 176
177 for (int j=0; j<6; j++) 177 for (int j=0; j<6; j++)
178 { 178 {
179 eProbaShellMap[j].clear(); 179 eProbaShellMap[j].clear();
180 pProbaShellMap[j].clear(); 180 pProbaShellMap[j].clear();
181 181
182 eDiffCrossSectionData[j].clear(); 182 eDiffCrossSectionData[j].clear();
183 pDiffCrossSectionData[j].clear(); 183 pDiffCrossSectionData[j].clear();
184 184
185 eNrjTransfData[j].clear(); 185 eNrjTransfData[j].clear();
186 pNrjTransfData[j].clear(); 186 pNrjTransfData[j].clear();
187 } 187 }
188 188
189 // 189 //
190 eTdummyVec.push_back(0.); 190 eTdummyVec.push_back(0.);
191 while(!eDiffCrossSection.eof()) 191 while(!eDiffCrossSection.eof())
192 { 192 {
193 G4double tDummy; 193 G4double tDummy;
194 G4double eDummy; 194 G4double eDummy;
195 eDiffCrossSection>>tDummy>>eDummy; 195 eDiffCrossSection>>tDummy>>eDummy;
196 if (tDummy != eTdummyVec.back()) eTdummyVe 196 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
197 197
198 G4double tmp; 198 G4double tmp;
199 for (int j=0; j<6; j++) 199 for (int j=0; j<6; j++)
200 { 200 {
201 eDiffCrossSection>> tmp; 201 eDiffCrossSection>> tmp;
202 202
203 eDiffCrossSectionData[j][tDummy][eDummy] 203 eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
204 204
205 if (fasterCode) 205 if (fasterCode)
206 { 206 {
207 eNrjTransfData[j][tDummy][eDiffCrossSe 207 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
208 eProbaShellMap[j][tDummy].push_back(eD 208 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
209 } 209 }
210 else 210 else
211 { 211 {
212 // SI - only if eof is not reached ! 212 // SI - only if eof is not reached !
213 if (!eDiffCrossSection.eof()) eDiffCrossSect 213 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
214 eVecm[tDummy].push_back(eDummy); 214 eVecm[tDummy].push_back(eDummy);
215 } 215 }
216 216
217 } 217 }
218 } 218 }
219 // 219 //
220 220
221 // *** PROTON 221 // *** PROTON
222 tableFile[proton] = fileProton; 222 tableFile[proton] = fileProton;
223 lowEnergyLimit[proton] = 50. * keV; 223 lowEnergyLimit[proton] = 50. * keV;
224 highEnergyLimit[proton] = 10. * GeV; 224 highEnergyLimit[proton] = 10. * GeV;
225 225
226 // Cross section 226 // Cross section
227 G4MicroElecCrossSectionDataSet* tableP = new 227 G4MicroElecCrossSectionDataSet* tableP = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
228 tableP->LoadData(fileProton); 228 tableP->LoadData(fileProton);
229 tableData[proton] = tableP; 229 tableData[proton] = tableP;
230 230
231 // Final state 231 // Final state
232 std::ostringstream pFullFileName; 232 std::ostringstream pFullFileName;
233 233
234 if (fasterCode) pFullFileName << path << "/m 234 if (fasterCode) pFullFileName << path << "/microelec/sigmadiff_cumulated_inelastic_p_Si.dat";
235 else pFullFileName << path << "/microelec/si 235 else pFullFileName << path << "/microelec/sigmadiff_inelastic_p_Si.dat";
236 236
237 std::ifstream pDiffCrossSection(pFullFileNam 237 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
238 238
239 if (!pDiffCrossSection) 239 if (!pDiffCrossSection)
240 { 240 {
241 if (fasterCode) G4Exception("G4MicroElecIn 241 if (fasterCode) G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
242 FatalException,"Missing data file:/micro 242 FatalException,"Missing data file:/microelec/sigmadiff_cumulated_inelastic_p_Si.dat");
243 243
244 else G4Exception("G4MicroElecInelasticMode 244 else G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
245 FatalException,"Missing data file:/micro 245 FatalException,"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
246 } 246 }
247 247
248 pTdummyVec.push_back(0.); 248 pTdummyVec.push_back(0.);
249 while(!pDiffCrossSection.eof()) 249 while(!pDiffCrossSection.eof())
250 { 250 {
251 G4double tDummy; 251 G4double tDummy;
252 G4double eDummy; 252 G4double eDummy;
253 pDiffCrossSection>>tDummy>>eDummy; 253 pDiffCrossSection>>tDummy>>eDummy;
254 if (tDummy != pTdummyVec.back()) pTdummyVe 254 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
255 for (int j=0; j<6; j++) 255 for (int j=0; j<6; j++)
256 { 256 {
257 pDiffCrossSection>>pDiffCrossSectionData 257 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
258 258
259 if (fasterCode) 259 if (fasterCode)
260 { 260 {
261 pNrjTransfData[j][tDummy][pDiffCrossSe 261 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
262 pProbaShellMap[j][tDummy].push_back(pD 262 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
263 } 263 }
264 else 264 else
265 { 265 {
266 // SI - only if eof is not reached ! 266 // SI - only if eof is not reached !
267 if (!pDiffCrossSection.eof()) pDiffCrossSect 267 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
268 pVecm[tDummy].push_back(eDummy); 268 pVecm[tDummy].push_back(eDummy);
269 } 269 }
270 } 270 }
271 } 271 }
272 272
273 if (particle==electronDef) 273 if (particle==electronDef)
274 { 274 {
275 SetLowEnergyLimit(lowEnergyLimit[electron] 275 SetLowEnergyLimit(lowEnergyLimit[electron]);
276 SetHighEnergyLimit(highEnergyLimit[electro 276 SetHighEnergyLimit(highEnergyLimit[electron]);
277 } 277 }
278 278
279 if (particle==protonDef) 279 if (particle==protonDef)
280 { 280 {
281 SetLowEnergyLimit(lowEnergyLimit[proton]); 281 SetLowEnergyLimit(lowEnergyLimit[proton]);
282 SetHighEnergyLimit(highEnergyLimit[proton] 282 SetHighEnergyLimit(highEnergyLimit[proton]);
283 } 283 }
284 284
285 if( verboseLevel>0 ) 285 if( verboseLevel>0 )
286 { 286 {
287 G4cout << "MicroElec Inelastic model is in 287 G4cout << "MicroElec Inelastic model is initialized " << G4endl
288 << "Energy range: " 288 << "Energy range: "
289 << LowEnergyLimit() / keV << " keV - " 289 << LowEnergyLimit() / keV << " keV - "
290 << HighEnergyLimit() / MeV << " MeV for " 290 << HighEnergyLimit() / MeV << " MeV for "
291 << particle->GetParticleName() 291 << particle->GetParticleName()
292 << " with mass (amu) " << particle->GetPD 292 << " with mass (amu) " << particle->GetPDGMass()/proton_mass_c2
293 << " and charge " << particle->GetPDGChar 293 << " and charge " << particle->GetPDGCharge()
294 << G4endl << G4endl ; 294 << G4endl << G4endl ;
295 } 295 }
296 296
297 // 297 //
298 fAtomDeexcitation = G4LossTableManager::Ins 298 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
299 299
300 if (isInitialised) { return; } 300 if (isInitialised) { return; }
301 fParticleChangeForGamma = GetParticleChangeF 301 fParticleChangeForGamma = GetParticleChangeForGamma();
302 isInitialised = true; 302 isInitialised = true;
303 } 303 }
304 304
305 //....oooOO0OOooo........oooOO0OOooo........oo 305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
306 306
307 G4double G4MicroElecInelasticModel::CrossSecti 307 G4double G4MicroElecInelasticModel::CrossSectionPerVolume(const G4Material* material,
308 308 const G4ParticleDefinition* particleDefinition,
309 309 G4double ekin,
310 310 G4double,
311 311 G4double)
312 { 312 {
313 if (verboseLevel > 3) 313 if (verboseLevel > 3)
314 G4cout << "Calling CrossSectionPerVolume() 314 G4cout << "Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" << G4endl;
315 315
316 G4double density = material->GetTotNbOfAtoms 316 G4double density = material->GetTotNbOfAtomsPerVolume();
317 317
318 // Calculate total cross section for model 318 // Calculate total cross section for model
319 G4double lowLim = 0; 319 G4double lowLim = 0;
320 G4double highLim = 0; 320 G4double highLim = 0;
321 G4double sigma=0; 321 G4double sigma=0;
322 322
323 const G4String& particleName = particleDefin 323 const G4String& particleName = particleDefinition->GetParticleName();
324 G4String nameLocal = particleName ; 324 G4String nameLocal = particleName ;
325 325
326 G4double Zeff2 = 1.0; 326 G4double Zeff2 = 1.0;
327 G4double Mion_c2 = particleDefinition->GetPD 327 G4double Mion_c2 = particleDefinition->GetPDGMass();
328 328
329 if (Mion_c2 > proton_mass_c2) 329 if (Mion_c2 > proton_mass_c2)
330 { 330 {
331 G4ionEffectiveCharge EffCharge ; 331 G4ionEffectiveCharge EffCharge ;
332 G4double Zeff = EffCharge.EffectiveCharge( 332 G4double Zeff = EffCharge.EffectiveCharge(particleDefinition, material,ekin);
333 Zeff2 = Zeff*Zeff; 333 Zeff2 = Zeff*Zeff;
334 334
335 if (verboseLevel > 3) 335 if (verboseLevel > 3)
336 G4cout << "Before scaling : " << G4endl 336 G4cout << "Before scaling : " << G4endl
337 << "Particle : " << nameLocal << ", mass 337 << "Particle : " << nameLocal << ", mass : " << Mion_c2/proton_mass_c2 << "*mp, charge " << Zeff
338 << ", Ekin (eV) = " << ekin/eV << G4endl 338 << ", Ekin (eV) = " << ekin/eV << G4endl ;
339 339
340 ekin *= proton_mass_c2/Mion_c2 ; 340 ekin *= proton_mass_c2/Mion_c2 ;
341 nameLocal = "proton" ; 341 nameLocal = "proton" ;
342 342
343 if (verboseLevel > 3) 343 if (verboseLevel > 3)
344 G4cout << "After scaling : " << G4endl 344 G4cout << "After scaling : " << G4endl
345 << "Particle : " << nameLocal << ", Eki 345 << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl ;
346 } 346 }
347 347
348 if (material == nistSi || material->GetBaseM 348 if (material == nistSi || material->GetBaseMaterial() == nistSi)
349 { 349 {
350 auto pos1 = lowEnergyLimit.find(nameLocal) 350 auto pos1 = lowEnergyLimit.find(nameLocal);
351 if (pos1 != lowEnergyLimit.end()) 351 if (pos1 != lowEnergyLimit.end())
352 { 352 {
353 lowLim = pos1->second; 353 lowLim = pos1->second;
354 } 354 }
355 355
356 auto pos2 = highEnergyLimit.find(nameLocal 356 auto pos2 = highEnergyLimit.find(nameLocal);
357 if (pos2 != highEnergyLimit.end()) 357 if (pos2 != highEnergyLimit.end())
358 { 358 {
359 highLim = pos2->second; 359 highLim = pos2->second;
360 } 360 }
361 361
362 if (ekin >= lowLim && ekin < highLim) 362 if (ekin >= lowLim && ekin < highLim)
363 { 363 {
364 auto pos = tableData.find(nameLocal); 364 auto pos = tableData.find(nameLocal);
365 if (pos != tableData.end()) 365 if (pos != tableData.end())
366 { 366 {
367 G4MicroElecCrossSectionDataSet* table 367 G4MicroElecCrossSectionDataSet* table = pos->second;
368 if (table != nullptr) 368 if (table != nullptr)
369 { 369 {
370 sigma = table->FindValue(ekin); 370 sigma = table->FindValue(ekin);
371 } 371 }
372 } 372 }
373 else 373 else
374 { 374 {
375 G4Exception("G4MicroElecInelasticModel 375 G4Exception("G4MicroElecInelasticModel::CrossSectionPerVolume","em0002",
376 FatalException,"Model not applicable t 376 FatalException,"Model not applicable to particle type.");
377 } 377 }
378 } 378 }
379 else 379 else
380 { 380 {
381 if (nameLocal!="e-") 381 if (nameLocal!="e-")
382 { 382 {
383 // G4cout << "Particle : " << nameLoca 383 // G4cout << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl;
384 // G4cout << "### Warning: particle en 384 // G4cout << "### Warning: particle energy out of bounds! ###" << G4endl;
385 } 385 }
386 } 386 }
387 387
388 if (verboseLevel > 3) 388 if (verboseLevel > 3)
389 { 389 {
390 G4cout << "---> Kinetic energy (eV)=" << 390 G4cout << "---> Kinetic energy (eV)=" << ekin/eV << G4endl;
391 G4cout << " - Cross section per Si atom 391 G4cout << " - Cross section per Si atom (cm^2)=" << sigma*Zeff2/cm2 << G4endl;
392 G4cout << " - Cross section per Si atom 392 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./cm) << G4endl;
393 } 393 }
394 394
395 } // if (SiMaterial) 395 } // if (SiMaterial)
396 return sigma*density*Zeff2; 396 return sigma*density*Zeff2;
397 } 397 }
398 398
399 //....oooOO0OOooo........oooOO0OOooo........oo 399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 400
401 void G4MicroElecInelasticModel::SampleSecondar 401 void G4MicroElecInelasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
402 402 const G4MaterialCutsCouple* couple,
403 403 const G4DynamicParticle* particle,
404 404 G4double,
405 405 G4double)
406 { 406 {
407 407
408 if (verboseLevel > 3) 408 if (verboseLevel > 3)
409 G4cout << "Calling SampleSecondaries() of 409 G4cout << "Calling SampleSecondaries() of G4MicroElecInelasticModel" << G4endl;
410 410
411 G4double lowLim = 0; 411 G4double lowLim = 0;
412 G4double highLim = 0; 412 G4double highLim = 0;
413 413
414 G4double ekin = particle->GetKineticEnergy() 414 G4double ekin = particle->GetKineticEnergy();
415 G4double k = ekin ; 415 G4double k = ekin ;
416 416
417 G4ParticleDefinition* PartDef = particle->Ge 417 G4ParticleDefinition* PartDef = particle->GetDefinition();
418 const G4String& particleName = PartDef->GetP 418 const G4String& particleName = PartDef->GetParticleName();
419 G4String nameLocal2 = particleName ; 419 G4String nameLocal2 = particleName ;
420 G4double particleMass = particle->GetDefinit 420 G4double particleMass = particle->GetDefinition()->GetPDGMass();
421 421
422 if (particleMass > proton_mass_c2) 422 if (particleMass > proton_mass_c2)
423 { 423 {
424 k *= proton_mass_c2/particleMass ; 424 k *= proton_mass_c2/particleMass ;
425 PartDef = G4Proton::ProtonDefinition(); 425 PartDef = G4Proton::ProtonDefinition();
426 nameLocal2 = "proton" ; 426 nameLocal2 = "proton" ;
427 } 427 }
428 428
429 auto pos1 = lowEnergyLimit.find(nameLocal2); 429 auto pos1 = lowEnergyLimit.find(nameLocal2);
430 if (pos1 != lowEnergyLimit.end()) 430 if (pos1 != lowEnergyLimit.end())
431 { 431 {
432 lowLim = pos1->second; 432 lowLim = pos1->second;
433 } 433 }
434 434
435 auto pos2 = highEnergyLimit.find(nameLocal2) 435 auto pos2 = highEnergyLimit.find(nameLocal2);
436 if (pos2 != highEnergyLimit.end()) 436 if (pos2 != highEnergyLimit.end())
437 { 437 {
438 highLim = pos2->second; 438 highLim = pos2->second;
439 } 439 }
440 440
441 if (k >= lowLim && k < highLim) 441 if (k >= lowLim && k < highLim)
442 { 442 {
443 G4ParticleMomentum primaryDirection = part 443 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
444 G4double totalEnergy = ekin + particleMass 444 G4double totalEnergy = ekin + particleMass;
445 G4double pSquare = ekin * (totalEnergy + p 445 G4double pSquare = ekin * (totalEnergy + particleMass);
446 G4double totalMomentum = std::sqrt(pSquare 446 G4double totalMomentum = std::sqrt(pSquare);
447 447
448 G4int Shell = 0; 448 G4int Shell = 0;
449 449
450 /* if (!fasterCode)*/ Shell = RandomSelect 450 /* if (!fasterCode)*/ Shell = RandomSelect(k,nameLocal2);
451 451
452 // SI: The following protection is necessa 452 // SI: The following protection is necessary to avoid infinite loops :
453 // sigmadiff_ionisation_e_born.dat has no 453 // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
454 // sigmadiff_cumulated_ionisation_e_born. 454 // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
455 // this is due to the fact that the max a 455 // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
456 // strictly above this value have non zer 456 // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
457 457
458 G4double bindingEnergy = SiStructure.Energ 458 G4double bindingEnergy = SiStructure.Energy(Shell);
459 459
460 if (verboseLevel > 3) 460 if (verboseLevel > 3)
461 { 461 {
462 G4cout << "---> Kinetic energy (eV)=" << 462 G4cout << "---> Kinetic energy (eV)=" << k/eV << G4endl ;
463 G4cout << "Shell: " << Shell << ", energ 463 G4cout << "Shell: " << Shell << ", energy: " << bindingEnergy/eV << G4endl;
464 } 464 }
465 465
466 // sample deexcitation 466 // sample deexcitation
467 std::size_t secNumberInit = 0; // need to 467 std::size_t secNumberInit = 0; // need to know at a certain point the energy of secondaries
468 std::size_t secNumberFinal = 0; // So I'll 468 std::size_t secNumberFinal = 0; // So I'll make the difference and then sum the energies
469 469
470 //SI: additional protection if tcs interpo 470 //SI: additional protection if tcs interpolation method is modified
471 if (k<bindingEnergy) return; 471 if (k<bindingEnergy) return;
472 472
473 G4int Z = 14; 473 G4int Z = 14;
474 474
475 if(fAtomDeexcitation && Shell > 2) { 475 if(fAtomDeexcitation && Shell > 2) {
476 476
477 G4AtomicShellEnumerator as = fKShell; 477 G4AtomicShellEnumerator as = fKShell;
478 478
479 if (Shell == 4) 479 if (Shell == 4)
480 { 480 {
481 as = G4AtomicShellEnumerator(1); 481 as = G4AtomicShellEnumerator(1);
482 } 482 }
483 else if (Shell == 3) 483 else if (Shell == 3)
484 { 484 {
485 as = G4AtomicShellEnumerator(3); 485 as = G4AtomicShellEnumerator(3);
486 } 486 }
487 487
488 const G4AtomicShell* shell = fAtomDeexci 488 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
489 secNumberInit = fvect->size(); 489 secNumberInit = fvect->size();
490 fAtomDeexcitation->GenerateParticles(fve 490 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
491 secNumberFinal = fvect->size(); 491 secNumberFinal = fvect->size();
492 } 492 }
493 493
494 G4double secondaryKinetic=-1000*eV; 494 G4double secondaryKinetic=-1000*eV;
495 495
496 if (!fasterCode) 496 if (!fasterCode)
497 { 497 {
498 secondaryKinetic = RandomizeEjectedElect 498 secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef,k,Shell);
499 } 499 }
500 else 500 else
501 { 501 {
502 secondaryKinetic = RandomizeEjectedElect 502 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(PartDef,k,Shell);
503 } 503 }
504 504
505 if (verboseLevel > 3) 505 if (verboseLevel > 3)
506 { 506 {
507 G4cout << "Ionisation process" << G4endl 507 G4cout << "Ionisation process" << G4endl;
508 G4cout << "Shell: " << Shell << " Kin. e 508 G4cout << "Shell: " << Shell << " Kin. energy (eV)=" << k/eV
509 << " Sec. energy (eV)=" << secondaryKine 509 << " Sec. energy (eV)=" << secondaryKinetic/eV << G4endl;
510 } 510 }
511 511
512 G4ThreeVector deltaDirection = 512 G4ThreeVector deltaDirection =
513 GetAngularDistribution()->SampleDirectionF 513 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
514 514 Z, Shell,
515 515 couple->GetMaterial());
516 516
517 if (particle->GetDefinition() == G4Electro 517 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
518 { 518 {
519 G4double deltaTotalMomentum = std::sqrt( 519 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
520 G4double finalPx = totalMomentum*primary 520 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
521 G4double finalPy = totalMomentum*primary 521 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
522 G4double finalPz = totalMomentum*primary 522 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
523 G4double finalMomentum = std::sqrt(final 523 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
524 finalPx /= finalMomentum; 524 finalPx /= finalMomentum;
525 finalPy /= finalMomentum; 525 finalPy /= finalMomentum;
526 finalPz /= finalMomentum; 526 finalPz /= finalMomentum;
527 527
528 G4ThreeVector direction; 528 G4ThreeVector direction;
529 direction.set(finalPx,finalPy,finalPz); 529 direction.set(finalPx,finalPy,finalPz);
530 530
531 fParticleChangeForGamma->ProposeMomentum 531 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
532 } 532 }
533 else 533 else
534 fParticleChangeForGamma->ProposeMomentum 534 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
535 535
536 // note that secondaryKinetic is the energ 536 // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
537 G4double deexSecEnergy = 0; 537 G4double deexSecEnergy = 0;
538 for (std::size_t j=secNumberInit; j < secN 538 for (std::size_t j=secNumberInit; j < secNumberFinal; ++j) {
539 deexSecEnergy = deexSecEnergy + (*fvect) 539 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
540 540
541 fParticleChangeForGamma->SetProposedKineti 541 fParticleChangeForGamma->SetProposedKineticEnergy(ekin-bindingEnergy-secondaryKinetic);
542 fParticleChangeForGamma->ProposeLocalEnerg 542 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy-deexSecEnergy);
543 543
544 if (secondaryKinetic>0) 544 if (secondaryKinetic>0)
545 { 545 {
546 G4DynamicParticle* dp = new G4DynamicPar 546 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
547 fvect->push_back(dp); 547 fvect->push_back(dp);
548 } 548 }
549 } 549 }
550 } 550 }
551 551
552 //....oooOO0OOooo........oooOO0OOooo........oo 552 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553 553
554 G4double G4MicroElecInelasticModel::RandomizeE 554 G4double G4MicroElecInelasticModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
555 555 G4double k, G4int shell)
556 { 556 {
557 if (particleDefinition == G4Electron::Electr 557 if (particleDefinition == G4Electron::ElectronDefinition())
558 { 558 {
559 G4double maximumEnergyTransfer=0.; 559 G4double maximumEnergyTransfer=0.;
560 if ((k+SiStructure.Energy(shell))/2. > k) 560 if ((k+SiStructure.Energy(shell))/2. > k) maximumEnergyTransfer=k;
561 else maximumEnergyTransfer = (k+SiStructur 561 else maximumEnergyTransfer = (k+SiStructure.Energy(shell))/2.;
562 562
563 G4double crossSectionMaximum = 0.; 563 G4double crossSectionMaximum = 0.;
564 564
565 G4double minEnergy = SiStructure.Energy(sh 565 G4double minEnergy = SiStructure.Energy(shell);
566 G4double maxEnergy = maximumEnergyTransfer 566 G4double maxEnergy = maximumEnergyTransfer;
567 G4int nEnergySteps = 100; 567 G4int nEnergySteps = 100;
568 568
569 G4double value(minEnergy); 569 G4double value(minEnergy);
570 G4double stpEnergy(std::pow(maxEnergy/valu 570 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
571 G4int step(nEnergySteps); 571 G4int step(nEnergySteps);
572 while (step>0) 572 while (step>0)
573 { 573 {
574 step--; 574 step--;
575 G4double differentialCrossSection = Diff 575 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
576 if(differentialCrossSection >= crossSect 576 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
577 value*=stpEnergy; 577 value*=stpEnergy;
578 } 578 }
579 579
580 G4double secondaryElectronKineticEnergy=0. 580 G4double secondaryElectronKineticEnergy=0.;
581 do 581 do
582 { 582 {
583 secondaryElectronKineticEnergy = G4Unifo 583 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
584 } while(G4UniformRand()*crossSectionMaximu 584 } while(G4UniformRand()*crossSectionMaximum >
585 DifferentialCrossSection(particleD 585 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
586 586
587 return secondaryElectronKineticEnergy; 587 return secondaryElectronKineticEnergy;
588 } 588 }
589 589
590 if (particleDefinition == G4Proton::ProtonDe 590 if (particleDefinition == G4Proton::ProtonDefinition())
591 { 591 {
592 G4double maximumEnergyTransfer = 4.* (elec 592 G4double maximumEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
593 G4double crossSectionMaximum = 0.; 593 G4double crossSectionMaximum = 0.;
594 594
595 G4double minEnergy = SiStructure.Energy(sh 595 G4double minEnergy = SiStructure.Energy(shell);
596 G4double maxEnergy = maximumEnergyTransfer 596 G4double maxEnergy = maximumEnergyTransfer;
597 G4int nEnergySteps = 100; 597 G4int nEnergySteps = 100;
598 598
599 G4double value(minEnergy); 599 G4double value(minEnergy);
600 G4double stpEnergy(std::pow(maxEnergy/valu 600 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
601 G4int step(nEnergySteps); 601 G4int step(nEnergySteps);
602 while (step>0) 602 while (step>0)
603 { 603 {
604 step--; 604 step--;
605 G4double differentialCrossSection = Diff 605 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
606 if(differentialCrossSection >= crossSect 606 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
607 value*=stpEnergy; 607 value*=stpEnergy;
608 } 608 }
609 609
610 G4double secondaryElectronKineticEnergy = 610 G4double secondaryElectronKineticEnergy = 0.;
611 do 611 do
612 { 612 {
613 secondaryElectronKineticEnergy = G4Unifo 613 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
614 614
615 } while(G4UniformRand()*crossSectionMaximu 615 } while(G4UniformRand()*crossSectionMaximum >
616 DifferentialCrossSection(particleD 616 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
617 return secondaryElectronKineticEnergy; 617 return secondaryElectronKineticEnergy;
618 } 618 }
619 619
620 return 0; 620 return 0;
621 } 621 }
622 622
623 //....oooOO0OOooo........oooOO0OOooo........oo 623 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
624 624
625 // The following section is not used anymore b 625 // The following section is not used anymore but is kept for memory
626 // GetAngularDistribution()->SampleDirectionFo 626 // GetAngularDistribution()->SampleDirectionForShell is used instead
627 627
628 /*void G4MicroElecInelasticModel::RandomizeEje 628 /*void G4MicroElecInelasticModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
629 G4double k, 629 G4double k,
630 G4double secKinetic, 630 G4double secKinetic,
631 G4double & cosTheta, 631 G4double & cosTheta,
632 G4double & phi ) 632 G4double & phi )
633 { 633 {
634 if (particleDefinition == G4Electron::Electro 634 if (particleDefinition == G4Electron::ElectronDefinition())
635 { 635 {
636 phi = twopi * G4UniformRand(); 636 phi = twopi * G4UniformRand();
637 G4double sin2O = (1.-secKinetic/k) / (1.+secK 637 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
638 cosTheta = std::sqrt(1.-sin2O); 638 cosTheta = std::sqrt(1.-sin2O);
639 } 639 }
640 640
641 if (particleDefinition == G4Proton::ProtonDef 641 if (particleDefinition == G4Proton::ProtonDefinition())
642 { 642 {
643 G4double maxSecKinetic = 4.* (electron_mass_c 643 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
644 phi = twopi * G4UniformRand(); 644 phi = twopi * G4UniformRand();
645 cosTheta = std::sqrt(secKinetic / maxSecKinet 645 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
646 } 646 }
647 647
648 else 648 else
649 { 649 {
650 G4double maxSecKinetic = 4.* (electron_mass_c 650 G4double maxSecKinetic = 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
651 phi = twopi * G4UniformRand(); 651 phi = twopi * G4UniformRand();
652 cosTheta = std::sqrt(secKinetic / maxSecKinet 652 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
653 } 653 }
654 } 654 }
655 */ 655 */
656 656
657 //....oooOO0OOooo........oooOO0OOooo........oo 657 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
658 658
659 G4double G4MicroElecInelasticModel::Differenti 659 G4double G4MicroElecInelasticModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
660 660 G4double k,
661 661 G4double energyTransfer,
662 662 G4int LevelIndex)
663 { 663 {
664 G4double sigma = 0.; 664 G4double sigma = 0.;
665 665
666 if (energyTransfer >= SiStructure.Energy(Lev 666 if (energyTransfer >= SiStructure.Energy(LevelIndex))
667 { 667 {
668 G4double valueT1 = 0; 668 G4double valueT1 = 0;
669 G4double valueT2 = 0; 669 G4double valueT2 = 0;
670 G4double valueE21 = 0; 670 G4double valueE21 = 0;
671 G4double valueE22 = 0; 671 G4double valueE22 = 0;
672 G4double valueE12 = 0; 672 G4double valueE12 = 0;
673 G4double valueE11 = 0; 673 G4double valueE11 = 0;
674 674
675 G4double xs11 = 0; 675 G4double xs11 = 0;
676 G4double xs12 = 0; 676 G4double xs12 = 0;
677 G4double xs21 = 0; 677 G4double xs21 = 0;
678 G4double xs22 = 0; 678 G4double xs22 = 0;
679 679
680 if (particleDefinition == G4Electron::Elec 680 if (particleDefinition == G4Electron::ElectronDefinition())
681 { 681 {
682 // k should be in eV and energy transfer 682 // k should be in eV and energy transfer eV also
683 auto t2 = std::upper_bound(eTdummyVec.be 683 auto t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
684 auto t1 = t2-1; 684 auto t1 = t2-1;
685 // The following condition avoids situat 685 // The following condition avoids situations where energyTransfer >last vector element
686 if (energyTransfer <= eVecm[(*t1)].back( 686 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
687 { 687 {
688 auto e12 = std::upper_bound(eVecm[(*t1 688 auto e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
689 auto e11 = e12-1; 689 auto e11 = e12-1;
690 auto e22 = std::upper_bound(eVecm[(*t2 690 auto e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
691 auto e21 = e22-1; 691 auto e21 = e22-1;
692 692
693 valueT1 =*t1; 693 valueT1 =*t1;
694 valueT2 =*t2; 694 valueT2 =*t2;
695 valueE21 =*e21; 695 valueE21 =*e21;
696 valueE22 =*e22; 696 valueE22 =*e22;
697 valueE12 =*e12; 697 valueE12 =*e12;
698 valueE11 =*e11; 698 valueE11 =*e11;
699 699
700 xs11 = eDiffCrossSectionData[LevelInde 700 xs11 = eDiffCrossSectionData[LevelIndex][valueT1][valueE11];
701 xs12 = eDiffCrossSectionData[LevelInde 701 xs12 = eDiffCrossSectionData[LevelIndex][valueT1][valueE12];
702 xs21 = eDiffCrossSectionData[LevelInde 702 xs21 = eDiffCrossSectionData[LevelIndex][valueT2][valueE21];
703 xs22 = eDiffCrossSectionData[LevelInde 703 xs22 = eDiffCrossSectionData[LevelIndex][valueT2][valueE22];
704 } 704 }
705 } 705 }
706 706
707 if (particleDefinition == G4Proton::Proton 707 if (particleDefinition == G4Proton::ProtonDefinition())
708 { 708 {
709 // k should be in eV and energy transfer 709 // k should be in eV and energy transfer eV also
710 auto t2 = std::upper_bound(pTdummyVec.be 710 auto t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
711 auto t1 = t2-1; 711 auto t1 = t2-1;
712 if (energyTransfer <= pVecm[(*t1)].back( 712 if (energyTransfer <= pVecm[(*t1)].back() && energyTransfer <= pVecm[(*t2)].back() )
713 { 713 {
714 auto e12 = std::upper_bound(pVecm[(*t1 714 auto e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
715 auto e11 = e12-1; 715 auto e11 = e12-1;
716 716
717 auto e22 = std::upper_bound(pVecm[(*t2 717 auto e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
718 auto e21 = e22-1; 718 auto e21 = e22-1;
719 719
720 valueT1 =*t1; 720 valueT1 =*t1;
721 valueT2 =*t2; 721 valueT2 =*t2;
722 valueE21 =*e21; 722 valueE21 =*e21;
723 valueE22 =*e22; 723 valueE22 =*e22;
724 valueE12 =*e12; 724 valueE12 =*e12;
725 valueE11 =*e11; 725 valueE11 =*e11;
726 726
727 xs11 = pDiffCrossSectionData[LevelInde 727 xs11 = pDiffCrossSectionData[LevelIndex][valueT1][valueE11];
728 xs12 = pDiffCrossSectionData[LevelInde 728 xs12 = pDiffCrossSectionData[LevelIndex][valueT1][valueE12];
729 xs21 = pDiffCrossSectionData[LevelInde 729 xs21 = pDiffCrossSectionData[LevelIndex][valueT2][valueE21];
730 xs22 = pDiffCrossSectionData[LevelInde 730 xs22 = pDiffCrossSectionData[LevelIndex][valueT2][valueE22];
731 } 731 }
732 } 732 }
733 733
734 sigma = QuadInterpolator( valueE11, va 734 sigma = QuadInterpolator( valueE11, valueE12,
735 valueE21, valueE22, 735 valueE21, valueE22,
736 xs11, xs12, 736 xs11, xs12,
737 xs21, xs22, 737 xs21, xs22,
738 valueT1, valueT2, 738 valueT1, valueT2,
739 k, energyTransfer); 739 k, energyTransfer);
740 } 740 }
741 return sigma; 741 return sigma;
742 } 742 }
743 743
744 //....oooOO0OOooo........oooOO0OOooo........oo 744 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
745 745
746 G4double G4MicroElecInelasticModel::Interpolat 746 G4double G4MicroElecInelasticModel::Interpolate(G4double e1,
747 747 G4double e2,
748 748 G4double e,
749 749 G4double xs1,
750 750 G4double xs2)
751 { 751 {
752 G4double value = 0.; 752 G4double value = 0.;
753 753
754 // Log-log interpolation by default 754 // Log-log interpolation by default
755 if (e1 != 0 && e2 != 0 && (std::log10(e2) - 755 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
756 && !fasterCode) 756 && !fasterCode)
757 { 757 {
758 G4double a = (std::log10(xs2)-std::log10(x 758 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
759 G4double b = std::log10(xs2) - a*std::log1 759 G4double b = std::log10(xs2) - a*std::log10(e2);
760 G4double sigma = a*std::log10(e) + b; 760 G4double sigma = a*std::log10(e) + b;
761 value = (std::pow(10.,sigma)); 761 value = (std::pow(10.,sigma));
762 762
763 } 763 }
764 764
765 // Switch to log-lin interpolation for faste 765 // Switch to log-lin interpolation for faster code
766 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 & 766 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
767 { 767 {
768 G4double d1 = std::log10(xs1); 768 G4double d1 = std::log10(xs1);
769 G4double d2 = std::log10(xs2); 769 G4double d2 = std::log10(xs2);
770 value = std::pow(10., (d1 + (d2 - d1) * (e 770 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
771 } 771 }
772 772
773 // Switch to lin-lin interpolation for faste 773 // Switch to lin-lin interpolation for faster code
774 // in case one of xs1 or xs2 (=cum proba) va 774 // in case one of xs1 or xs2 (=cum proba) value is zero
775 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) 775 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0)) // && fasterCode)
776 { 776 {
777 G4double d1 = xs1; 777 G4double d1 = xs1;
778 G4double d2 = xs2; 778 G4double d2 = xs2;
779 value = (d1 + (d2 - d1) * (e - e1) / (e2 - 779 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
780 } 780 }
781 781
782 return value; 782 return value;
783 } 783 }
784 784
785 //....oooOO0OOooo........oooOO0OOooo........oo 785 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
786 786
787 G4double G4MicroElecInelasticModel::QuadInterp 787 G4double G4MicroElecInelasticModel::QuadInterpolator(G4double e11, G4double e12,
788 788 G4double e21, G4double e22,
789 789 G4double xs11, G4double xs12,
790 790 G4double xs21, G4double xs22,
791 791 G4double t1, G4double t2,
792 792 G4double t, G4double e)
793 { 793 {
794 G4double interpolatedvalue1 = Interpolate(e1 794 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
795 G4double interpolatedvalue2 = Interpolate(e2 795 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
796 G4double value = Interpolate(t1, t2, t, inte 796 G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
797 return value; 797 return value;
798 } 798 }
799 799
800 //....oooOO0OOooo........oooOO0OOooo........oo 800 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
801 801
802 G4int G4MicroElecInelasticModel::RandomSelect( 802 G4int G4MicroElecInelasticModel::RandomSelect(G4double k, const G4String& particle )
803 { 803 {
804 G4int level = 0; 804 G4int level = 0;
805 805
806 auto pos = tableData.find(particle); 806 auto pos = tableData.find(particle);
807 if (pos != tableData.cend()) 807 if (pos != tableData.cend())
808 { 808 {
809 G4MicroElecCrossSectionDataSet* table = po 809 G4MicroElecCrossSectionDataSet* table = pos->second;
810 810
811 if (table != nullptr) 811 if (table != nullptr)
812 { 812 {
813 G4double* valuesBuffer = new G4double[ta 813 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
814 const G4int n = (G4int)table->NumberOfCo 814 const G4int n = (G4int)table->NumberOfComponents();
815 G4int i(n); 815 G4int i(n);
816 G4double value = 0.; 816 G4double value = 0.;
817 817
818 while (i>0) 818 while (i>0)
819 { 819 {
820 --i; 820 --i;
821 valuesBuffer[i] = table->GetComponent( 821 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
822 value += valuesBuffer[i]; 822 value += valuesBuffer[i];
823 } 823 }
824 824
825 value *= G4UniformRand(); 825 value *= G4UniformRand();
826 826
827 i = n; 827 i = n;
828 828
829 while (i > 0) 829 while (i > 0)
830 { 830 {
831 --i; 831 --i;
832 832
833 if (valuesBuffer[i] > value) 833 if (valuesBuffer[i] > value)
834 { 834 {
835 delete[] valuesBuffer; 835 delete[] valuesBuffer;
836 return i; 836 return i;
837 } 837 }
838 value -= valuesBuffer[i]; 838 value -= valuesBuffer[i];
839 } 839 }
840 840
841 if (valuesBuffer) delete[] valuesBuffer; 841 if (valuesBuffer) delete[] valuesBuffer;
842 842
843 } 843 }
844 } 844 }
845 else 845 else
846 { 846 {
847 G4Exception("G4MicroElecInelasticModel::Ra 847 G4Exception("G4MicroElecInelasticModel::RandomSelect","em0002",FatalException,"Model not applicable to particle type.");
848 } 848 }
849 849
850 return level; 850 return level;
851 } 851 }
852 852
853 //....oooOO0OOooo........oooOO0OOooo........oo 853 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
854 854
855 G4double G4MicroElecInelasticModel::RandomizeE 855 G4double G4MicroElecInelasticModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
856 856 G4double k,
857 857 G4int shell)
858 { 858 {
859 859
860 G4double secondaryElectronKineticEnergy = 0. 860 G4double secondaryElectronKineticEnergy = 0.;
861 G4double random = G4UniformRand(); 861 G4double random = G4UniformRand();
862 secondaryElectronKineticEnergy = TransferedE 862 secondaryElectronKineticEnergy = TransferedEnergy(particleDefinition,
863 863 k / eV,
864 864 shell,
865 865 random) * eV
866 - SiStructure.Energy(shell); 866 - SiStructure.Energy(shell);
867 867
868 if (secondaryElectronKineticEnergy < 0.) 868 if (secondaryElectronKineticEnergy < 0.)
869 return 0.; 869 return 0.;
870 // 870 //
871 return secondaryElectronKineticEnergy; 871 return secondaryElectronKineticEnergy;
872 } 872 }
873 873
874 //....oooOO0OOooo........oooOO0OOooo........oo 874 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
875 875
876 G4double G4MicroElecInelasticModel::Transfered 876 G4double G4MicroElecInelasticModel::TransferedEnergy(G4ParticleDefinition* particleDefinition,
877 877 G4double k,
878 878 G4int ionizationLevelIndex,
879 879 G4double random)
880 { 880 {
881 G4double nrj = 0.; 881 G4double nrj = 0.;
882 882
883 G4double valueK1 = 0; 883 G4double valueK1 = 0;
884 G4double valueK2 = 0; 884 G4double valueK2 = 0;
885 G4double valuePROB21 = 0; 885 G4double valuePROB21 = 0;
886 G4double valuePROB22 = 0; 886 G4double valuePROB22 = 0;
887 G4double valuePROB12 = 0; 887 G4double valuePROB12 = 0;
888 G4double valuePROB11 = 0; 888 G4double valuePROB11 = 0;
889 889
890 G4double nrjTransf11 = 0; 890 G4double nrjTransf11 = 0;
891 G4double nrjTransf12 = 0; 891 G4double nrjTransf12 = 0;
892 G4double nrjTransf21 = 0; 892 G4double nrjTransf21 = 0;
893 G4double nrjTransf22 = 0; 893 G4double nrjTransf22 = 0;
894 894
895 G4double maximumEnergyTransfer1 = 0; 895 G4double maximumEnergyTransfer1 = 0;
896 G4double maximumEnergyTransfer2 = 0; 896 G4double maximumEnergyTransfer2 = 0;
897 G4double maximumEnergyTransferP = 4.* (elect 897 G4double maximumEnergyTransferP = 4.* (electron_mass_c2 / proton_mass_c2) * k;
898 G4double bindingEnergy = SiStructure.Energy( 898 G4double bindingEnergy = SiStructure.Energy(ionizationLevelIndex)*1e6;
899 899
900 if (particleDefinition == G4Electron::Electr 900 if (particleDefinition == G4Electron::ElectronDefinition())
901 { 901 {
902 // k should be in eV 902 // k should be in eV
903 auto k2 = std::upper_bound(eTdummyVec.begi 903 auto k2 = std::upper_bound(eTdummyVec.begin(),
904 eTdummyVec.end(), 904 eTdummyVec.end(),
905 k); 905 k);
906 auto k1 = k2 - 1; 906 auto k1 = k2 - 1;
907 907
908 /* 908 /*
909 G4cout << "----> k=" << k 909 G4cout << "----> k=" << k
910 << " " << *k1 910 << " " << *k1
911 << " " << *k2 911 << " " << *k2
912 << " " << random 912 << " " << random
913 << " " << ionizationLevelIndex 913 << " " << ionizationLevelIndex
914 << " " << eProbaShellMap[ionizationLevelI 914 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
915 << " " << eProbaShellMap[ionizationLevelI 915 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
916 << G4endl; 916 << G4endl;
917 */ 917 */
918 918
919 // SI : the following condition avoids sit 919 // SI : the following condition avoids situations where random >last vector element
920 if (random <= eProbaShellMap[ionizationLev 920 if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
921 && random <= eProbaShellMap[ionization 921 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
922 { 922 {
923 auto prob12 = 923 auto prob12 =
924 std::upper_bound(eProbaShellMap[ioni 924 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
925 eProbaShellMap[ioni 925 eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
926 random); 926 random);
927 auto prob11 = prob12 - 1; 927 auto prob11 = prob12 - 1;
928 auto prob22 = 928 auto prob22 =
929 std::upper_bound(eProbaShellMap[ioni 929 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
930 eProbaShellMap[ioni 930 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
931 random); 931 random);
932 auto prob21 = prob22 - 1; 932 auto prob21 = prob22 - 1;
933 933
934 valueK1 = *k1; 934 valueK1 = *k1;
935 valueK2 = *k2; 935 valueK2 = *k2;
936 valuePROB21 = *prob21; 936 valuePROB21 = *prob21;
937 valuePROB22 = *prob22; 937 valuePROB22 = *prob22;
938 valuePROB12 = *prob12; 938 valuePROB12 = *prob12;
939 valuePROB11 = *prob11; 939 valuePROB11 = *prob11;
940 940
941 // The following condition avoid getting 941 // The following condition avoid getting transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
942 if(valuePROB11 == 0) nrjTransf11 = bindi 942 if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
943 else nrjTransf11 = eNrjTransfData[ioniza 943 else nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
944 if(valuePROB12 == 1) 944 if(valuePROB12 == 1)
945 { 945 {
946 if ((valueK1+bindingEnergy)/2. > valueK1) ma 946 if ((valueK1+bindingEnergy)/2. > valueK1) maximumEnergyTransfer1=valueK1;
947 else maximumEnergyTransfer1 = (valueK1+b 947 else maximumEnergyTransfer1 = (valueK1+bindingEnergy)/2.;
948 948
949 nrjTransf12 = maximumEnergyTransfer1; 949 nrjTransf12 = maximumEnergyTransfer1;
950 } 950 }
951 else nrjTransf12 = eNrjTransfData[ioniza 951 else nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
952 952
953 if(valuePROB21 == 0) nrjTransf21 = bindi 953 if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
954 else nrjTransf21 = eNrjTransfData[ioniza 954 else nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
955 if(valuePROB22 == 1) 955 if(valuePROB22 == 1)
956 { 956 {
957 if ((valueK2+bindingEnergy)/2. > valueK2) ma 957 if ((valueK2+bindingEnergy)/2. > valueK2) maximumEnergyTransfer2=valueK2;
958 else maximumEnergyTransfer2 = (valueK2+b 958 else maximumEnergyTransfer2 = (valueK2+bindingEnergy)/2.;
959 959
960 nrjTransf22 = maximumEnergyTransfer2; 960 nrjTransf22 = maximumEnergyTransfer2;
961 } 961 }
962 else 962 else
963 nrjTransf22 = eNrjTransfData[ionizationLevel 963 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
964 } 964 }
965 // Avoids cases where cum xs is zero for k 965 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
966 if (random > eProbaShellMap[ionizationLeve 966 if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
967 { 967 {
968 auto prob22 = 968 auto prob22 =
969 std::upper_bound(eProbaShellMap[ioni 969 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
970 eProbaShellMap[ioni 970 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
971 random); 971 random);
972 972
973 auto prob21 = prob22 - 1; 973 auto prob21 = prob22 - 1;
974 974
975 valueK1 = *k1; 975 valueK1 = *k1;
976 valueK2 = *k2; 976 valueK2 = *k2;
977 valuePROB21 = *prob21; 977 valuePROB21 = *prob21;
978 valuePROB22 = *prob22; 978 valuePROB22 = *prob22;
979 979
980 nrjTransf21 = eNrjTransfData[ionizationL 980 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
981 nrjTransf22 = eNrjTransfData[ionizationL 981 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
982 982
983 G4double interpolatedvalue2 = Interpolat 983 G4double interpolatedvalue2 = Interpolate(valuePROB21,
984 984 valuePROB22,
985 985 random,
986 986 nrjTransf21,
987 987 nrjTransf22);
988 988
989 // zeros are explicitly set 989 // zeros are explicitly set
990 G4double value = Interpolate(valueK1, va 990 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
991 return value; 991 return value;
992 } 992 }
993 } 993 }
994 // 994 //
995 else if (particleDefinition == G4Proton::Pro 995 else if (particleDefinition == G4Proton::ProtonDefinition())
996 { 996 {
997 // k should be in eV 997 // k should be in eV
998 auto k2 = std::upper_bound(pTdummyVec.begi 998 auto k2 = std::upper_bound(pTdummyVec.begin(),
999 pTdummyVec.end(), 999 pTdummyVec.end(),
1000 k); 1000 k);
1001 1001
1002 auto k1 = k2 - 1; 1002 auto k1 = k2 - 1;
1003 1003
1004 /* 1004 /*
1005 G4cout << "----> k=" << k 1005 G4cout << "----> k=" << k
1006 << " " << *k1 1006 << " " << *k1
1007 << " " << *k2 1007 << " " << *k2
1008 << " " << random 1008 << " " << random
1009 << " " << ionizationLevelIndex 1009 << " " << ionizationLevelIndex
1010 << " " << pProbaShellMap[ionizationLevel 1010 << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1011 << " " << pProbaShellMap[ionizationLevel 1011 << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1012 << G4endl; 1012 << G4endl;
1013 */ 1013 */
1014 1014
1015 // SI : the following condition avoids si 1015 // SI : the following condition avoids situations where random > last vector element,
1016 // for eg. when the last element is 1016 // for eg. when the last element is zero
1017 if (random <= pProbaShellMap[ionizationLe 1017 if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1018 && random <= pProbaShellMap[ionizatio 1018 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1019 { 1019 {
1020 auto prob12 = 1020 auto prob12 =
1021 std::upper_bound(pProbaShellMap[ion 1021 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1022 pProbaShellMap[ion 1022 pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1023 random); 1023 random);
1024 1024
1025 auto prob11 = prob12 - 1; 1025 auto prob11 = prob12 - 1;
1026 1026
1027 auto prob22 = 1027 auto prob22 =
1028 std::upper_bound(pProbaShellMap[ionizationL 1028 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1029 pProbaShellMap[ionizationLevelIndex][( 1029 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1030 random); 1030 random);
1031 1031
1032 auto prob21 = prob22 - 1; 1032 auto prob21 = prob22 - 1;
1033 1033
1034 valueK1 = *k1; 1034 valueK1 = *k1;
1035 valueK2 = *k2; 1035 valueK2 = *k2;
1036 valuePROB21 = *prob21; 1036 valuePROB21 = *prob21;
1037 valuePROB22 = *prob22; 1037 valuePROB22 = *prob22;
1038 valuePROB12 = *prob12; 1038 valuePROB12 = *prob12;
1039 valuePROB11 = *prob11; 1039 valuePROB11 = *prob11;
1040 1040
1041 // The following condition avoid gettin 1041 // The following condition avoid getting transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
1042 if(valuePROB11 == 0) nrjTransf11 = bind 1042 if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
1043 else nrjTransf11 = pNrjTransfData[ioniz 1043 else nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1044 if(valuePROB12 == 1) nrjTransf12 = maxi 1044 if(valuePROB12 == 1) nrjTransf12 = maximumEnergyTransferP;
1045 else nrjTransf12 = pNrjTransfData[ioniz 1045 else nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1046 if(valuePROB21 == 0) nrjTransf21 = bind 1046 if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
1047 else nrjTransf21 = pNrjTransfData[ioniz 1047 else nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1048 if(valuePROB22 == 1) nrjTransf22 = maxi 1048 if(valuePROB22 == 1) nrjTransf22 = maximumEnergyTransferP;
1049 else nrjTransf22 = pNrjTransfData[ioniz 1049 else nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1050 1050
1051 } 1051 }
1052 1052
1053 // Avoids cases where cum xs is zero for 1053 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1054 if (random > pProbaShellMap[ionizationLev 1054 if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1055 { 1055 {
1056 auto prob22 = 1056 auto prob22 =
1057 std::upper_bound(pProbaShellMap[ion 1057 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1058 pProbaShellMap[ion 1058 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1059 random); 1059 random);
1060 1060
1061 auto prob21 = prob22 - 1; 1061 auto prob21 = prob22 - 1;
1062 1062
1063 valueK1 = *k1; 1063 valueK1 = *k1;
1064 valueK2 = *k2; 1064 valueK2 = *k2;
1065 valuePROB21 = *prob21; 1065 valuePROB21 = *prob21;
1066 valuePROB22 = *prob22; 1066 valuePROB22 = *prob22;
1067 1067
1068 nrjTransf21 = pNrjTransfData[ionization 1068 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1069 nrjTransf22 = pNrjTransfData[ionization 1069 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1070 1070
1071 G4double interpolatedvalue2 = Interpola 1071 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1072 1072 valuePROB22,
1073 1073 random,
1074 1074 nrjTransf21,
1075 1075 nrjTransf22);
1076 1076
1077 // zeros are explicitly set 1077 // zeros are explicitly set
1078 G4double value = Interpolate(valueK1, v 1078 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1079 1079
1080 return value; 1080 return value;
1081 } 1081 }
1082 } 1082 }
1083 // End electron and proton cases 1083 // End electron and proton cases
1084 1084
1085 G4double nrjTransfProduct = nrjTransf11 * n 1085 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1086 * nrjTransf22; 1086 * nrjTransf22;
1087 1087
1088 if (nrjTransfProduct != 0.) 1088 if (nrjTransfProduct != 0.)
1089 { 1089 {
1090 nrj = QuadInterpolator(valuePROB11, 1090 nrj = QuadInterpolator(valuePROB11,
1091 valuePROB12, 1091 valuePROB12,
1092 valuePROB21, 1092 valuePROB21,
1093 valuePROB22, 1093 valuePROB22,
1094 nrjTransf11, 1094 nrjTransf11,
1095 nrjTransf12, 1095 nrjTransf12,
1096 nrjTransf21, 1096 nrjTransf21,
1097 nrjTransf22, 1097 nrjTransf22,
1098 valueK1, 1098 valueK1,
1099 valueK2, 1099 valueK2,
1100 k, 1100 k,
1101 random); 1101 random);
1102 } 1102 }
1103 1103
1104 return nrj; 1104 return nrj;
1105 } 1105 }
1106 1106
1107 1107
1108 1108