| 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 // G4MicroElecInelasticModel_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
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........oooOO0OOooo........oooOO0OOooo....
92
93 using namespace std;
94
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96
97 G4MicroElecInelasticModel_new::G4MicroElecInelasticModel_new(
98 const G4ParticleDefinition*, const G4String& nam)
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 openings, sampling of atoms
108 // 4 = entering in methods
109
110 if( verboseLevel>0 )
111 {
112 G4cout << "MicroElec inelastic model is constructed " << G4endl;
113 }
114
115 //Mark this model as "applicable" for atomic deexcitation
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........oooOO0OOooo........oooOO0OOooo....
129
130 G4MicroElecInelasticModel_new::~G4MicroElecInelasticModel_new()
131 {
132 // Cross section
133 // (0)
134 for (auto const& p : tableTCS) {
135 MapData* tableData = p.second;
136 for (auto const& pos : *tableData) { delete pos.second; }
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 : tableMaterialsStructures) {
191 delete p.second;
192 }
193 }
194
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196
197 void G4MicroElecInelasticModel_new::Initialise(const G4ParticleDefinition* particle,
198 const G4DataVector& /*cuts*/)
199 {
200 if (isInitialised) { return; }
201
202 if (verboseLevel > 3)
203 G4cout << "Calling G4MicroElecInelasticModel_new::Initialise()" << G4endl;
204
205 const char* path = G4FindDataDir("G4LEDATA");
206 if (!path)
207 {
208 G4Exception("G4MicroElecElasticModel_new::Initialise","em0006",FatalException,"G4LEDATA environment variable not set.");
209 return;
210 }
211
212 G4String modelName = "mermin";
213 G4cout << "****************************" << G4endl;
214 G4cout << modelName << " model loaded !" << G4endl;
215
216 // Energy limits
217 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
218 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
219 G4String electron = electronDef->GetParticleName();
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 = G4ProductionCutsTable::GetProductionCutsTable();
230 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
231
232 for (G4int i = 0; i < numOfCouples; ++i) {
233 const G4Material* material = theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
234 G4cout << "Material " << i + 1 << " / " << numOfCouples << " : " << material->GetName() << G4endl;
235 if (material->GetName() == "Vacuum") continue;
236 G4String mat = material->GetName().substr(3, material->GetName().size());
237 MapData* tableData = new MapData;
238 currentMaterialStructure = new G4MicroElecMaterialStructure(mat);
239
240 tableMaterialsStructures[mat] = currentMaterialStructure;
241 if (particle == electronDef) {
242 //TCS
243 G4String fileElectron("Inelastic/" + modelName + "_sigma_inelastic_e-_" + mat);
244 G4cout << fileElectron << G4endl;
245 G4MicroElecCrossSectionDataSet_new* tableE = new G4MicroElecCrossSectionDataSet_new(new G4LogLogInterpolation, MeV, scaleFactor);
246 tableE->LoadData(fileElectron);
247 tableData->insert(make_pair(electron, tableE));
248
249 // DCS
250 std::ostringstream eFullFileName;
251 if (fasterCode) {
252 eFullFileName << path << "/microelec/Inelastic/cumulated_" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat";
253 G4cout << "Faster code = true" << G4endl;
254 G4cout << "Inelastic/cumulated_" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat" << G4endl;
255 }
256 else {
257 eFullFileName << path << "/microelec/Inelastic/" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat";
258 G4cout << "Faster code = false" << G4endl;
259 G4cout << "Inelastic/" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat" << G4endl;
260 }
261
262 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
263 if (!eDiffCrossSection)
264 {
265 std::stringstream ss;
266 ss << "Missing data " << eFullFileName.str().c_str();
267 std::string sortieString = ss.str();
268
269 if (fasterCode) G4Exception("G4MicroElecInelasticModel_new::Initialise", "em0003",
270 FatalException, sortieString.c_str());
271 else {
272 G4Exception("G4MicroElecInelasticModel_new::Initialise", "em0003",
273 FatalException, "Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
274 }
275 }
276
277 // Clear the arrays for re-initialization case (MT mode)
278 // Octobre 22nd, 2014 - Melanie Raine
279 //Creating vectors of maps for DCS and Cumulated DCS for the current material.
280 //Each vector is storing one map for each shell.
281 G4bool isUsed1 = false;
282 vector<TriDimensionMap>* eDiffCrossSectionData =
283 new vector<TriDimensionMap>; //Storage of [IncidentEnergy, TransfEnergy, DCS values], used in slower code
284 vector<TriDimensionMap>* eNrjTransfData =
285 new vector<TriDimensionMap>; //Storage of possible transfer energies by shell
286 vector<VecMap>* eProbaShellMap = new vector<VecMap>; //Storage of the vectors containing all cumulated DCS values for an initial energy, by shell
287 vector<G4double>* eTdummyVec = new vector<G4double>; //Storage of incident energies for interpolation
288 VecMap* eVecm = new VecMap; //Transfered energy map for slower code
289
290 for (G4int j = 0; j < currentMaterialStructure->NumberOfLevels(); ++j) //Filling the map vectors with an empty map for each shell
291 {
292 eDiffCrossSectionData->push_back(TriDimensionMap());
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()) eTdummyVec->push_back(tDummy);
304
305 G4double tmp; //probability
306 for (G4int j = 0; j < currentMaterialStructure->NumberOfLevels(); ++j)
307 {
308 eDiffCrossSection >> tmp;
309 (*eDiffCrossSectionData)[j][tDummy][eDummy] = tmp;
310
311 if (fasterCode)
312 {
313 (*eNrjTransfData)[j][tDummy][(*eDiffCrossSectionData)[j][tDummy][eDummy]] = eDummy;
314 (*eProbaShellMap)[j][tDummy].push_back((*eDiffCrossSectionData)[j][tDummy][eDummy]);
315 }
316 else { // SI - only if eof is not reached !
317 (*eDiffCrossSectionData)[j][tDummy][eDummy] *= scaleFactor;
318 (*eVecm)[tDummy].push_back(eDummy);
319 }
320 }
321 }
322 //
323 //G4cout << "add to material vector" << G4endl;
324
325 //Filing maps for the current material into the master maps
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 + "_sigma_inelastic_p_" + mat); G4cout << fileProton << G4endl;
352 G4MicroElecCrossSectionDataSet_new* tableP = new G4MicroElecCrossSectionDataSet_new(new G4LogLogInterpolation, MeV, scaleFactor);
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/Inelastic/cumulated_" + modelName + "_sigmadiff_inelastic_p_" + mat + ".dat";
360 G4cout << "Faster code = true" << G4endl;
361 G4cout << "Inelastic/cumulated_" + modelName + "_sigmadiff_inelastic_p_" + mat + ".dat" << G4endl;
362 }
363 else {
364 pFullFileName << path << "/microelec/Inelastic/" + modelName + "_sigmadiff_inelastic_p_" + mat + ".dat";
365 G4cout << "Faster code = false" << G4endl;
366 G4cout << "Inelastic/" + modelName + "_sigmadiff_inelastic_e-_" + mat + ".dat" << G4endl;
367 }
368
369 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
370 if (!pDiffCrossSection)
371 {
372 if (fasterCode) G4Exception("G4MicroElecInelasticModel_new::Initialise", "em0003",
373 FatalException, "Missing data file:/microelec/sigmadiff_cumulated_inelastic_p_Si.dat");
374 else {
375 G4Exception("G4MicroElecInelasticModel_new::Initialise", "em0003",
376 FatalException, "Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
377 }
378 }
379
380 //
381 // Clear the arrays for re-initialization case (MT mode)
382 // Octobre 22nd, 2014 - Melanie Raine
383 //Creating vectors of maps for DCS and Cumulated DCS for the current material.
384 //Each vector is storing one map for each shell.
385
386 G4bool isUsed1 = false;
387 vector<TriDimensionMap>* pDiffCrossSectionData =
388 new vector<TriDimensionMap>; //Storage of [IncidentEnergy, TransfEnergy, DCS values], used in slower code
389 vector<TriDimensionMap>* pNrjTransfData =
390 new vector<TriDimensionMap>; //Storage of possible transfer energies by shell
391 vector<VecMap>* pProbaShellMap =
392 new vector<VecMap>; //Storage of the vectors containing all cumulated DCS values for an initial energy, by shell
393 vector<G4double>* pTdummyVec =
394 new vector<G4double>; //Storage of incident energies for interpolation
395 VecMap* pVecm = new VecMap; //Transfered energy map for slower code
396
397 for (G4int j = 0; j < currentMaterialStructure->NumberOfLevels(); ++j)
398 //Filling the map vectors with an empty map for each shell
399 {
400 pDiffCrossSectionData->push_back(TriDimensionMap());
401 pNrjTransfData->push_back(TriDimensionMap());
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()) pTdummyVec->push_back(tDummy);
412
413 G4double tmp; //probability
414 for (G4int j = 0; j < currentMaterialStructure->NumberOfLevels(); j++)
415 {
416 pDiffCrossSection >> tmp;
417 (*pDiffCrossSectionData)[j][tDummy][eDummy] = tmp;
418 // ArrayofMaps[j] -> fill with 3DMap(incidentEnergy,
419 // 2Dmap (transferedEnergy,proba=tmp) ) -> fill map for shell j
420 // with proba for transfered energy eDummy
421
422 if (fasterCode)
423 {
424 (*pNrjTransfData)[j][tDummy][(*pDiffCrossSectionData)[j][tDummy][eDummy]] = eDummy;
425 (*pProbaShellMap)[j][tDummy].push_back((*pDiffCrossSectionData)[j][tDummy][eDummy]);
426 }
427 else { // SI - only if eof is not reached !
428 (*pDiffCrossSectionData)[j][tDummy][eDummy] *= scaleFactor;
429 (*pVecm)[tDummy].push_back(eDummy);
430 }
431 }
432 }
433
434 //Filing maps for the current material into the master maps
435 if (fasterCode) {
436 isUsed1 = true;
437 pNrjTransStorage[mat] = pNrjTransfData;
438 pProbaShellStorage[mat] = pProbaShellMap;
439 }
440 else {
441 pDiffDatatable[mat] = pDiffCrossSectionData;
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[electron]);
460 SetHighEnergyLimit(highEnergyLimit[electron]);
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 initialized " << G4endl
472 << "Energy range: "
473 << LowEnergyLimit() / keV << " keV - "
474 << HighEnergyLimit() / MeV << " MeV for "
475 << particle->GetParticleName()
476 << " with mass (amu) " << particle->GetPDGMass()/proton_mass_c2
477 << " and charge " << particle->GetPDGCharge()
478 << G4endl << G4endl ;
479 }
480
481 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
482
483 fParticleChangeForGamma = GetParticleChangeForGamma();
484 isInitialised = true;
485 }
486
487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
488
489 G4double G4MicroElecInelasticModel_new::CrossSectionPerVolume(const G4Material* material,
490 const G4ParticleDefinition* particleDefinition,
491 G4double ekin,
492 G4double,
493 G4double)
494 {
495 if (verboseLevel > 3) G4cout << "Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" << G4endl;
496
497 G4double density = material->GetTotNbOfAtomsPerVolume();
498 currentMaterial = material->GetName().substr(3, material->GetName().size());
499
500 MapStructure::iterator structPos;
501 structPos = tableMaterialsStructures.find(currentMaterial);
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 found!";
511 G4Exception("G4MicroElecInelasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, str);
512 return 0;
513 }
514 else if(structPos == tableMaterialsStructures.end())
515 {
516 G4String str = "Material ";
517 str += currentMaterial + " Structure not found!";
518 G4Exception("G4MicroElecInelasticModel_new::ComputeCrossSectionPerVolume", "em0002", FatalException, str);
519 return 0;
520 }
521 else {
522 MapData* tableData = tablepos->second;
523 currentMaterialStructure = structPos->second;
524
525 G4double sigma = 0;
526
527 const G4String& particleName = particleDefinition->GetParticleName();
528 G4String nameLocal = particleName;
529 G4int pdg = particleDefinition->GetPDGEncoding();
530 G4int Z = particleDefinition->GetAtomicNumber();
531
532 G4double Zeff = 1.0, Zeff2 = Zeff*Zeff;
533 G4double Mion_c2 = particleDefinition->GetPDGMass();
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->GetInelasticModelLowLimit(pdg);
542 G4double highLim = currentMaterialStructure->GetInelasticModelHighLimit(pdg);
543
544 if (ekin >= lowLim && ekin < highLim)
545 {
546 std::map< G4String, G4MicroElecCrossSectionDataSet_new*, std::less<G4String> >::iterator pos;
547 pos = tableData->find(nameLocal); //find particle type
548
549 if (pos != tableData->end())
550 {
551 G4MicroElecCrossSectionDataSet_new* table = pos->second;
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 < currentMaterialStructure->NumberOfLevels(); i++) {
559 Zeff = BKZ(ekin / (proton_mass_c2 / Mion_c2), Mion_c2 / c_squared, Z, currentMaterialStructure->Energy(i)); // il faut garder le vrai ekin car le calcul à l'interieur de la methode convertie l'énergie en vitesse
560 Zeff2 = Zeff*Zeff;
561 sigma += Zeff2*table->FindShellValue(ekin, i);
562 // il faut utiliser le ekin mis à l'echelle pour chercher la bonne
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_new::CrossSectionPerVolume",
575 "em0002", FatalException,
576 "Model not applicable to particle type.");
577 }
578 }
579 else
580 {
581 return 1 / DBL_MAX;
582 }
583
584 if (verboseLevel > 3)
585 {
586 G4cout << "---> Kinetic energy (eV)=" << ekin / eV << G4endl;
587 G4cout << " - Cross section per Si atom (cm^2)=" << sigma / cm2 << G4endl;
588 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density / (1. / cm) << G4endl;
589 }
590
591 return (sigma)*density;
592 }
593 }
594
595 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
596
597 void G4MicroElecInelasticModel_new::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
598 const G4MaterialCutsCouple* couple,
599 const G4DynamicParticle* particle,
600 G4double,
601 G4double)
602 {
603
604 if (verboseLevel > 3)
605 G4cout << "Calling SampleSecondaries() of G4MicroElecInelasticModel" << G4endl;
606
607 G4int pdg = particle->GetParticleDefinition()->GetPDGEncoding();
608 G4double lowLim = currentMaterialStructure->GetInelasticModelLowLimit(pdg);
609 G4double highLim = currentMaterialStructure->GetInelasticModelHighLimit(pdg);
610
611 G4double ekin = particle->GetKineticEnergy();
612 G4double k = ekin ;
613
614 G4ParticleDefinition* PartDef = particle->GetDefinition();
615 const G4String& particleName = PartDef->GetParticleName();
616 G4String nameLocal2 = particleName ;
617 G4double particleMass = particle->GetDefinition()->GetPDGMass();
618 G4double originalMass = particleMass; // a passer en argument dans samplesecondaryenergy pour évaluer correctement Qmax
619 G4int originalZ = particle->GetDefinition()->GetAtomicNumber();
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 = particle->GetMomentumDirection();
631 G4double totalEnergy = ekin + particleMass;
632 G4double pSquare = ekin * (totalEnergy + particleMass);
633 G4double totalMomentum = std::sqrt(pSquare);
634
635 G4int Shell = 1;
636
637 Shell = RandomSelect(k,nameLocal2,originalMass, originalZ);
638
639 G4double bindingEnergy = currentMaterialStructure->Energy(Shell);
640 G4double limitEnergy = currentMaterialStructure->GetLimitEnergy(Shell);
641
642 G4bool weaklyBound = currentMaterialStructure->IsShellWeaklyBound(Shell);
643
644 if (verboseLevel > 3)
645 {
646 G4cout << "---> Kinetic energy (eV)=" << k/eV << G4endl ;
647 G4cout << "Shell: " << Shell << ", energy: " << bindingEnergy/eV << G4endl;
648 }
649
650
651 if (k<limitEnergy) {
652 if (weaklyBound && k > currentMaterialStructure->GetEnergyGap()) {
653 limitEnergy = currentMaterialStructure->GetEnergyGap();
654 }
655 else return; }
656
657 // sample deexcitation
658
659 std::size_t secNumberInit = 0; // need to know at a certain point the energy of secondaries
660 std::size_t secNumberFinal = 0; // So I'll make the difference and then sum the energies
661
662 // G4cout << currentMaterial << G4endl;
663 G4int Z = currentMaterialStructure->GetZ(Shell);
664 G4int shellEnum = currentMaterialStructure->GetEADL_Enumerator(Shell);
665 if (currentMaterialStructure->IsShellWeaklyBound(Shell)) { shellEnum = -1; }
666
667 if(fAtomDeexcitation && shellEnum >=0)
668 {
669 // G4cout << "enter if deex and shell 0" << G4endl;
670 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellEnum);
671 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
672 secNumberInit = fvect->size();
673 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
674 secNumberFinal = fvect->size();
675 }
676
677 G4double secondaryKinetic=-1000*eV;
678 SEFromFermiLevel = false;
679 if (!fasterCode)
680 {
681 secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef, k, Shell, originalMass, originalZ);
682 }
683 else
684 {
685 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(PartDef, k, Shell) ;
686 }
687
688 if (verboseLevel > 3)
689 {
690 G4cout << "Ionisation process" << G4endl;
691 G4cout << "Shell: " << Shell << " Kin. energy (eV)=" << k/eV
692 << " Sec. energy (eV)=" << secondaryKinetic/eV << G4endl;
693 }
694 G4ThreeVector deltaDirection =
695 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
696 Z, Shell,
697 couple->GetMaterial());
698
699 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
700 {
701 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
702
703 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
704 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
705 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
706 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
707 finalPx /= finalMomentum;
708 finalPy /= finalMomentum;
709 finalPz /= finalMomentum;
710
711 G4ThreeVector direction;
712 direction.set(finalPx,finalPy,finalPz);
713
714 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
715 }
716 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
717
718 // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
719 G4double deexSecEnergy = 0;
720 for (std::size_t j=secNumberInit; j < secNumberFinal; ++j) {
721 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
722 }
723 // correction CI 12/01/2023 limit energy = gap
724 //if (SEFromFermiLevel) limitEnergy = currentMaterialStructure->GetEnergyGap();
725 //fParticleChangeForGamma->SetProposedKineticEnergy(ekin - secondaryKinetic - limitEnergy); //Ef = Ei-(Q-El)-El = Ei-Q
726 //fParticleChangeForGamma->ProposeLocalEnergyDeposit(limitEnergy - deexSecEnergy);
727
728 // correction CI 09/03/2022 limit energy = gap
729 //if (!SEFromFermiLevel && weaklyBound) limitEnergy += currentMaterialStructure->GetInitialEnergy();
730 //if (!SEFromFermiLevel && weaklyBound) limitEnergy += currentMaterialStructure->GetEnergyGap();
731 fParticleChangeForGamma->SetProposedKineticEnergy(ekin - secondaryKinetic-limitEnergy); //Ef = Ei-(Q-El)-El = Ei-Q
732 fParticleChangeForGamma->ProposeLocalEnergyDeposit(limitEnergy-deexSecEnergy);
733
734 if (secondaryKinetic>0)
735 {
736 G4DynamicParticle* dp = new G4DynamicParticle(G4Electron::Electron(), deltaDirection, secondaryKinetic); //Esec = Q-El
737 fvect->push_back(dp);
738 }
739 }
740 }
741
742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
743
744 G4double G4MicroElecInelasticModel_new::RandomizeEjectedElectronEnergy(
745 const G4ParticleDefinition* particleDefinition,
746 G4double k, G4int shell, G4double originalMass, G4int)
747 {
748 G4double secondaryElectronKineticEnergy=0.;
749 if (particleDefinition == G4Electron::ElectronDefinition())
750 {
751 G4double maximumEnergyTransfer=k;
752 G4double crossSectionMaximum = 0.;
753 G4double minEnergy = currentMaterialStructure->GetLimitEnergy(shell);
754 G4double maxEnergy = maximumEnergyTransfer;
755 G4int nEnergySteps = 100;
756
757 G4double value(minEnergy);
758 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
759 G4int step(nEnergySteps);
760 while (step>0)
761 {
762 --step;
763 G4double differentialCrossSection =
764 DifferentialCrossSection(particleDefinition, k, value, shell);
765 crossSectionMaximum = std::max(crossSectionMaximum, differentialCrossSection);
766 value*=stpEnergy;
767 }
768
769 do
770 {
771 secondaryElectronKineticEnergy = G4UniformRand() *
772 (maximumEnergyTransfer-currentMaterialStructure->GetLimitEnergy(shell));
773 } while(G4UniformRand()*crossSectionMaximum >
774 DifferentialCrossSection(particleDefinition, k,
775 (secondaryElectronKineticEnergy+currentMaterialStructure->GetLimitEnergy(shell)),shell));
776 // added 12/01/2023
777 return secondaryElectronKineticEnergy;
778 }
779 else if (particleDefinition == G4Proton::ProtonDefinition())
780 {
781 G4double maximumEnergyTransfer =
782 ComputeElasticQmax(k/(proton_mass_c2/originalMass),
783 currentMaterialStructure->Energy(shell),
784 originalMass/c_squared, electron_mass_c2/c_squared);
785
786 G4double crossSectionMaximum = 0.;
787
788 G4double minEnergy = currentMaterialStructure->GetLimitEnergy(shell);
789 G4double maxEnergy = maximumEnergyTransfer;
790 G4int nEnergySteps = 100;
791
792 G4double value(minEnergy);
793 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
794 G4int step(nEnergySteps);
795
796 while (step>0)
797 {
798 --step;
799 G4double differentialCrossSection =
800 DifferentialCrossSection(particleDefinition, k, value, shell);
801 crossSectionMaximum = std::max(crossSectionMaximum, differentialCrossSection);
802 value*=stpEnergy;
803 }
804
805 G4double energyTransfer = 0.;
806 do
807 {
808 energyTransfer = G4UniformRand() * maximumEnergyTransfer;
809 } while(G4UniformRand()*crossSectionMaximum >
810 DifferentialCrossSection(particleDefinition, k,energyTransfer,shell));
811
812 secondaryElectronKineticEnergy =
813 energyTransfer-currentMaterialStructure->GetLimitEnergy(shell);
814
815 }
816 return std::max(secondaryElectronKineticEnergy, 0.0);
817 }
818
819 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
820
821 G4double G4MicroElecInelasticModel_new::RandomizeEjectedElectronEnergyFromCumulatedDcs(
822 const G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
823 {
824 G4double secondaryElectronKineticEnergy = 0.;
825 G4double random = G4UniformRand();
826 G4bool weaklyBound = currentMaterialStructure->IsShellWeaklyBound(shell);
827 G4double transf = TransferedEnergy(particleDefinition, k, shell, random);
828 if (!weaklyBound) {
829 secondaryElectronKineticEnergy = transf - currentMaterialStructure->GetLimitEnergy(shell); //shell energy for core electrons,
830 if(secondaryElectronKineticEnergy <= 0.) {
831 secondaryElectronKineticEnergy = 0.0;
832 }
833 }
834 else {
835 secondaryElectronKineticEnergy = transf - currentMaterialStructure->GetLimitEnergy(shell);
836 // for weaklybound electrons = gap + average energy in the energy band
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........oooOO0OOooo........oooOO0OOooo......
847
848 G4double G4MicroElecInelasticModel_new::TransferedEnergy(
849 const G4ParticleDefinition* particleDefinition,
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.* (electron_mass_c2 / proton_mass_c2) * k;
869 G4double bindingEnergy = currentMaterialStructure->GetLimitEnergy(ionizationLevelIndex);
870
871 if (particleDefinition == G4Electron::ElectronDefinition())
872 {
873 dataDiffCSMap::iterator iterator_Nrj;
874 iterator_Nrj = eNrjTransStorage.find(currentMaterial);
875
876 dataProbaShellMap::iterator iterator_Proba;
877 iterator_Proba = eProbaShellStorage.find(currentMaterial);
878
879 incidentEnergyMap::iterator iterator_Tdummy;
880 iterator_Tdummy = eIncidentEnergyStorage.find(currentMaterial);
881
882 if(iterator_Nrj == eNrjTransStorage.end() || iterator_Proba == eProbaShellStorage.end() ||
883 iterator_Tdummy == eIncidentEnergyStorage.end())
884 {
885 G4String str = "Material ";
886 str += currentMaterial + " not found!";
887 G4Exception("G4MicroElecInelasticModel_new::TransferedEnergy", "em0002",
888 FatalException, str);
889 }
890 else {
891 vector<TriDimensionMap>* eNrjTransfData = iterator_Nrj->second; //Storage of possible transfer energies
892 vector<VecMap>* eProbaShellMap = iterator_Proba->second; //Storage of probabilities for energy transfer
893 vector<G4double>* eTdummyVec = iterator_Tdummy->second; //Incident energies for interpolation
894
895 // k should be in eV auto, :std::vector<double>::iterator
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 situations where random >last vector element
902 if (random <= (*eProbaShellMap)[ionizationLevelIndex][(*k1)].back()
903 && random <= (*eProbaShellMap)[ionizationLevelIndex][(*k2)].back())
904 {
905 auto prob12 =
906 std::upper_bound((*eProbaShellMap)[ionizationLevelIndex][(*k1)].begin(),
907 (*eProbaShellMap)[ionizationLevelIndex][(*k1)].end(),
908 random);
909
910 auto prob11 = prob12 - 1;
911
912 auto prob22 =
913 std::upper_bound((*eProbaShellMap)[ionizationLevelIndex][(*k2)].begin(),
914 (*eProbaShellMap)[ionizationLevelIndex][(*k2)].end(),
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 transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
927 if (valuePROB11 == 0) nrjTransf11 = bindingEnergy;
928 else nrjTransf11 = (*eNrjTransfData)[ionizationLevelIndex][valueK1][valuePROB11];
929 if (valuePROB12 == 1)
930 {
931 if ((valueK1 + bindingEnergy) / 2. > valueK1)
932 maximumEnergyTransfer1 = valueK1;
933 else
934 maximumEnergyTransfer1 = (valueK1 + bindingEnergy) / 2.;
935
936 nrjTransf12 = maximumEnergyTransfer1;
937 }
938 else
939 nrjTransf12 = (*eNrjTransfData)[ionizationLevelIndex][valueK1][valuePROB12];
940
941 if (valuePROB21 == 0) nrjTransf21 = bindingEnergy;
942 else nrjTransf21 = (*eNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB21];
943 if (valuePROB22 == 1)
944 {
945 if ((valueK2 + bindingEnergy) / 2. > valueK2) maximumEnergyTransfer2 = valueK2;
946 else maximumEnergyTransfer2 = (valueK2 + bindingEnergy) / 2.;
947
948 nrjTransf22 = maximumEnergyTransfer2;
949 }
950 else nrjTransf22 = (*eNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB22];
951
952 }
953 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
954 if (random > (*eProbaShellMap)[ionizationLevelIndex][(*k1)].back())
955 {
956 auto prob22 =
957 std::upper_bound((*eProbaShellMap)[ionizationLevelIndex][(*k2)].begin(),
958 (*eProbaShellMap)[ionizationLevelIndex][(*k2)].end(),
959 random);
960 auto prob21 = prob22 - 1;
961
962 valueK1 = *k1;
963 valueK2 = *k2;
964 valuePROB21 = *prob21;
965 valuePROB22 = *prob22;
966
967 nrjTransf21 = (*eNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB21];
968 nrjTransf22 = (*eNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB22];
969
970 G4double interpolatedvalue2 = Interpolate(valuePROB21,
971 valuePROB22,
972 random,
973 nrjTransf21,
974 nrjTransf22);
975
976 // zeros are explicitly set
977 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
978
979 return value;
980 }
981 }
982 }
983 else if (particleDefinition == G4Proton::ProtonDefinition())
984 {
985 // k should be in eV
986 dataDiffCSMap::iterator iterator_Nrj;
987 iterator_Nrj = pNrjTransStorage.find(currentMaterial);
988
989 dataProbaShellMap::iterator iterator_Proba;
990 iterator_Proba = pProbaShellStorage.find(currentMaterial);
991
992 incidentEnergyMap::iterator iterator_Tdummy;
993 iterator_Tdummy = pIncidentEnergyStorage.find(currentMaterial);
994
995 if (iterator_Nrj == pNrjTransStorage.end() || iterator_Proba == pProbaShellStorage.end() ||
996 iterator_Tdummy == pIncidentEnergyStorage.end())
997 {
998 G4String str = "Material ";
999 str += currentMaterial + " not found!";
1000 G4Exception("G4MicroElecInelasticModel_new::TransferedEnergy", "em0002",
1001 FatalException, str);
1002 }
1003 else
1004 {
1005 vector<TriDimensionMap>* pNrjTransfData = iterator_Nrj->second; //Storage of possible transfer energies
1006 vector<VecMap>* pProbaShellMap = iterator_Proba->second; //Storage of probabilities for energy transfer
1007 vector<G4double>* pTdummyVec = iterator_Tdummy->second; //Incident energies for interpolation
1008
1009 auto k2 = std::upper_bound(pTdummyVec->begin(),
1010 pTdummyVec->end(),
1011 k);
1012
1013 auto k1 = k2 - 1;
1014
1015 // SI : the following condition avoids situations where random > last vector element,
1016 // for eg. when the last element is zero
1017 if (random <= (*pProbaShellMap)[ionizationLevelIndex][(*k1)].back()
1018 && random <= (*pProbaShellMap)[ionizationLevelIndex][(*k2)].back())
1019 {
1020 auto prob12 =
1021 std::upper_bound((*pProbaShellMap)[ionizationLevelIndex][(*k1)].begin(),
1022 (*pProbaShellMap)[ionizationLevelIndex][(*k1)].end(),
1023 random);
1024 auto prob11 = prob12 - 1;
1025 auto prob22 =
1026 std::upper_bound((*pProbaShellMap)[ionizationLevelIndex][(*k2)].begin(),
1027 (*pProbaShellMap)[ionizationLevelIndex][(*k2)].end(),
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 getting transfered energy < binding energy
1039 // and forces cumxs = 1 for maximum energy transfer.
1040 if (valuePROB11 == 0) nrjTransf11 = bindingEnergy;
1041 else nrjTransf11 = (*pNrjTransfData)[ionizationLevelIndex][valueK1][valuePROB11];
1042
1043 if (valuePROB12 == 1) nrjTransf12 = maximumEnergyTransferP;
1044 else nrjTransf12 = (*pNrjTransfData)[ionizationLevelIndex][valueK1][valuePROB12];
1045
1046 if (valuePROB21 == 0) nrjTransf21 = bindingEnergy;
1047 else nrjTransf21 = (*pNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB21];
1048
1049 if (valuePROB22 == 1) nrjTransf22 = maximumEnergyTransferP;
1050 else nrjTransf22 = (*pNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB22];
1051
1052 }
1053
1054 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1055 if (random > (*pProbaShellMap)[ionizationLevelIndex][(*k1)].back())
1056 {
1057 auto prob22 =
1058 std::upper_bound((*pProbaShellMap)[ionizationLevelIndex][(*k2)].begin(),
1059 (*pProbaShellMap)[ionizationLevelIndex][(*k2)].end(),
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)[ionizationLevelIndex][valueK2][valuePROB21];
1070 nrjTransf22 = (*pNrjTransfData)[ionizationLevelIndex][valueK2][valuePROB22];
1071
1072 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1073 valuePROB22,
1074 random,
1075 nrjTransf21,
1076 nrjTransf22);
1077
1078 // zeros are explicitly set
1079 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1080 return value;
1081 }
1082 }
1083 }
1084 // End electron and proton cases
1085
1086 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
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........oooOO0OOooo........oooOO0OOooo......
1108
1109 G4double G4MicroElecInelasticModel_new::DifferentialCrossSection(
1110 const G4ParticleDefinition * particleDefinition,
1111 G4double k,
1112 G4double energyTransfer,
1113 G4int LevelIndex)
1114 {
1115 G4double sigma = 0.;
1116
1117 if (energyTransfer >= currentMaterialStructure->GetLimitEnergy(LevelIndex))
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::ElectronDefinition())
1132 {
1133
1134 dataDiffCSMap::iterator iterator_Proba;
1135 iterator_Proba = eDiffDatatable.find(currentMaterial);
1136
1137 incidentEnergyMap::iterator iterator_Nrj;
1138 iterator_Nrj = eIncidentEnergyStorage.find(currentMaterial);
1139
1140 TranfEnergyMap::iterator iterator_TransfNrj;
1141 iterator_TransfNrj = eVecmStorage.find(currentMaterial);
1142
1143 if (iterator_Proba != eDiffDatatable.end() && iterator_Nrj != eIncidentEnergyStorage.end()
1144 && iterator_TransfNrj!= eVecmStorage.end())
1145 {
1146 vector<TriDimensionMap>* eDiffCrossSectionData = (iterator_Proba->second);
1147 vector<G4double>* eTdummyVec = iterator_Nrj->second; //Incident energies for interpolation
1148 VecMap* eVecm = iterator_TransfNrj->second;
1149
1150 // k should be in eV and energy transfer eV also
1151 auto t2 = std::upper_bound(eTdummyVec->begin(), eTdummyVec->end(), k);
1152 auto t1 = t2 - 1;
1153 // SI : the following condition avoids situations where energyTransfer >last vector element
1154 if (energyTransfer <= (*eVecm)[(*t1)].back() && energyTransfer <= (*eVecm)[(*t2)].back())
1155 {
1156 auto e12 = std::upper_bound((*eVecm)[(*t1)].begin(), (*eVecm)[(*t1)].end(), energyTransfer);
1157 auto e11 = e12 - 1;
1158 auto e22 = std::upper_bound((*eVecm)[(*t2)].begin(), (*eVecm)[(*t2)].end(), energyTransfer);
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)[LevelIndex][valueT1][valueE11];
1169 xs12 = (*eDiffCrossSectionData)[LevelIndex][valueT1][valueE12];
1170 xs21 = (*eDiffCrossSectionData)[LevelIndex][valueT2][valueE21];
1171 xs22 = (*eDiffCrossSectionData)[LevelIndex][valueT2][valueE22];
1172 }
1173 }
1174 else {
1175 G4String str = "Material ";
1176 str += currentMaterial + " not found!";
1177 G4Exception("G4MicroElecDielectricModels::DifferentialCrossSection", "em0002", FatalException, str);
1178 }
1179 }
1180
1181 if (particleDefinition == G4Proton::ProtonDefinition())
1182 {
1183 dataDiffCSMap::iterator iterator_Proba;
1184 iterator_Proba = pDiffDatatable.find(currentMaterial);
1185
1186 incidentEnergyMap::iterator iterator_Nrj;
1187 iterator_Nrj = pIncidentEnergyStorage.find(currentMaterial);
1188
1189 TranfEnergyMap::iterator iterator_TransfNrj;
1190 iterator_TransfNrj = pVecmStorage.find(currentMaterial);
1191
1192 if (iterator_Proba != pDiffDatatable.end() && iterator_Nrj != pIncidentEnergyStorage.end()
1193 && iterator_TransfNrj != pVecmStorage.end())
1194 {
1195 vector<TriDimensionMap>* pDiffCrossSectionData = (iterator_Proba->second);
1196 vector<G4double>* pTdummyVec = iterator_Nrj->second; //Incident energies for interpolation
1197 VecMap* pVecm = iterator_TransfNrj->second;
1198
1199 // k should be in eV and energy transfer eV also
1200 auto t2 =
1201 std::upper_bound(pTdummyVec->begin(), pTdummyVec->end(), k);
1202 auto t1 = t2 - 1;
1203 if (energyTransfer <= (*pVecm)[(*t1)].back() && energyTransfer <= (*pVecm)[(*t2)].back())
1204 {
1205 auto e12 = std::upper_bound((*pVecm)[(*t1)].begin(), (*pVecm)[(*t1)].end(), energyTransfer);
1206 auto e11 = e12 - 1;
1207 auto e22 = std::upper_bound((*pVecm)[(*t2)].begin(), (*pVecm)[(*t2)].end(), energyTransfer);
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)[LevelIndex][valueT1][valueE11];
1218 xs12 = (*pDiffCrossSectionData)[LevelIndex][valueT1][valueE12];
1219 xs21 = (*pDiffCrossSectionData)[LevelIndex][valueT2][valueE21];
1220 xs22 = (*pDiffCrossSectionData)[LevelIndex][valueT2][valueE22];
1221 }
1222 }
1223 else {
1224 G4String str = "Material ";
1225 str += currentMaterial + " not found!";
1226 G4Exception("G4MicroElecDielectricModels::DifferentialCrossSection", "em0002", FatalException, str);
1227 }
1228 }
1229
1230 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
1231 if (xsProduct != 0.)
1232 {
1233 sigma = QuadInterpolator( valueE11, valueE12,
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........oooOO0OOooo........oooOO0OOooo......
1247
1248
1249 G4double G4MicroElecInelasticModel_new::Interpolate(G4double e1,
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 > 0. && !fasterCode)
1261 {
1262 G4double a = std::log(xs2/xs1)/ std::log(e2/e1);
1263 G4double b = std::log(xs2) - a * std::log(e2);
1264 G4double sigma = a * std::log(e) + b;
1265 value = (std::exp(sigma));
1266 }
1267
1268 // Switch to log-lin interpolation for faster code
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 - e1) / (e2 - e1)));
1274 }
1275
1276 // Switch to lin-lin interpolation for faster code
1277 // in case one of xs1 or xs2 (=cum proba) value is zero
1278 else
1279 {
1280 G4double d1 = xs1;
1281 G4double d2 = xs2;
1282 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
1283 }
1284
1285 return value;
1286 }
1287
1288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1289
1290 G4double G4MicroElecInelasticModel_new::QuadInterpolator(G4double e11, G4double e12,
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(e11, e12, e, xs11, xs12);
1298 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
1299 G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
1300 return value;
1301 }
1302
1303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1304
1305 G4int G4MicroElecInelasticModel_new::RandomSelect(G4double k, const G4String& particle, G4double originalMass, G4int originalZ )
1306 {
1307 G4int level = 0;
1308
1309 TCSMap::iterator tablepos;
1310 tablepos = tableTCS.find(currentMaterial);
1311 if (tablepos == tableTCS.end()) {
1312 G4Exception("G4MicroElecInelasticModel_new::RandomSelect","em0002",FatalException,"Model not applicable to material");
1313 return level;
1314 }
1315
1316 MapData* tableData = tablepos->second;
1317
1318 std::map< G4String,G4MicroElecCrossSectionDataSet_new*,std::less<G4String> >::iterator pos;
1319 pos = tableData->find(particle);
1320
1321 std::vector<G4double> Zeff(currentMaterialStructure->NumberOfLevels(), 1.0);
1322 if(originalMass>proton_mass_c2) {
1323 for(G4int nl=0;nl<currentMaterialStructure->NumberOfLevels();nl++) {
1324 Zeff[nl] = BKZ(k/(proton_mass_c2/originalMass), originalMass/c_squared, originalZ, currentMaterialStructure->Energy(nl));
1325 }
1326 }
1327
1328 if (pos != tableData->end())
1329 {
1330 G4MicroElecCrossSectionDataSet_new* table = pos->second;
1331
1332 if (table != 0)
1333 {
1334 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1335 const G4int n = (G4int)table->NumberOfComponents();
1336 G4int i = (G4int)n;
1337 G4double value = 0.;
1338
1339 while (i>0)
1340 {
1341 --i;
1342 valuesBuffer[i] = table->GetComponent(i)->FindValue(k)*Zeff[i]*Zeff[i];
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_new::RandomSelect","em0002",FatalException,"Model not applicable to particle type.");
1368 }
1369
1370 return level;
1371 }
1372
1373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1374
1375 G4double G4MicroElecInelasticModel_new::ComputeRelativistVelocity(G4double E, G4double mass) {
1376 G4double x = E/mass;
1377 return c_light*std::sqrt(x*(x + 2.0))/(x + 1.0);
1378 }
1379
1380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1381
1382 G4double G4MicroElecInelasticModel_new::ComputeElasticQmax(G4double T1i, G4double T2i, G4double M1, G4double M2) {
1383 G4double v1i = ComputeRelativistVelocity(T1i, M1);
1384 G4double v2i = ComputeRelativistVelocity(T2i, M2);
1385
1386 G4double v2f = 2*M1/(M1+M2)*v1i + (M2-M1)/(M1+M2)*-1*v2i;
1387 G4double vtransfer2a = v2f*v2f-v2i*v2i;
1388
1389 v2f = 2*M1/(M1+M2)*v1i + (M2-M1)/(M1+M2)*v2i;
1390 G4double vtransfer2b = v2f*v2f-v2i*v2i;
1391
1392 G4double vtransfer2 = std::max(vtransfer2a, vtransfer2b);
1393 return 0.5*M2*vtransfer2;
1394 }
1395
1396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1397
1398 G4double G4MicroElecInelasticModel_new::stepFunc(G4double x) {
1399 return (x < 0.) ? 1.0 : 0.0;
1400 }
1401
1402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1403
1404 G4double G4MicroElecInelasticModel_new::vrkreussler(G4double v, G4double vF)
1405 {
1406 G4double r = vF*( std::pow(v/vF+1., 3.) - fabs(std::pow(v/vF-1., 3.))
1407 + 4.*(v/vF)*(v/vF) ) + stepFunc(v/vF-1.) * (3./2.*v/vF -
1408 4.*(v/vF)*(v/vF) + 3.*std::pow(v/vF, 3.)
1409 - 0.5*std::pow(v/vF, 5.));
1410 return r/(10.*v/vF);
1411 }
1412
1413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1414
1415 G4double G4MicroElecInelasticModel_new::BKZ(G4double Ep, G4double mp, G4int Zp, G4double Eplasmon)
1416 {
1417 // need atomic unit conversion
1418 G4double hbar = hbar_Planck, hbar2 = hbar*hbar, me = electron_mass_c2/c_squared, Ry = me*elm_coupling*elm_coupling/(2*hbar2);
1419 G4double hartree = 2*Ry, a0 = Bohr_radius, velocity = a0*hartree/hbar;
1420 G4double vp = ComputeRelativistVelocity(Ep, mp);
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./3.);
1427 G4double c = 0.9;
1428 G4double vr = vrkreussler(vp /*in u.a*/, vF /*in u.a*/);
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/7.);
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.+std::pow(2.*l0*vF,2.)));
1440 }
1441
1442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1443