| 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_new.cc, 2011/08/2 27 // G4MicroElecInelasticModel_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 75
76 76
77 #include "globals.hh" 77 #include "globals.hh"
78 #include "G4MicroElecInelasticModel_new.hh" << 78 #include "G4MicroElecInelasticModel_new.hh"
79 #include "G4PhysicalConstants.hh" 79 #include "G4PhysicalConstants.hh"
80 #include "G4SystemOfUnits.hh" 80 #include "G4SystemOfUnits.hh"
81 #include "G4ios.hh" 81 #include "G4ios.hh"
82 #include "G4UnitsTable.hh" 82 #include "G4UnitsTable.hh"
83 #include "G4UAtomicDeexcitation.hh" 83 #include "G4UAtomicDeexcitation.hh"
84 #include "G4LossTableManager.hh" 84 #include "G4LossTableManager.hh"
85 #include "G4ionEffectiveCharge.hh" 85 #include "G4ionEffectiveCharge.hh"
86 #include "G4MicroElecMaterialStructure.hh" 86 #include "G4MicroElecMaterialStructure.hh"
87 #include "G4DeltaAngle.hh" 87 #include "G4DeltaAngle.hh"
88 88
89 #include <sstream> 89 #include <sstream>
90 90
91 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 92
93 using namespace std; 93 using namespace std;
94 94
95 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 96
97 G4MicroElecInelasticModel_new::G4MicroElecInel 97 G4MicroElecInelasticModel_new::G4MicroElecInelasticModel_new(
98 const G4ParticleDefinition*, const G4St 98 const G4ParticleDefinition*, const G4String& nam)
99 :G4VEmModel(nam),isInitialised(false) 99 :G4VEmModel(nam),isInitialised(false)
100 { 100 {
101 101
102 verboseLevel= 0; 102 verboseLevel= 0;
103 // Verbosity scale: 103 // Verbosity scale:
104 // 0 = nothing 104 // 0 = nothing
105 // 1 = warning for energy non-conservation 105 // 1 = warning for energy non-conservation
106 // 2 = details of energy budget 106 // 2 = details of energy budget
107 // 3 = calculation of cross sections, file o 107 // 3 = calculation of cross sections, file openings, sampling of atoms
108 // 4 = entering in methods 108 // 4 = entering in methods
109 109
110 if( verboseLevel>0 ) 110 if( verboseLevel>0 )
111 { 111 {
112 G4cout << "MicroElec inelastic model is co 112 G4cout << "MicroElec inelastic model is constructed " << G4endl;
113 } 113 }
114 114
115 //Mark this model as "applicable" for atomic 115 //Mark this model as "applicable" for atomic deexcitation
116 SetDeexcitationFlag(true); 116 SetDeexcitationFlag(true);
117 fAtomDeexcitation = nullptr; 117 fAtomDeexcitation = nullptr;
118 fParticleChangeForGamma = nullptr; 118 fParticleChangeForGamma = nullptr;
119 119
120 // default generator 120 // default generator
121 SetAngularDistribution(new G4DeltaAngle()); 121 SetAngularDistribution(new G4DeltaAngle());
122 122
123 // Selection of computation method 123 // Selection of computation method
124 fasterCode = true; 124 fasterCode = true;
125 SEFromFermiLevel = false; 125 SEFromFermiLevel = false;
126 } 126 }
127 127
128 //....oooOO0OOooo........oooOO0OOooo........oo 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 129
130 G4MicroElecInelasticModel_new::~G4MicroElecIne 130 G4MicroElecInelasticModel_new::~G4MicroElecInelasticModel_new()
131 { 131 {
132 // Cross section 132 // Cross section
133 // (0) 133 // (0)
134 for (auto const& p : tableTCS) { << 134 TCSMap::iterator pos2;
135 MapData* tableData = p.second; << 135 for (pos2 = tableTCS.begin(); pos2 != tableTCS.end(); ++pos2) {
136 for (auto const& pos : *tableData) { delet << 136 MapData* tableData = pos2->second;
>> 137 for (auto pos = tableData->begin(); pos != tableData->end(); ++pos)
>> 138 {
>> 139 G4MicroElecCrossSectionDataSet_new* table = pos->second;
>> 140 delete table;
>> 141 }
137 delete tableData; 142 delete tableData;
138 } 143 }
139 tableTCS.clear(); 144 tableTCS.clear();
140 << 145
>> 146 dataDiffCSMap::iterator iterator_proba;
141 // (1) 147 // (1)
142 for (auto const & obj : eDiffDatatable) { << 148 for (iterator_proba = eNrjTransStorage.begin(); iterator_proba != eNrjTransStorage.end(); ++iterator_proba) {
143 auto ptr = obj.second; << 149 vector<TriDimensionMap>* eNrjTransfData = iterator_proba->second;
144 ptr->clear(); << 150 eNrjTransfData->clear();
145 delete ptr; << 151 delete eNrjTransfData;
146 } 152 }
>> 153 eNrjTransStorage.clear();
147 154
148 for (auto const & obj : pDiffDatatable) { << 155 for (iterator_proba = pNrjTransStorage.begin(); iterator_proba != pNrjTransStorage.end(); ++iterator_proba) {
149 auto ptr = obj.second; << 156 vector<TriDimensionMap>* pNrjTransfData = iterator_proba->second;
150 ptr->clear(); << 157 pNrjTransfData->clear();
151 delete ptr; << 158 delete pNrjTransfData;
152 } 159 }
153 << 160 pNrjTransStorage.clear();
>> 161
154 // (2) 162 // (2)
155 for (auto const & obj : eNrjTransStorage) { << 163 for (iterator_proba = eDiffDatatable.begin(); iterator_proba != eDiffDatatable.end(); ++iterator_proba) {
156 delete obj.second; << 164 vector<TriDimensionMap>* eDiffCrossSectionData = iterator_proba->second;
157 } << 165 eDiffCrossSectionData->clear();
158 for (auto const & obj : pNrjTransStorage) { << 166 delete eDiffCrossSectionData;
159 delete obj.second; <<
160 } 167 }
>> 168 eDiffDatatable.clear();
161 169
>> 170 for (iterator_proba = pDiffDatatable.begin(); iterator_proba != pDiffDatatable.end(); ++iterator_proba) {
>> 171 vector<TriDimensionMap>* pDiffCrossSectionData = iterator_proba->second;
>> 172 pDiffCrossSectionData->clear();
>> 173 delete pDiffCrossSectionData;
>> 174 }
>> 175 pDiffDatatable.clear();
>> 176
162 // (3) 177 // (3)
163 for (auto const& p : eProbaShellStorage) { << 178 dataProbaShellMap::iterator iterator_probaShell;
164 delete p.second; << 179
>> 180 for (iterator_probaShell = eProbaShellStorage.begin(); iterator_probaShell != eProbaShellStorage.end(); ++iterator_probaShell) {
>> 181 vector<VecMap>* eProbaShellMap = iterator_probaShell->second;
>> 182 eProbaShellMap->clear();
>> 183 delete eProbaShellMap;
165 } 184 }
166 << 185 eProbaShellStorage.clear();
167 for (auto const& p : pProbaShellStorage) { << 186
168 delete p.second; << 187 for (iterator_probaShell = pProbaShellStorage.begin(); iterator_probaShell != pProbaShellStorage.end(); ++iterator_probaShell) {
>> 188 vector<VecMap>* pProbaShellMap = iterator_probaShell->second;
>> 189 pProbaShellMap->clear();
>> 190 delete pProbaShellMap;
169 } 191 }
>> 192 pProbaShellStorage.clear();
170 193
171 // (4) 194 // (4)
172 for (auto const& p : eIncidentEnergyStorage) << 195 TranfEnergyMap::iterator iterator_nrjtransf;
173 delete p.second; << 196 for (iterator_nrjtransf = eVecmStorage.begin(); iterator_nrjtransf != eVecmStorage.end(); ++iterator_nrjtransf) {
>> 197 VecMap* eVecm = iterator_nrjtransf->second;
>> 198 eVecm->clear();
>> 199 delete eVecm;
174 } 200 }
175 << 201 eVecmStorage.clear();
176 for (auto const& p : pIncidentEnergyStorage) << 202 for (iterator_nrjtransf = pVecmStorage.begin(); iterator_nrjtransf != pVecmStorage.end(); ++iterator_nrjtransf) {
177 delete p.second; << 203 VecMap* pVecm = iterator_nrjtransf->second;
>> 204 pVecm->clear();
>> 205 delete pVecm;
178 } 206 }
179 << 207 pVecmStorage.clear();
>> 208
180 // (5) 209 // (5)
181 for (auto const& p : eVecmStorage) { << 210 incidentEnergyMap::iterator iterator_energy;
182 delete p.second; << 211 for (iterator_energy = eIncidentEnergyStorage.begin(); iterator_energy != eIncidentEnergyStorage.end(); ++iterator_energy) {
>> 212 std::vector<G4double>* eTdummyVec = iterator_energy->second;
>> 213 eTdummyVec->clear();
>> 214 delete eTdummyVec;
183 } 215 }
184 << 216 eIncidentEnergyStorage.clear();
185 for (auto const& p : pVecmStorage) { << 217
186 delete p.second; << 218 for (iterator_energy = pIncidentEnergyStorage.begin(); iterator_energy != pIncidentEnergyStorage.end(); ++iterator_energy) {
>> 219 std::vector<G4double>* pTdummyVec = iterator_energy->second;
>> 220 pTdummyVec->clear();
>> 221 delete pTdummyVec;
187 } 222 }
188 << 223 pIncidentEnergyStorage.clear();
>> 224
189 // (6) 225 // (6)
190 for (auto const& p : tableMaterialsStructure << 226 MapStructure::iterator iterator_matStructure;
191 delete p.second; << 227 for (iterator_matStructure = tableMaterialsStructures.begin();
>> 228 iterator_matStructure != tableMaterialsStructures.end(); ++iterator_matStructure) {
>> 229 currentMaterialStructure = iterator_matStructure->second;
>> 230 delete currentMaterialStructure;
192 } 231 }
>> 232 tableMaterialsStructures.clear();
>> 233 currentMaterialStructure = nullptr;
193 } 234 }
194 235
195 //....oooOO0OOooo........oooOO0OOooo........oo 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196 237
197 void G4MicroElecInelasticModel_new::Initialise 238 void G4MicroElecInelasticModel_new::Initialise(const G4ParticleDefinition* particle,
198 const G4DataVector& /*cuts*/) 239 const G4DataVector& /*cuts*/)
199 { 240 {
200 if (isInitialised) { return; } 241 if (isInitialised) { return; }
201 242
202 if (verboseLevel > 3) 243 if (verboseLevel > 3)
203 G4cout << "Calling G4MicroElecInelasticMod 244 G4cout << "Calling G4MicroElecInelasticModel_new::Initialise()" << G4endl;
204 245
205 const char* path = G4FindDataDir("G4LEDATA") 246 const char* path = G4FindDataDir("G4LEDATA");
206 if (!path) 247 if (!path)
207 { 248 {
208 G4Exception("G4MicroElecElasticModel_new 249 G4Exception("G4MicroElecElasticModel_new::Initialise","em0006",FatalException,"G4LEDATA environment variable not set.");
209 return; 250 return;
210 } 251 }
211 252
212 G4String modelName = "mermin"; 253 G4String modelName = "mermin";
213 G4cout << "****************************" << 254 G4cout << "****************************" << G4endl;
214 G4cout << modelName << " model loaded !" << 255 G4cout << modelName << " model loaded !" << G4endl;
215 256
216 // Energy limits 257 // Energy limits
217 G4ParticleDefinition* electronDef = G4Electr 258 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
218 G4ParticleDefinition* protonDef = G4Proton:: 259 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
219 G4String electron = electronDef->GetParticle 260 G4String electron = electronDef->GetParticleName();
220 G4String proton = protonDef->GetParticleName 261 G4String proton = protonDef->GetParticleName();
221 262
222 G4double scaleFactor = 1.0; 263 G4double scaleFactor = 1.0;
223 264
224 // *** ELECTRON 265 // *** ELECTRON
225 lowEnergyLimit[electron] = 2 * eV; 266 lowEnergyLimit[electron] = 2 * eV;
226 highEnergyLimit[electron] = 10.0 * MeV; 267 highEnergyLimit[electron] = 10.0 * MeV;
227 268
228 // Cross section 269 // Cross section
229 G4ProductionCutsTable* theCoupleTable = G4Pr 270 G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
230 G4int numOfCouples = (G4int)theCoupleTable-> 271 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
231 272
232 for (G4int i = 0; i < numOfCouples; ++i) { 273 for (G4int i = 0; i < numOfCouples; ++i) {
233 const G4Material* material = theCoupleTabl 274 const G4Material* material = theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
234 G4cout << "Material " << i + 1 << " / " << 275 G4cout << "Material " << i + 1 << " / " << numOfCouples << " : " << material->GetName() << G4endl;
235 if (material->GetName() == "Vacuum") conti 276 if (material->GetName() == "Vacuum") continue;
236 G4String mat = material->GetName().substr( 277 G4String mat = material->GetName().substr(3, material->GetName().size());
237 MapData* tableData = new MapData; 278 MapData* tableData = new MapData;
238 currentMaterialStructure = new G4MicroElec 279 currentMaterialStructure = new G4MicroElecMaterialStructure(mat);
239 280
240 tableMaterialsStructures[mat] = currentMat 281 tableMaterialsStructures[mat] = currentMaterialStructure;
241 if (particle == electronDef) { 282 if (particle == electronDef) {
242 //TCS 283 //TCS
243 G4String fileElectron("Inelastic/" + mod 284 G4String fileElectron("Inelastic/" + modelName + "_sigma_inelastic_e-_" + mat);
244 G4cout << fileElectron << G4endl; 285 G4cout << fileElectron << G4endl;
245 G4MicroElecCrossSectionDataSet_new* tabl 286 G4MicroElecCrossSectionDataSet_new* tableE = new G4MicroElecCrossSectionDataSet_new(new G4LogLogInterpolation, MeV, scaleFactor);
246 tableE->LoadData(fileElectron); 287 tableE->LoadData(fileElectron);
247 tableData->insert(make_pair(electron, ta 288 tableData->insert(make_pair(electron, tableE));
248 289
249 // DCS 290 // DCS
250 std::ostringstream eFullFileName; 291 std::ostringstream eFullFileName;
251 if (fasterCode) { 292 if (fasterCode) {
252 eFullFileName << path << "/microelec/Inelast 293 eFullFileName << path << "/microelec/Inelastic/cumulated_" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat";
253 G4cout << "Faster code = true" << G4endl; 294 G4cout << "Faster code = true" << G4endl;
254 G4cout << "Inelastic/cumulated_" + modelName 295 G4cout << "Inelastic/cumulated_" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat" << G4endl;
255 } 296 }
256 else { 297 else {
257 eFullFileName << path << "/microelec/Inelast 298 eFullFileName << path << "/microelec/Inelastic/" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat";
258 G4cout << "Faster code = false" << G4endl; 299 G4cout << "Faster code = false" << G4endl;
259 G4cout << "Inelastic/" + modelName + "_sigma 300 G4cout << "Inelastic/" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat" << G4endl;
260 } 301 }
261 302
262 std::ifstream eDiffCrossSection(eFullFil 303 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
263 if (!eDiffCrossSection) 304 if (!eDiffCrossSection)
264 { 305 {
265 std::stringstream ss; 306 std::stringstream ss;
266 ss << "Missing data " << eFullFileName.str 307 ss << "Missing data " << eFullFileName.str().c_str();
267 std::string sortieString = ss.str(); 308 std::string sortieString = ss.str();
268 309
269 if (fasterCode) G4Exception("G4MicroElecIn 310 if (fasterCode) G4Exception("G4MicroElecInelasticModel_new::Initialise", "em0003",
270 FatalException, sortieString.c_s 311 FatalException, sortieString.c_str());
>> 312
271 else { 313 else {
272 G4Exception("G4MicroElecInelasticModel_n 314 G4Exception("G4MicroElecInelasticModel_new::Initialise", "em0003",
273 FatalException, "Missing data file:/micr 315 FatalException, "Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
274 } 316 }
275 } 317 }
276 318
277 // Clear the arrays for re-initializatio 319 // Clear the arrays for re-initialization case (MT mode)
278 // Octobre 22nd, 2014 - Melanie Raine 320 // Octobre 22nd, 2014 - Melanie Raine
279 //Creating vectors of maps for DCS and C 321 //Creating vectors of maps for DCS and Cumulated DCS for the current material.
280 //Each vector is storing one map for eac 322 //Each vector is storing one map for each shell.
281 G4bool isUsed1 = false; <<
282 vector<TriDimensionMap>* eDiffCrossSecti 323 vector<TriDimensionMap>* eDiffCrossSectionData =
283 new vector<TriDimensionMap>; //Storage of [I 324 new vector<TriDimensionMap>; //Storage of [IncidentEnergy, TransfEnergy, DCS values], used in slower code
284 vector<TriDimensionMap>* eNrjTransfData 325 vector<TriDimensionMap>* eNrjTransfData =
285 new vector<TriDimensionMap>; //Storage of po 326 new vector<TriDimensionMap>; //Storage of possible transfer energies by shell
286 vector<VecMap>* eProbaShellMap = new vec 327 vector<VecMap>* eProbaShellMap = new vector<VecMap>; //Storage of the vectors containing all cumulated DCS values for an initial energy, by shell
287 vector<G4double>* eTdummyVec = new vecto 328 vector<G4double>* eTdummyVec = new vector<G4double>; //Storage of incident energies for interpolation
288 VecMap* eVecm = new VecMap; //Transfered 329 VecMap* eVecm = new VecMap; //Transfered energy map for slower code
289 330
290 for (G4int j = 0; j < currentMaterialStr << 331 for (int j = 0; j < currentMaterialStructure->NumberOfLevels(); j++) //Filling the map vectors with an empty map for each shell
291 { 332 {
292 eDiffCrossSectionData->push_back(TriDimens 333 eDiffCrossSectionData->push_back(TriDimensionMap());
293 eNrjTransfData->push_back(TriDimensionMap( 334 eNrjTransfData->push_back(TriDimensionMap());
294 eProbaShellMap->push_back(VecMap()); 335 eProbaShellMap->push_back(VecMap());
295 } 336 }
296 337
297 eTdummyVec->push_back(0.); 338 eTdummyVec->push_back(0.);
298 while (!eDiffCrossSection.eof()) 339 while (!eDiffCrossSection.eof())
299 { 340 {
300 G4double tDummy; //incident energy 341 G4double tDummy; //incident energy
301 G4double eDummy; //transfered energy 342 G4double eDummy; //transfered energy
302 eDiffCrossSection >> tDummy >> eDummy; 343 eDiffCrossSection >> tDummy >> eDummy;
303 if (tDummy != eTdummyVec->back()) eTdummyV 344 if (tDummy != eTdummyVec->back()) eTdummyVec->push_back(tDummy);
304 345
305 G4double tmp; //probability 346 G4double tmp; //probability
306 for (G4int j = 0; j < currentMaterialStruc << 347 for (int j = 0; j < currentMaterialStructure->NumberOfLevels(); j++)
307 { 348 {
308 eDiffCrossSection >> tmp; 349 eDiffCrossSection >> tmp;
309 (*eDiffCrossSectionData)[j][tDummy][eD 350 (*eDiffCrossSectionData)[j][tDummy][eDummy] = tmp;
310 351
311 if (fasterCode) 352 if (fasterCode)
312 { 353 {
313 (*eNrjTransfData)[j][tDummy][(*eDiffCros 354 (*eNrjTransfData)[j][tDummy][(*eDiffCrossSectionData)[j][tDummy][eDummy]] = eDummy;
314 (*eProbaShellMap)[j][tDummy].push_back(( 355 (*eProbaShellMap)[j][tDummy].push_back((*eDiffCrossSectionData)[j][tDummy][eDummy]);
315 } 356 }
316 else { // SI - only if eof is not rea 357 else { // SI - only if eof is not reached !
317 (*eDiffCrossSectionData)[j][tDummy][eDummy << 358 if (!eDiffCrossSection.eof()) (*eDiffCrossSectionData)[j][tDummy][eDummy] *= scaleFactor;
318 (*eVecm)[tDummy].push_back(eDummy); 359 (*eVecm)[tDummy].push_back(eDummy);
319 } 360 }
320 } 361 }
321 } 362 }
322 // 363 //
323 //G4cout << "add to material vector" << << 364 G4cout << "add to material vector" << G4endl;
324 365
325 //Filing maps for the current material i 366 //Filing maps for the current material into the master maps
326 if (fasterCode) { 367 if (fasterCode) {
327 isUsed1 = true; <<
328 eNrjTransStorage[mat] = eNrjTransfData; 368 eNrjTransStorage[mat] = eNrjTransfData;
329 eProbaShellStorage[mat] = eProbaShellMap; 369 eProbaShellStorage[mat] = eProbaShellMap;
330 } 370 }
331 else { 371 else {
332 eDiffDatatable[mat] = eDiffCrossSectionData; 372 eDiffDatatable[mat] = eDiffCrossSectionData;
333 eVecmStorage[mat] = eVecm; 373 eVecmStorage[mat] = eVecm;
334 } 374 }
335 eIncidentEnergyStorage[mat] = eTdummyVec 375 eIncidentEnergyStorage[mat] = eTdummyVec;
336 376
337 //Cleanup support vectors 377 //Cleanup support vectors
338 if (!isUsed1) { << 378 // delete eProbaShellMap;
339 delete eProbaShellMap; << 379 // delete eDiffCrossSectionData;
340 delete eNrjTransfData; << 380 // delete eNrjTransfData;
341 } else { <<
342 delete eDiffCrossSectionData; <<
343 delete eVecm; <<
344 } <<
345 } 381 }
346 382
347 // *** PROTON 383 // *** PROTON
348 if (particle == protonDef) 384 if (particle == protonDef)
349 { 385 {
350 // Cross section 386 // Cross section
351 G4String fileProton("Inelastic/" + modelName 387 G4String fileProton("Inelastic/" + modelName + "_sigma_inelastic_p_" + mat); G4cout << fileProton << G4endl;
352 G4MicroElecCrossSectionDataSet_new* tableP = 388 G4MicroElecCrossSectionDataSet_new* tableP = new G4MicroElecCrossSectionDataSet_new(new G4LogLogInterpolation, MeV, scaleFactor);
353 tableP->LoadData(fileProton); 389 tableP->LoadData(fileProton);
354 tableData->insert(make_pair(proton, tableP)) 390 tableData->insert(make_pair(proton, tableP));
355 391
356 // DCS 392 // DCS
357 std::ostringstream pFullFileName; 393 std::ostringstream pFullFileName;
358 if (fasterCode) { 394 if (fasterCode) {
359 pFullFileName << path << "/microelec/Inela 395 pFullFileName << path << "/microelec/Inelastic/cumulated_" + modelName + "_sigmadiff_inelastic_p_" + mat + ".dat";
360 G4cout << "Faster code = true" << G4endl; 396 G4cout << "Faster code = true" << G4endl;
361 G4cout << "Inelastic/cumulated_" + modelNa 397 G4cout << "Inelastic/cumulated_" + modelName + "_sigmadiff_inelastic_p_" + mat + ".dat" << G4endl;
362 } 398 }
363 else { 399 else {
364 pFullFileName << path << "/microelec/Inela 400 pFullFileName << path << "/microelec/Inelastic/" + modelName + "_sigmadiff_inelastic_p_" + mat + ".dat";
365 G4cout << "Faster code = false" << G4endl; 401 G4cout << "Faster code = false" << G4endl;
366 G4cout << "Inelastic/" + modelName + "_sig 402 G4cout << "Inelastic/" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat" << G4endl;
367 } 403 }
368 404
369 std::ifstream pDiffCrossSection(pFullFileNam 405 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
370 if (!pDiffCrossSection) 406 if (!pDiffCrossSection)
371 { 407 {
372 if (fasterCode) G4Exception("G4MicroElec 408 if (fasterCode) G4Exception("G4MicroElecInelasticModel_new::Initialise", "em0003",
373 FatalException, "Missing data file:/ 409 FatalException, "Missing data file:/microelec/sigmadiff_cumulated_inelastic_p_Si.dat");
374 else { 410 else {
375 G4Exception("G4MicroElecInelasticModel 411 G4Exception("G4MicroElecInelasticModel_new::Initialise", "em0003",
376 FatalException, "Missing data file:/mi 412 FatalException, "Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
377 } 413 }
378 } 414 }
379 415
380 // 416 //
381 // Clear the arrays for re-initialization ca 417 // Clear the arrays for re-initialization case (MT mode)
382 // Octobre 22nd, 2014 - Melanie Raine 418 // Octobre 22nd, 2014 - Melanie Raine
383 //Creating vectors of maps for DCS and Cumul 419 //Creating vectors of maps for DCS and Cumulated DCS for the current material.
384 //Each vector is storing one map for each sh 420 //Each vector is storing one map for each shell.
385 421
386 G4bool isUsed1 = false; <<
387 vector<TriDimensionMap>* pDiffCrossSectionDa 422 vector<TriDimensionMap>* pDiffCrossSectionData =
388 new vector<TriDimensionMap>; //Storage of 423 new vector<TriDimensionMap>; //Storage of [IncidentEnergy, TransfEnergy, DCS values], used in slower code
389 vector<TriDimensionMap>* pNrjTransfData = 424 vector<TriDimensionMap>* pNrjTransfData =
390 new vector<TriDimensionMap>; //Storage of 425 new vector<TriDimensionMap>; //Storage of possible transfer energies by shell
391 vector<VecMap>* pProbaShellMap = 426 vector<VecMap>* pProbaShellMap =
392 new vector<VecMap>; //Storage of the vecto 427 new vector<VecMap>; //Storage of the vectors containing all cumulated DCS values for an initial energy, by shell
393 vector<G4double>* pTdummyVec = 428 vector<G4double>* pTdummyVec =
394 new vector<G4double>; //Storage of inciden 429 new vector<G4double>; //Storage of incident energies for interpolation
395 VecMap* pVecm = new VecMap; //Transfered ene << 430 VecMap* eVecm = new VecMap; //Transfered energy map for slower code
396 431
397 for (G4int j = 0; j < currentMaterialStructu << 432 for (int j = 0; j < currentMaterialStructure->NumberOfLevels(); ++j)
398 //Filling the map vectors with an empty ma 433 //Filling the map vectors with an empty map for each shell
399 { 434 {
400 pDiffCrossSectionData->push_back(TriDime 435 pDiffCrossSectionData->push_back(TriDimensionMap());
401 pNrjTransfData->push_back(TriDimensionMa 436 pNrjTransfData->push_back(TriDimensionMap());
402 pProbaShellMap->push_back(VecMap()); 437 pProbaShellMap->push_back(VecMap());
403 } 438 }
404 439
405 pTdummyVec->push_back(0.); 440 pTdummyVec->push_back(0.);
406 while (!pDiffCrossSection.eof()) 441 while (!pDiffCrossSection.eof())
407 { 442 {
408 G4double tDummy; //incident energy 443 G4double tDummy; //incident energy
409 G4double eDummy; //transfered energy 444 G4double eDummy; //transfered energy
410 pDiffCrossSection >> tDummy >> eDummy; 445 pDiffCrossSection >> tDummy >> eDummy;
411 if (tDummy != pTdummyVec->back()) pTdumm 446 if (tDummy != pTdummyVec->back()) pTdummyVec->push_back(tDummy);
412 447
413 G4double tmp; //probability 448 G4double tmp; //probability
414 for (G4int j = 0; j < currentMaterialStr << 449 for (int j = 0; j < currentMaterialStructure->NumberOfLevels(); j++)
415 { 450 {
416 pDiffCrossSection >> tmp; 451 pDiffCrossSection >> tmp;
417 (*pDiffCrossSectionData)[j][tDummy][eDummy 452 (*pDiffCrossSectionData)[j][tDummy][eDummy] = tmp;
418 // ArrayofMaps[j] -> fill with 3DMap(incid 453 // ArrayofMaps[j] -> fill with 3DMap(incidentEnergy,
419 // 2Dmap (transferedEnergy,proba=tmp) ) -> 454 // 2Dmap (transferedEnergy,proba=tmp) ) -> fill map for shell j
420 // with proba for transfered energy eDummy 455 // with proba for transfered energy eDummy
421 456
422 if (fasterCode) 457 if (fasterCode)
423 { 458 {
424 (*pNrjTransfData)[j][tDummy][(*pDiffCr 459 (*pNrjTransfData)[j][tDummy][(*pDiffCrossSectionData)[j][tDummy][eDummy]] = eDummy;
425 (*pProbaShellMap)[j][tDummy].push_back 460 (*pProbaShellMap)[j][tDummy].push_back((*pDiffCrossSectionData)[j][tDummy][eDummy]);
426 } 461 }
427 else { // SI - only if eof is not reached 462 else { // SI - only if eof is not reached !
428 (*pDiffCrossSectionData)[j][tDummy][eDum << 463 if (!pDiffCrossSection.eof()) (*pDiffCrossSectionData)[j][tDummy][eDummy] *= scaleFactor;
429 (*pVecm)[tDummy].push_back(eDummy); << 464 (*eVecm)[tDummy].push_back(eDummy);
430 } 465 }
431 } 466 }
432 } 467 }
433 468
434 //Filing maps for the current material into 469 //Filing maps for the current material into the master maps
435 if (fasterCode) { 470 if (fasterCode) {
436 isUsed1 = true; <<
437 pNrjTransStorage[mat] = pNrjTransfData; 471 pNrjTransStorage[mat] = pNrjTransfData;
438 pProbaShellStorage[mat] = pProbaShellMap; 472 pProbaShellStorage[mat] = pProbaShellMap;
439 } 473 }
440 else { 474 else {
441 pDiffDatatable[mat] = pDiffCrossSectionDat 475 pDiffDatatable[mat] = pDiffCrossSectionData;
442 pVecmStorage[mat] = pVecm; << 476 pVecmStorage[mat] = eVecm;
443 } 477 }
444 pIncidentEnergyStorage[mat] = pTdummyVec; 478 pIncidentEnergyStorage[mat] = pTdummyVec;
445 << 479
446 //Cleanup support vectors 480 //Cleanup support vectors
447 if (!isUsed1) { << 481 // delete pNrjTransfData;
448 delete pProbaShellMap; << 482 // delete eVecm;
449 delete pNrjTransfData; << 483 // delete pDiffCrossSectionData;
450 } else { << 484 // delete pProbaShellMap;
451 delete pDiffCrossSectionData; <<
452 delete pVecm; <<
453 } <<
454 } 485 }
455 tableTCS[mat] = tableData; 486 tableTCS[mat] = tableData;
456 } << 487 }
457 if (particle==electronDef) 488 if (particle==electronDef)
458 { 489 {
459 SetLowEnergyLimit(lowEnergyLimit[electro 490 SetLowEnergyLimit(lowEnergyLimit[electron]);
460 SetHighEnergyLimit(highEnergyLimit[elect 491 SetHighEnergyLimit(highEnergyLimit[electron]);
461 } 492 }
462 493
463 if (particle==protonDef) 494 if (particle==protonDef)
464 { 495 {
465 SetLowEnergyLimit(100*eV); 496 SetLowEnergyLimit(100*eV);
466 SetHighEnergyLimit(10*MeV); 497 SetHighEnergyLimit(10*MeV);
467 } 498 }
468 499
469 if( verboseLevel>1 ) 500 if( verboseLevel>1 )
470 { 501 {
471 G4cout << "MicroElec Inelastic model is 502 G4cout << "MicroElec Inelastic model is initialized " << G4endl
472 << "Energy range: " 503 << "Energy range: "
473 << LowEnergyLimit() / keV << " keV - " 504 << LowEnergyLimit() / keV << " keV - "
474 << HighEnergyLimit() / MeV << " MeV for 505 << HighEnergyLimit() / MeV << " MeV for "
475 << particle->GetParticleName() 506 << particle->GetParticleName()
476 << " with mass (amu) " << particle->Get 507 << " with mass (amu) " << particle->GetPDGMass()/proton_mass_c2
477 << " and charge " << particle->GetPDGCh 508 << " and charge " << particle->GetPDGCharge()
478 << G4endl << G4endl ; 509 << G4endl << G4endl ;
479 } 510 }
480 511
481 fAtomDeexcitation = G4LossTableManager::Ins 512 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
482 << 513
483 fParticleChangeForGamma = GetParticleChangeF 514 fParticleChangeForGamma = GetParticleChangeForGamma();
484 isInitialised = true; 515 isInitialised = true;
485 } 516 }
486 517
487 //....oooOO0OOooo........oooOO0OOooo........oo 518 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
488 519
489 G4double G4MicroElecInelasticModel_new::CrossS 520 G4double G4MicroElecInelasticModel_new::CrossSectionPerVolume(const G4Material* material,
490 const G4ParticleDefinition* par 521 const G4ParticleDefinition* particleDefinition,
491 G4double ekin, 522 G4double ekin,
492 G4double, 523 G4double,
493 G4double) 524 G4double)
494 { 525 {
495 if (verboseLevel > 3) G4cout << "Calling Cro 526 if (verboseLevel > 3) G4cout << "Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" << G4endl;
496 527
497 G4double density = material->GetTotNbOfAtoms 528 G4double density = material->GetTotNbOfAtomsPerVolume();
498 currentMaterial = material->GetName().substr 529 currentMaterial = material->GetName().substr(3, material->GetName().size());
499 530
500 MapStructure::iterator structPos; 531 MapStructure::iterator structPos;
501 structPos = tableMaterialsStructures.find(cu 532 structPos = tableMaterialsStructures.find(currentMaterial);
502 533
503 // Calculate total cross section for model 534 // Calculate total cross section for model
504 TCSMap::iterator tablepos; 535 TCSMap::iterator tablepos;
505 tablepos = tableTCS.find(currentMaterial); 536 tablepos = tableTCS.find(currentMaterial);
506 537
507 if (tablepos == tableTCS.end() ) 538 if (tablepos == tableTCS.end() )
508 { 539 {
509 G4String str = "Material "; 540 G4String str = "Material ";
510 str += currentMaterial + " TCS Table not 541 str += currentMaterial + " TCS Table not found!";
511 G4Exception("G4MicroElecInelasticModel_n 542 G4Exception("G4MicroElecInelasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, str);
512 return 0; 543 return 0;
513 } 544 }
514 else if(structPos == tableMaterialsStructure 545 else if(structPos == tableMaterialsStructures.end())
515 { 546 {
516 G4String str = "Material "; 547 G4String str = "Material ";
517 str += currentMaterial + " Structure not 548 str += currentMaterial + " Structure not found!";
518 G4Exception("G4MicroElecInelasticModel_n 549 G4Exception("G4MicroElecInelasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, str);
519 return 0; 550 return 0;
520 } 551 }
521 else { 552 else {
522 MapData* tableData = tablepos->second; 553 MapData* tableData = tablepos->second;
523 currentMaterialStructure = structPos->seco 554 currentMaterialStructure = structPos->second;
524 555
525 G4double sigma = 0; 556 G4double sigma = 0;
526 557
527 const G4String& particleName = particleDef 558 const G4String& particleName = particleDefinition->GetParticleName();
528 G4String nameLocal = particleName; 559 G4String nameLocal = particleName;
529 G4int pdg = particleDefinition->GetPDGEnco 560 G4int pdg = particleDefinition->GetPDGEncoding();
530 G4int Z = particleDefinition->GetAtomicNum 561 G4int Z = particleDefinition->GetAtomicNumber();
531 562
532 G4double Zeff = 1.0, Zeff2 = Zeff*Zeff; 563 G4double Zeff = 1.0, Zeff2 = Zeff*Zeff;
533 G4double Mion_c2 = particleDefinition->Get 564 G4double Mion_c2 = particleDefinition->GetPDGMass();
534 565
535 if (Mion_c2 > proton_mass_c2) 566 if (Mion_c2 > proton_mass_c2)
536 { 567 {
537 ekin *= proton_mass_c2 / Mion_c2; 568 ekin *= proton_mass_c2 / Mion_c2;
538 nameLocal = "proton"; 569 nameLocal = "proton";
539 } 570 }
540 571
541 G4double lowLim = currentMaterialStructure 572 G4double lowLim = currentMaterialStructure->GetInelasticModelLowLimit(pdg);
542 G4double highLim = currentMaterialStructur 573 G4double highLim = currentMaterialStructure->GetInelasticModelHighLimit(pdg);
543 574
544 if (ekin >= lowLim && ekin < highLim) 575 if (ekin >= lowLim && ekin < highLim)
545 { 576 {
546 std::map< G4String, G4MicroElecCrossSectionD 577 std::map< G4String, G4MicroElecCrossSectionDataSet_new*, std::less<G4String> >::iterator pos;
547 pos = tableData->find(nameLocal); //find par 578 pos = tableData->find(nameLocal); //find particle type
548 579
549 if (pos != tableData->end()) 580 if (pos != tableData->end())
550 { 581 {
551 G4MicroElecCrossSectionDataSet_new* tabl 582 G4MicroElecCrossSectionDataSet_new* table = pos->second;
552 if (table != 0) 583 if (table != 0)
553 { 584 {
554 sigma = table->FindValue(ekin); 585 sigma = table->FindValue(ekin);
555 586
556 if (Mion_c2 > proton_mass_c2) { 587 if (Mion_c2 > proton_mass_c2) {
557 sigma = 0.; 588 sigma = 0.;
558 for (G4int i = 0; i < currentMaterialStr 589 for (G4int i = 0; i < currentMaterialStructure->NumberOfLevels(); i++) {
559 Zeff = BKZ(ekin / (proton_mass_c2 / Mi 590 Zeff = BKZ(ekin / (proton_mass_c2 / Mion_c2), Mion_c2 / c_squared, Z, currentMaterialStructure->Energy(i)); // il faut garder le vrai ekin car le calcul à l'interieur de la methode convertie l'énergie en vitesse
560 Zeff2 = Zeff*Zeff; 591 Zeff2 = Zeff*Zeff;
561 sigma += Zeff2*table->FindShellValue(e 592 sigma += Zeff2*table->FindShellValue(ekin, i);
562 // il faut utiliser le ekin mis à l'echelle p 593 // il faut utiliser le ekin mis à l'echelle pour chercher la bonne
563 // valeur dans les tables proton 594 // valeur dans les tables proton
564 595
565 } 596 }
566 } 597 }
567 else { 598 else {
568 sigma = table->FindValue(ekin); 599 sigma = table->FindValue(ekin);
569 } 600 }
570 } 601 }
571 } 602 }
572 else 603 else
573 { 604 {
574 G4Exception("G4MicroElecInelasticModel_n 605 G4Exception("G4MicroElecInelasticModel_new::CrossSectionPerVolume",
575 "em0002", FatalExcepti 606 "em0002", FatalException,
576 "Model not applicable 607 "Model not applicable to particle type.");
577 } 608 }
578 } 609 }
579 else 610 else
580 { 611 {
581 return 1 / DBL_MAX; 612 return 1 / DBL_MAX;
582 } 613 }
583 614
584 if (verboseLevel > 3) 615 if (verboseLevel > 3)
585 { 616 {
586 G4cout << "---> Kinetic energy (eV)=" << eki 617 G4cout << "---> Kinetic energy (eV)=" << ekin / eV << G4endl;
587 G4cout << " - Cross section per Si atom (cm^ 618 G4cout << " - Cross section per Si atom (cm^2)=" << sigma / cm2 << G4endl;
588 G4cout << " - Cross section per Si atom (cm^ 619 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density / (1. / cm) << G4endl;
589 } 620 }
590 621
591 return (sigma)*density; << 622 return (sigma)*density;}
592 } << 623
593 } 624 }
594 625
595 //....oooOO0OOooo........oooOO0OOooo........oo 626 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
596 627
597 void G4MicroElecInelasticModel_new::SampleSeco 628 void G4MicroElecInelasticModel_new::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
598 629 const G4MaterialCutsCouple* couple,
599 630 const G4DynamicParticle* particle,
600 631 G4double,
601 632 G4double)
602 { 633 {
603 634
604 if (verboseLevel > 3) 635 if (verboseLevel > 3)
605 G4cout << "Calling SampleSecondaries() of 636 G4cout << "Calling SampleSecondaries() of G4MicroElecInelasticModel" << G4endl;
606 637
607 G4int pdg = particle->GetParticleDefinition( 638 G4int pdg = particle->GetParticleDefinition()->GetPDGEncoding();
608 G4double lowLim = currentMaterialStructure-> 639 G4double lowLim = currentMaterialStructure->GetInelasticModelLowLimit(pdg);
609 G4double highLim = currentMaterialStructure- 640 G4double highLim = currentMaterialStructure->GetInelasticModelHighLimit(pdg);
610 641
611 G4double ekin = particle->GetKineticEnergy() 642 G4double ekin = particle->GetKineticEnergy();
612 G4double k = ekin ; 643 G4double k = ekin ;
613 644
614 G4ParticleDefinition* PartDef = particle->Ge 645 G4ParticleDefinition* PartDef = particle->GetDefinition();
615 const G4String& particleName = PartDef->GetP 646 const G4String& particleName = PartDef->GetParticleName();
616 G4String nameLocal2 = particleName ; 647 G4String nameLocal2 = particleName ;
617 G4double particleMass = particle->GetDefinit 648 G4double particleMass = particle->GetDefinition()->GetPDGMass();
618 G4double originalMass = particleMass; // a p 649 G4double originalMass = particleMass; // a passer en argument dans samplesecondaryenergy pour évaluer correctement Qmax
619 G4int originalZ = particle->GetDefinition()- 650 G4int originalZ = particle->GetDefinition()->GetAtomicNumber();
620 651
621 if (particleMass > proton_mass_c2) 652 if (particleMass > proton_mass_c2)
622 { 653 {
623 k *= proton_mass_c2/particleMass ; 654 k *= proton_mass_c2/particleMass ;
624 PartDef = G4Proton::ProtonDefinition(); 655 PartDef = G4Proton::ProtonDefinition();
625 nameLocal2 = "proton" ; 656 nameLocal2 = "proton" ;
626 } 657 }
627 658
628 if (k >= lowLim && k < highLim) 659 if (k >= lowLim && k < highLim)
629 { 660 {
630 G4ParticleMomentum primaryDirection = pa 661 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
631 G4double totalEnergy = ekin + particleMa 662 G4double totalEnergy = ekin + particleMass;
632 G4double pSquare = ekin * (totalEnergy + 663 G4double pSquare = ekin * (totalEnergy + particleMass);
633 G4double totalMomentum = std::sqrt(pSqua 664 G4double totalMomentum = std::sqrt(pSquare);
634 665
635 G4int Shell = 1; 666 G4int Shell = 1;
636 667
637 Shell = RandomSelect(k,nameLocal2,origin 668 Shell = RandomSelect(k,nameLocal2,originalMass, originalZ);
638 669
639 G4double bindingEnergy = currentMaterial 670 G4double bindingEnergy = currentMaterialStructure->Energy(Shell);
640 G4double limitEnergy = currentMaterialSt 671 G4double limitEnergy = currentMaterialStructure->GetLimitEnergy(Shell);
641 <<
642 G4bool weaklyBound = currentMaterialStruct <<
643 672
644 if (verboseLevel > 3) 673 if (verboseLevel > 3)
645 { 674 {
646 G4cout << "---> Kinetic energy (eV)=" << k 675 G4cout << "---> Kinetic energy (eV)=" << k/eV << G4endl ;
647 G4cout << "Shell: " << Shell << ", energy: 676 G4cout << "Shell: " << Shell << ", energy: " << bindingEnergy/eV << G4endl;
648 } 677 }
649 678
650 <<
651 if (k<limitEnergy) { <<
652 if (weaklyBound && k > currentMaterialStru <<
653 limitEnergy = currentMaterialStructure-> <<
654 } <<
655 else return; } <<
656 <<
657 // sample deexcitation 679 // sample deexcitation
658 680
659 std::size_t secNumberInit = 0; // need 681 std::size_t secNumberInit = 0; // need to know at a certain point the energy of secondaries
660 std::size_t secNumberFinal = 0; // So I' 682 std::size_t secNumberFinal = 0; // So I'll make the difference and then sum the energies
661 << 683
>> 684 //SI: additional protection if tcs interpolation method is modified
>> 685 //if (k<bindingEnergy) return;
>> 686 if (k<limitEnergy) return;
662 // G4cout << currentMaterial << G4endl; 687 // G4cout << currentMaterial << G4endl;
663 G4int Z = currentMaterialStructure->GetZ 688 G4int Z = currentMaterialStructure->GetZ(Shell);
664 G4int shellEnum = currentMaterialStructu 689 G4int shellEnum = currentMaterialStructure->GetEADL_Enumerator(Shell);
665 if (currentMaterialStructure->IsShellWea 690 if (currentMaterialStructure->IsShellWeaklyBound(Shell)) { shellEnum = -1; }
666 691
667 if(fAtomDeexcitation && shellEnum >=0) 692 if(fAtomDeexcitation && shellEnum >=0)
668 { 693 {
669 // G4cout << "enter if deex and shell 0 694 // G4cout << "enter if deex and shell 0" << G4endl;
670 G4AtomicShellEnumerator as = G4AtomicShell 695 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellEnum);
671 const G4AtomicShell* shell = fAtomDeexcita 696 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
672 secNumberInit = fvect->size(); 697 secNumberInit = fvect->size();
673 fAtomDeexcitation->GenerateParticles(fvect 698 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
674 secNumberFinal = fvect->size(); 699 secNumberFinal = fvect->size();
675 } 700 }
676 701
677 G4double secondaryKinetic=-1000*eV; 702 G4double secondaryKinetic=-1000*eV;
678 SEFromFermiLevel = false; 703 SEFromFermiLevel = false;
679 if (!fasterCode) 704 if (!fasterCode)
680 { 705 {
681 secondaryKinetic = RandomizeEjectedElectro 706 secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef, k, Shell, originalMass, originalZ);
682 } 707 }
683 else 708 else
684 { 709 {
685 secondaryKinetic = RandomizeEjectedElectro 710 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(PartDef, k, Shell) ;
686 } 711 }
687 712
688 if (verboseLevel > 3) 713 if (verboseLevel > 3)
689 { 714 {
690 G4cout << "Ionisation process" << G4endl; 715 G4cout << "Ionisation process" << G4endl;
691 G4cout << "Shell: " << Shell << " Kin. ene 716 G4cout << "Shell: " << Shell << " Kin. energy (eV)=" << k/eV
692 << " Sec. energy (eV)=" << secondaryKinet 717 << " Sec. energy (eV)=" << secondaryKinetic/eV << G4endl;
693 } 718 }
694 G4ThreeVector deltaDirection = 719 G4ThreeVector deltaDirection =
695 GetAngularDistribution()->SampleDirectionFor 720 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
696 Z, Shell, 721 Z, Shell,
697 couple->GetMaterial()); 722 couple->GetMaterial());
698 723
699 if (particle->GetDefinition() == G4Elect 724 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
700 { 725 {
701 G4double deltaTotalMomentum = std::sqrt(se 726 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
702 727
703 G4double finalPx = totalMomentum*primaryDi 728 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
704 G4double finalPy = totalMomentum*primaryDi 729 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
705 G4double finalPz = totalMomentum*primaryDi 730 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
706 G4double finalMomentum = std::sqrt(finalPx 731 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
707 finalPx /= finalMomentum; 732 finalPx /= finalMomentum;
708 finalPy /= finalMomentum; 733 finalPy /= finalMomentum;
709 finalPz /= finalMomentum; 734 finalPz /= finalMomentum;
710 735
711 G4ThreeVector direction; 736 G4ThreeVector direction;
712 direction.set(finalPx,finalPy,finalPz); 737 direction.set(finalPx,finalPy,finalPz);
713 738
714 fParticleChangeForGamma->ProposeMomentumDi 739 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
715 } 740 }
716 else fParticleChangeForGamma->ProposeMom 741 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
717 742
718 // note that secondaryKinetic is the ene 743 // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
719 G4double deexSecEnergy = 0; 744 G4double deexSecEnergy = 0;
720 for (std::size_t j=secNumberInit; j < se 745 for (std::size_t j=secNumberInit; j < secNumberFinal; ++j) {
721 deexSecEnergy = deexSecEnergy + (*fvec 746 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
722 } 747 }
723 // correction CI 12/01/2023 limit energy << 748 if (SEFromFermiLevel) limitEnergy = currentMaterialStructure->GetEnergyGap();
724 //if (SEFromFermiLevel) limitEnergy = curr << 749 fParticleChangeForGamma->SetProposedKineticEnergy(ekin - secondaryKinetic - limitEnergy); //Ef = Ei-(Q-El)-El = Ei-Q
725 //fParticleChangeForGamma->SetProposedKi << 750 fParticleChangeForGamma->ProposeLocalEnergyDeposit(limitEnergy - deexSecEnergy);
726 //fParticleChangeForGamma->ProposeLocalE << 751
727 <<
728 // correction CI 09/03/2022 limit energy <<
729 //if (!SEFromFermiLevel && weaklyBound) li <<
730 //if (!SEFromFermiLevel && weaklyBound) li <<
731 fParticleChangeForGamma->SetProposedKineti <<
732 fParticleChangeForGamma->ProposeLocalEnerg <<
733 <<
734 if (secondaryKinetic>0) 752 if (secondaryKinetic>0)
735 { 753 {
736 G4DynamicParticle* dp = new G4DynamicParti 754 G4DynamicParticle* dp = new G4DynamicParticle(G4Electron::Electron(), deltaDirection, secondaryKinetic); //Esec = Q-El
737 fvect->push_back(dp); 755 fvect->push_back(dp);
738 } 756 }
739 } 757 }
740 } 758 }
741 759
742 //....oooOO0OOooo........oooOO0OOooo........oo 760 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
743 761
744 G4double G4MicroElecInelasticModel_new::Random 762 G4double G4MicroElecInelasticModel_new::RandomizeEjectedElectronEnergy(
745 const G4ParticleDefinition* particleD 763 const G4ParticleDefinition* particleDefinition,
746 G4double k, G4int shell, G4double originalM 764 G4double k, G4int shell, G4double originalMass, G4int)
747 { 765 {
748 G4double secondaryElectronKineticEnergy=0.; 766 G4double secondaryElectronKineticEnergy=0.;
749 if (particleDefinition == G4Electron::Electr 767 if (particleDefinition == G4Electron::ElectronDefinition())
750 { 768 {
751 G4double maximumEnergyTransfer=k; 769 G4double maximumEnergyTransfer=k;
752 G4double crossSectionMaximum = 0.; 770 G4double crossSectionMaximum = 0.;
753 G4double minEnergy = currentMaterialStru 771 G4double minEnergy = currentMaterialStructure->GetLimitEnergy(shell);
754 G4double maxEnergy = maximumEnergyTransf 772 G4double maxEnergy = maximumEnergyTransfer;
755 G4int nEnergySteps = 100; 773 G4int nEnergySteps = 100;
756 774
757 G4double value(minEnergy); 775 G4double value(minEnergy);
758 G4double stpEnergy(std::pow(maxEnergy/va 776 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
759 G4int step(nEnergySteps); 777 G4int step(nEnergySteps);
760 while (step>0) 778 while (step>0)
761 { 779 {
762 --step; 780 --step;
763 G4double differentialCrossSection = 781 G4double differentialCrossSection =
764 DifferentialCrossSection(particleDefinit 782 DifferentialCrossSection(particleDefinition, k, value, shell);
765 crossSectionMaximum = std::max(crossSectio 783 crossSectionMaximum = std::max(crossSectionMaximum, differentialCrossSection);
766 value*=stpEnergy; 784 value*=stpEnergy;
767 } 785 }
768 786
769 do 787 do
770 { 788 {
771 secondaryElectronKineticEnergy = G4Uniform 789 secondaryElectronKineticEnergy = G4UniformRand() *
772 (maximumEnergyTransfer-currentMaterialSt 790 (maximumEnergyTransfer-currentMaterialStructure->GetLimitEnergy(shell));
773 } while(G4UniformRand()*crossSectionMaximum 791 } while(G4UniformRand()*crossSectionMaximum >
774 DifferentialCrossSection(particleDefinitio 792 DifferentialCrossSection(particleDefinition, k,
775 (secondaryElectronKineticEnergy+currentM 793 (secondaryElectronKineticEnergy+currentMaterialStructure->GetLimitEnergy(shell)),shell));
776 // added 12/01/2023 << 794 }
777 return secondaryElectronKineticEnergy; <<
778 } <<
779 else if (particleDefinition == G4Proton::Pro 795 else if (particleDefinition == G4Proton::ProtonDefinition())
780 { 796 {
781 G4double maximumEnergyTransfer = 797 G4double maximumEnergyTransfer =
782 ComputeElasticQmax(k/(proton_mass_c2/origina 798 ComputeElasticQmax(k/(proton_mass_c2/originalMass),
783 currentMaterialStructure->Energy(shel 799 currentMaterialStructure->Energy(shell),
784 originalMass/c_squared, electron_mass 800 originalMass/c_squared, electron_mass_c2/c_squared);
785 801
786 G4double crossSectionMaximum = 0.; 802 G4double crossSectionMaximum = 0.;
787 803
788 G4double minEnergy = currentMaterialStru 804 G4double minEnergy = currentMaterialStructure->GetLimitEnergy(shell);
789 G4double maxEnergy = maximumEnergyTransf 805 G4double maxEnergy = maximumEnergyTransfer;
790 G4int nEnergySteps = 100; 806 G4int nEnergySteps = 100;
791 807
792 G4double value(minEnergy); 808 G4double value(minEnergy);
793 G4double stpEnergy(std::pow(maxEnergy/va 809 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
794 G4int step(nEnergySteps); 810 G4int step(nEnergySteps);
795 811
796 while (step>0) 812 while (step>0)
797 { 813 {
798 --step; 814 --step;
799 G4double differentialCrossSection = 815 G4double differentialCrossSection =
800 DifferentialCrossSection(particleDefinit 816 DifferentialCrossSection(particleDefinition, k, value, shell);
801 crossSectionMaximum = std::max(crossSectio 817 crossSectionMaximum = std::max(crossSectionMaximum, differentialCrossSection);
802 value*=stpEnergy; 818 value*=stpEnergy;
803 } 819 }
804 820
805 G4double energyTransfer = 0.; 821 G4double energyTransfer = 0.;
806 do 822 do
807 { 823 {
808 energyTransfer = G4UniformRand() * maximum 824 energyTransfer = G4UniformRand() * maximumEnergyTransfer;
809 } while(G4UniformRand()*crossSectionMaximum 825 } while(G4UniformRand()*crossSectionMaximum >
810 DifferentialCrossSection(particleDefinitio 826 DifferentialCrossSection(particleDefinition, k,energyTransfer,shell));
811 827
812 secondaryElectronKineticEnergy = 828 secondaryElectronKineticEnergy =
813 energyTransfer-currentMaterialStructure->Get 829 energyTransfer-currentMaterialStructure->GetLimitEnergy(shell);
814 830
815 } 831 }
816 return std::max(secondaryElectronKineticEner 832 return std::max(secondaryElectronKineticEnergy, 0.0);
817 } 833 }
818 834
819 //....oooOO0OOooo........oooOO0OOooo........oo 835 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
820 836
821 G4double G4MicroElecInelasticModel_new::Random 837 G4double G4MicroElecInelasticModel_new::RandomizeEjectedElectronEnergyFromCumulatedDcs(
822 const G4ParticleDefinition* particleD 838 const G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
823 { 839 {
824 G4double secondaryElectronKineticEnergy = 0. 840 G4double secondaryElectronKineticEnergy = 0.;
825 G4double random = G4UniformRand(); 841 G4double random = G4UniformRand();
826 G4bool weaklyBound = currentMaterialStructur << 842
827 G4double transf = TransferedEnergy(particleD << 843 secondaryElectronKineticEnergy = TransferedEnergy(particleDefinition, k, shell, random)
828 if (!weaklyBound) { << 844 - currentMaterialStructure->GetLimitEnergy(shell) ;
829 secondaryElectronKineticEnergy = transf - << 845
830 if(secondaryElectronKineticEnergy <= 0.) { << 846 if (isnan(secondaryElectronKineticEnergy)) { secondaryElectronKineticEnergy = k - currentMaterialStructure->GetLimitEnergy(shell); }
831 secondaryElectronKineticEnergy = 0.0; << 847
832 } << 848 if (secondaryElectronKineticEnergy < 0.) {
833 } << 849 secondaryElectronKineticEnergy = k - currentMaterialStructure->GetEnergyGap();
834 else { << 850 SEFromFermiLevel = true;
835 secondaryElectronKineticEnergy = transf - << 851 }
836 // for weaklybound electrons = gap + aver << 852 return secondaryElectronKineticEnergy;
837 if (secondaryElectronKineticEnergy <= 0.) <<
838 secondaryElectronKineticEnergy = 0.0; <<
839 SEFromFermiLevel = true; <<
840 } <<
841 } <<
842 //corrections CI 07/02/2022 - added <<
843 return secondaryElectronKineticEnergy; <<
844 } 853 }
845 854
846 //....oooOO0OOooo........oooOO0OOooo........oo 855 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
847 856
848 G4double G4MicroElecInelasticModel_new::Transf 857 G4double G4MicroElecInelasticModel_new::TransferedEnergy(
849 const G4ParticleDefinition* particleDefinit 858 const G4ParticleDefinition* particleDefinition,
850 G4double k, 859 G4double k,
851 G4int ionizationLevelIndex, 860 G4int ionizationLevelIndex,
852 G4double random) 861 G4double random)
853 { 862 {
854 G4double nrj = 0.; 863 G4double nrj = 0.;
855 G4double valueK1 = 0; 864 G4double valueK1 = 0;
856 G4double valueK2 = 0; 865 G4double valueK2 = 0;
857 G4double valuePROB21 = 0; 866 G4double valuePROB21 = 0;
858 G4double valuePROB22 = 0; 867 G4double valuePROB22 = 0;
859 G4double valuePROB12 = 0; 868 G4double valuePROB12 = 0;
860 G4double valuePROB11 = 0; 869 G4double valuePROB11 = 0;
861 G4double nrjTransf11 = 0; 870 G4double nrjTransf11 = 0;
862 G4double nrjTransf12 = 0; 871 G4double nrjTransf12 = 0;
863 G4double nrjTransf21 = 0; 872 G4double nrjTransf21 = 0;
864 G4double nrjTransf22 = 0; 873 G4double nrjTransf22 = 0;
865 874
866 G4double maximumEnergyTransfer1 = 0; 875 G4double maximumEnergyTransfer1 = 0;
867 G4double maximumEnergyTransfer2 = 0; 876 G4double maximumEnergyTransfer2 = 0;
868 G4double maximumEnergyTransferP = 4.* (elect 877 G4double maximumEnergyTransferP = 4.* (electron_mass_c2 / proton_mass_c2) * k;
869 G4double bindingEnergy = currentMaterialStru 878 G4double bindingEnergy = currentMaterialStructure->GetLimitEnergy(ionizationLevelIndex);
870 879
871 if (particleDefinition == G4Electron::Electr 880 if (particleDefinition == G4Electron::ElectronDefinition())
872 { 881 {
873 dataDiffCSMap::iterator iterator_Nrj; 882 dataDiffCSMap::iterator iterator_Nrj;
874 iterator_Nrj = eNrjTransStorage.find(cur 883 iterator_Nrj = eNrjTransStorage.find(currentMaterial);
875 884
876 dataProbaShellMap::iterator iterator_Pro 885 dataProbaShellMap::iterator iterator_Proba;
877 iterator_Proba = eProbaShellStorage.find 886 iterator_Proba = eProbaShellStorage.find(currentMaterial);
878 887
879 incidentEnergyMap::iterator iterator_Tdu 888 incidentEnergyMap::iterator iterator_Tdummy;
880 iterator_Tdummy = eIncidentEnergyStorage 889 iterator_Tdummy = eIncidentEnergyStorage.find(currentMaterial);
881 890
882 if(iterator_Nrj == eNrjTransStorage.end( 891 if(iterator_Nrj == eNrjTransStorage.end() || iterator_Proba == eProbaShellStorage.end() ||
883 iterator_Tdummy == eIncidentEnergyStorage.e 892 iterator_Tdummy == eIncidentEnergyStorage.end())
884 { 893 {
885 G4String str = "Material "; 894 G4String str = "Material ";
886 str += currentMaterial + " not found!"; 895 str += currentMaterial + " not found!";
887 G4Exception("G4MicroElecInelasticModel_new 896 G4Exception("G4MicroElecInelasticModel_new::TransferedEnergy", "em0002",
888 FatalException, str); 897 FatalException, str);
889 } 898 }
890 else { 899 else {
891 vector<TriDimensionMap>* eNrjTransfData = it 900 vector<TriDimensionMap>* eNrjTransfData = iterator_Nrj->second; //Storage of possible transfer energies
892 vector<VecMap>* eProbaShellMap = iterator_Pr 901 vector<VecMap>* eProbaShellMap = iterator_Proba->second; //Storage of probabilities for energy transfer
893 vector<G4double>* eTdummyVec = iterator_Tdum 902 vector<G4double>* eTdummyVec = iterator_Tdummy->second; //Incident energies for interpolation
894 903
895 // k should be in eV auto, :std::vector<doub << 904 // k should be in eV
896 auto k2 = std::upper_bound(eTdummyVec->begin 905 auto k2 = std::upper_bound(eTdummyVec->begin(),
897 eTdummyVec->end(), 906 eTdummyVec->end(),
898 k); 907 k);
899 auto k1 = k2 - 1; 908 auto k1 = k2 - 1;
900 909
901 // SI : the following condition avoids situa 910 // SI : the following condition avoids situations where random >last vector element
902 if (random <= (*eProbaShellMap)[ionizationLe 911 if (random <= (*eProbaShellMap)[ionizationLevelIndex][(*k1)].back()
903 && random <= (*eProbaShellMap)[ionizatio 912 && random <= (*eProbaShellMap)[ionizationLevelIndex][(*k2)].back())
904 { 913 {
905 auto prob12 = 914 auto prob12 =
906 std::upper_bound((*eProbaShellMap)[ion 915 std::upper_bound((*eProbaShellMap)[ionizationLevelIndex][(*k1)].begin(),
907 (*eProbaShellMap)[ionizationLevel 916 (*eProbaShellMap)[ionizationLevelIndex][(*k1)].end(),
908 random); 917 random);
909 918
910 auto prob11 = prob12 - 1; 919 auto prob11 = prob12 - 1;
911 920
912 auto prob22 = 921 auto prob22 =
913 std::upper_bound((*eProbaShellMap)[ion 922 std::upper_bound((*eProbaShellMap)[ionizationLevelIndex][(*k2)].begin(),
914 (*eProbaShellMap)[ionizationLevel 923 (*eProbaShellMap)[ionizationLevelIndex][(*k2)].end(),
915 random); 924 random);
916 925
917 auto prob21 = prob22 - 1; 926 auto prob21 = prob22 - 1;
918 927
919 valueK1 = *k1; 928 valueK1 = *k1;
920 valueK2 = *k2; 929 valueK2 = *k2;
921 valuePROB21 = *prob21; 930 valuePROB21 = *prob21;
922 valuePROB22 = *prob22; 931 valuePROB22 = *prob22;
923 valuePROB12 = *prob12; 932 valuePROB12 = *prob12;
924 valuePROB11 = *prob11; 933 valuePROB11 = *prob11;
925 934
926 // The following condition avoid getting 935 // The following condition avoid getting transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
927 if (valuePROB11 == 0) nrjTransf11 = bind 936 if (valuePROB11 == 0) nrjTransf11 = bindingEnergy;
928 else nrjTransf11 = (*eNrjTransfData)[ion 937 else nrjTransf11 = (*eNrjTransfData)[ionizationLevelIndex][valueK1][valuePROB11];
929 if (valuePROB12 == 1) 938 if (valuePROB12 == 1)
930 { 939 {
931 if ((valueK1 + bindingEnergy) / 2. > value 940 if ((valueK1 + bindingEnergy) / 2. > valueK1)
932 maximumEnergyTransfer1 = valueK1; 941 maximumEnergyTransfer1 = valueK1;
933 else 942 else
934 maximumEnergyTransfer1 = (valueK1 + bind 943 maximumEnergyTransfer1 = (valueK1 + bindingEnergy) / 2.;
935 944
936 nrjTransf12 = maximumEnergyTransfer1; 945 nrjTransf12 = maximumEnergyTransfer1;
937 } 946 }
938 else 947 else
939 nrjTransf12 = (*eNrjTransfData)[ioniza 948 nrjTransf12 = (*eNrjTransfData)[ionizationLevelIndex][valueK1][valuePROB12];
940 949
941 if (valuePROB21 == 0) nrjTransf21 = bind 950 if (valuePROB21 == 0) nrjTransf21 = bindingEnergy;
942 else nrjTransf21 = (*eNrjTransfData)[ion 951 else nrjTransf21 = (*eNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB21];
943 if (valuePROB22 == 1) 952 if (valuePROB22 == 1)
944 { 953 {
945 if ((valueK2 + bindingEnergy) / 2. > value 954 if ((valueK2 + bindingEnergy) / 2. > valueK2) maximumEnergyTransfer2 = valueK2;
946 else maximumEnergyTransfer2 = (valueK2 + b 955 else maximumEnergyTransfer2 = (valueK2 + bindingEnergy) / 2.;
947 956
948 nrjTransf22 = maximumEnergyTransfer2; 957 nrjTransf22 = maximumEnergyTransfer2;
949 } 958 }
950 else nrjTransf22 = (*eNrjTransfData)[ion 959 else nrjTransf22 = (*eNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB22];
951 960
952 } 961 }
953 // Avoids cases where cum xs is zero for k1 962 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
954 if (random > (*eProbaShellMap)[ionizationLev 963 if (random > (*eProbaShellMap)[ionizationLevelIndex][(*k1)].back())
955 { 964 {
956 auto prob22 = 965 auto prob22 =
957 std::upper_bound((*eProbaShellMap)[ion 966 std::upper_bound((*eProbaShellMap)[ionizationLevelIndex][(*k2)].begin(),
958 (*eProbaShellMap)[ionizationLevel 967 (*eProbaShellMap)[ionizationLevelIndex][(*k2)].end(),
959 random); 968 random);
960 auto prob21 = prob22 - 1; 969 auto prob21 = prob22 - 1;
961 970
962 valueK1 = *k1; 971 valueK1 = *k1;
963 valueK2 = *k2; 972 valueK2 = *k2;
964 valuePROB21 = *prob21; 973 valuePROB21 = *prob21;
965 valuePROB22 = *prob22; 974 valuePROB22 = *prob22;
966 975
967 nrjTransf21 = (*eNrjTransfData)[ionizati 976 nrjTransf21 = (*eNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB21];
968 nrjTransf22 = (*eNrjTransfData)[ionizati 977 nrjTransf22 = (*eNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB22];
969 978
970 G4double interpolatedvalue2 = Interpolat 979 G4double interpolatedvalue2 = Interpolate(valuePROB21,
971 valuePROB22, 980 valuePROB22,
972 random, 981 random,
973 nrjTransf21, 982 nrjTransf21,
974 nrjTransf22); 983 nrjTransf22);
975 984
976 // zeros are explicitly set 985 // zeros are explicitly set
977 G4double value = Interpolate(valueK1, va 986 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
978 987
979 return value; 988 return value;
980 } 989 }
981 } 990 }
982 } 991 }
983 else if (particleDefinition == G4Proton::Pro 992 else if (particleDefinition == G4Proton::ProtonDefinition())
984 { 993 {
985 // k should be in eV 994 // k should be in eV
986 dataDiffCSMap::iterator iterator_Nrj; 995 dataDiffCSMap::iterator iterator_Nrj;
987 iterator_Nrj = pNrjTransStorage.find(cur 996 iterator_Nrj = pNrjTransStorage.find(currentMaterial);
988 997
989 dataProbaShellMap::iterator iterator_Pro 998 dataProbaShellMap::iterator iterator_Proba;
990 iterator_Proba = pProbaShellStorage.find 999 iterator_Proba = pProbaShellStorage.find(currentMaterial);
991 1000
992 incidentEnergyMap::iterator iterator_Tdu 1001 incidentEnergyMap::iterator iterator_Tdummy;
993 iterator_Tdummy = pIncidentEnergyStorage 1002 iterator_Tdummy = pIncidentEnergyStorage.find(currentMaterial);
994 1003
995 if (iterator_Nrj == pNrjTransStorage.end() << 1004 if (iterator_Nrj == pNrjTransStorage.end() || iterator_Proba == pProbaShellStorage.end() ||
996 iterator_Tdummy == pIncidentEnergyStorag << 1005 iterator_Tdummy == pIncidentEnergyStorage.end())
997 { << 1006 {
998 G4String str = "Material "; << 1007 G4String str = "Material ";
999 str += currentMaterial + " not found!"; << 1008 str += currentMaterial + " not found!";
1000 G4Exception("G4MicroElecInelasticModel_ << 1009 G4Exception("G4MicroElecInelasticModel_new::TransferedEnergy", "em0002",
1001 FatalException, str); << 1010 FatalException, str);
1002 } << 1011 }
1003 else 1012 else
1004 { 1013 {
1005 vector<TriDimensionMap>* pNrjTransfData = 1014 vector<TriDimensionMap>* pNrjTransfData = iterator_Nrj->second; //Storage of possible transfer energies
1006 vector<VecMap>* pProbaShellMap = iterator 1015 vector<VecMap>* pProbaShellMap = iterator_Proba->second; //Storage of probabilities for energy transfer
1007 vector<G4double>* pTdummyVec = iterator_T 1016 vector<G4double>* pTdummyVec = iterator_Tdummy->second; //Incident energies for interpolation
1008 1017
1009 auto k2 = std::upper_bound(pTdummyVec->be 1018 auto k2 = std::upper_bound(pTdummyVec->begin(),
1010 pTdummyVec->end(), 1019 pTdummyVec->end(),
1011 k); 1020 k);
1012 1021
1013 auto k1 = k2 - 1; 1022 auto k1 = k2 - 1;
1014 1023
1015 // SI : the following condition avoids si 1024 // SI : the following condition avoids situations where random > last vector element,
1016 // for eg. when the last element is 1025 // for eg. when the last element is zero
1017 if (random <= (*pProbaShellMap)[ionizatio 1026 if (random <= (*pProbaShellMap)[ionizationLevelIndex][(*k1)].back()
1018 && random <= (*pProbaShellMap)[ioniza 1027 && random <= (*pProbaShellMap)[ionizationLevelIndex][(*k2)].back())
1019 { 1028 {
1020 auto prob12 = 1029 auto prob12 =
1021 std::upper_bound((*pProbaShellMap)[ioniza 1030 std::upper_bound((*pProbaShellMap)[ionizationLevelIndex][(*k1)].begin(),
1022 (*pProbaShellMap)[ionizationLevelInd 1031 (*pProbaShellMap)[ionizationLevelIndex][(*k1)].end(),
1023 random); 1032 random);
1024 auto prob11 = prob12 - 1; 1033 auto prob11 = prob12 - 1;
1025 auto prob22 = 1034 auto prob22 =
1026 std::upper_bound((*pProbaShellMap)[ioniza 1035 std::upper_bound((*pProbaShellMap)[ionizationLevelIndex][(*k2)].begin(),
1027 (*pProbaShellMap)[ionizationLevelInd 1036 (*pProbaShellMap)[ionizationLevelIndex][(*k2)].end(),
1028 random); 1037 random);
1029 auto prob21 = prob22 - 1; 1038 auto prob21 = prob22 - 1;
1030 1039
1031 valueK1 = *k1; 1040 valueK1 = *k1;
1032 valueK2 = *k2; 1041 valueK2 = *k2;
1033 valuePROB21 = *prob21; 1042 valuePROB21 = *prob21;
1034 valuePROB22 = *prob22; 1043 valuePROB22 = *prob22;
1035 valuePROB12 = *prob12; 1044 valuePROB12 = *prob12;
1036 valuePROB11 = *prob11; 1045 valuePROB11 = *prob11;
1037 1046
1038 // The following condition avoid gett 1047 // The following condition avoid getting transfered energy < binding energy
1039 // and forces cumxs = 1 for maximum e 1048 // and forces cumxs = 1 for maximum energy transfer.
1040 if (valuePROB11 == 0) nrjTransf11 = b 1049 if (valuePROB11 == 0) nrjTransf11 = bindingEnergy;
1041 else nrjTransf11 = (*pNrjTransfData)[ 1050 else nrjTransf11 = (*pNrjTransfData)[ionizationLevelIndex][valueK1][valuePROB11];
1042 1051
1043 if (valuePROB12 == 1) nrjTransf12 = m 1052 if (valuePROB12 == 1) nrjTransf12 = maximumEnergyTransferP;
1044 else nrjTransf12 = (*pNrjTransfData)[ 1053 else nrjTransf12 = (*pNrjTransfData)[ionizationLevelIndex][valueK1][valuePROB12];
1045 1054
1046 if (valuePROB21 == 0) nrjTransf21 = b 1055 if (valuePROB21 == 0) nrjTransf21 = bindingEnergy;
1047 else nrjTransf21 = (*pNrjTransfData)[ 1056 else nrjTransf21 = (*pNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB21];
1048 1057
1049 if (valuePROB22 == 1) nrjTransf22 = m 1058 if (valuePROB22 == 1) nrjTransf22 = maximumEnergyTransferP;
1050 else nrjTransf22 = (*pNrjTransfData)[ 1059 else nrjTransf22 = (*pNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB22];
1051 1060
1052 } 1061 }
1053 1062
1054 // Avoids cases where cum xs is zero for 1063 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1055 if (random > (*pProbaShellMap)[ionization 1064 if (random > (*pProbaShellMap)[ionizationLevelIndex][(*k1)].back())
1056 { 1065 {
1057 auto prob22 = 1066 auto prob22 =
1058 std::upper_bound((*pProbaShellMap)[ioniza 1067 std::upper_bound((*pProbaShellMap)[ionizationLevelIndex][(*k2)].begin(),
1059 (*pProbaShellMap)[ionizationLevelInd 1068 (*pProbaShellMap)[ionizationLevelIndex][(*k2)].end(),
1060 random); 1069 random);
1061 1070
1062 auto prob21 = prob22 - 1; 1071 auto prob21 = prob22 - 1;
1063 1072
1064 valueK1 = *k1; 1073 valueK1 = *k1;
1065 valueK2 = *k2; 1074 valueK2 = *k2;
1066 valuePROB21 = *prob21; 1075 valuePROB21 = *prob21;
1067 valuePROB22 = *prob22; 1076 valuePROB22 = *prob22;
1068 1077
1069 nrjTransf21 = (*pNrjTransfData)[ioniz 1078 nrjTransf21 = (*pNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB21];
1070 nrjTransf22 = (*pNrjTransfData)[ioniz 1079 nrjTransf22 = (*pNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB22];
1071 1080
1072 G4double interpolatedvalue2 = Interpo 1081 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1073 valuePROB22, 1082 valuePROB22,
1074 random, 1083 random,
1075 nrjTransf21, 1084 nrjTransf21,
1076 nrjTransf22); 1085 nrjTransf22);
1077 1086
1078 // zeros are explicitly set 1087 // zeros are explicitly set
1079 G4double value = Interpolate(valueK1, << 1088 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1080 return value; 1089 return value;
1081 } 1090 }
1082 } 1091 }
1083 } 1092 }
1084 // End electron and proton cases 1093 // End electron and proton cases
1085 1094
1086 G4double nrjTransfProduct = nrjTransf11 * n 1095 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
1087 1096
1088 if (nrjTransfProduct != 0.) 1097 if (nrjTransfProduct != 0.)
1089 { 1098 {
1090 nrj = QuadInterpolator(valuePROB11, 1099 nrj = QuadInterpolator(valuePROB11,
1091 valuePROB12, 1100 valuePROB12,
1092 valuePROB21, 1101 valuePROB21,
1093 valuePROB22, 1102 valuePROB22,
1094 nrjTransf11, 1103 nrjTransf11,
1095 nrjTransf12, 1104 nrjTransf12,
1096 nrjTransf21, 1105 nrjTransf21,
1097 nrjTransf22, 1106 nrjTransf22,
1098 valueK1, 1107 valueK1,
1099 valueK2, 1108 valueK2,
1100 k, 1109 k,
1101 random); 1110 random);
1102 } 1111 }
1103 1112
1104 return nrj; 1113 return nrj;
1105 } 1114 }
1106 1115
1107 //....oooOO0OOooo........oooOO0OOooo........o 1116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1108 1117
1109 G4double G4MicroElecInelasticModel_new::Diffe 1118 G4double G4MicroElecInelasticModel_new::DifferentialCrossSection(
1110 const G4ParticleDefinition * particl 1119 const G4ParticleDefinition * particleDefinition,
1111 G4double k, 1120 G4double k,
1112 G4double energyTransfer, 1121 G4double energyTransfer,
1113 G4int LevelIndex) 1122 G4int LevelIndex)
1114 { 1123 {
1115 G4double sigma = 0.; 1124 G4double sigma = 0.;
1116 1125
1117 if (energyTransfer >= currentMaterialStruct 1126 if (energyTransfer >= currentMaterialStructure->GetLimitEnergy(LevelIndex))
1118 { 1127 {
1119 G4double valueT1 = 0; 1128 G4double valueT1 = 0;
1120 G4double valueT2 = 0; 1129 G4double valueT2 = 0;
1121 G4double valueE21 = 0; 1130 G4double valueE21 = 0;
1122 G4double valueE22 = 0; 1131 G4double valueE22 = 0;
1123 G4double valueE12 = 0; 1132 G4double valueE12 = 0;
1124 G4double valueE11 = 0; 1133 G4double valueE11 = 0;
1125 1134
1126 G4double xs11 = 0; 1135 G4double xs11 = 0;
1127 G4double xs12 = 0; 1136 G4double xs12 = 0;
1128 G4double xs21 = 0; 1137 G4double xs21 = 0;
1129 G4double xs22 = 0; 1138 G4double xs22 = 0;
1130 1139
1131 if (particleDefinition == G4Electron::E 1140 if (particleDefinition == G4Electron::ElectronDefinition())
1132 { 1141 {
1133 1142
1134 dataDiffCSMap::iterator iterator_Proba; 1143 dataDiffCSMap::iterator iterator_Proba;
1135 iterator_Proba = eDiffDatatable.find(curr 1144 iterator_Proba = eDiffDatatable.find(currentMaterial);
1136 1145
1137 incidentEnergyMap::iterator iterator_Nrj; 1146 incidentEnergyMap::iterator iterator_Nrj;
1138 iterator_Nrj = eIncidentEnergyStorage.fin 1147 iterator_Nrj = eIncidentEnergyStorage.find(currentMaterial);
1139 1148
1140 TranfEnergyMap::iterator iterator_TransfN 1149 TranfEnergyMap::iterator iterator_TransfNrj;
1141 iterator_TransfNrj = eVecmStorage.find(cu 1150 iterator_TransfNrj = eVecmStorage.find(currentMaterial);
1142 1151
1143 if (iterator_Proba != eDiffDatatable.end( 1152 if (iterator_Proba != eDiffDatatable.end() && iterator_Nrj != eIncidentEnergyStorage.end()
1144 && iterator_TransfNrj!= eVecmStorage. 1153 && iterator_TransfNrj!= eVecmStorage.end())
1145 { 1154 {
1146 vector<TriDimensionMap>* eDiffCrossSe 1155 vector<TriDimensionMap>* eDiffCrossSectionData = (iterator_Proba->second);
1147 vector<G4double>* eTdummyVec = iterat 1156 vector<G4double>* eTdummyVec = iterator_Nrj->second; //Incident energies for interpolation
1148 VecMap* eVecm = iterator_TransfNrj->s 1157 VecMap* eVecm = iterator_TransfNrj->second;
1149 1158
1150 // k should be in eV and energy trans 1159 // k should be in eV and energy transfer eV also
1151 auto t2 = std::upper_bound(eTdummyVec 1160 auto t2 = std::upper_bound(eTdummyVec->begin(), eTdummyVec->end(), k);
1152 auto t1 = t2 - 1; 1161 auto t1 = t2 - 1;
1153 // SI : the following condition avoid 1162 // SI : the following condition avoids situations where energyTransfer >last vector element
1154 if (energyTransfer <= (*eVecm)[(*t1)] 1163 if (energyTransfer <= (*eVecm)[(*t1)].back() && energyTransfer <= (*eVecm)[(*t2)].back())
1155 { 1164 {
1156 auto e12 = std::upper_bound((*eVecm)[(* 1165 auto e12 = std::upper_bound((*eVecm)[(*t1)].begin(), (*eVecm)[(*t1)].end(), energyTransfer);
1157 auto e11 = e12 - 1; 1166 auto e11 = e12 - 1;
1158 auto e22 = std::upper_bound((*eVecm)[(* 1167 auto e22 = std::upper_bound((*eVecm)[(*t2)].begin(), (*eVecm)[(*t2)].end(), energyTransfer);
1159 auto e21 = e22 - 1; 1168 auto e21 = e22 - 1;
1160 1169
1161 valueT1 = *t1; 1170 valueT1 = *t1;
1162 valueT2 = *t2; 1171 valueT2 = *t2;
1163 valueE21 = *e21; 1172 valueE21 = *e21;
1164 valueE22 = *e22; 1173 valueE22 = *e22;
1165 valueE12 = *e12; 1174 valueE12 = *e12;
1166 valueE11 = *e11; 1175 valueE11 = *e11;
1167 1176
1168 xs11 = (*eDiffCrossSectionData)[LevelIn 1177 xs11 = (*eDiffCrossSectionData)[LevelIndex][valueT1][valueE11];
1169 xs12 = (*eDiffCrossSectionData)[LevelIn 1178 xs12 = (*eDiffCrossSectionData)[LevelIndex][valueT1][valueE12];
1170 xs21 = (*eDiffCrossSectionData)[LevelIn 1179 xs21 = (*eDiffCrossSectionData)[LevelIndex][valueT2][valueE21];
1171 xs22 = (*eDiffCrossSectionData)[LevelIn 1180 xs22 = (*eDiffCrossSectionData)[LevelIndex][valueT2][valueE22];
1172 } 1181 }
1173 } 1182 }
1174 else { 1183 else {
1175 G4String str = "Material "; 1184 G4String str = "Material ";
1176 str += currentMaterial + " not found!"; 1185 str += currentMaterial + " not found!";
1177 G4Exception("G4MicroElecDielectricModel 1186 G4Exception("G4MicroElecDielectricModels::DifferentialCrossSection", "em0002", FatalException, str);
1178 } 1187 }
1179 } 1188 }
1180 1189
1181 if (particleDefinition == G4Proton::Pro 1190 if (particleDefinition == G4Proton::ProtonDefinition())
1182 { 1191 {
1183 dataDiffCSMap::iterator iterator_Proba; 1192 dataDiffCSMap::iterator iterator_Proba;
1184 iterator_Proba = pDiffDatatable.find(curr 1193 iterator_Proba = pDiffDatatable.find(currentMaterial);
1185 1194
1186 incidentEnergyMap::iterator iterator_Nrj; 1195 incidentEnergyMap::iterator iterator_Nrj;
1187 iterator_Nrj = pIncidentEnergyStorage.fin 1196 iterator_Nrj = pIncidentEnergyStorage.find(currentMaterial);
1188 1197
1189 TranfEnergyMap::iterator iterator_TransfN 1198 TranfEnergyMap::iterator iterator_TransfNrj;
1190 iterator_TransfNrj = pVecmStorage.find(cu 1199 iterator_TransfNrj = pVecmStorage.find(currentMaterial);
1191 1200
1192 if (iterator_Proba != pDiffDatatable.end( 1201 if (iterator_Proba != pDiffDatatable.end() && iterator_Nrj != pIncidentEnergyStorage.end()
1193 && iterator_TransfNrj != pVecmStorage 1202 && iterator_TransfNrj != pVecmStorage.end())
1194 { 1203 {
1195 vector<TriDimensionMap>* pDiffCrossSe 1204 vector<TriDimensionMap>* pDiffCrossSectionData = (iterator_Proba->second);
1196 vector<G4double>* pTdummyVec = iterat 1205 vector<G4double>* pTdummyVec = iterator_Nrj->second; //Incident energies for interpolation
1197 VecMap* pVecm = iterator_TransfNrj->s 1206 VecMap* pVecm = iterator_TransfNrj->second;
1198 1207
1199 // k should be in eV and energy trans 1208 // k should be in eV and energy transfer eV also
1200 auto t2 = 1209 auto t2 =
1201 std::upper_bound(pTdummyVec->begin(), pTd 1210 std::upper_bound(pTdummyVec->begin(), pTdummyVec->end(), k);
1202 auto t1 = t2 - 1; 1211 auto t1 = t2 - 1;
1203 if (energyTransfer <= (*pVecm)[(*t1)] 1212 if (energyTransfer <= (*pVecm)[(*t1)].back() && energyTransfer <= (*pVecm)[(*t2)].back())
1204 { 1213 {
1205 auto e12 = std::upper_bound((*pVecm)[(* 1214 auto e12 = std::upper_bound((*pVecm)[(*t1)].begin(), (*pVecm)[(*t1)].end(), energyTransfer);
1206 auto e11 = e12 - 1; 1215 auto e11 = e12 - 1;
1207 auto e22 = std::upper_bound((*pVecm)[(* 1216 auto e22 = std::upper_bound((*pVecm)[(*t2)].begin(), (*pVecm)[(*t2)].end(), energyTransfer);
1208 auto e21 = e22 - 1; 1217 auto e21 = e22 - 1;
1209 1218
1210 valueT1 = *t1; 1219 valueT1 = *t1;
1211 valueT2 = *t2; 1220 valueT2 = *t2;
1212 valueE21 = *e21; 1221 valueE21 = *e21;
1213 valueE22 = *e22; 1222 valueE22 = *e22;
1214 valueE12 = *e12; 1223 valueE12 = *e12;
1215 valueE11 = *e11; 1224 valueE11 = *e11;
1216 1225
1217 xs11 = (*pDiffCrossSectionData)[LevelIn 1226 xs11 = (*pDiffCrossSectionData)[LevelIndex][valueT1][valueE11];
1218 xs12 = (*pDiffCrossSectionData)[LevelIn 1227 xs12 = (*pDiffCrossSectionData)[LevelIndex][valueT1][valueE12];
1219 xs21 = (*pDiffCrossSectionData)[LevelIn 1228 xs21 = (*pDiffCrossSectionData)[LevelIndex][valueT2][valueE21];
1220 xs22 = (*pDiffCrossSectionData)[LevelIn 1229 xs22 = (*pDiffCrossSectionData)[LevelIndex][valueT2][valueE22];
1221 } 1230 }
1222 } 1231 }
1223 else { 1232 else {
1224 G4String str = "Material "; 1233 G4String str = "Material ";
1225 str += currentMaterial + " not found!"; 1234 str += currentMaterial + " not found!";
1226 G4Exception("G4MicroElecDielectricModel 1235 G4Exception("G4MicroElecDielectricModels::DifferentialCrossSection", "em0002", FatalException, str);
1227 } 1236 }
1228 } 1237 }
1229 1238
1230 G4double xsProduct = xs11 * xs12 * xs21 1239 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
1231 if (xsProduct != 0.) 1240 if (xsProduct != 0.)
1232 { 1241 {
1233 sigma = QuadInterpolator( valueE11, v 1242 sigma = QuadInterpolator( valueE11, valueE12,
1234 valueE21, valueE22, 1243 valueE21, valueE22,
1235 xs11, xs12, 1244 xs11, xs12,
1236 xs21, xs22, 1245 xs21, xs22,
1237 valueT1, valueT2, 1246 valueT1, valueT2,
1238 k, energyTransfer); 1247 k, energyTransfer);
1239 } 1248 }
1240 1249
1241 } 1250 }
1242 1251
1243 return sigma; 1252 return sigma;
1244 } 1253 }
1245 1254
1246 //....oooOO0OOooo........oooOO0OOooo........o 1255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1247 1256
1248 1257
1249 G4double G4MicroElecInelasticModel_new::Inter 1258 G4double G4MicroElecInelasticModel_new::Interpolate(G4double e1,
1250 G4double e2, 1259 G4double e2,
1251 G4double e, 1260 G4double e,
1252 G4double xs1, 1261 G4double xs1,
1253 G4double xs2) 1262 G4double xs2)
1254 { 1263 {
1255 G4double value = 0.; 1264 G4double value = 0.;
1256 if (e == e1 || e1 == e2 ) { return xs1; } <<
1257 if (e == e2) { return xs2; } <<
1258 1265
1259 // Log-log interpolation by default 1266 // Log-log interpolation by default
1260 if (e1 > 0. && e2 > 0. && xs1 > 0. && xs2 > << 1267 if (e1 != 0 && e2 != 0 && (e2-e1) != 0 && !fasterCode)
1261 { 1268 {
1262 G4double a = std::log(xs2/xs1)/ std::lo 1269 G4double a = std::log(xs2/xs1)/ std::log(e2/e1);
1263 G4double b = std::log(xs2) - a * std::l 1270 G4double b = std::log(xs2) - a * std::log(e2);
1264 G4double sigma = a * std::log(e) + b; 1271 G4double sigma = a * std::log(e) + b;
1265 value = (std::exp(sigma)); 1272 value = (std::exp(sigma));
1266 } 1273 }
1267 1274
1268 // Switch to log-lin interpolation for fast 1275 // Switch to log-lin interpolation for faster code
1269 else if (xs1 > 0. && xs2 > 0. && fasterCode << 1276 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
1270 { 1277 {
1271 G4double d1 = std::log(xs1); 1278 G4double d1 = std::log(xs1);
1272 G4double d2 = std::log(xs2); 1279 G4double d2 = std::log(xs2);
1273 value = std::exp((d1 + (d2 - d1) * (e - 1280 value = std::exp((d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
1274 } 1281 }
1275 1282
1276 // Switch to lin-lin interpolation for fast 1283 // Switch to lin-lin interpolation for faster code
1277 // in case one of xs1 or xs2 (=cum proba) v 1284 // in case one of xs1 or xs2 (=cum proba) value is zero
1278 else << 1285 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
1279 { 1286 {
1280 G4double d1 = xs1; 1287 G4double d1 = xs1;
1281 G4double d2 = xs2; 1288 G4double d2 = xs2;
1282 value = (d1 + (d2 - d1) * (e - e1) / (e 1289 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
1283 } 1290 }
1284 1291
1285 return value; 1292 return value;
1286 } 1293 }
1287 1294
1288 //....oooOO0OOooo........oooOO0OOooo........o 1295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1289 1296
1290 G4double G4MicroElecInelasticModel_new::QuadI 1297 G4double G4MicroElecInelasticModel_new::QuadInterpolator(G4double e11, G4double e12,
1291 G4double e21, G4double e22, 1298 G4double e21, G4double e22,
1292 G4double xs11, G4double xs12, 1299 G4double xs11, G4double xs12,
1293 G4double xs21, G4double xs22, 1300 G4double xs21, G4double xs22,
1294 G4double t1, G4double t2, 1301 G4double t1, G4double t2,
1295 G4double t, G4double e) 1302 G4double t, G4double e)
1296 { 1303 {
1297 G4double interpolatedvalue1 = Interpolate(e 1304 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
1298 G4double interpolatedvalue2 = Interpolate(e 1305 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
1299 G4double value = Interpolate(t1, t2, t, int 1306 G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
1300 return value; 1307 return value;
1301 } 1308 }
1302 1309
1303 //....oooOO0OOooo........oooOO0OOooo........o 1310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1304 1311
1305 G4int G4MicroElecInelasticModel_new::RandomSe 1312 G4int G4MicroElecInelasticModel_new::RandomSelect(G4double k, const G4String& particle, G4double originalMass, G4int originalZ )
1306 { 1313 {
1307 G4int level = 0; 1314 G4int level = 0;
1308 1315
1309 TCSMap::iterator tablepos; 1316 TCSMap::iterator tablepos;
1310 tablepos = tableTCS.find(currentMaterial); 1317 tablepos = tableTCS.find(currentMaterial);
1311 if (tablepos == tableTCS.end()) { <<
1312 G4Exception("G4MicroElecInelasticModel_ne <<
1313 return level; <<
1314 } <<
1315 <<
1316 MapData* tableData = tablepos->second; 1318 MapData* tableData = tablepos->second;
1317 1319
1318 std::map< G4String,G4MicroElecCrossSectionD 1320 std::map< G4String,G4MicroElecCrossSectionDataSet_new*,std::less<G4String> >::iterator pos;
1319 pos = tableData->find(particle); 1321 pos = tableData->find(particle);
1320 1322
1321 std::vector<G4double> Zeff(currentMaterialS 1323 std::vector<G4double> Zeff(currentMaterialStructure->NumberOfLevels(), 1.0);
1322 if(originalMass>proton_mass_c2) { 1324 if(originalMass>proton_mass_c2) {
1323 for(G4int nl=0;nl<currentMaterialStructur 1325 for(G4int nl=0;nl<currentMaterialStructure->NumberOfLevels();nl++) {
1324 Zeff[nl] = BKZ(k/(proton_mass_c2/origin 1326 Zeff[nl] = BKZ(k/(proton_mass_c2/originalMass), originalMass/c_squared, originalZ, currentMaterialStructure->Energy(nl));
1325 } 1327 }
1326 } 1328 }
1327 1329
1328 if (pos != tableData->end()) 1330 if (pos != tableData->end())
1329 { 1331 {
1330 G4MicroElecCrossSectionDataSet_new* tab 1332 G4MicroElecCrossSectionDataSet_new* table = pos->second;
1331 1333
1332 if (table != 0) 1334 if (table != 0)
1333 { 1335 {
1334 G4double* valuesBuffer = new G4double[tab 1336 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1335 const G4int n = (G4int)table->NumberOfCom 1337 const G4int n = (G4int)table->NumberOfComponents();
1336 G4int i = (G4int)n; 1338 G4int i = (G4int)n;
1337 G4double value = 0.; 1339 G4double value = 0.;
1338 1340
1339 while (i>0) 1341 while (i>0)
1340 { 1342 {
1341 --i; 1343 --i;
1342 valuesBuffer[i] = table->GetComponent 1344 valuesBuffer[i] = table->GetComponent(i)->FindValue(k)*Zeff[i]*Zeff[i];
1343 value += valuesBuffer[i]; 1345 value += valuesBuffer[i];
1344 } 1346 }
1345 value *= G4UniformRand(); 1347 value *= G4UniformRand();
1346 1348
1347 i = n; 1349 i = n;
1348 1350
1349 while (i > 0) 1351 while (i > 0)
1350 { 1352 {
1351 --i; 1353 --i;
1352 1354
1353 if (valuesBuffer[i] > value) 1355 if (valuesBuffer[i] > value)
1354 { 1356 {
1355 delete[] valuesBuffer; 1357 delete[] valuesBuffer;
1356 return i; 1358 return i;
1357 } 1359 }
1358 value -= valuesBuffer[i]; 1360 value -= valuesBuffer[i];
1359 } 1361 }
1360 1362
1361 if (valuesBuffer) delete[] valuesBuffer; 1363 if (valuesBuffer) delete[] valuesBuffer;
1362 1364
1363 } 1365 }
1364 } 1366 }
1365 else 1367 else
1366 { 1368 {
1367 G4Exception("G4MicroElecInelasticModel_ 1369 G4Exception("G4MicroElecInelasticModel_new::RandomSelect","em0002",FatalException,"Model not applicable to particle type.");
1368 } 1370 }
1369 1371
1370 return level; 1372 return level;
1371 } 1373 }
1372 1374
1373 //....oooOO0OOooo........oooOO0OOooo........o 1375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1374 1376
1375 G4double G4MicroElecInelasticModel_new::Compu 1377 G4double G4MicroElecInelasticModel_new::ComputeRelativistVelocity(G4double E, G4double mass) {
1376 G4double x = E/mass; 1378 G4double x = E/mass;
1377 return c_light*std::sqrt(x*(x + 2.0))/(x + 1379 return c_light*std::sqrt(x*(x + 2.0))/(x + 1.0);
1378 } 1380 }
1379 1381
1380 //....oooOO0OOooo........oooOO0OOooo........o 1382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1381 1383
1382 G4double G4MicroElecInelasticModel_new::Compu 1384 G4double G4MicroElecInelasticModel_new::ComputeElasticQmax(G4double T1i, G4double T2i, G4double M1, G4double M2) {
1383 G4double v1i = ComputeRelativistVelocity(T1 1385 G4double v1i = ComputeRelativistVelocity(T1i, M1);
1384 G4double v2i = ComputeRelativistVelocity(T2 1386 G4double v2i = ComputeRelativistVelocity(T2i, M2);
1385 1387
1386 G4double v2f = 2*M1/(M1+M2)*v1i + (M2-M1)/( 1388 G4double v2f = 2*M1/(M1+M2)*v1i + (M2-M1)/(M1+M2)*-1*v2i;
1387 G4double vtransfer2a = v2f*v2f-v2i*v2i; 1389 G4double vtransfer2a = v2f*v2f-v2i*v2i;
1388 1390
1389 v2f = 2*M1/(M1+M2)*v1i + (M2-M1)/(M1+M2)*v2 1391 v2f = 2*M1/(M1+M2)*v1i + (M2-M1)/(M1+M2)*v2i;
1390 G4double vtransfer2b = v2f*v2f-v2i*v2i; 1392 G4double vtransfer2b = v2f*v2f-v2i*v2i;
1391 1393
1392 G4double vtransfer2 = std::max(vtransfer2a, 1394 G4double vtransfer2 = std::max(vtransfer2a, vtransfer2b);
1393 return 0.5*M2*vtransfer2; 1395 return 0.5*M2*vtransfer2;
1394 } 1396 }
1395 1397
1396 //....oooOO0OOooo........oooOO0OOooo........o 1398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1397 1399
1398 G4double G4MicroElecInelasticModel_new::stepF 1400 G4double G4MicroElecInelasticModel_new::stepFunc(G4double x) {
1399 return (x < 0.) ? 1.0 : 0.0; 1401 return (x < 0.) ? 1.0 : 0.0;
1400 } 1402 }
1401 1403
1402 //....oooOO0OOooo........oooOO0OOooo........o 1404 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1403 1405
1404 G4double G4MicroElecInelasticModel_new::vrkre 1406 G4double G4MicroElecInelasticModel_new::vrkreussler(G4double v, G4double vF)
1405 { 1407 {
1406 G4double r = vF*( std::pow(v/vF+1., 3.) - f 1408 G4double r = vF*( std::pow(v/vF+1., 3.) - fabs(std::pow(v/vF-1., 3.))
1407 + 4.*(v/vF)*(v/vF) ) + stepFunc(v/vF- 1409 + 4.*(v/vF)*(v/vF) ) + stepFunc(v/vF-1.) * (3./2.*v/vF -
1408 4.*(v/vF)*(v/vF) + 3.*std::po 1410 4.*(v/vF)*(v/vF) + 3.*std::pow(v/vF, 3.)
1409 - 0.5*std::pow(v/vF, 5.)); 1411 - 0.5*std::pow(v/vF, 5.));
1410 return r/(10.*v/vF); 1412 return r/(10.*v/vF);
1411 } 1413 }
1412 1414
1413 //....oooOO0OOooo........oooOO0OOooo........o 1415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1414 1416
1415 G4double G4MicroElecInelasticModel_new::BKZ(G 1417 G4double G4MicroElecInelasticModel_new::BKZ(G4double Ep, G4double mp, G4int Zp, G4double Eplasmon)
1416 { 1418 {
1417 // need atomic unit conversion 1419 // need atomic unit conversion
1418 G4double hbar = hbar_Planck, hbar2 = hbar*h 1420 G4double hbar = hbar_Planck, hbar2 = hbar*hbar, me = electron_mass_c2/c_squared, Ry = me*elm_coupling*elm_coupling/(2*hbar2);
1419 G4double hartree = 2*Ry, a0 = Bohr_radius, 1421 G4double hartree = 2*Ry, a0 = Bohr_radius, velocity = a0*hartree/hbar;
1420 G4double vp = ComputeRelativistVelocity(Ep, 1422 G4double vp = ComputeRelativistVelocity(Ep, mp);
1421 1423
1422 vp /= velocity; 1424 vp /= velocity;
1423 1425
1424 G4double wp = Eplasmon/hartree; 1426 G4double wp = Eplasmon/hartree;
1425 G4double a = std::pow(4./9./CLHEP::pi, 1./3 1427 G4double a = std::pow(4./9./CLHEP::pi, 1./3.);
1426 G4double vF = std::pow(wp*wp/(3.*a*a*a), 1. 1428 G4double vF = std::pow(wp*wp/(3.*a*a*a), 1./3.);
1427 G4double c = 0.9; 1429 G4double c = 0.9;
1428 G4double vr = vrkreussler(vp /*in u.a*/, vF 1430 G4double vr = vrkreussler(vp /*in u.a*/, vF /*in u.a*/);
1429 G4double yr = vr/std::pow(Zp, 2./3.); 1431 G4double yr = vr/std::pow(Zp, 2./3.);
1430 G4double q = 0.; 1432 G4double q = 0.;
1431 if(Zp==2) q = 1-exp(-c*vr/(Zp-5./16.)); 1433 if(Zp==2) q = 1-exp(-c*vr/(Zp-5./16.));
1432 else q = 1.-exp(-c*(yr-0.07)); 1434 else q = 1.-exp(-c*(yr-0.07));
1433 G4double Neq = Zp*(1.-q); 1435 G4double Neq = Zp*(1.-q);
1434 G4double l0 = 0.; 1436 G4double l0 = 0.;
1435 if(Neq<=2) l0 = 3./(Zp-0.3*(Neq-1))/2.; 1437 if(Neq<=2) l0 = 3./(Zp-0.3*(Neq-1))/2.;
1436 else l0 = 0.48*std::pow(Neq, 2./3.)/(Zp-Neq 1438 else l0 = 0.48*std::pow(Neq, 2./3.)/(Zp-Neq/7.);
1437 if(Zp==2) c = 1.0; 1439 if(Zp==2) c = 1.0;
1438 else c = 3./2.; 1440 else c = 3./2.;
1439 return Zp*(q + c*(1.-q)/vF/vF/2.0 * log(1.+ 1441 return Zp*(q + c*(1.-q)/vF/vF/2.0 * log(1.+std::pow(2.*l0*vF,2.)));
1440 } 1442 }
1441 1443
1442 //....oooOO0OOooo........oooOO0OOooo........o 1444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1443 1445