| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // G4MicroElecElasticModel_new.cc, 2011/08/29
28 // 2020/05/20 P. Caro
29 // Q. Gibaru is with CEA [a]
30 // M. Raine and D. Lambert a
31 //
32 // A part of this work has been funded by the
33 // [a] CEA, DAM, DIF - 91297 ARPAJON, France
34 // [b] ONERA - DPHY, 2 avenue E.Belin, 31055 T
35 // [c] CNES, 18 av.E.Belin, 31401 Toulouse CED
36 //
37 // Based on the following publications
38 // - A.Valentin, M. Raine,
39 // Inelastic cross-sections of low energy e
40 // for the simulation of heavy ion trac
41 // NSS Conf. Record 2010, pp. 80-85
42 // https://doi.org/10.1109/NSSMIC.
43 //
44 // - A.Valentin, M. Raine, M.Gaillardin,
45 // Geant4 physics processes for microdo
46 // very low energy electromagnetic mode
47 // https://doi.org/10.1016/j.nimb.
48 // NIM B, vol. 288, pp. 66-73, 2012, pa
49 // heavy ions in Si, NIM B, vol. 287, p
50 // https://doi.org/10.1016/j.nimb.
51 //
52 // - M. Raine, M. Gaillardin, P. Paillet
53 // Geant4 physics processes for silicon
54 // Improvements and extension of the en
55 // NIM B, vol. 325, pp. 97-100, 2014
56 // https://doi.org/10.1016/j.nimb.
57 //
58 // - J. Pierron, C. Inguimbert, M. Belhaj
59 // Electron emission yield for low ener
60 // Monte Carlo simulation and experimen
61 // Journal of Applied Physics 121 (2017
62 // https://doi.org/10.1063/1.498
63 //
64 // - P. Caron,
65 // Study of Electron-Induced Single-Eve
66 // PHD, 16th October 2019
67 //
68 // - Q.Gibaru, C.Inguimbert, P.Caron, M.Raine
69 // Geant4 physics processes for microdo
70 // Extension of MicroElec to very low e
71 // NIM B, 2020, in review.
72 //
73 //
74 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
84
85 using namespace std;
86
87 //....oooOO0OOooo........oooOO0OOooo........oo
88
89 G4MicroElecElasticModel_new::G4MicroElecElasti
90 const G4String& nam)
91 :G4VEmModel(nam), isInitialised(false)
92 {
93 killBelowEnergy = 0.1*eV; // Minimum e-
94 lowEnergyLimit = 0.1 * eV;
95 lowEnergyLimitOfModel = 10 * eV; // The mode
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 o
106 // 4 = entering in methods
107
108 if( verboseLevel>0 )
109 {
110 G4cout << "MicroElec Elastic model is co
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........oo
125
126 G4MicroElecElasticModel_new::~G4MicroElecElast
127 {
128 // For total cross section
129 TCSMap::iterator pos2;
130 for (pos2 = tableTCS.begin(); pos2 != tableT
131 MapData* tableData = pos2->second;
132 std::map< G4String, G4MicroElecCrossSectio
133 for (pos = tableData->begin(); pos != tabl
134 {
135 G4MicroElecCrossSectionDataSet_new* table =
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
145 TriDimensionMap* eDiffCrossSectionData = i
146 eDiffCrossSectionData->clear();
147 delete eDiffCrossSectionData;
148 }
149
150 energyMap::iterator iterator_energy;
151 for (iterator_energy = eIncidentEnergyStorag
152 std::vector<G4double>* eTdummyVec = iterat
153 eTdummyVec->clear();
154 delete eTdummyVec;
155 }
156
157 ProbaMap::iterator iterator_proba;
158 for (iterator_proba = eProbaStorage.begin();
159 VecMap* eVecm = iterator_proba->second;
160 eVecm->clear();
161 delete eVecm;
162 }
163
164 }
165
166 //....oooOO0OOooo........oooOO0OOooo........oo
167
168 void G4MicroElecElasticModel_new::Initialise(c
169 const G4DataVector& /*cuts*/)
170 {
171 if (isOkToBeInitialised == true && isInitial
172
173 if (verboseLevel > -1)
174 G4cout << "Calling G4MicroElecElasticMod
175 // Energy limits
176 // Reading of data files
177
178 G4double scaleFactor = 1e-18 * cm * cm;
179
180 G4ProductionCutsTable* theCoupleTable =
181 G4ProductionCutsTable::GetProductionCuts
182 G4int numOfCouples = (G4int)theCoupleTable
183
184 for (G4int i = 0; i < numOfCouples; ++i) {
185 const G4Material* material =
186 theCoupleTable->GetMaterialCutsCouple(i)->Ge
187
188 //theCoupleTable->GetMaterialCutsCouple(
189
190 G4cout << "MicroElasticModel, Material "
191 if (material->GetName() == "Vacuum") con
192
193 G4String matName = material->GetName().s
194 G4cout<< matName<< G4endl;
195
196 currentMaterialStructure = new G4MicroEl
197 lowEnergyLimitTable[matName]=currentMate
198 highEnergyLimitTable[matName]=currentMat
199 workFunctionTable[matName] = currentMate
200
201 delete currentMaterialStructure;
202
203 G4cout << "Reading TCS file" << G4endl;
204 G4String fileElectron = "Elastic/elsepa_
205 G4cout << "Elastic Total Cross file : "
206
207 G4ParticleDefinition* electronDef = G4El
208 G4String electron = electronDef->GetPart
209
210 // For total cross section
211 MapData* tableData = new MapData();
212
213 G4MicroElecCrossSectionDataSet_new* tabl
214 tableE->LoadData(fileElectron);
215 tableData->insert(make_pair(electron, ta
216 tableTCS[matName] = tableData; //Storage
217
218 // For final state
219 const char* path = G4FindDataDir("G4LEDA
220 if (!path)
221 {
222 G4Exception("G4MicroElecElasticModel_new::
223 return;
224 }
225
226 //Reading DCS file
227 std::ostringstream eFullFileName;
228 eFullFileName << path << "/microelec/Ela
229 G4cout << "Elastic Cumulated Diff Cross
230 std::ifstream eDiffCrossSection(eFullFil
231
232 if (!eDiffCrossSection)
233 G4Exception("G4MicroElecElasticModel_new::In
234
235 // October 21th, 2014 - Melanie Raine
236 // Added clear for MT
237 // Diff Cross Sections in cumulated mode
238 TriDimensionMap* eDiffCrossSectionData =
239 std::vector<G4double>* eTdummyVec = new
240 VecMap* eProbVec = new VecMap; //Probabi
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); //addin
254 (*eProbVec)[tDummy].push_back(0.); //a
255 }
256
257 eDiffCrossSection >> (*eDiffCrossSectionDa
258
259 if (eProb != (*eProbVec)[tDummy].back()) {
260 (*eProbVec)[tDummy].push_back(eProb); //
261 }
262
263 }
264
265 //Filling maps for the material
266 thetaDataStorage[matName] = eDiffCrossSe
267 eIncidentEnergyStorage[matName] = eTdumm
268 eProbaStorage[matName] = eProbVec;
269 }
270 // End final state
271
272 if (verboseLevel > 2)
273 G4cout << "Loaded cross section files fo
274
275 if (verboseLevel > 0)
276 {
277 G4cout << "MicroElec Elastic model is initia
278 << "Energy range: "
279 << LowEnergyLimit() / eV << " eV - "
280 << HighEnergyLimit() / MeV << " MeV"
281 << G4endl; // system("pause"); linux
282 }
283
284 if (isInitialised) { return; }
285 fParticleChangeForGamma = GetParticleChang
286 isInitialised = true;
287 }
288 }
289
290 //....oooOO0OOooo........oooOO0OOooo........oo
291
292 G4double G4MicroElecElasticModel_new::CrossSec
293 const G4ParticleDefinition* p,
294 G4double ekin,
295 G4double,
296 G4double)
297 {
298 if (verboseLevel > 3)
299 G4cout << "Calling CrossSectionPerVolume()
300
301 isOkToBeInitialised = true;
302 currentMaterialName = material->GetName().su
303 const G4DataVector cuts;
304 Initialise(p, cuts);
305 // Calculate total cross section for model
306 MapEnergy::iterator lowEPos;
307 lowEPos = lowEnergyLimitTable.find(currentMa
308
309 MapEnergy::iterator highEPos;
310 highEPos = highEnergyLimitTable.find(current
311
312 MapEnergy::iterator killEPos;
313 killEPos = workFunctionTable.find(currentMat
314
315 if (lowEPos == lowEnergyLimitTable.end() ||
316 {
317 G4String str = "Material ";
318 str += currentMaterialName + " not found
319 G4Exception("G4MicroElecElasticModel_new
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"
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 pot
347 prefactor = 2.2;// Facteur pour modifier
348
349 return AcousticCrossSectionPerVolume(ekin,
350 }
351
352 else if (currentMaterialName == "ALUMINUM_OXID
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 screeni
360 Eac = 2.1622471654789847e-18, //C deforma
361 prefactor = 1;
362 return AcousticCrossSectionPerVolume(ekin,
363 }
364 //Elastic
365 else {
366 acousticModelEnabled = false;
367
368 G4double density = material->GetTotNbOfAto
369 const G4String& particleName = p->GetParti
370
371 TCSMap::iterator tablepos;
372 tablepos = tableTCS.find(currentMaterialNa
373
374 if (tablepos != tableTCS.end())
375 {
376 MapData* tableData = tablepos->second;
377
378 if (ekin >= lowEnergyLimit && ekin < highEne
379 {
380 std::map< G4String, G4MicroElecCrossSect
381 pos = tableData->find(particleName);
382
383 if (pos != tableData->end())
384 {
385 G4MicroElecCrossSectionDataSet_new* table
386 if (table != 0)
387 {
388 sigma = table->FindValue(ekin);
389 }
390 }
391 else
392 {
393 G4Exception("G4MicroElecElasticModel_new::
394 }
395 }
396 else return 1 / DBL_MAX;
397 }
398 else
399 {
400 G4String str = "Material ";
401 str += currentMaterialName + " TCS table not
402 G4Exception("G4MicroElecElasticModel_new::Co
403 }
404
405 if (verboseLevel > 3)
406 {
407 G4cout << "---> Kinetic energy(eV)=" << ekin
408 G4cout << " - Cross section per Si atom (cm^
409 G4cout << " - Cross section per Si atom (cm^
410 }
411
412 // Hsing-YinChangaAndrewAlvaradoaTreyWeber
413 if (currentMaterialName == "BORON_NITRIDE"
414 sigma = sigma * tanh(0.5 * pow(ekin /
415 }
416 return sigma*density;
417 }
418 }
419
420 //....oooOO0OOooo........oooOO0OOooo........oo
421
422 G4double G4MicroElecElasticModel_new::Acoustic
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(p
438
439 // Parametres SiO2
440 G4double T = 300,
441 Ebz = (std::pow(h, 2) * std::pow(kbz, 2))
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,
449 }
450
451 else if (E > Ebz) //Screened relationship
452 {
453 Pac = ((2 * pi * m0 * (2 * nbz + 1)) / (
454 }
455 else //Linear interpolation
456 {
457 G4double fEbz = ((2 * pi * m0 * (2 * nbz
458 G4double fEbz4 = ((pi * kb * T) / (h * s
459 G4double alpha = ((fEbz - fEbz4) / (Ebz
460 Pac = alpha * E + (fEbz - alpha * Ebz);
461 }
462
463 G4double MFP = (std::sqrt(2 * E / m0) / (pre
464
465 return 1 / MFP;
466 }
467
468
469 //....oooOO0OOooo........oooOO0OOooo........oo
470
471 void G4MicroElecElasticModel_new::SampleSecond
472 const G4MaterialCutsCouple* /*coup
473 const G4DynamicParticle* aDynamicE
474 G4double,
475 G4double)
476 {
477
478 if (verboseLevel > 3)
479 G4cout << "Calling SampleSecondaries() of
480
481 G4double electronEnergy0 = aDynamicElectron-
482
483 if (electronEnergy0 < killBelowEnergy)
484 {
485 fParticleChangeForGamma->SetProposedKinet
486 fParticleChangeForGamma->ProposeTrackStat
487 fParticleChangeForGamma->ProposeLocalEner
488 return;
489 }
490
491 if (electronEnergy0 < highEnergyLimit)
492 {
493 G4double cosTheta = 0;
494 if (acousticModelEnabled)
495 {
496 cosTheta = 1 - 2 * G4UniformRand(); //Isotr
497 }
498 else if (electronEnergy0 >= lowEnergyLimi
499 {
500 cosTheta = RandomizeCosTheta(electronEnergy
501 }
502
503 G4double phi = 2. * pi * G4UniformRand();
504
505 G4ThreeVector zVers = aDynamicElectron->G
506 G4ThreeVector xVers = zVers.orthogonal();
507 G4ThreeVector yVers = zVers.cross(xVers);
508
509 G4double xDir = std::sqrt(1. - cosTheta*c
510 G4double yDir = xDir;
511 xDir *= std::cos(phi);
512 yDir *= std::sin(phi);
513
514 G4ThreeVector zPrimeVers((xDir*xVers + yD
515
516 fParticleChangeForGamma->ProposeMomentumD
517 fParticleChangeForGamma->SetProposedKinet
518 }
519 }
520
521 //....oooOO0OOooo........oooOO0OOooo........oo
522
523 G4double G4MicroElecElasticModel_new::DamageE
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,
534 epsilon_d=0.01014*std::pow(Z2,(-7./3.))*(T/e
535 g_epsilon_d= epsilon_d+0.40244*std::pow(epsi
536
537 E_nu=1./(1.+ k_d*g_epsilon_d);
538
539 return E_nu;
540 }
541
542 //....oooOO0OOooo........oooOO0OOooo........oo
543
544 G4double G4MicroElecElasticModel_new::Theta
545 (G4ParticleDefinition * particleDefinition,
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::Electr
561 {
562 ThetaMap::iterator iterator_angle;
563 iterator_angle = thetaDataStorage.find(cur
564
565 energyMap::iterator iterator_energy;
566 iterator_energy = eIncidentEnergyStorage.f
567
568 ProbaMap::iterator iterator_proba;
569 iterator_proba = eProbaStorage.find(curren
570
571 if (iterator_angle != thetaDataStorage.end
572 {
573 TriDimensionMap* eDiffCrossSectionData = ite
574 std::vector<G4double>* eTdummyVec = iterator
575 VecMap* eVecm = iterator_proba->second;
576
577 auto t2 = std::upper_bound(eTdummyVec->begin
578 auto t1 = t2 - 1;
579 auto e12 = std::upper_bound((*eVecm)[(
580 auto e11 = e12 - 1;
581 auto e22 = std::upper_bound((*eVecm)[(*t2)].
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][val
592 xs12 = (*eDiffCrossSectionData)[valueT1][val
593 xs21 = (*eDiffCrossSectionData)[valueT2][val
594 xs22 = (*eDiffCrossSectionData)[valueT2][val
595 }
596 else
597 {
598 G4String str = "Material ";
599 str += currentMaterialName + " not found!";
600 G4Exception("G4MicroElecElasticModel_new::Co
601 }
602
603 }
604
605 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)
606
607 theta = QuadInterpolator( valueE11, valueE1
608 valueE21, valueE22,
609 xs11, xs12,
610 xs21, xs22,
611 valueT1, valueT2,
612 k, integrDiff );
613
614 return theta;
615 }
616
617 //....oooOO0OOooo........oooOO0OOooo........oo
618
619 G4double G4MicroElecElasticModel_new::LinLogIn
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 - e
628 return value;
629 }
630
631 //....oooOO0OOooo........oooOO0OOooo........oo
632
633 G4double G4MicroElecElasticModel_new::LinLinIn
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)/ (
642 return value;
643 }
644
645 //....oooOO0OOooo........oooOO0OOooo........oo
646
647 G4double G4MicroElecElasticModel_new::LogLogIn
648 G4double e2,
649 G4double e,
650 G4double xs1,
651 G4double xs2)
652 {
653 G4double a = (std::log10(xs2)-std::log10(xs1
654 G4double b = std::log10(xs2) - a*std::log10(
655 G4double sigma = a*std::log10(e) + b;
656 G4double value = (std::pow(10.,sigma));
657 return value;
658 }
659
660 //....oooOO0OOooo........oooOO0OOooo........oo
661
662 G4double G4MicroElecElasticModel_new::QuadInte
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 = LinLinInterpol
673 G4double interpolatedvalue2 = LinLinInterpol
674 G4double value = LinLinInterpolate(t1, t2, t
675
676 return value;
677 }
678
679 //....oooOO0OOooo........oooOO0OOooo........oo
680
681 G4double G4MicroElecElasticModel_new::Randomiz
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
690
691 cosTheta= std::cos(theta*pi/180.);
692
693 return cosTheta;
694 }
695
696 //....oooOO0OOooo........oooOO0OOooo........oo
697
698
699 void G4MicroElecElasticModel_new::SetKillBelow
700 {
701 killBelowEnergy = threshold;
702
703 if (threshold < 5*CLHEP::eV)
704 {
705 G4Exception ("*** WARNING : the G4MicroE
706 threshold = 5*CLHEP::eV;
707 }
708 }
709
710 //....oooOO0OOooo........oooOO0OOooo........oo
711