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