| 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" <<
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" <<
62 53
63 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 55
65 using namespace std; 56 using namespace std;
66 57
67 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 59
69 G4MicroElecInelasticModel::G4MicroElecInelasti 60 G4MicroElecInelasticModel::G4MicroElecInelasticModel(const G4ParticleDefinition*,
70 << 61 const G4String& nam)
71 :G4VEmModel(nam),isInitialised(false) << 62 :G4VEmModel(nam),fAtomDeexcitation(0),isInitialised(false)
72 { 63 {
73 nistSi = G4NistManager::Instance()->FindOrBu 64 nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
74 65
75 verboseLevel= 0; 66 verboseLevel= 0;
76 // Verbosity scale: 67 // Verbosity scale:
77 // 0 = nothing << 68 // 0 = nothing
78 // 1 = warning for energy non-conservation << 69 // 1 = warning for energy non-conservation
79 // 2 = details of energy budget 70 // 2 = details of energy budget
80 // 3 = calculation of cross sections, file o 71 // 3 = calculation of cross sections, file openings, sampling of atoms
81 // 4 = entering in methods 72 // 4 = entering in methods
82 73
83 if( verboseLevel>0 ) << 74 if( verboseLevel>0 )
84 { << 75 {
85 G4cout << "MicroElec inelastic model is co 76 G4cout << "MicroElec inelastic model is constructed " << G4endl;
86 } 77 }
87 << 78
88 //Mark this model as "applicable" for atomic 79 //Mark this model as "applicable" for atomic deexcitation
89 SetDeexcitationFlag(true); 80 SetDeexcitationFlag(true);
90 fAtomDeexcitation = nullptr; << 81 fParticleChangeForGamma = 0;
91 fParticleChangeForGamma = nullptr; <<
92 <<
93 // default generator <<
94 SetAngularDistribution(new G4DeltaAngle()); <<
95 <<
96 // Selection of computation method <<
97 fasterCode = true; //false; <<
98 } 82 }
99 83
100 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 85
102 G4MicroElecInelasticModel::~G4MicroElecInelast 86 G4MicroElecInelasticModel::~G4MicroElecInelasticModel()
103 { << 87 {
104 // Cross section 88 // Cross section
105 for (auto & pos : tableData) << 89
>> 90 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 91 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
106 { 92 {
107 G4MicroElecCrossSectionDataSet* table = po << 93 G4MicroElecCrossSectionDataSet* table = pos->second;
108 delete table; 94 delete table;
109 } 95 }
110 96
111 // Final state 97 // Final state
>> 98
112 eVecm.clear(); 99 eVecm.clear();
113 pVecm.clear(); 100 pVecm.clear();
114 << 101
115 } 102 }
116 103
117 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 105
119 void G4MicroElecInelasticModel::Initialise(con 106 void G4MicroElecInelasticModel::Initialise(const G4ParticleDefinition* particle,
120 con << 107 const G4DataVector& /*cuts*/)
121 { 108 {
122 << 109
123 if (verboseLevel > 3) 110 if (verboseLevel > 3)
124 G4cout << "Calling G4MicroElecInelasticMod 111 G4cout << "Calling G4MicroElecInelasticModel::Initialise()" << G4endl;
125 << 112
126 // Energy limits 113 // Energy limits
127 << 114
128 G4String fileElectron("microelec/sigma_inela 115 G4String fileElectron("microelec/sigma_inelastic_e_Si");
129 G4String fileProton("microelec/sigma_inelast 116 G4String fileProton("microelec/sigma_inelastic_p_Si");
130 << 117
131 G4ParticleDefinition* electronDef = G4Electr 118 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
132 G4ParticleDefinition* protonDef = G4Proton:: 119 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
133 G4String electron = electronDef->GetParticle << 120
134 G4String proton = protonDef->GetParticleName << 121 G4String electron;
>> 122 G4String proton;
135 123
136 G4double scaleFactor = 1e-18 * cm *cm; 124 G4double scaleFactor = 1e-18 * cm *cm;
137 125
138 const char* path = G4FindDataDir("G4LEDATA") << 126 char *path = getenv("G4LEDATA");
139 127
140 // *** ELECTRON 128 // *** ELECTRON
141 tableFile[electron] = fileElectron; << 129 electron = electronDef->GetParticleName();
142 lowEnergyLimit[electron] = 16.7 * eV; <<
143 highEnergyLimit[electron] = 100.0 * MeV; <<
144 <<
145 // Cross section <<
146 G4MicroElecCrossSectionDataSet* tableE = new <<
147 tableE->LoadData(fileElectron); <<
148 <<
149 tableData[electron] = tableE; <<
150 <<
151 // Final state <<
152 <<
153 std::ostringstream eFullFileName; <<
154 <<
155 if (fasterCode) eFullFileName << path << "/m <<
156 else eFullFileName << path << "/microelec/si <<
157 <<
158 std::ifstream eDiffCrossSection(eFullFileNam <<
159 <<
160 if (!eDiffCrossSection) <<
161 { <<
162 if (fasterCode) G4Exception("G4MicroElecI <<
163 FatalException,"Missing data file:/microe <<
164 130
165 else G4Exception("G4MicroElecInelasticMod << 131 tableFile[electron] = fileElectron;
166 FatalException,"Missing data file:/microe <<
167 } <<
168 132
169 // Clear the arrays for re-initialization ca << 133 lowEnergyLimit[electron] = 16.7 * eV;
170 // Octobre 22nd, 2014 - Melanie Raine << 134 highEnergyLimit[electron] = 100.0 * MeV;
171 eTdummyVec.clear(); <<
172 pTdummyVec.clear(); <<
173 <<
174 eVecm.clear(); <<
175 pVecm.clear(); <<
176 <<
177 for (int j=0; j<6; j++) <<
178 { <<
179 eProbaShellMap[j].clear(); <<
180 pProbaShellMap[j].clear(); <<
181 135
182 eDiffCrossSectionData[j].clear(); << 136 // Cross section
183 pDiffCrossSectionData[j].clear(); <<
184 <<
185 eNrjTransfData[j].clear(); <<
186 pNrjTransfData[j].clear(); <<
187 } <<
188 <<
189 // <<
190 eTdummyVec.push_back(0.); <<
191 while(!eDiffCrossSection.eof()) <<
192 { <<
193 G4double tDummy; <<
194 G4double eDummy; <<
195 eDiffCrossSection>>tDummy>>eDummy; <<
196 if (tDummy != eTdummyVec.back()) eTdummyVe <<
197 137
198 G4double tmp; << 138 G4MicroElecCrossSectionDataSet* tableE = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
199 for (int j=0; j<6; j++) << 139 tableE->LoadData(fileElectron);
200 { <<
201 eDiffCrossSection>> tmp; <<
202 140
203 eDiffCrossSectionData[j][tDummy][eDummy] << 141 tableData[electron] = tableE;
>> 142
>> 143 // Final state
>> 144
>> 145 std::ostringstream eFullFileName;
>> 146 eFullFileName << path << "/microelec/sigmadiff_inelastic_e_Si.dat";
>> 147 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
>> 148
>> 149 if (!eDiffCrossSection)
>> 150 {
>> 151 G4Exception("G4MicroElecInelasticModel::Initialise","em0003",FatalException,"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
>> 152 }
204 153
205 if (fasterCode) << 154 eTdummyVec.push_back(0.);
206 { << 155 while(!eDiffCrossSection.eof())
207 eNrjTransfData[j][tDummy][eDiffCrossSe << 156 {
208 eProbaShellMap[j][tDummy].push_back(eD << 157 double tDummy;
209 } << 158 double eDummy;
210 else << 159 eDiffCrossSection>>tDummy>>eDummy;
>> 160 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
>> 161 for (int j=0; j<6; j++)
211 { 162 {
212 // SI - only if eof is not reached ! << 163 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
213 if (!eDiffCrossSection.eof()) eDiffCrossSect << 164
214 eVecm[tDummy].push_back(eDummy); << 165 // SI - only if eof is not reached !
>> 166 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
>> 167
>> 168 eVecm[tDummy].push_back(eDummy);
>> 169
215 } 170 }
216 <<
217 } 171 }
218 } << 172 //
219 // <<
220 <<
221 // *** PROTON <<
222 tableFile[proton] = fileProton; <<
223 lowEnergyLimit[proton] = 50. * keV; <<
224 highEnergyLimit[proton] = 10. * GeV; <<
225 <<
226 // Cross section <<
227 G4MicroElecCrossSectionDataSet* tableP = new <<
228 tableP->LoadData(fileProton); <<
229 tableData[proton] = tableP; <<
230 <<
231 // Final state <<
232 std::ostringstream pFullFileName; <<
233 <<
234 if (fasterCode) pFullFileName << path << "/m <<
235 else pFullFileName << path << "/microelec/si <<
236 173
237 std::ifstream pDiffCrossSection(pFullFileNam << 174 // *** PROTON
238 <<
239 if (!pDiffCrossSection) <<
240 { <<
241 if (fasterCode) G4Exception("G4MicroElecIn <<
242 FatalException,"Missing data file:/micro <<
243 175
244 else G4Exception("G4MicroElecInelasticMode << 176 proton = protonDef->GetParticleName();
245 FatalException,"Missing data file:/micro << 177
246 } << 178 tableFile[proton] = fileProton;
247 << 179
248 pTdummyVec.push_back(0.); << 180 lowEnergyLimit[proton] = 50. * keV;
249 while(!pDiffCrossSection.eof()) << 181 highEnergyLimit[proton] = 10. * GeV;
250 { << 182
251 G4double tDummy; << 183 // Cross section
252 G4double eDummy; << 184
253 pDiffCrossSection>>tDummy>>eDummy; << 185 G4MicroElecCrossSectionDataSet* tableP = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
254 if (tDummy != pTdummyVec.back()) pTdummyVe << 186 tableP->LoadData(fileProton);
255 for (int j=0; j<6; j++) <<
256 { <<
257 pDiffCrossSection>>pDiffCrossSectionData <<
258 187
259 if (fasterCode) << 188 tableData[proton] = tableP;
260 { << 189
261 pNrjTransfData[j][tDummy][pDiffCrossSe << 190 // Final state
262 pProbaShellMap[j][tDummy].push_back(pD << 191
263 } << 192 std::ostringstream pFullFileName;
264 else << 193 pFullFileName << path << "/microelec/sigmadiff_inelastic_p_Si.dat";
>> 194 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
>> 195
>> 196 if (!pDiffCrossSection)
>> 197 {
>> 198 G4Exception("G4MicroElecInelasticModel::Initialise","em0003",FatalException,"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
>> 199 }
>> 200
>> 201 pTdummyVec.push_back(0.);
>> 202 while(!pDiffCrossSection.eof())
>> 203 {
>> 204 double tDummy;
>> 205 double eDummy;
>> 206 pDiffCrossSection>>tDummy>>eDummy;
>> 207 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
>> 208 for (int j=0; j<6; j++)
265 { 209 {
266 // SI - only if eof is not reached ! << 210 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
267 if (!pDiffCrossSection.eof()) pDiffCrossSect << 211
268 pVecm[tDummy].push_back(eDummy); << 212 // SI - only if eof is not reached !
>> 213 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
>> 214
>> 215 pVecm[tDummy].push_back(eDummy);
269 } 216 }
270 } 217 }
271 } <<
272 218
273 if (particle==electronDef) << 219
>> 220 if (particle==electronDef)
274 { 221 {
275 SetLowEnergyLimit(lowEnergyLimit[electron] 222 SetLowEnergyLimit(lowEnergyLimit[electron]);
276 SetHighEnergyLimit(highEnergyLimit[electro 223 SetHighEnergyLimit(highEnergyLimit[electron]);
277 } 224 }
278 << 225
279 if (particle==protonDef) << 226 if (particle==protonDef)
280 { 227 {
281 SetLowEnergyLimit(lowEnergyLimit[proton]); 228 SetLowEnergyLimit(lowEnergyLimit[proton]);
282 SetHighEnergyLimit(highEnergyLimit[proton] 229 SetHighEnergyLimit(highEnergyLimit[proton]);
283 } 230 }
284 << 231
285 if( verboseLevel>0 ) << 232 if( verboseLevel>0 )
286 { << 233 {
287 G4cout << "MicroElec Inelastic model is in 234 G4cout << "MicroElec Inelastic model is initialized " << G4endl
288 << "Energy range: " << 235 << "Energy range: "
289 << LowEnergyLimit() / keV << " keV - " << 236 << LowEnergyLimit() / eV << " eV - "
290 << HighEnergyLimit() / MeV << " MeV for " << 237 << HighEnergyLimit() / keV << " keV for "
291 << particle->GetParticleName() << 238 << particle->GetParticleName()
292 << " with mass (amu) " << particle->GetPD 239 << " with mass (amu) " << particle->GetPDGMass()/proton_mass_c2
293 << " and charge " << particle->GetPDGChar 240 << " and charge " << particle->GetPDGCharge()
294 << G4endl << G4endl ; << 241 << G4endl << G4endl ;
295 } 242 }
296 243
297 // 244 //
298 fAtomDeexcitation = G4LossTableManager::Ins << 245
299 << 246 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
>> 247
300 if (isInitialised) { return; } 248 if (isInitialised) { return; }
301 fParticleChangeForGamma = GetParticleChangeF 249 fParticleChangeForGamma = GetParticleChangeForGamma();
302 isInitialised = true; 250 isInitialised = true;
>> 251
303 } 252 }
304 253
305 //....oooOO0OOooo........oooOO0OOooo........oo 254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
306 255
307 G4double G4MicroElecInelasticModel::CrossSecti 256 G4double G4MicroElecInelasticModel::CrossSectionPerVolume(const G4Material* material,
308 << 257 const G4ParticleDefinition* particleDefinition,
309 << 258 G4double ekin,
310 << 259 G4double,
311 << 260 G4double)
312 { 261 {
313 if (verboseLevel > 3) 262 if (verboseLevel > 3)
314 G4cout << "Calling CrossSectionPerVolume() 263 G4cout << "Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" << G4endl;
315 << 264
316 G4double density = material->GetTotNbOfAtoms 265 G4double density = material->GetTotNbOfAtomsPerVolume();
>> 266
>> 267 /* if (
>> 268 particleDefinition != G4Proton::ProtonDefinition()
>> 269 &&
>> 270 particleDefinition != G4Electron::ElectronDefinition()
>> 271 &&
>> 272 particleDefinition != G4GenericIon::GenericIonDefinition()
>> 273 )
>> 274
>> 275 return 0;*/
317 276
318 // Calculate total cross section for model << 277 // Calculate total cross section for model
>> 278
319 G4double lowLim = 0; 279 G4double lowLim = 0;
320 G4double highLim = 0; 280 G4double highLim = 0;
321 G4double sigma=0; 281 G4double sigma=0;
322 << 282
323 const G4String& particleName = particleDefin 283 const G4String& particleName = particleDefinition->GetParticleName();
324 G4String nameLocal = particleName ; 284 G4String nameLocal = particleName ;
325 << 285
326 G4double Zeff2 = 1.0; 286 G4double Zeff2 = 1.0;
327 G4double Mion_c2 = particleDefinition->GetPD 287 G4double Mion_c2 = particleDefinition->GetPDGMass();
328 << 288
329 if (Mion_c2 > proton_mass_c2) 289 if (Mion_c2 > proton_mass_c2)
330 { 290 {
331 G4ionEffectiveCharge EffCharge ; << 291 G4ionEffectiveCharge EffCharge ;
332 G4double Zeff = EffCharge.EffectiveCharge( 292 G4double Zeff = EffCharge.EffectiveCharge(particleDefinition, material,ekin);
333 Zeff2 = Zeff*Zeff; 293 Zeff2 = Zeff*Zeff;
334 << 294
335 if (verboseLevel > 3) << 295 if (verboseLevel > 3)
336 G4cout << "Before scaling : " << G4endl << 296 G4cout << "Before scaling : " << G4endl
337 << "Particle : " << nameLocal << ", mass << 297 << "Particle : " << nameLocal << ", mass : " << Mion_c2/proton_mass_c2 << "*mp, charge " << Zeff
338 << ", Ekin (eV) = " << ekin/eV << G4endl << 298 << ", Ekin (eV) = " << ekin/eV << G4endl ;
339 << 299
340 ekin *= proton_mass_c2/Mion_c2 ; 300 ekin *= proton_mass_c2/Mion_c2 ;
341 nameLocal = "proton" ; 301 nameLocal = "proton" ;
342 << 302
343 if (verboseLevel > 3) << 303 if (verboseLevel > 3)
344 G4cout << "After scaling : " << G4endl << 304 G4cout << "After scaling : " << G4endl
345 << "Particle : " << nameLocal << ", Eki << 305 << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl ;
346 } 306 }
347 << 307
348 if (material == nistSi || material->GetBaseM 308 if (material == nistSi || material->GetBaseMaterial() == nistSi)
349 { << 309 {
350 auto pos1 = lowEnergyLimit.find(nameLocal) << 310
>> 311 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
>> 312 pos1 = lowEnergyLimit.find(nameLocal);
351 if (pos1 != lowEnergyLimit.end()) 313 if (pos1 != lowEnergyLimit.end())
352 { 314 {
353 lowLim = pos1->second; 315 lowLim = pos1->second;
354 } 316 }
355 << 317
356 auto pos2 = highEnergyLimit.find(nameLocal << 318 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
>> 319 pos2 = highEnergyLimit.find(nameLocal);
357 if (pos2 != highEnergyLimit.end()) 320 if (pos2 != highEnergyLimit.end())
358 { 321 {
359 highLim = pos2->second; 322 highLim = pos2->second;
360 } 323 }
361 << 324
362 if (ekin >= lowLim && ekin < highLim) 325 if (ekin >= lowLim && ekin < highLim)
363 { 326 {
364 auto pos = tableData.find(nameLocal); << 327 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 328 pos = tableData.find(nameLocal);
>> 329
365 if (pos != tableData.end()) 330 if (pos != tableData.end())
366 { 331 {
367 G4MicroElecCrossSectionDataSet* table 332 G4MicroElecCrossSectionDataSet* table = pos->second;
368 if (table != nullptr) << 333 if (table != 0)
369 { 334 {
370 sigma = table->FindValue(ekin); << 335 sigma = table->FindValue(ekin);
371 } 336 }
372 } 337 }
373 else 338 else
374 { 339 {
375 G4Exception("G4MicroElecInelasticModel << 340 G4Exception("G4MicroElecInelasticModel::CrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
376 FatalException,"Model not applicable t <<
377 } 341 }
378 } 342 }
379 else << 343 else
380 { 344 {
381 if (nameLocal!="e-") << 345 if (nameLocal!="e-")
382 { << 346 {
383 // G4cout << "Particle : " << nameLoca << 347 // G4cout << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl;
384 // G4cout << "### Warning: particle en << 348 // G4cout << "### Warning: particle energy out of bounds! ###" << G4endl;
385 } << 349 }
386 } 350 }
387 << 351
388 if (verboseLevel > 3) 352 if (verboseLevel > 3)
389 { 353 {
390 G4cout << "---> Kinetic energy (eV)=" << 354 G4cout << "---> Kinetic energy (eV)=" << ekin/eV << G4endl;
391 G4cout << " - Cross section per Si atom 355 G4cout << " - Cross section per Si atom (cm^2)=" << sigma*Zeff2/cm2 << G4endl;
392 G4cout << " - Cross section per Si atom 356 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./cm) << G4endl;
393 } << 357 }
394 << 358
395 } // if (SiMaterial) 359 } // if (SiMaterial)
396 return sigma*density*Zeff2; << 360 return sigma*density*Zeff2;
>> 361
>> 362
397 } 363 }
398 364
399 //....oooOO0OOooo........oooOO0OOooo........oo 365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 366
401 void G4MicroElecInelasticModel::SampleSecondar 367 void G4MicroElecInelasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
402 << 368 const G4MaterialCutsCouple* /*couple*/,
403 << 369 const G4DynamicParticle* particle,
404 << 370 G4double,
405 << 371 G4double)
406 { 372 {
407 << 373
408 if (verboseLevel > 3) 374 if (verboseLevel > 3)
409 G4cout << "Calling SampleSecondaries() of 375 G4cout << "Calling SampleSecondaries() of G4MicroElecInelasticModel" << G4endl;
410 << 376
411 G4double lowLim = 0; 377 G4double lowLim = 0;
412 G4double highLim = 0; 378 G4double highLim = 0;
413 << 379
414 G4double ekin = particle->GetKineticEnergy() 380 G4double ekin = particle->GetKineticEnergy();
415 G4double k = ekin ; 381 G4double k = ekin ;
416 << 382
417 G4ParticleDefinition* PartDef = particle->Ge 383 G4ParticleDefinition* PartDef = particle->GetDefinition();
418 const G4String& particleName = PartDef->GetP 384 const G4String& particleName = PartDef->GetParticleName();
419 G4String nameLocal2 = particleName ; 385 G4String nameLocal2 = particleName ;
420 G4double particleMass = particle->GetDefinit 386 G4double particleMass = particle->GetDefinition()->GetPDGMass();
421 << 387
422 if (particleMass > proton_mass_c2) 388 if (particleMass > proton_mass_c2)
423 { 389 {
424 k *= proton_mass_c2/particleMass ; 390 k *= proton_mass_c2/particleMass ;
425 PartDef = G4Proton::ProtonDefinition(); 391 PartDef = G4Proton::ProtonDefinition();
426 nameLocal2 = "proton" ; 392 nameLocal2 = "proton" ;
427 } 393 }
428 << 394
429 auto pos1 = lowEnergyLimit.find(nameLocal2); << 395 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
>> 396 pos1 = lowEnergyLimit.find(nameLocal2);
>> 397
430 if (pos1 != lowEnergyLimit.end()) 398 if (pos1 != lowEnergyLimit.end())
431 { 399 {
432 lowLim = pos1->second; 400 lowLim = pos1->second;
433 } 401 }
434 << 402
435 auto pos2 = highEnergyLimit.find(nameLocal2) << 403 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
>> 404 pos2 = highEnergyLimit.find(nameLocal2);
>> 405
436 if (pos2 != highEnergyLimit.end()) 406 if (pos2 != highEnergyLimit.end())
437 { 407 {
438 highLim = pos2->second; 408 highLim = pos2->second;
439 } 409 }
440 << 410
441 if (k >= lowLim && k < highLim) 411 if (k >= lowLim && k < highLim)
442 { 412 {
443 G4ParticleMomentum primaryDirection = part 413 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
444 G4double totalEnergy = ekin + particleMass 414 G4double totalEnergy = ekin + particleMass;
445 G4double pSquare = ekin * (totalEnergy + p 415 G4double pSquare = ekin * (totalEnergy + particleMass);
446 G4double totalMomentum = std::sqrt(pSquare 416 G4double totalMomentum = std::sqrt(pSquare);
447 << 417
448 G4int Shell = 0; << 418 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 419 G4double bindingEnergy = SiStructure.Energy(Shell);
459 <<
460 if (verboseLevel > 3) 420 if (verboseLevel > 3)
461 { 421 {
462 G4cout << "---> Kinetic energy (eV)=" << << 422 G4cout << "---> Kinetic energy (eV)=" << k/eV << G4endl ;
463 G4cout << "Shell: " << Shell << ", energ << 423 G4cout << "Shell: " << Shell << ", energy: " << bindingEnergy/eV << G4endl;
464 } 424 }
465 << 425
466 // sample deexcitation << 426 // sample deexcitation
467 std::size_t secNumberInit = 0; // need to << 427
468 std::size_t secNumberFinal = 0; // So I'll << 428 G4int secNumberInit = 0; // need to know at a certain point the energy of secondaries
469 << 429 G4int secNumberFinal = 0; // So I'll make the difference and then sum the energies
470 //SI: additional protection if tcs interpo << 430
471 if (k<bindingEnergy) return; <<
472 <<
473 G4int Z = 14; <<
474 <<
475 if(fAtomDeexcitation && Shell > 2) { 431 if(fAtomDeexcitation && Shell > 2) {
476 << 432 G4int Z = 14;
477 G4AtomicShellEnumerator as = fKShell; 433 G4AtomicShellEnumerator as = fKShell;
478 << 434
479 if (Shell == 4) << 435 if (Shell == 4)
480 { << 436 {
481 as = G4AtomicShellEnumerator(1); << 437 as = G4AtomicShellEnumerator(1);
482 } << 438 }
483 else if (Shell == 3) 439 else if (Shell == 3)
484 { << 440 {
485 as = G4AtomicShellEnumerator(3); << 441 as = G4AtomicShellEnumerator(3);
486 } << 442 }
487 << 443
488 const G4AtomicShell* shell = fAtomDeexci << 444 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
489 secNumberInit = fvect->size(); 445 secNumberInit = fvect->size();
490 fAtomDeexcitation->GenerateParticles(fve 446 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
491 secNumberFinal = fvect->size(); 447 secNumberFinal = fvect->size();
492 } 448 }
493 <<
494 G4double secondaryKinetic=-1000*eV; <<
495 449
496 if (!fasterCode) << 450 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef,k,Shell);
497 { << 451
498 secondaryKinetic = RandomizeEjectedElect <<
499 } <<
500 else <<
501 { <<
502 secondaryKinetic = RandomizeEjectedElect <<
503 } <<
504 <<
505 if (verboseLevel > 3) 452 if (verboseLevel > 3)
506 { 453 {
507 G4cout << "Ionisation process" << G4endl << 454 G4cout << "Ionisation process" << G4endl;
508 G4cout << "Shell: " << Shell << " Kin. e << 455 G4cout << "Shell: " << Shell << " Kin. energy (eV)=" << k/eV
509 << " Sec. energy (eV)=" << secondaryKine << 456 << " Sec. energy (eV)=" << secondaryKinetic/eV << G4endl;
510 } << 457 }
511 << 458
512 G4ThreeVector deltaDirection = << 459 G4double cosTheta = 0.;
513 GetAngularDistribution()->SampleDirectionF << 460 G4double phi = 0.;
514 << 461 RandomizeEjectedElectronDirection(PartDef, k, secondaryKinetic, cosTheta, phi);
515 << 462
516 << 463 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
517 if (particle->GetDefinition() == G4Electro << 464 G4double dirX = sinTheta*std::cos(phi);
518 { << 465 G4double dirY = sinTheta*std::sin(phi);
519 G4double deltaTotalMomentum = std::sqrt( << 466 G4double dirZ = cosTheta;
>> 467 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
>> 468 deltaDirection.rotateUz(primaryDirection);
>> 469
>> 470 //if (particle->GetDefinition() == G4Electron::ElectronDefinition())
>> 471 //{
>> 472 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
>> 473
520 G4double finalPx = totalMomentum*primary 474 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
521 G4double finalPy = totalMomentum*primary 475 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
522 G4double finalPz = totalMomentum*primary 476 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
523 G4double finalMomentum = std::sqrt(final 477 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
524 finalPx /= finalMomentum; 478 finalPx /= finalMomentum;
525 finalPy /= finalMomentum; 479 finalPy /= finalMomentum;
526 finalPz /= finalMomentum; 480 finalPz /= finalMomentum;
527 << 481
528 G4ThreeVector direction; 482 G4ThreeVector direction;
529 direction.set(finalPx,finalPy,finalPz); 483 direction.set(finalPx,finalPy,finalPz);
530 <<
531 fParticleChangeForGamma->ProposeMomentum <<
532 } <<
533 else <<
534 fParticleChangeForGamma->ProposeMomentum <<
535 484
>> 485 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
>> 486 //}
>> 487 //else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
>> 488
536 // note that secondaryKinetic is the energ 489 // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
537 G4double deexSecEnergy = 0; 490 G4double deexSecEnergy = 0;
538 for (std::size_t j=secNumberInit; j < secN << 491 for (G4int j=secNumberInit; j < secNumberFinal; j++) {
539 deexSecEnergy = deexSecEnergy + (*fvect) << 492 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
540 << 493
541 fParticleChangeForGamma->SetProposedKineti 494 fParticleChangeForGamma->SetProposedKineticEnergy(ekin-bindingEnergy-secondaryKinetic);
542 fParticleChangeForGamma->ProposeLocalEnerg 495 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy-deexSecEnergy);
543 << 496
544 if (secondaryKinetic>0) << 497 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
545 { << 498 fvect->push_back(dp);
546 G4DynamicParticle* dp = new G4DynamicPar << 499
547 fvect->push_back(dp); <<
548 } <<
549 } 500 }
>> 501
550 } 502 }
551 503
552 //....oooOO0OOooo........oooOO0OOooo........oo 504 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553 505
554 G4double G4MicroElecInelasticModel::RandomizeE << 506 G4double G4MicroElecInelasticModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
555 << 507 G4double k, G4int shell)
556 { 508 {
557 if (particleDefinition == G4Electron::Electr << 509 if (particleDefinition == G4Electron::ElectronDefinition())
558 { 510 {
559 G4double maximumEnergyTransfer=0.; 511 G4double maximumEnergyTransfer=0.;
560 if ((k+SiStructure.Energy(shell))/2. > k) 512 if ((k+SiStructure.Energy(shell))/2. > k) maximumEnergyTransfer=k;
561 else maximumEnergyTransfer = (k+SiStructur 513 else maximumEnergyTransfer = (k+SiStructure.Energy(shell))/2.;
562 << 514
563 G4double crossSectionMaximum = 0.; 515 G4double crossSectionMaximum = 0.;
564 << 516
565 G4double minEnergy = SiStructure.Energy(sh 517 G4double minEnergy = SiStructure.Energy(shell);
566 G4double maxEnergy = maximumEnergyTransfer 518 G4double maxEnergy = maximumEnergyTransfer;
567 G4int nEnergySteps = 100; 519 G4int nEnergySteps = 100;
568 << 520
569 G4double value(minEnergy); 521 G4double value(minEnergy);
570 G4double stpEnergy(std::pow(maxEnergy/valu 522 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
571 G4int step(nEnergySteps); 523 G4int step(nEnergySteps);
572 while (step>0) 524 while (step>0)
573 { 525 {
574 step--; 526 step--;
575 G4double differentialCrossSection = Diff 527 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
576 if(differentialCrossSection >= crossSect 528 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
577 value*=stpEnergy; 529 value*=stpEnergy;
578 } 530 }
579 << 531
>> 532
580 G4double secondaryElectronKineticEnergy=0. 533 G4double secondaryElectronKineticEnergy=0.;
581 do << 534 do
582 { 535 {
583 secondaryElectronKineticEnergy = G4Unifo 536 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
584 } while(G4UniformRand()*crossSectionMaximu 537 } while(G4UniformRand()*crossSectionMaximum >
585 DifferentialCrossSection(particleD << 538 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
586 << 539
587 return secondaryElectronKineticEnergy; << 540 return secondaryElectronKineticEnergy;
>> 541
588 } 542 }
589 << 543
590 if (particleDefinition == G4Proton::ProtonDe << 544 if (particleDefinition == G4Proton::ProtonDefinition())
591 { 545 {
592 G4double maximumEnergyTransfer = 4.* (elec 546 G4double maximumEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
593 G4double crossSectionMaximum = 0.; 547 G4double crossSectionMaximum = 0.;
594 << 548
595 G4double minEnergy = SiStructure.Energy(sh 549 G4double minEnergy = SiStructure.Energy(shell);
596 G4double maxEnergy = maximumEnergyTransfer << 550 G4double maxEnergy = maximumEnergyTransfer;
597 G4int nEnergySteps = 100; 551 G4int nEnergySteps = 100;
598 << 552
599 G4double value(minEnergy); 553 G4double value(minEnergy);
600 G4double stpEnergy(std::pow(maxEnergy/valu 554 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
601 G4int step(nEnergySteps); 555 G4int step(nEnergySteps);
602 while (step>0) 556 while (step>0)
603 { 557 {
604 step--; << 558 step--;
605 G4double differentialCrossSection = Diff 559 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
606 if(differentialCrossSection >= crossSect 560 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
607 value*=stpEnergy; 561 value*=stpEnergy;
608 } 562 }
609 <<
610 G4double secondaryElectronKineticEnergy = 563 G4double secondaryElectronKineticEnergy = 0.;
611 do 564 do
612 { 565 {
613 secondaryElectronKineticEnergy = G4Unifo 566 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
614 << 567
615 } while(G4UniformRand()*crossSectionMaximu << 568 } while(G4UniformRand()*crossSectionMaximum >=
616 DifferentialCrossSection(particleD << 569 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
617 return secondaryElectronKineticEnergy; 570 return secondaryElectronKineticEnergy;
618 } 571 }
619 << 572
620 return 0; 573 return 0;
621 } 574 }
622 575
623 //....oooOO0OOooo........oooOO0OOooo........oo 576 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
624 577
625 // The following section is not used anymore b << 578 void G4MicroElecInelasticModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
626 // GetAngularDistribution()->SampleDirectionFo << 579 G4double k,
627 << 580 G4double secKinetic,
628 /*void G4MicroElecInelasticModel::RandomizeEje << 581 G4double & cosTheta,
629 G4double k, << 582 G4double & phi )
630 G4double secKinetic, << 583 {
631 G4double & cosTheta, << 584 if (particleDefinition == G4Electron::ElectronDefinition())
632 G4double & phi ) << 585 {
633 { << 586 phi = twopi * G4UniformRand();
634 if (particleDefinition == G4Electron::Electro << 587 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
635 { << 588 cosTheta = std::sqrt(1.-sin2O);
636 phi = twopi * G4UniformRand(); << 589 }
637 G4double sin2O = (1.-secKinetic/k) / (1.+secK <<
638 cosTheta = std::sqrt(1.-sin2O); <<
639 } <<
640 <<
641 if (particleDefinition == G4Proton::ProtonDef <<
642 { <<
643 G4double maxSecKinetic = 4.* (electron_mass_c <<
644 phi = twopi * G4UniformRand(); <<
645 cosTheta = std::sqrt(secKinetic / maxSecKinet <<
646 } <<
647 590
648 else << 591 if (particleDefinition == G4Proton::ProtonDefinition())
649 { << 592 {
650 G4double maxSecKinetic = 4.* (electron_mass_c << 593 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
651 phi = twopi * G4UniformRand(); << 594 phi = twopi * G4UniformRand();
652 cosTheta = std::sqrt(secKinetic / maxSecKinet << 595 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
653 } << 596 }
654 } << 597
655 */ << 598 else
>> 599 {
>> 600 G4double maxSecKinetic = 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
>> 601 phi = twopi * G4UniformRand();
>> 602 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
>> 603 }
>> 604 }
656 605
657 //....oooOO0OOooo........oooOO0OOooo........oo 606 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
658 607
659 G4double G4MicroElecInelasticModel::Differenti << 608 double G4MicroElecInelasticModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
660 << 609 G4double k,
661 << 610 G4double energyTransfer,
662 << 611 G4int LevelIndex)
663 { 612 {
664 G4double sigma = 0.; 613 G4double sigma = 0.;
665 << 614
666 if (energyTransfer >= SiStructure.Energy(Lev 615 if (energyTransfer >= SiStructure.Energy(LevelIndex))
667 { 616 {
668 G4double valueT1 = 0; 617 G4double valueT1 = 0;
669 G4double valueT2 = 0; 618 G4double valueT2 = 0;
670 G4double valueE21 = 0; 619 G4double valueE21 = 0;
671 G4double valueE22 = 0; 620 G4double valueE22 = 0;
672 G4double valueE12 = 0; 621 G4double valueE12 = 0;
673 G4double valueE11 = 0; 622 G4double valueE11 = 0;
674 << 623
675 G4double xs11 = 0; << 624 G4double xs11 = 0;
676 G4double xs12 = 0; << 625 G4double xs12 = 0;
677 G4double xs21 = 0; << 626 G4double xs21 = 0;
678 G4double xs22 = 0; << 627 G4double xs22 = 0;
679 << 628
680 if (particleDefinition == G4Electron::Elec << 629 if (particleDefinition == G4Electron::ElectronDefinition())
681 { 630 {
682 // k should be in eV and energy transfer << 631 // k should be in eV and energy transfer eV also
683 auto t2 = std::upper_bound(eTdummyVec.be << 632
684 auto t1 = t2-1; << 633 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
685 // The following condition avoids situat << 634 std::vector<double>::iterator t1 = t2-1;
>> 635 // SI : the following condition avoids situations where energyTransfer >last vector element
686 if (energyTransfer <= eVecm[(*t1)].back( 636 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
687 { 637 {
688 auto e12 = std::upper_bound(eVecm[(*t1 << 638 std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
689 auto e11 = e12-1; << 639 std::vector<double>::iterator e11 = e12-1;
690 auto e22 = std::upper_bound(eVecm[(*t2 << 640
691 auto e21 = e22-1; << 641 std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
692 << 642 std::vector<double>::iterator e21 = e22-1;
>> 643
693 valueT1 =*t1; 644 valueT1 =*t1;
694 valueT2 =*t2; 645 valueT2 =*t2;
695 valueE21 =*e21; 646 valueE21 =*e21;
696 valueE22 =*e22; 647 valueE22 =*e22;
697 valueE12 =*e12; 648 valueE12 =*e12;
698 valueE11 =*e11; 649 valueE11 =*e11;
699 << 650
700 xs11 = eDiffCrossSectionData[LevelInde 651 xs11 = eDiffCrossSectionData[LevelIndex][valueT1][valueE11];
701 xs12 = eDiffCrossSectionData[LevelInde 652 xs12 = eDiffCrossSectionData[LevelIndex][valueT1][valueE12];
702 xs21 = eDiffCrossSectionData[LevelInde 653 xs21 = eDiffCrossSectionData[LevelIndex][valueT2][valueE21];
703 xs22 = eDiffCrossSectionData[LevelInde 654 xs22 = eDiffCrossSectionData[LevelIndex][valueT2][valueE22];
704 } << 655 }
>> 656
705 } 657 }
706 << 658
707 if (particleDefinition == G4Proton::Proton << 659 if (particleDefinition == G4Proton::ProtonDefinition())
708 { << 660 {
709 // k should be in eV and energy transfer 661 // k should be in eV and energy transfer eV also
710 auto t2 = std::upper_bound(pTdummyVec.be << 662 std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
711 auto t1 = t2-1; << 663 std::vector<double>::iterator t1 = t2-1;
712 if (energyTransfer <= pVecm[(*t1)].back( 664 if (energyTransfer <= pVecm[(*t1)].back() && energyTransfer <= pVecm[(*t2)].back() )
713 { 665 {
714 auto e12 = std::upper_bound(pVecm[(*t1 << 666 std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
715 auto e11 = e12-1; << 667 std::vector<double>::iterator e11 = e12-1;
716 << 668
717 auto e22 = std::upper_bound(pVecm[(*t2 << 669 std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
718 auto e21 = e22-1; << 670 std::vector<double>::iterator e21 = e22-1;
719 << 671
720 valueT1 =*t1; 672 valueT1 =*t1;
721 valueT2 =*t2; 673 valueT2 =*t2;
722 valueE21 =*e21; 674 valueE21 =*e21;
723 valueE22 =*e22; 675 valueE22 =*e22;
724 valueE12 =*e12; 676 valueE12 =*e12;
725 valueE11 =*e11; 677 valueE11 =*e11;
726 << 678
727 xs11 = pDiffCrossSectionData[LevelInde << 679 xs11 = pDiffCrossSectionData[LevelIndex][valueT1][valueE11];
728 xs12 = pDiffCrossSectionData[LevelInde 680 xs12 = pDiffCrossSectionData[LevelIndex][valueT1][valueE12];
729 xs21 = pDiffCrossSectionData[LevelInde 681 xs21 = pDiffCrossSectionData[LevelIndex][valueT2][valueE21];
730 xs22 = pDiffCrossSectionData[LevelInde 682 xs22 = pDiffCrossSectionData[LevelIndex][valueT2][valueE22];
731 } << 683 }
732 } << 684 }
733 << 685
734 sigma = QuadInterpolator( valueE11, va << 686 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
735 valueE21, valueE22, << 687 if (xsProduct != 0.)
736 xs11, xs12, << 688 {
737 xs21, xs22, << 689 sigma = QuadInterpolator( valueE11, valueE12,
738 valueT1, valueT2, << 690 valueE21, valueE22,
739 k, energyTransfer); << 691 xs11, xs12,
740 } << 692 xs21, xs22,
>> 693 valueT1, valueT2,
>> 694 k, energyTransfer);
>> 695 }
>> 696
>> 697 }
>> 698
741 return sigma; 699 return sigma;
742 } 700 }
743 701
744 //....oooOO0OOooo........oooOO0OOooo........oo 702 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
745 703
746 G4double G4MicroElecInelasticModel::Interpolat << 704 G4double G4MicroElecInelasticModel::LogLogInterpolate(G4double e1,
747 << 705 G4double e2,
748 << 706 G4double e,
749 << 707 G4double xs1,
750 << 708 G4double xs2)
751 { << 709 {
752 G4double value = 0.; << 710 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
753 << 711 G4double b = std::log10(xs2) - a*std::log10(e2);
754 // Log-log interpolation by default << 712 G4double sigma = a*std::log10(e) + b;
755 if (e1 != 0 && e2 != 0 && (std::log10(e2) - << 713 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; 714 return value;
783 } 715 }
784 716
785 //....oooOO0OOooo........oooOO0OOooo........oo 717 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
786 718
787 G4double G4MicroElecInelasticModel::QuadInterp << 719 G4double G4MicroElecInelasticModel::QuadInterpolator(G4double e11, G4double e12,
788 << 720 G4double e21, G4double e22,
789 << 721 G4double xs11, G4double xs12,
790 << 722 G4double xs21, G4double xs22,
791 << 723 G4double t1, G4double t2,
792 << 724 G4double t, G4double e)
793 { << 725 {
794 G4double interpolatedvalue1 = Interpolate(e1 << 726 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
795 G4double interpolatedvalue2 = Interpolate(e2 << 727 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
796 G4double value = Interpolate(t1, t2, t, inte << 728 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
797 return value; 729 return value;
798 } 730 }
799 731
800 //....oooOO0OOooo........oooOO0OOooo........oo 732 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
801 733
802 G4int G4MicroElecInelasticModel::RandomSelect( 734 G4int G4MicroElecInelasticModel::RandomSelect(G4double k, const G4String& particle )
803 { << 735 {
804 G4int level = 0; 736 G4int level = 0;
805 << 737
806 auto pos = tableData.find(particle); << 738 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
807 if (pos != tableData.cend()) << 739 pos = tableData.find(particle);
>> 740
>> 741 if (pos != tableData.end())
808 { 742 {
809 G4MicroElecCrossSectionDataSet* table = po 743 G4MicroElecCrossSectionDataSet* table = pos->second;
810 << 744
811 if (table != nullptr) << 745 if (table != 0)
812 { 746 {
813 G4double* valuesBuffer = new G4double[ta 747 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
814 const G4int n = (G4int)table->NumberOfCo << 748 const size_t n(table->NumberOfComponents());
815 G4int i(n); << 749 size_t i(n);
816 G4double value = 0.; 750 G4double value = 0.;
817 << 751
818 while (i>0) 752 while (i>0)
819 { << 753 {
820 --i; << 754 i--;
821 valuesBuffer[i] = table->GetComponent( << 755 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
822 value += valuesBuffer[i]; << 756 value += valuesBuffer[i];
823 } 757 }
824 << 758
825 value *= G4UniformRand(); 759 value *= G4UniformRand();
826 << 760
827 i = n; 761 i = n;
828 << 762
829 while (i > 0) 763 while (i > 0)
830 { 764 {
831 --i; << 765 i--;
832 << 766
833 if (valuesBuffer[i] > value) << 767 if (valuesBuffer[i] > value)
834 { 768 {
835 delete[] valuesBuffer; 769 delete[] valuesBuffer;
836 return i; 770 return i;
837 } 771 }
838 value -= valuesBuffer[i]; << 772 value -= valuesBuffer[i];
839 } 773 }
840 << 774
841 if (valuesBuffer) delete[] valuesBuffer; 775 if (valuesBuffer) delete[] valuesBuffer;
842 << 776
843 } 777 }
844 } 778 }
845 else 779 else
846 { 780 {
847 G4Exception("G4MicroElecInelasticModel::Ra 781 G4Exception("G4MicroElecInelasticModel::RandomSelect","em0002",FatalException,"Model not applicable to particle type.");
848 } 782 }
849 << 783
850 return level; 784 return level;
851 } 785 }
852 786
853 //....oooOO0OOooo........oooOO0OOooo........oo 787 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
854 <<
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 788
1107 789
1108 790