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