| 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 // G4MicroElecElasticModel_new.cc, 2011/08/29 27 // G4MicroElecElasticModel_new.cc, 2011/08/29 A.Valentin, M. Raine are with CEA [a]
28 // 2020/05/20 P. Caro 28 // 2020/05/20 P. Caron, C. Inguimbert are with ONERA [b]
29 // Q. Gibaru is with CEA [a] 29 // Q. Gibaru is with CEA [a], ONERA [b] and CNES [c]
30 // M. Raine and D. Lambert a 30 // M. Raine and D. Lambert are with CEA [a]
31 // 31 //
32 // A part of this work has been funded by the 32 // A part of this work has been funded by the French space agency(CNES[c])
33 // [a] CEA, DAM, DIF - 91297 ARPAJON, France 33 // [a] CEA, DAM, DIF - 91297 ARPAJON, France
34 // [b] ONERA - DPHY, 2 avenue E.Belin, 31055 T 34 // [b] ONERA - DPHY, 2 avenue E.Belin, 31055 Toulouse, France
35 // [c] CNES, 18 av.E.Belin, 31401 Toulouse CED 35 // [c] CNES, 18 av.E.Belin, 31401 Toulouse CEDEX, France
36 // 36 //
37 // Based on the following publications 37 // Based on the following publications
38 // - A.Valentin, M. Raine, 38 // - A.Valentin, M. Raine,
39 // Inelastic cross-sections of low energy e 39 // Inelastic cross-sections of low energy electrons in silicon
40 // for the simulation of heavy ion trac 40 // for the simulation of heavy ion tracks with the Geant4-DNA toolkit,
41 // NSS Conf. Record 2010, pp. 80-85 41 // NSS Conf. Record 2010, pp. 80-85
42 // https://doi.org/10.1109/NSSMIC. 42 // https://doi.org/10.1109/NSSMIC.2010.5873720
43 // 43 //
44 // - A.Valentin, M. Raine, M.Gaillardin, 44 // - A.Valentin, M. Raine, M.Gaillardin, P.Paillet
45 // Geant4 physics processes for microdo 45 // Geant4 physics processes for microdosimetry simulation:
46 // very low energy electromagnetic mode 46 // very low energy electromagnetic models for electrons in Silicon,
47 // https://doi.org/10.1016/j.nimb. 47 // https://doi.org/10.1016/j.nimb.2012.06.007
48 // NIM B, vol. 288, pp. 66-73, 2012, pa 48 // NIM B, vol. 288, pp. 66-73, 2012, part A
49 // heavy ions in Si, NIM B, vol. 287, p 49 // heavy ions in Si, NIM B, vol. 287, pp. 124-129, 2012, part B
50 // https://doi.org/10.1016/j.nimb. 50 // https://doi.org/10.1016/j.nimb.2012.07.028
51 // 51 //
52 // - M. Raine, M. Gaillardin, P. Paillet 52 // - M. Raine, M. Gaillardin, P. Paillet
53 // Geant4 physics processes for silicon 53 // Geant4 physics processes for silicon microdosimetry simulation:
54 // Improvements and extension of the en 54 // Improvements and extension of the energy-range validity up to 10 GeV/nucleon
55 // NIM B, vol. 325, pp. 97-100, 2014 55 // NIM B, vol. 325, pp. 97-100, 2014
56 // https://doi.org/10.1016/j.nimb. 56 // https://doi.org/10.1016/j.nimb.2014.01.014
57 // 57 //
58 // - J. Pierron, C. Inguimbert, M. Belhaj 58 // - J. Pierron, C. Inguimbert, M. Belhaj, T. Gineste, J. Puech, M. Raine
59 // Electron emission yield for low ener 59 // Electron emission yield for low energy electrons:
60 // Monte Carlo simulation and experimen 60 // Monte Carlo simulation and experimental comparison for Al, Ag, and Si
61 // Journal of Applied Physics 121 (2017 61 // Journal of Applied Physics 121 (2017) 215107.
62 // https://doi.org/10.1063/1.498 62 // https://doi.org/10.1063/1.4984761
63 // 63 //
64 // - P. Caron, 64 // - P. Caron,
65 // Study of Electron-Induced Single-Eve 65 // Study of Electron-Induced Single-Event Upset in Integrated Memory Devices
66 // PHD, 16th October 2019 66 // PHD, 16th October 2019
67 // 67 //
68 // - Q.Gibaru, C.Inguimbert, P.Caron, M.Raine 68 // - Q.Gibaru, C.Inguimbert, P.Caron, M.Raine, D.Lambert, J.Puech,
69 // Geant4 physics processes for microdo 69 // Geant4 physics processes for microdosimetry and secondary electron emission simulation :
70 // Extension of MicroElec to very low e 70 // Extension of MicroElec to very low energies and new materials
71 // NIM B, 2020, in review. 71 // NIM B, 2020, in review.
72 // 72 //
73 // 73 //
74 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 #include "G4MicroElecElasticModel_new.hh" 75 #include "G4MicroElecElasticModel_new.hh"
76 #include "G4PhysicalConstants.hh" 76 #include "G4PhysicalConstants.hh"
77 #include "G4SystemOfUnits.hh" 77 #include "G4SystemOfUnits.hh"
78 #include "G4Exp.hh" 78 #include "G4Exp.hh"
79 #include "G4Material.hh" 79 #include "G4Material.hh"
80 #include "G4String.hh" 80 #include "G4String.hh"
81 81
82 82
83 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 84
85 using namespace std; 85 using namespace std;
86 86
87 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 88
89 G4MicroElecElasticModel_new::G4MicroElecElasti 89 G4MicroElecElasticModel_new::G4MicroElecElasticModel_new(const G4ParticleDefinition*,
90 const G4String& nam) 90 const G4String& nam)
91 :G4VEmModel(nam), isInitialised(false) 91 :G4VEmModel(nam), isInitialised(false)
92 { 92 {
93 killBelowEnergy = 0.1*eV; // Minimum e- 93 killBelowEnergy = 0.1*eV; // Minimum e- energy for energy loss by excitation
94 lowEnergyLimit = 0.1 * eV; 94 lowEnergyLimit = 0.1 * eV;
95 lowEnergyLimitOfModel = 10 * eV; // The mode 95 lowEnergyLimitOfModel = 10 * eV; // The model lower energy is 10 eV
96 highEnergyLimit = 500. * keV; 96 highEnergyLimit = 500. * keV;
97 SetLowEnergyLimit(lowEnergyLimit); 97 SetLowEnergyLimit(lowEnergyLimit);
98 SetHighEnergyLimit(highEnergyLimit); 98 SetHighEnergyLimit(highEnergyLimit);
99 99
100 verboseLevel= 0; 100 verboseLevel= 0;
101 // Verbosity scale: 101 // Verbosity scale:
102 // 0 = nothing 102 // 0 = nothing
103 // 1 = warning for energy non-conservation 103 // 1 = warning for energy non-conservation
104 // 2 = details of energy budget 104 // 2 = details of energy budget
105 // 3 = calculation of cross sections, file o 105 // 3 = calculation of cross sections, file openings, sampling of atoms
106 // 4 = entering in methods 106 // 4 = entering in methods
107 107
108 if( verboseLevel>0 ) 108 if( verboseLevel>0 )
109 { 109 {
110 G4cout << "MicroElec Elastic model is co 110 G4cout << "MicroElec Elastic model is constructed " << G4endl
111 << "Energy range: " 111 << "Energy range: "
112 << lowEnergyLimit / eV << " eV - " 112 << lowEnergyLimit / eV << " eV - "
113 << highEnergyLimit / MeV << " MeV" 113 << highEnergyLimit / MeV << " MeV"
114 << G4endl; 114 << G4endl;
115 } 115 }
116 fParticleChangeForGamma = 0; 116 fParticleChangeForGamma = 0;
117 117
118 killElectron = false; 118 killElectron = false;
119 acousticModelEnabled = false; 119 acousticModelEnabled = false;
120 currentMaterialName = ""; 120 currentMaterialName = "";
121 isOkToBeInitialised = false; 121 isOkToBeInitialised = false;
122 } 122 }
123 123
124 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125 125
126 G4MicroElecElasticModel_new::~G4MicroElecElast 126 G4MicroElecElasticModel_new::~G4MicroElecElasticModel_new()
127 { 127 {
128 // For total cross section 128 // For total cross section
129 TCSMap::iterator pos2; 129 TCSMap::iterator pos2;
130 for (pos2 = tableTCS.begin(); pos2 != tableT 130 for (pos2 = tableTCS.begin(); pos2 != tableTCS.end(); ++pos2) {
131 MapData* tableData = pos2->second; 131 MapData* tableData = pos2->second;
132 std::map< G4String, G4MicroElecCrossSectio 132 std::map< G4String, G4MicroElecCrossSectionDataSet_new*, std::less<G4String> >::iterator pos;
133 for (pos = tableData->begin(); pos != tabl 133 for (pos = tableData->begin(); pos != tableData->end(); ++pos)
134 { 134 {
135 G4MicroElecCrossSectionDataSet_new* table = 135 G4MicroElecCrossSectionDataSet_new* table = pos->second;
136 delete table; 136 delete table;
137 } 137 }
138 delete tableData; 138 delete tableData;
139 } 139 }
140 140
141 //Clearing DCS maps 141 //Clearing DCS maps
142 142
143 ThetaMap::iterator iterator_angle; 143 ThetaMap::iterator iterator_angle;
144 for (iterator_angle = thetaDataStorage.begin 144 for (iterator_angle = thetaDataStorage.begin(); iterator_angle != thetaDataStorage.end(); ++iterator_angle) {
145 TriDimensionMap* eDiffCrossSectionData = i 145 TriDimensionMap* eDiffCrossSectionData = iterator_angle->second;
146 eDiffCrossSectionData->clear(); 146 eDiffCrossSectionData->clear();
147 delete eDiffCrossSectionData; 147 delete eDiffCrossSectionData;
148 } 148 }
149 149
150 energyMap::iterator iterator_energy; 150 energyMap::iterator iterator_energy;
151 for (iterator_energy = eIncidentEnergyStorag 151 for (iterator_energy = eIncidentEnergyStorage.begin(); iterator_energy != eIncidentEnergyStorage.end(); ++iterator_energy) {
152 std::vector<G4double>* eTdummyVec = iterat 152 std::vector<G4double>* eTdummyVec = iterator_energy->second;
153 eTdummyVec->clear(); 153 eTdummyVec->clear();
154 delete eTdummyVec; 154 delete eTdummyVec;
155 } 155 }
156 156
157 ProbaMap::iterator iterator_proba; 157 ProbaMap::iterator iterator_proba;
158 for (iterator_proba = eProbaStorage.begin(); 158 for (iterator_proba = eProbaStorage.begin(); iterator_proba != eProbaStorage.end(); ++iterator_proba) {
159 VecMap* eVecm = iterator_proba->second; 159 VecMap* eVecm = iterator_proba->second;
160 eVecm->clear(); 160 eVecm->clear();
161 delete eVecm; 161 delete eVecm;
162 } 162 }
163 163
164 } 164 }
165 165
166 //....oooOO0OOooo........oooOO0OOooo........oo 166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
167 167
168 void G4MicroElecElasticModel_new::Initialise(c 168 void G4MicroElecElasticModel_new::Initialise(const G4ParticleDefinition* /*particle*/,
169 const G4DataVector& /*cuts*/) 169 const G4DataVector& /*cuts*/)
170 { 170 {
171 if (isOkToBeInitialised == true && isInitial 171 if (isOkToBeInitialised == true && isInitialised == false) {
172 172
173 if (verboseLevel > -1) 173 if (verboseLevel > -1)
174 G4cout << "Calling G4MicroElecElasticMod 174 G4cout << "Calling G4MicroElecElasticModel_new::Initialise()" << G4endl;
175 // Energy limits 175 // Energy limits
176 // Reading of data files 176 // Reading of data files
177 177
178 G4double scaleFactor = 1e-18 * cm * cm; 178 G4double scaleFactor = 1e-18 * cm * cm;
179 179
180 G4ProductionCutsTable* theCoupleTable = 180 G4ProductionCutsTable* theCoupleTable =
181 G4ProductionCutsTable::GetProductionCuts 181 G4ProductionCutsTable::GetProductionCutsTable();
182 G4int numOfCouples = (G4int)theCoupleTable 182 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
183 183
184 for (G4int i = 0; i < numOfCouples; ++i) { 184 for (G4int i = 0; i < numOfCouples; ++i) {
185 const G4Material* material = 185 const G4Material* material =
186 theCoupleTable->GetMaterialCutsCouple(i)->Ge 186 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
187 187
188 //theCoupleTable->GetMaterialCutsCouple( 188 //theCoupleTable->GetMaterialCutsCouple(i)->;
189 189
190 G4cout << "MicroElasticModel, Material " 190 G4cout << "MicroElasticModel, Material " << i + 1 << " / " << numOfCouples << " : " << material->GetName() << G4endl;
191 if (material->GetName() == "Vacuum") con 191 if (material->GetName() == "Vacuum") continue;
192 192
193 G4String matName = material->GetName().s 193 G4String matName = material->GetName().substr(3, material->GetName().size());
194 G4cout<< matName<< G4endl; 194 G4cout<< matName<< G4endl;
195 195
196 currentMaterialStructure = new G4MicroEl 196 currentMaterialStructure = new G4MicroElecMaterialStructure(matName);
197 lowEnergyLimitTable[matName]=currentMate 197 lowEnergyLimitTable[matName]=currentMaterialStructure->GetElasticModelLowLimit();
198 highEnergyLimitTable[matName]=currentMat 198 highEnergyLimitTable[matName]=currentMaterialStructure->GetElasticModelHighLimit();
199 workFunctionTable[matName] = currentMate 199 workFunctionTable[matName] = currentMaterialStructure->GetWorkFunction();
200 200
201 delete currentMaterialStructure; 201 delete currentMaterialStructure;
202 202
203 G4cout << "Reading TCS file" << G4endl; 203 G4cout << "Reading TCS file" << G4endl;
204 G4String fileElectron = "Elastic/elsepa_ 204 G4String fileElectron = "Elastic/elsepa_elastic_cross_e_" + matName;
205 G4cout << "Elastic Total Cross file : " 205 G4cout << "Elastic Total Cross file : " << fileElectron << G4endl;
206 206
207 G4ParticleDefinition* electronDef = G4El 207 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
208 G4String electron = electronDef->GetPart 208 G4String electron = electronDef->GetParticleName();
209 209
210 // For total cross section 210 // For total cross section
211 MapData* tableData = new MapData(); 211 MapData* tableData = new MapData();
212 212
213 G4MicroElecCrossSectionDataSet_new* tabl 213 G4MicroElecCrossSectionDataSet_new* tableE = new G4MicroElecCrossSectionDataSet_new(new G4LogLogInterpolation, eV, scaleFactor);
214 tableE->LoadData(fileElectron); 214 tableE->LoadData(fileElectron);
215 tableData->insert(make_pair(electron, ta 215 tableData->insert(make_pair(electron, tableE));
216 tableTCS[matName] = tableData; //Storage 216 tableTCS[matName] = tableData; //Storage of TCS
217 217
218 // For final state 218 // For final state
219 const char* path = G4FindDataDir("G4LEDA 219 const char* path = G4FindDataDir("G4LEDATA");
220 if (!path) 220 if (!path)
221 { 221 {
222 G4Exception("G4MicroElecElasticModel_new:: 222 G4Exception("G4MicroElecElasticModel_new::Initialise","em0006",FatalException,"G4LEDATA environment variable not set.");
223 return; 223 return;
224 } 224 }
225 225
226 //Reading DCS file 226 //Reading DCS file
227 std::ostringstream eFullFileName; 227 std::ostringstream eFullFileName;
228 eFullFileName << path << "/microelec/Ela 228 eFullFileName << path << "/microelec/Elastic/elsepa_elastic_cumulated_diffcross_e_" + matName + ".dat";
229 G4cout << "Elastic Cumulated Diff Cross 229 G4cout << "Elastic Cumulated Diff Cross : " << eFullFileName.str().c_str() << G4endl;
230 std::ifstream eDiffCrossSection(eFullFil 230 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
231 231
232 if (!eDiffCrossSection) 232 if (!eDiffCrossSection)
233 G4Exception("G4MicroElecElasticModel_new::In 233 G4Exception("G4MicroElecElasticModel_new::Initialise", "em0003", FatalException, "Missing data file: /microelec/sigmadiff_cumulated_elastic_e_Si.dat");
234 234
235 // October 21th, 2014 - Melanie Raine 235 // October 21th, 2014 - Melanie Raine
236 // Added clear for MT 236 // Added clear for MT
237 // Diff Cross Sections in cumulated mode 237 // Diff Cross Sections in cumulated mode
238 TriDimensionMap* eDiffCrossSectionData = 238 TriDimensionMap* eDiffCrossSectionData = new TriDimensionMap(); //Angles
239 std::vector<G4double>* eTdummyVec = new 239 std::vector<G4double>* eTdummyVec = new std::vector<G4double>; //Incident energy vector
240 VecMap* eProbVec = new VecMap; //Probabi 240 VecMap* eProbVec = new VecMap; //Probabilities
241 241
242 eTdummyVec->push_back(0.); 242 eTdummyVec->push_back(0.);
243 243
244 while (!eDiffCrossSection.eof()) 244 while (!eDiffCrossSection.eof())
245 { 245 {
246 G4double tDummy; //incident energy 246 G4double tDummy; //incident energy
247 G4double eProb; //Proba 247 G4double eProb; //Proba
248 eDiffCrossSection >> tDummy >> eProb; 248 eDiffCrossSection >> tDummy >> eProb;
249 249
250 // SI : mandatory eVecm initialization 250 // SI : mandatory eVecm initialization
251 if (tDummy != eTdummyVec->back()) 251 if (tDummy != eTdummyVec->back())
252 { 252 {
253 eTdummyVec->push_back(tDummy); //addin 253 eTdummyVec->push_back(tDummy); //adding values for incident energy points
254 (*eProbVec)[tDummy].push_back(0.); //a 254 (*eProbVec)[tDummy].push_back(0.); //adding probability for the first angle, equal to 0
255 } 255 }
256 256
257 eDiffCrossSection >> (*eDiffCrossSectionDa 257 eDiffCrossSection >> (*eDiffCrossSectionData)[tDummy][eProb]; //adding Angle Value to map
258 258
259 if (eProb != (*eProbVec)[tDummy].back()) { 259 if (eProb != (*eProbVec)[tDummy].back()) {
260 (*eProbVec)[tDummy].push_back(eProb); // 260 (*eProbVec)[tDummy].push_back(eProb); //Adding cumulated proba to map
261 } 261 }
262 262
263 } 263 }
264 264
265 //Filling maps for the material 265 //Filling maps for the material
266 thetaDataStorage[matName] = eDiffCrossSe 266 thetaDataStorage[matName] = eDiffCrossSectionData;
267 eIncidentEnergyStorage[matName] = eTdumm 267 eIncidentEnergyStorage[matName] = eTdummyVec;
268 eProbaStorage[matName] = eProbVec; 268 eProbaStorage[matName] = eProbVec;
269 } 269 }
270 // End final state 270 // End final state
271 271
272 if (verboseLevel > 2) 272 if (verboseLevel > 2)
273 G4cout << "Loaded cross section files fo 273 G4cout << "Loaded cross section files for MicroElec Elastic model" << G4endl;
274 274
275 if (verboseLevel > 0) 275 if (verboseLevel > 0)
276 { 276 {
277 G4cout << "MicroElec Elastic model is initia 277 G4cout << "MicroElec Elastic model is initialized " << G4endl
278 << "Energy range: " 278 << "Energy range: "
279 << LowEnergyLimit() / eV << " eV - " 279 << LowEnergyLimit() / eV << " eV - "
280 << HighEnergyLimit() / MeV << " MeV" 280 << HighEnergyLimit() / MeV << " MeV"
281 << G4endl; // system("pause"); linux 281 << G4endl; // system("pause"); linux doesn't like
282 } 282 }
283 283
284 if (isInitialised) { return; } 284 if (isInitialised) { return; }
285 fParticleChangeForGamma = GetParticleChang 285 fParticleChangeForGamma = GetParticleChangeForGamma();
286 isInitialised = true; 286 isInitialised = true;
287 } 287 }
288 } 288 }
289 289
290 //....oooOO0OOooo........oooOO0OOooo........oo 290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291 291
292 G4double G4MicroElecElasticModel_new::CrossSec 292 G4double G4MicroElecElasticModel_new::CrossSectionPerVolume(const G4Material* material,
293 const G4ParticleDefinition* p, 293 const G4ParticleDefinition* p,
294 G4double ekin, 294 G4double ekin,
295 G4double, 295 G4double,
296 G4double) 296 G4double)
297 { 297 {
298 if (verboseLevel > 3) 298 if (verboseLevel > 3)
299 G4cout << "Calling CrossSectionPerVolume() 299 G4cout << "Calling CrossSectionPerVolume() of G4MicroElecElasticModel" << G4endl;
300 300
301 isOkToBeInitialised = true; 301 isOkToBeInitialised = true;
302 currentMaterialName = material->GetName().su 302 currentMaterialName = material->GetName().substr(3, material->GetName().size());
303 const G4DataVector cuts; 303 const G4DataVector cuts;
304 Initialise(p, cuts); 304 Initialise(p, cuts);
305 // Calculate total cross section for model 305 // Calculate total cross section for model
306 MapEnergy::iterator lowEPos; 306 MapEnergy::iterator lowEPos;
307 lowEPos = lowEnergyLimitTable.find(currentMa 307 lowEPos = lowEnergyLimitTable.find(currentMaterialName);
308 308
309 MapEnergy::iterator highEPos; 309 MapEnergy::iterator highEPos;
310 highEPos = highEnergyLimitTable.find(current 310 highEPos = highEnergyLimitTable.find(currentMaterialName);
311 311
312 MapEnergy::iterator killEPos; 312 MapEnergy::iterator killEPos;
313 killEPos = workFunctionTable.find(currentMat 313 killEPos = workFunctionTable.find(currentMaterialName);
314 314
315 if (lowEPos == lowEnergyLimitTable.end() || 315 if (lowEPos == lowEnergyLimitTable.end() || highEPos == highEnergyLimitTable.end() || killEPos == workFunctionTable.end())
316 { 316 {
317 G4String str = "Material "; 317 G4String str = "Material ";
318 str += currentMaterialName + " not found 318 str += currentMaterialName + " not found!";
319 G4Exception("G4MicroElecElasticModel_new 319 G4Exception("G4MicroElecElasticModel_new::EnergyLimits", "em0002", FatalException, str);
320 return 0; 320 return 0;
321 } 321 }
322 else { 322 else {
323 // G4cout << "normal elastic " << G4endl 323 // G4cout << "normal elastic " << G4endl;
324 lowEnergyLimit = lowEPos->second; 324 lowEnergyLimit = lowEPos->second;
325 highEnergyLimit = highEPos->second; 325 highEnergyLimit = highEPos->second;
326 killBelowEnergy = killEPos->second; 326 killBelowEnergy = killEPos->second;
327 327
328 } 328 }
329 329
330 if (ekin < killBelowEnergy) { 330 if (ekin < killBelowEnergy) {
331 331
332 return DBL_MAX; } 332 return DBL_MAX; }
333 333
334 G4double sigma=0; 334 G4double sigma=0;
335 335
336 //Phonon for SiO2 336 //Phonon for SiO2
337 if (currentMaterialName == "SILICON_DIOXIDE" 337 if (currentMaterialName == "SILICON_DIOXIDE" && ekin < 100 * eV) {
338 acousticModelEnabled = true; 338 acousticModelEnabled = true;
339 339
340 //Values for SiO2 340 //Values for SiO2
341 G4double kbz = 11.54e9, 341 G4double kbz = 11.54e9,
342 rho = 2.2 * 1000, // [g/cm3] * 1000 342 rho = 2.2 * 1000, // [g/cm3] * 1000
343 cs = 3560, //Sound speed 343 cs = 3560, //Sound speed
344 Ebz = 5.1 * 1.6e-19, 344 Ebz = 5.1 * 1.6e-19,
345 Aac = 17 * Ebz, //A screening parameter 345 Aac = 17 * Ebz, //A screening parameter
346 Eac = 3.5 * 1.6e-19, //C deformation pot 346 Eac = 3.5 * 1.6e-19, //C deformation potential
347 prefactor = 2.2;// Facteur pour modifier 347 prefactor = 2.2;// Facteur pour modifier les MFP
348 348
349 return AcousticCrossSectionPerVolume(ekin, 349 return AcousticCrossSectionPerVolume(ekin, kbz, rho, cs, Aac, Eac, prefactor);
350 } 350 }
351 351
352 else if (currentMaterialName == "ALUMINUM_OXID 352 else if (currentMaterialName == "ALUMINUM_OXIDE" && ekin < 20 * eV) {
353 acousticModelEnabled = true; 353 acousticModelEnabled = true;
354 354
355 //Values for Al2O3 355 //Values for Al2O3
356 G4double kbz = 8871930614.247564, 356 G4double kbz = 8871930614.247564,
357 rho = 3.97 * 1000, // [g/cm3] * 1000 357 rho = 3.97 * 1000, // [g/cm3] * 1000
358 cs = 233329.07733059773, //Sound speed 358 cs = 233329.07733059773, //Sound speed
359 Aac = 2.9912494342262614e-19, //A screeni 359 Aac = 2.9912494342262614e-19, //A screening parameter
360 Eac = 2.1622471654789847e-18, //C deforma 360 Eac = 2.1622471654789847e-18, //C deformation potential
361 prefactor = 1; 361 prefactor = 1;
362 return AcousticCrossSectionPerVolume(ekin, 362 return AcousticCrossSectionPerVolume(ekin, kbz, rho, cs, Aac, Eac, prefactor);
363 } 363 }
364 //Elastic 364 //Elastic
365 else { 365 else {
366 acousticModelEnabled = false; 366 acousticModelEnabled = false;
367 367
368 G4double density = material->GetTotNbOfAto 368 G4double density = material->GetTotNbOfAtomsPerVolume();
369 const G4String& particleName = p->GetParti 369 const G4String& particleName = p->GetParticleName();
370 370
371 TCSMap::iterator tablepos; 371 TCSMap::iterator tablepos;
372 tablepos = tableTCS.find(currentMaterialNa 372 tablepos = tableTCS.find(currentMaterialName);
373 373
374 if (tablepos != tableTCS.end()) 374 if (tablepos != tableTCS.end())
375 { 375 {
376 MapData* tableData = tablepos->second; 376 MapData* tableData = tablepos->second;
377 377
378 if (ekin >= lowEnergyLimit && ekin < highEne 378 if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
379 { 379 {
380 std::map< G4String, G4MicroElecCrossSect 380 std::map< G4String, G4MicroElecCrossSectionDataSet_new*, std::less<G4String> >::iterator pos;
381 pos = tableData->find(particleName); 381 pos = tableData->find(particleName);
382 382
383 if (pos != tableData->end()) 383 if (pos != tableData->end())
384 { 384 {
385 G4MicroElecCrossSectionDataSet_new* table 385 G4MicroElecCrossSectionDataSet_new* table = pos->second;
386 if (table != 0) 386 if (table != 0)
387 { 387 {
388 sigma = table->FindValue(ekin); 388 sigma = table->FindValue(ekin);
389 } 389 }
390 } 390 }
391 else 391 else
392 { 392 {
393 G4Exception("G4MicroElecElasticModel_new:: 393 G4Exception("G4MicroElecElasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, "Model not applicable to particle type.");
394 } 394 }
395 } 395 }
396 else return 1 / DBL_MAX; 396 else return 1 / DBL_MAX;
397 } 397 }
398 else 398 else
399 { 399 {
400 G4String str = "Material "; 400 G4String str = "Material ";
401 str += currentMaterialName + " TCS table not 401 str += currentMaterialName + " TCS table not found!";
402 G4Exception("G4MicroElecElasticModel_new::Co 402 G4Exception("G4MicroElecElasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, str);
403 } 403 }
404 404
405 if (verboseLevel > 3) 405 if (verboseLevel > 3)
406 { 406 {
407 G4cout << "---> Kinetic energy(eV)=" << ekin 407 G4cout << "---> Kinetic energy(eV)=" << ekin / eV << G4endl;
408 G4cout << " - Cross section per Si atom (cm^ 408 G4cout << " - Cross section per Si atom (cm^2)=" << sigma / cm / cm << G4endl;
409 G4cout << " - Cross section per Si atom (cm^ 409 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density / (1. / cm) << G4endl;
410 } 410 }
411 411
412 // Hsing-YinChangaAndrewAlvaradoaTreyWeber 412 // Hsing-YinChangaAndrewAlvaradoaTreyWeberaJaimeMarianab Monte Carlo modeling of low-energy electron-induced secondary electron emission yields in micro-architected boron nitride surfaces - ScienceDirect, (n.d.). https://www.sciencedirect.com/science/article/pii/S0168583X19304069 (accessed April 1, 2022).
413 if (currentMaterialName == "BORON_NITRIDE" 413 if (currentMaterialName == "BORON_NITRIDE") {
414 sigma = sigma * tanh(0.5 * pow(ekin / 414 sigma = sigma * tanh(0.5 * pow(ekin / 5.2e-6, 2));
415 } 415 }
416 return sigma*density; 416 return sigma*density;
417 } 417 }
418 } 418 }
419 419
420 //....oooOO0OOooo........oooOO0OOooo........oo 420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
421 421
422 G4double G4MicroElecElasticModel_new::Acoustic 422 G4double G4MicroElecElasticModel_new::AcousticCrossSectionPerVolume(G4double ekin,
423 G4double kbz, 423 G4double kbz,
424 G4double rho, 424 G4double rho,
425 G4double cs, 425 G4double cs,
426 G4double Aac, 426 G4double Aac,
427 G4double Eac, 427 G4double Eac,
428 G4double prefactor) 428 G4double prefactor)
429 { 429 {
430 430
431 G4double e = 1.6e-19, 431 G4double e = 1.6e-19,
432 m0 = 9.10938356e-31, 432 m0 = 9.10938356e-31,
433 h = 1.0546e-34, 433 h = 1.0546e-34,
434 kb = 1.38e-23; 434 kb = 1.38e-23;
435 435
436 G4double E = (ekin / eV) * e; 436 G4double E = (ekin / eV) * e;
437 G4double D = (2 / (std::sqrt(2) * std::pow(p 437 G4double D = (2 / (std::sqrt(2) * std::pow(pi, 2) * std::pow(h, 3))) * (1 + 2 * E) * std::pow(m0, 1.5) * std::sqrt(E);
438 438
439 // Parametres SiO2 439 // Parametres SiO2
440 G4double T = 300, 440 G4double T = 300,
441 Ebz = (std::pow(h, 2) * std::pow(kbz, 2)) 441 Ebz = (std::pow(h, 2) * std::pow(kbz, 2)) / (2 * m0),
442 hwbz = cs * kbz * h, 442 hwbz = cs * kbz * h,
443 nbz = 1.0 / (exp(hwbz / (kb * T)) - 1), 443 nbz = 1.0 / (exp(hwbz / (kb * T)) - 1),
444 Pac; 444 Pac;
445 445
446 if (E < Ebz / 4.0) 446 if (E < Ebz / 4.0)
447 { 447 {
448 Pac = ((pi * kb * T) / (h * std::pow(cs, 448 Pac = ((pi * kb * T) / (h * std::pow(cs, 2) * rho)) * (std::pow(Eac, 2) * D) / (1 + (E / Aac));
449 } 449 }
450 450
451 else if (E > Ebz) //Screened relationship 451 else if (E > Ebz) //Screened relationship
452 { 452 {
453 Pac = ((2 * pi * m0 * (2 * nbz + 1)) / ( 453 Pac = ((2 * pi * m0 * (2 * nbz + 1)) / (h * rho * hwbz)) * std::pow(Eac, 2) * D * E * 2 * std::pow((Aac / E), 2) * (((-E / Aac) / (1 + (E / Aac))) + log(1 + (E / Aac)));
454 } 454 }
455 else //Linear interpolation 455 else //Linear interpolation
456 { 456 {
457 G4double fEbz = ((2 * pi * m0 * (2 * nbz 457 G4double fEbz = ((2 * pi * m0 * (2 * nbz + 1)) / (h * rho * hwbz)) * std::pow(Eac, 2) * D * Ebz * 2 * std::pow((Aac / Ebz), 2) * (((-Ebz / Aac) / (1 + (Ebz / Aac))) + log(1 + (Ebz / Aac)));
458 G4double fEbz4 = ((pi * kb * T) / (h * s 458 G4double fEbz4 = ((pi * kb * T) / (h * std::pow(cs, 2) * rho)) * (std::pow(Eac, 2) * D) / (1 + ((Ebz / 4) / Aac));
459 G4double alpha = ((fEbz - fEbz4) / (Ebz 459 G4double alpha = ((fEbz - fEbz4) / (Ebz - (Ebz / 4)));
460 Pac = alpha * E + (fEbz - alpha * Ebz); 460 Pac = alpha * E + (fEbz - alpha * Ebz);
461 } 461 }
462 462
463 G4double MFP = (std::sqrt(2 * E / m0) / (pre 463 G4double MFP = (std::sqrt(2 * E / m0) / (prefactor * Pac)) * m;
464 464
465 return 1 / MFP; 465 return 1 / MFP;
466 } 466 }
467 467
468 468
469 //....oooOO0OOooo........oooOO0OOooo........oo 469 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
470 470
471 void G4MicroElecElasticModel_new::SampleSecond 471 void G4MicroElecElasticModel_new::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
472 const G4MaterialCutsCouple* /*coup 472 const G4MaterialCutsCouple* /*couple*/,
473 const G4DynamicParticle* aDynamicE 473 const G4DynamicParticle* aDynamicElectron,
474 G4double, 474 G4double,
475 G4double) 475 G4double)
476 { 476 {
477 477
478 if (verboseLevel > 3) 478 if (verboseLevel > 3)
479 G4cout << "Calling SampleSecondaries() of 479 G4cout << "Calling SampleSecondaries() of G4MicroElecElasticModel" << G4endl;
480 480
481 G4double electronEnergy0 = aDynamicElectron- 481 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
482 482
483 if (electronEnergy0 < killBelowEnergy) 483 if (electronEnergy0 < killBelowEnergy)
484 { 484 {
485 fParticleChangeForGamma->SetProposedKinet 485 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
486 fParticleChangeForGamma->ProposeTrackStat 486 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
487 fParticleChangeForGamma->ProposeLocalEner 487 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
488 return; 488 return;
489 } 489 }
490 490
491 if (electronEnergy0 < highEnergyLimit) 491 if (electronEnergy0 < highEnergyLimit)
492 { 492 {
493 G4double cosTheta = 0; 493 G4double cosTheta = 0;
494 if (acousticModelEnabled) 494 if (acousticModelEnabled)
495 { 495 {
496 cosTheta = 1 - 2 * G4UniformRand(); //Isotr 496 cosTheta = 1 - 2 * G4UniformRand(); //Isotrope
497 } 497 }
498 else if (electronEnergy0 >= lowEnergyLimi 498 else if (electronEnergy0 >= lowEnergyLimit)
499 { 499 {
500 cosTheta = RandomizeCosTheta(electronEnergy 500 cosTheta = RandomizeCosTheta(electronEnergy0);
501 } 501 }
502 502
503 G4double phi = 2. * pi * G4UniformRand(); 503 G4double phi = 2. * pi * G4UniformRand();
504 504
505 G4ThreeVector zVers = aDynamicElectron->G 505 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
506 G4ThreeVector xVers = zVers.orthogonal(); 506 G4ThreeVector xVers = zVers.orthogonal();
507 G4ThreeVector yVers = zVers.cross(xVers); 507 G4ThreeVector yVers = zVers.cross(xVers);
508 508
509 G4double xDir = std::sqrt(1. - cosTheta*c 509 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
510 G4double yDir = xDir; 510 G4double yDir = xDir;
511 xDir *= std::cos(phi); 511 xDir *= std::cos(phi);
512 yDir *= std::sin(phi); 512 yDir *= std::sin(phi);
513 513
514 G4ThreeVector zPrimeVers((xDir*xVers + yD 514 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
515 515
516 fParticleChangeForGamma->ProposeMomentumD 516 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
517 fParticleChangeForGamma->SetProposedKinet 517 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
518 } 518 }
519 } 519 }
520 520
521 //....oooOO0OOooo........oooOO0OOooo........oo 521 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
522 522
523 G4double G4MicroElecElasticModel_new::DamageE 523 G4double G4MicroElecElasticModel_new::DamageEnergy(G4double T,G4double A, G4double Z)
524 { 524 {
525 //.................. T in eV!!!!!!!!!!!!! 525 //.................. T in eV!!!!!!!!!!!!!
526 G4double Z2= Z; 526 G4double Z2= Z;
527 G4double M2= A; 527 G4double M2= A;
528 G4double k_d; 528 G4double k_d;
529 G4double epsilon_d; 529 G4double epsilon_d;
530 G4double g_epsilon_d; 530 G4double g_epsilon_d;
531 G4double E_nu; 531 G4double E_nu;
532 532
533 k_d=0.1334*std::pow(Z2,(2./3.))*std::pow(M2, 533 k_d=0.1334*std::pow(Z2,(2./3.))*std::pow(M2,(-1./2.));
534 epsilon_d=0.01014*std::pow(Z2,(-7./3.))*(T/e 534 epsilon_d=0.01014*std::pow(Z2,(-7./3.))*(T/eV);
535 g_epsilon_d= epsilon_d+0.40244*std::pow(epsi 535 g_epsilon_d= epsilon_d+0.40244*std::pow(epsilon_d,(3./4.))+3.4008*std::pow(epsilon_d,(1./6.));
536 536
537 E_nu=1./(1.+ k_d*g_epsilon_d); 537 E_nu=1./(1.+ k_d*g_epsilon_d);
538 538
539 return E_nu; 539 return E_nu;
540 } 540 }
541 541
542 //....oooOO0OOooo........oooOO0OOooo........oo 542 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
543 543
544 G4double G4MicroElecElasticModel_new::Theta 544 G4double G4MicroElecElasticModel_new::Theta
545 (G4ParticleDefinition * particleDefinition, 545 (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)
546 { 546 {
547 547
548 G4double theta = 0.; 548 G4double theta = 0.;
549 G4double valueT1 = 0; 549 G4double valueT1 = 0;
550 G4double valueT2 = 0; 550 G4double valueT2 = 0;
551 G4double valueE21 = 0; 551 G4double valueE21 = 0;
552 G4double valueE22 = 0; 552 G4double valueE22 = 0;
553 G4double valueE12 = 0; 553 G4double valueE12 = 0;
554 G4double valueE11 = 0; 554 G4double valueE11 = 0;
555 G4double xs11 = 0; 555 G4double xs11 = 0;
556 G4double xs12 = 0; 556 G4double xs12 = 0;
557 G4double xs21 = 0; 557 G4double xs21 = 0;
558 G4double xs22 = 0; 558 G4double xs22 = 0;
559 559
560 if (particleDefinition == G4Electron::Electr 560 if (particleDefinition == G4Electron::ElectronDefinition())
561 { 561 {
562 ThetaMap::iterator iterator_angle; 562 ThetaMap::iterator iterator_angle;
563 iterator_angle = thetaDataStorage.find(cur 563 iterator_angle = thetaDataStorage.find(currentMaterialName);
564 564
565 energyMap::iterator iterator_energy; 565 energyMap::iterator iterator_energy;
566 iterator_energy = eIncidentEnergyStorage.f 566 iterator_energy = eIncidentEnergyStorage.find(currentMaterialName);
567 567
568 ProbaMap::iterator iterator_proba; 568 ProbaMap::iterator iterator_proba;
569 iterator_proba = eProbaStorage.find(curren 569 iterator_proba = eProbaStorage.find(currentMaterialName);
570 570
571 if (iterator_angle != thetaDataStorage.end 571 if (iterator_angle != thetaDataStorage.end() && iterator_energy != eIncidentEnergyStorage.end() && iterator_proba != eProbaStorage.end())
572 { 572 {
573 TriDimensionMap* eDiffCrossSectionData = ite 573 TriDimensionMap* eDiffCrossSectionData = iterator_angle->second; //Theta points
574 std::vector<G4double>* eTdummyVec = iterator 574 std::vector<G4double>* eTdummyVec = iterator_energy->second;
575 VecMap* eVecm = iterator_proba->second; 575 VecMap* eVecm = iterator_proba->second;
576 576
577 auto t2 = std::upper_bound(eTdummyVec->begin 577 auto t2 = std::upper_bound(eTdummyVec->begin(), eTdummyVec->end(), k);
578 auto t1 = t2 - 1; 578 auto t1 = t2 - 1;
579 auto e12 = std::upper_bound((*eVecm)[( 579 auto e12 = std::upper_bound((*eVecm)[(*t1)].begin(), (*eVecm)[(*t1)].end(), integrDiff);
580 auto e11 = e12 - 1; 580 auto e11 = e12 - 1;
581 auto e22 = std::upper_bound((*eVecm)[(*t2)]. 581 auto e22 = std::upper_bound((*eVecm)[(*t2)].begin(), (*eVecm)[(*t2)].end(), integrDiff);
582 auto e21 = e22 - 1; 582 auto e21 = e22 - 1;
583 583
584 valueT1 = *t1; 584 valueT1 = *t1;
585 valueT2 = *t2; 585 valueT2 = *t2;
586 valueE21 = *e21; 586 valueE21 = *e21;
587 valueE22 = *e22; 587 valueE22 = *e22;
588 valueE12 = *e12; 588 valueE12 = *e12;
589 valueE11 = *e11; 589 valueE11 = *e11;
590 590
591 xs11 = (*eDiffCrossSectionData)[valueT1][val 591 xs11 = (*eDiffCrossSectionData)[valueT1][valueE11];
592 xs12 = (*eDiffCrossSectionData)[valueT1][val 592 xs12 = (*eDiffCrossSectionData)[valueT1][valueE12];
593 xs21 = (*eDiffCrossSectionData)[valueT2][val 593 xs21 = (*eDiffCrossSectionData)[valueT2][valueE21];
594 xs22 = (*eDiffCrossSectionData)[valueT2][val 594 xs22 = (*eDiffCrossSectionData)[valueT2][valueE22];
595 } 595 }
596 else 596 else
597 { 597 {
598 G4String str = "Material "; 598 G4String str = "Material ";
599 str += currentMaterialName + " not found!"; 599 str += currentMaterialName + " not found!";
600 G4Exception("G4MicroElecElasticModel_new::Co 600 G4Exception("G4MicroElecElasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, str);
601 } 601 }
602 602
603 } 603 }
604 604
605 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) 605 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
606 606
607 theta = QuadInterpolator( valueE11, valueE1 607 theta = QuadInterpolator( valueE11, valueE12,
608 valueE21, valueE22, 608 valueE21, valueE22,
609 xs11, xs12, 609 xs11, xs12,
610 xs21, xs22, 610 xs21, xs22,
611 valueT1, valueT2, 611 valueT1, valueT2,
612 k, integrDiff ); 612 k, integrDiff );
613 613
614 return theta; 614 return theta;
615 } 615 }
616 616
617 //....oooOO0OOooo........oooOO0OOooo........oo 617 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
618 618
619 G4double G4MicroElecElasticModel_new::LinLogIn 619 G4double G4MicroElecElasticModel_new::LinLogInterpolate(G4double e1,
620 G4double e2, 620 G4double e2,
621 G4double e, 621 G4double e,
622 G4double xs1, 622 G4double xs1,
623 G4double xs2) 623 G4double xs2)
624 { 624 {
625 G4double d1 = std::log(xs1); 625 G4double d1 = std::log(xs1);
626 G4double d2 = std::log(xs2); 626 G4double d2 = std::log(xs2);
627 G4double value = G4Exp(d1 + (d2 - d1)*(e - e 627 G4double value = G4Exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
628 return value; 628 return value;
629 } 629 }
630 630
631 //....oooOO0OOooo........oooOO0OOooo........oo 631 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
632 632
633 G4double G4MicroElecElasticModel_new::LinLinIn 633 G4double G4MicroElecElasticModel_new::LinLinInterpolate(G4double e1,
634 G4double e2, 634 G4double e2,
635 G4double e, 635 G4double e,
636 G4double xs1, 636 G4double xs1,
637 G4double xs2) 637 G4double xs2)
638 { 638 {
639 G4double d1 = xs1; 639 G4double d1 = xs1;
640 G4double d2 = xs2; 640 G4double d2 = xs2;
641 G4double value = (d1 + (d2 - d1)*(e - e1)/ ( 641 G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
642 return value; 642 return value;
643 } 643 }
644 644
645 //....oooOO0OOooo........oooOO0OOooo........oo 645 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646 646
647 G4double G4MicroElecElasticModel_new::LogLogIn 647 G4double G4MicroElecElasticModel_new::LogLogInterpolate(G4double e1,
648 G4double e2, 648 G4double e2,
649 G4double e, 649 G4double e,
650 G4double xs1, 650 G4double xs1,
651 G4double xs2) 651 G4double xs2)
652 { 652 {
653 G4double a = (std::log10(xs2)-std::log10(xs1 653 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
654 G4double b = std::log10(xs2) - a*std::log10( 654 G4double b = std::log10(xs2) - a*std::log10(e2);
655 G4double sigma = a*std::log10(e) + b; 655 G4double sigma = a*std::log10(e) + b;
656 G4double value = (std::pow(10.,sigma)); 656 G4double value = (std::pow(10.,sigma));
657 return value; 657 return value;
658 } 658 }
659 659
660 //....oooOO0OOooo........oooOO0OOooo........oo 660 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
661 661
662 G4double G4MicroElecElasticModel_new::QuadInte 662 G4double G4MicroElecElasticModel_new::QuadInterpolator(G4double e11, G4double e12,
663 G4double e21, G4double e22, 663 G4double e21, G4double e22,
664 G4double xs11, G4double xs12, 664 G4double xs11, G4double xs12,
665 G4double xs21, G4double xs22, 665 G4double xs21, G4double xs22,
666 G4double t1, G4double t2, 666 G4double t1, G4double t2,
667 G4double t, G4double e) 667 G4double t, G4double e)
668 { 668 {
669 669
670 670
671 // Lin-Lin 671 // Lin-Lin
672 G4double interpolatedvalue1 = LinLinInterpol 672 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
673 G4double interpolatedvalue2 = LinLinInterpol 673 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
674 G4double value = LinLinInterpolate(t1, t2, t 674 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
675 675
676 return value; 676 return value;
677 } 677 }
678 678
679 //....oooOO0OOooo........oooOO0OOooo........oo 679 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
680 680
681 G4double G4MicroElecElasticModel_new::Randomiz 681 G4double G4MicroElecElasticModel_new::RandomizeCosTheta(G4double k)
682 { 682 {
683 G4double integrdiff=0; 683 G4double integrdiff=0;
684 G4double uniformRand=G4UniformRand(); 684 G4double uniformRand=G4UniformRand();
685 integrdiff = uniformRand; 685 integrdiff = uniformRand;
686 686
687 G4double theta=0.; 687 G4double theta=0.;
688 G4double cosTheta=0.; 688 G4double cosTheta=0.;
689 theta = Theta(G4Electron::ElectronDefinition 689 theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
690 690
691 cosTheta= std::cos(theta*pi/180.); 691 cosTheta= std::cos(theta*pi/180.);
692 692
693 return cosTheta; 693 return cosTheta;
694 } 694 }
695 695
696 //....oooOO0OOooo........oooOO0OOooo........oo 696 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
697 697
698 698
699 void G4MicroElecElasticModel_new::SetKillBelow 699 void G4MicroElecElasticModel_new::SetKillBelowThreshold (G4double threshold)
700 { 700 {
701 killBelowEnergy = threshold; 701 killBelowEnergy = threshold;
702 702
703 if (threshold < 5*CLHEP::eV) 703 if (threshold < 5*CLHEP::eV)
704 { 704 {
705 G4Exception ("*** WARNING : the G4MicroE 705 G4Exception ("*** WARNING : the G4MicroElecElasticModel class is not validated below 5 eV !","",JustWarning,"") ;
706 threshold = 5*CLHEP::eV; 706 threshold = 5*CLHEP::eV;
707 } 707 }
708 } 708 }
709 709
710 //....oooOO0OOooo........oooOO0OOooo........oo 710 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
711 711