Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4MicroElecElasticModel_new.cc
Go to the documentation of this file.
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
85using namespace std;
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88
89G4MicroElecElasticModel_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}
122
123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124
125G4MicroElecElasticModel_new::~G4MicroElecElasticModel_new()
126{
127 // For total cross section
128 TCSMap::iterator pos2;
129 for (pos2 = tableTCS.begin(); pos2 != tableTCS.end(); ++pos2) {
130 MapData* tableData = pos2->second;
131 std::map< G4String, G4MicroElecCrossSectionDataSet_new*, std::less<G4String> >::iterator pos;
132 for (pos = tableData->begin(); pos != tableData->end(); ++pos)
133 {
134 G4MicroElecCrossSectionDataSet_new* table = pos->second;
135 delete table;
136 }
137 delete tableData;
138 }
139
140 //Clearing DCS maps
141
142 ThetaMap::iterator iterator_angle;
143 for (iterator_angle = thetaDataStorage.begin(); iterator_angle != thetaDataStorage.end(); ++iterator_angle) {
144 TriDimensionMap* eDiffCrossSectionData = iterator_angle->second;
145 eDiffCrossSectionData->clear();
146 delete eDiffCrossSectionData;
147 }
148
149 energyMap::iterator iterator_energy;
150 for (iterator_energy = eIncidentEnergyStorage.begin(); iterator_energy != eIncidentEnergyStorage.end(); ++iterator_energy) {
151 std::vector<G4double>* eTdummyVec = iterator_energy->second;
152 eTdummyVec->clear();
153 delete eTdummyVec;
154 }
155
156 ProbaMap::iterator iterator_proba;
157 for (iterator_proba = eProbaStorage.begin(); iterator_proba != eProbaStorage.end(); ++iterator_proba) {
158 VecMap* eVecm = iterator_proba->second;
159 eVecm->clear();
160 delete eVecm;
161 }
162
163}
164
165//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166
167void G4MicroElecElasticModel_new::Initialise(const G4ParticleDefinition* /*particle*/,
168 const G4DataVector& /*cuts*/)
169{
170 if (verboseLevel > -1)
171 G4cout << "Calling G4MicroElecElasticModel_new::Initialise()" << G4endl;
172 // Energy limits
173 // Reading of data files
174
175 G4double scaleFactor = 1e-18 * cm * cm;
176
177 G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
178 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
179
180 for (G4int i = 0; i < numOfCouples; ++i) {
181 const G4Material* material = theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
182
183 G4cout << "MicroElasticModel, Material " << i + 1 << " / " << numOfCouples << " : " << material->GetName() << G4endl;
184 if (material->GetName() == "Vacuum") continue;
185
186 G4String matName = material->GetName().substr(3, material->GetName().size());
187 G4cout<< matName<< G4endl;
188
189 currentMaterialStructure = new G4MicroElecMaterialStructure(matName);
190 lowEnergyLimitTable[matName]=currentMaterialStructure->GetElasticModelLowLimit();
191 highEnergyLimitTable[matName]=currentMaterialStructure->GetElasticModelHighLimit();
192 workFunctionTable[matName] = currentMaterialStructure->GetWorkFunction();
193
194 delete currentMaterialStructure;
195
196 G4cout << "Reading TCS file" << G4endl;
197 G4String fileElectron = "Elastic/elsepa_elastic_cross_e_" + matName;
198 G4cout << "Elastic Total Cross file : " << fileElectron << G4endl;
199
200 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
201 G4String electron = electronDef->GetParticleName();
202
203 // For total cross section
204 MapData* tableData = new MapData();
205
206 G4MicroElecCrossSectionDataSet_new* tableE = new G4MicroElecCrossSectionDataSet_new(new G4LogLogInterpolation, eV, scaleFactor);
207 tableE->LoadData(fileElectron);
208 tableData->insert(make_pair(electron, tableE));
209 tableTCS[matName] = tableData; //Storage of TCS
210
211 // For final state
212 const char* path = G4FindDataDir("G4LEDATA");
213 if (!path)
214 {
215 G4Exception("G4MicroElecElasticModel_new::Initialise","em0006",FatalException,"G4LEDATA environment variable not set.");
216 return;
217 }
218
219 //Reading DCS file
220 std::ostringstream eFullFileName;
221 eFullFileName << path << "/microelec/Elastic/elsepa_elastic_cumulated_diffcross_e_" + matName + ".dat";
222 G4cout << "Elastic Cumulated Diff Cross : " << eFullFileName.str().c_str() << G4endl;
223 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
224
225 if (!eDiffCrossSection)
226 G4Exception("G4MicroElecElasticModel_new::Initialise", "em0003", FatalException, "Missing data file: /microelec/sigmadiff_cumulated_elastic_e_Si.dat");
227
228 // October 21th, 2014 - Melanie Raine
229 // Added clear for MT
230 // Diff Cross Sections in cumulated mode
231 TriDimensionMap* eDiffCrossSectionData = new TriDimensionMap(); //Angles
232 std::vector<G4double>* eTdummyVec = new std::vector<G4double>; //Incident energy vector
233 VecMap* eProbVec = new VecMap; //Probabilities
234
235 eTdummyVec->push_back(0.);
236
237 while (!eDiffCrossSection.eof())
238 {
239 G4double tDummy; //incident energy
240 G4double eProb; //Proba
241 eDiffCrossSection >> tDummy >> eProb;
242
243 // SI : mandatory eVecm initialization
244 if (tDummy != eTdummyVec->back())
245 {
246 eTdummyVec->push_back(tDummy); //adding values for incident energy points
247 (*eProbVec)[tDummy].push_back(0.); //adding probability for the first angle, equal to 0
248 }
249
250 eDiffCrossSection >> (*eDiffCrossSectionData)[tDummy][eProb]; //adding Angle Value to map
251
252 if (eProb != (*eProbVec)[tDummy].back()) {
253 (*eProbVec)[tDummy].push_back(eProb); //Adding cumulated proba to map
254 }
255
256 }
257
258 //Filling maps for the material
259 thetaDataStorage[matName] = eDiffCrossSectionData;
260 eIncidentEnergyStorage[matName] = eTdummyVec;
261 eProbaStorage[matName] = eProbVec;
262 }
263 // End final state
264
265 if (verboseLevel > 2)
266 G4cout << "Loaded cross section files for MicroElec Elastic model" << G4endl;
267
268 if (verboseLevel > 0) {
269 G4cout << "MicroElec Elastic model is initialized " << G4endl
270 << "Energy range: "
271 << LowEnergyLimit() / eV << " eV - "
272 << HighEnergyLimit() / MeV << " MeV"
273 << G4endl; // system("pause"); linux doesn't like
274 }
275
276 if (isInitialised) { return; }
277 fParticleChangeForGamma = GetParticleChangeForGamma();
278 isInitialised = true;
279}
280
281//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282
283G4double G4MicroElecElasticModel_new::CrossSectionPerVolume(const G4Material* material,
284 const G4ParticleDefinition* p,
285 G4double ekin,
286 G4double,
287 G4double)
288{
289 if (verboseLevel > 3)
290 G4cout << "Calling CrossSectionPerVolume() of G4MicroElecElasticModel" << G4endl;
291
292 currentMaterialName = material->GetName().substr(3, material->GetName().size());
293 const G4DataVector cuts;
294 // Calculate total cross section for model
295 MapEnergy::iterator lowEPos;
296 lowEPos = lowEnergyLimitTable.find(currentMaterialName);
297
298 MapEnergy::iterator highEPos;
299 highEPos = highEnergyLimitTable.find(currentMaterialName);
300
301 MapEnergy::iterator killEPos;
302 killEPos = workFunctionTable.find(currentMaterialName);
303
304 if (lowEPos == lowEnergyLimitTable.end() || highEPos == highEnergyLimitTable.end() || killEPos == workFunctionTable.end())
305 {
306 G4String str = "Material ";
307 str += currentMaterialName + " not found!";
308 G4Exception("G4MicroElecElasticModel_new::EnergyLimits", "em0002", FatalException, str);
309 return 0;
310 }
311 else {
312 // G4cout << "normal elastic " << G4endl;
313 lowEnergyLimit = lowEPos->second;
314 highEnergyLimit = highEPos->second;
315 killBelowEnergy = killEPos->second;
316 }
317
318 if (ekin < killBelowEnergy) { return DBL_MAX; }
319
320 G4double sigma=0;
321
322 //Phonon for SiO2
323 if (currentMaterialName == "SILICON_DIOXIDE" && ekin < 100 * eV) {
324 acousticModelEnabled = true;
325
326 //Values for SiO2
327 G4double kbz = 11.54e9,
328 rho = 2.2 * 1000, // [g/cm3] * 1000
329 cs = 3560, //Sound speed
330 Ebz = 5.1 * 1.6e-19,
331 Aac = 17 * Ebz, //A screening parameter
332 Eac = 3.5 * 1.6e-19, //C deformation potential
333 prefactor = 2.2;// Facteur pour modifier les MFP
334
335 return AcousticCrossSectionPerVolume(ekin, kbz, rho, cs, Aac, Eac, prefactor);
336 }
337
338 else if (currentMaterialName == "ALUMINUM_OXIDE" && ekin < 20 * eV) {
339 acousticModelEnabled = true;
340
341 //Values for Al2O3
342 G4double kbz = 8871930614.247564,
343 rho = 3.97 * 1000, // [g/cm3] * 1000
344 cs = 233329.07733059773, //Sound speed
345 Aac = 2.9912494342262614e-19, //A screening parameter
346 Eac = 2.1622471654789847e-18, //C deformation potential
347 prefactor = 1;
348 return AcousticCrossSectionPerVolume(ekin, kbz, rho, cs, Aac, Eac, prefactor);
349 }
350 //Elastic
351 else {
352 acousticModelEnabled = false;
353
354 G4double density = material->GetTotNbOfAtomsPerVolume();
355 const G4String& particleName = p->GetParticleName();
356
357 TCSMap::iterator tablepos;
358 tablepos = tableTCS.find(currentMaterialName);
359
360 if (tablepos != tableTCS.end())
361 {
362 MapData* tableData = tablepos->second;
363
364 if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
365 {
366 std::map< G4String, G4MicroElecCrossSectionDataSet_new*, std::less<G4String> >::iterator pos;
367 pos = tableData->find(particleName);
368
369 if (pos != tableData->end()){
370 G4MicroElecCrossSectionDataSet_new* table = pos->second;
371 if (table != 0)
372 {
373 sigma = table->FindValue(ekin);
374 }
375 }
376 else
377 {
378 G4Exception("G4MicroElecElasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, "Model not applicable to particle type.");
379 }
380 }
381 else return 1 / DBL_MAX;
382 }
383 else {
384 G4String str = "Material ";
385 str += currentMaterialName + " TCS table not found!";
386 G4Exception("G4MicroElecElasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, str);
387 }
388
389 if (verboseLevel > 3) {
390 G4cout << "---> Kinetic energy(eV)=" << ekin / eV << G4endl;
391 G4cout << " - Cross section per Si atom (cm^2)=" << sigma / cm / cm << G4endl;
392 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density / (1. / cm) << G4endl;
393 }
394
395 // 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).
396 if (currentMaterialName == "BORON_NITRIDE") {
397 sigma = sigma * tanh(0.5 * pow(ekin / 5.2e-6, 2));
398 }
399 return sigma*density;
400 }
401}
402
403//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
404
405G4double G4MicroElecElasticModel_new::AcousticCrossSectionPerVolume(G4double ekin,
406 G4double kbz,
407 G4double rho,
408 G4double cs,
409 G4double Aac,
410 G4double Eac,
411 G4double prefactor)
412{
413
414 G4double e = 1.6e-19,
415 m0 = 9.10938356e-31,
416 h = 1.0546e-34,
417 kb = 1.38e-23;
418
419 G4double E = (ekin / eV) * e;
420 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);
421
422 // Parametres SiO2
423 G4double T = 300,
424 Ebz = (std::pow(h, 2) * std::pow(kbz, 2)) / (2 * m0),
425 hwbz = cs * kbz * h,
426 nbz = 1.0 / (exp(hwbz / (kb * T)) - 1),
427 Pac;
428
429 if (E < Ebz / 4.0)
430 {
431 Pac = ((pi * kb * T) / (h * std::pow(cs, 2) * rho)) * (std::pow(Eac, 2) * D) / (1 + (E / Aac));
432 }
433
434 else if (E > Ebz) //Screened relationship
435 {
436 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)));
437 }
438 else //Linear interpolation
439 {
440 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)));
441 G4double fEbz4 = ((pi * kb * T) / (h * std::pow(cs, 2) * rho)) * (std::pow(Eac, 2) * D) / (1 + ((Ebz / 4) / Aac));
442 G4double alpha = ((fEbz - fEbz4) / (Ebz - (Ebz / 4)));
443 Pac = alpha * E + (fEbz - alpha * Ebz);
444 }
445
446 G4double MFP = (std::sqrt(2 * E / m0) / (prefactor * Pac)) * m;
447
448 return 1 / MFP;
449}
450
451
452//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
453
454void G4MicroElecElasticModel_new::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
455 const G4MaterialCutsCouple* /*couple*/,
456 const G4DynamicParticle* aDynamicElectron,
457 G4double,
458 G4double)
459{
460 if (verboseLevel > 3)
461 G4cout << "Calling SampleSecondaries() of G4MicroElecElasticModel" << G4endl;
462
463 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
464
465 if (electronEnergy0 < killBelowEnergy)
466 {
467 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
468 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
469 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
470 return;
471 }
472
473 if (electronEnergy0 < highEnergyLimit)
474 {
475 G4double cosTheta = 0;
476 if (acousticModelEnabled)
477 {
478 cosTheta = 1 - 2 * G4UniformRand(); //Isotrope
479 }
480 else if (electronEnergy0 >= lowEnergyLimit)
481 {
482 cosTheta = RandomizeCosTheta(electronEnergy0);
483 }
484 G4double phi = 2. * pi * G4UniformRand();
485
486 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
487 G4ThreeVector xVers = zVers.orthogonal();
488 G4ThreeVector yVers = zVers.cross(xVers);
489
490 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
491 G4double yDir = xDir;
492 xDir *= std::cos(phi);
493 yDir *= std::sin(phi);
494
495 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
496
497 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
498 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
499 }
500}
501
502//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
503
504G4double G4MicroElecElasticModel_new::DamageEnergy(G4double T,G4double A, G4double Z)
505{
506 //.................. T in eV!!!!!!!!!!!!!
507 G4double Z2= Z;
508 G4double M2= A;
509 G4double k_d;
510 G4double epsilon_d;
511 G4double g_epsilon_d;
512 G4double E_nu;
513
514 k_d=0.1334*std::pow(Z2,(2./3.))*std::pow(M2,(-1./2.));
515 epsilon_d=0.01014*std::pow(Z2,(-7./3.))*(T/eV);
516 g_epsilon_d= epsilon_d+0.40244*std::pow(epsilon_d,(3./4.))+3.4008*std::pow(epsilon_d,(1./6.));
517
518 E_nu=1./(1.+ k_d*g_epsilon_d);
519
520 return E_nu;
521}
522
523//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
524
525G4double G4MicroElecElasticModel_new::Theta
526 (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)
527{
528
529 G4double theta = 0.;
530 G4double valueT1 = 0;
531 G4double valueT2 = 0;
532 G4double valueE21 = 0;
533 G4double valueE22 = 0;
534 G4double valueE12 = 0;
535 G4double valueE11 = 0;
536 G4double xs11 = 0;
537 G4double xs12 = 0;
538 G4double xs21 = 0;
539 G4double xs22 = 0;
540
541 if (particleDefinition == G4Electron::ElectronDefinition())
542 {
543 ThetaMap::iterator iterator_angle;
544 iterator_angle = thetaDataStorage.find(currentMaterialName);
545
546 energyMap::iterator iterator_energy;
547 iterator_energy = eIncidentEnergyStorage.find(currentMaterialName);
548
549 ProbaMap::iterator iterator_proba;
550 iterator_proba = eProbaStorage.find(currentMaterialName);
551
552 if (iterator_angle != thetaDataStorage.end() && iterator_energy != eIncidentEnergyStorage.end() && iterator_proba != eProbaStorage.end())
553 {
554 TriDimensionMap* eDiffCrossSectionData = iterator_angle->second; //Theta points
555 std::vector<G4double>* eTdummyVec = iterator_energy->second;
556 VecMap* eVecm = iterator_proba->second;
557
558 auto t2 = std::upper_bound(eTdummyVec->begin(), eTdummyVec->end(), k);
559 auto t1 = t2 - 1;
560 auto e12 = std::upper_bound((*eVecm)[(*t1)].begin(), (*eVecm)[(*t1)].end(), integrDiff);
561 auto e11 = e12 - 1;
562 auto e22 = std::upper_bound((*eVecm)[(*t2)].begin(), (*eVecm)[(*t2)].end(), integrDiff);
563 auto e21 = e22 - 1;
564
565 valueT1 = *t1;
566 valueT2 = *t2;
567 valueE21 = *e21;
568 valueE22 = *e22;
569 valueE12 = *e12;
570 valueE11 = *e11;
571
572 xs11 = (*eDiffCrossSectionData)[valueT1][valueE11];
573 xs12 = (*eDiffCrossSectionData)[valueT1][valueE12];
574 xs21 = (*eDiffCrossSectionData)[valueT2][valueE21];
575 xs22 = (*eDiffCrossSectionData)[valueT2][valueE22];
576 }
577 else
578 {
579 G4String str = "Material ";
580 str += currentMaterialName + " not found!";
581 G4Exception("G4MicroElecElasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, str);
582 }
583 }
584
585 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
586
587 theta = QuadInterpolator( valueE11, valueE12,
588 valueE21, valueE22,
589 xs11, xs12,
590 xs21, xs22,
591 valueT1, valueT2,
592 k, integrDiff );
593
594 return theta;
595}
596
597//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
598
599G4double G4MicroElecElasticModel_new::LinLogInterpolate(G4double e1,
600 G4double e2,
601 G4double e,
602 G4double xs1,
603 G4double xs2)
604{
605 G4double d1 = std::log(xs1);
606 G4double d2 = std::log(xs2);
607 G4double value = G4Exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
608 return value;
609}
610
611//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
612
613G4double G4MicroElecElasticModel_new::LinLinInterpolate(G4double e1,
614 G4double e2,
615 G4double e,
616 G4double xs1,
617 G4double xs2)
618{
619 G4double d1 = xs1;
620 G4double d2 = xs2;
621 G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
622 return value;
623}
624
625//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
626
627G4double G4MicroElecElasticModel_new::LogLogInterpolate(G4double e1,
628 G4double e2,
629 G4double e,
630 G4double xs1,
631 G4double xs2)
632{
633 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
634 G4double b = std::log10(xs2) - a*std::log10(e2);
635 G4double sigma = a*std::log10(e) + b;
636 G4double value = (std::pow(10.,sigma));
637 return value;
638}
639
640//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
641
642G4double G4MicroElecElasticModel_new::QuadInterpolator(G4double e11, G4double e12,
643 G4double e21, G4double e22,
644 G4double xs11, G4double xs12,
645 G4double xs21, G4double xs22,
646 G4double t1, G4double t2,
647 G4double t, G4double e)
648{
649
650
651 // Lin-Lin
652 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
653 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
654 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
655
656 return value;
657}
658
659//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
660
661G4double G4MicroElecElasticModel_new::RandomizeCosTheta(G4double k)
662{
663 G4double integrdiff=0;
664 G4double uniformRand=G4UniformRand();
665 integrdiff = uniformRand;
666
667 G4double theta=0.;
668 G4double cosTheta=0.;
669 theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
670
671 cosTheta= std::cos(theta*pi/180.);
672
673 return cosTheta;
674}
675
676//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
677
678
679void G4MicroElecElasticModel_new::SetKillBelowThreshold (G4double threshold)
680{
681 killBelowEnergy = threshold;
682
683 if (threshold < 5*CLHEP::eV)
684 {
685 G4Exception ("*** WARNING : the G4MicroElecElasticModel class is not validated below 5 eV !","",JustWarning,"") ;
686 threshold = 5*CLHEP::eV;
687 }
688}
689
690//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
D
G4double D(G4double temp)
Definition G4DNAElectronHoleRecombination.cc:93
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition G4FindDataDir.cc:94
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:132
G4MicroElecElasticModel_new.hh
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
fStopAndKill
@ fStopAndKill
Definition G4TrackStatus.hh:46
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
A
const G4double A[17]
Definition G4WaterStopping.cc:53
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::orthogonal
Hep3Vector orthogonal() const
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition G4Electron.cc:86
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition G4MaterialCutsCouple.hh:166
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition G4Material.hh:197
G4Material::GetName
const G4String & GetName() const
Definition G4Material.hh:171
G4MicroElecCrossSectionDataSet_new
Definition G4MicroElecCrossSectionDataSet_new.hh:37
G4MicroElecCrossSectionDataSet_new::LoadData
G4bool LoadData(const G4String &argFileName) override
Definition G4MicroElecCrossSectionDataSet_new.cc:59
G4MicroElecElasticModel_new::AcousticCrossSectionPerVolume
G4double AcousticCrossSectionPerVolume(G4double ekin, G4double kbz, G4double rho, G4double cs, G4double Aac, G4double Eac, G4double prefactor)
Definition G4MicroElecElasticModel_new.cc:405
G4MicroElecElasticModel_new::SetKillBelowThreshold
void SetKillBelowThreshold(G4double threshold)
Definition G4MicroElecElasticModel_new.cc:679
G4MicroElecElasticModel_new::~G4MicroElecElasticModel_new
~G4MicroElecElasticModel_new() override
Definition G4MicroElecElasticModel_new.cc:125
G4MicroElecElasticModel_new::G4MicroElecElasticModel_new
G4MicroElecElasticModel_new(const G4ParticleDefinition *p=0, const G4String &nam="MicroElecElasticModel")
Definition G4MicroElecElasticModel_new.cc:89
G4MicroElecElasticModel_new::SampleSecondaries
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
Definition G4MicroElecElasticModel_new.cc:454
G4MicroElecElasticModel_new::Initialise
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition G4MicroElecElasticModel_new.cc:167
G4MicroElecElasticModel_new::DamageEnergy
G4double DamageEnergy(G4double T, G4double A, G4double Z)
Definition G4MicroElecElasticModel_new.cc:504
G4MicroElecElasticModel_new::CrossSectionPerVolume
G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax) override
Definition G4MicroElecElasticModel_new.cc:283
G4MicroElecElasticModel_new::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition G4MicroElecElasticModel_new.hh:123
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition G4ProductionCutsTable.hh:256
G4ProductionCutsTable::GetTableSize
std::size_t GetTableSize() const
Definition G4ProductionCutsTable.hh:249
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition G4ProductionCutsTable.cc:58
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition G4VEmModel.hh:734
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition G4VEmModel.cc:121
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition G4VEmModel.hh:643
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition G4VEmModel.hh:636
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition G4VEmModel.hh:741
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition G4VEmModel.cc:67
DBL_MAX
#define DBL_MAX
Definition templates.hh:62