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