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