Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4MicroElecInelasticModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// G4MicroElecInelasticModel.cc, 2011/08/29 A.Valentin, M. Raine
28//
29// Based on the following publications
30//
31// - Inelastic cross-sections of low energy electrons in silicon
32// for the simulation of heavy ion tracks with theGeant4-DNA toolkit,
33// NSS Conf. Record 2010, pp. 80-85.
34// - Geant4 physics processes for microdosimetry simulation:
35// very low energy electromagnetic models for electrons in Si,
36// NIM B, vol. 288, pp. 66 - 73, 2012.
37// - Geant4 physics processes for microdosimetry simulation:
38// very low energy electromagnetic models for protons and
39// heavy ions in Si, NIM B, vol. 287, pp. 124 - 129, 2012.
40//
41//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42
44
45#include "globals.hh"
47#include "G4SystemOfUnits.hh"
48#include "G4ios.hh"
49#include "G4UnitsTable.hh"
52#include "G4LossTableManager.hh"
55#include "G4Electron.hh"
56#include "G4Proton.hh"
57#include "G4GenericIon.hh"
59#include "G4NistManager.hh"
61#include "G4DeltaAngle.hh"
62
63//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64
65using namespace std;
66
67//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68
70 const G4String& nam)
71:G4VEmModel(nam),isInitialised(false)
72{
74
75 verboseLevel= 0;
76 // Verbosity scale:
77 // 0 = nothing
78 // 1 = warning for energy non-conservation
79 // 2 = details of energy budget
80 // 3 = calculation of cross sections, file openings, sampling of atoms
81 // 4 = entering in methods
82
83 if( verboseLevel>0 )
84 {
85 G4cout << "MicroElec inelastic model is constructed " << G4endl;
86 }
87
88 //Mark this model as "applicable" for atomic deexcitation
90 fAtomDeexcitation = nullptr;
92
93 // default generator
95
96 // Selection of computation method
97 fasterCode = true; //false;
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
103{
104 // Cross section
105 for (auto & pos : tableData)
106 {
107 G4MicroElecCrossSectionDataSet* table = pos.second;
108 delete table;
109 }
110
111 // Final state
112 eVecm.clear();
113 pVecm.clear();
114
115}
116
117//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118
120 const G4DataVector& /*cuts*/)
121{
122
123 if (verboseLevel > 3)
124 G4cout << "Calling G4MicroElecInelasticModel::Initialise()" << G4endl;
125
126 // Energy limits
127
128 G4String fileElectron("microelec/sigma_inelastic_e_Si");
129 G4String fileProton("microelec/sigma_inelastic_p_Si");
130
133 G4String electron = electronDef->GetParticleName();
134 G4String proton = protonDef->GetParticleName();;
135
136 G4double scaleFactor = 1e-18 * cm *cm;
137
138 const char* path = G4FindDataDir("G4LEDATA");
139
140 // *** ELECTRON
141 tableFile[electron] = fileElectron;
142 lowEnergyLimit[electron] = 16.7 * eV;
143 highEnergyLimit[electron] = 100.0 * MeV;
144
145 // Cross section
147 tableE->LoadData(fileElectron);
148
149 tableData[electron] = tableE;
150
151 // Final state
152
153 std::ostringstream eFullFileName;
154
155 if (fasterCode) eFullFileName << path << "/microelec/sigmadiff_cumulated_inelastic_e_Si.dat";
156 else eFullFileName << path << "/microelec/sigmadiff_inelastic_e_Si.dat";
157
158 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
159
160 if (!eDiffCrossSection)
161 {
162 if (fasterCode) G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
163 FatalException,"Missing data file:/microelec/sigmadiff_cumulated_inelastic_e_Si.dat");
164
165 else G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
166 FatalException,"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
167 }
168
169 // Clear the arrays for re-initialization case (MT mode)
170 // Octobre 22nd, 2014 - Melanie Raine
171 eTdummyVec.clear();
172 pTdummyVec.clear();
173
174 eVecm.clear();
175 pVecm.clear();
176
177 for (int j=0; j<6; j++)
178 {
179 eProbaShellMap[j].clear();
180 pProbaShellMap[j].clear();
181
182 eDiffCrossSectionData[j].clear();
183 pDiffCrossSectionData[j].clear();
184
185 eNrjTransfData[j].clear();
186 pNrjTransfData[j].clear();
187 }
188
189 //
190 eTdummyVec.push_back(0.);
191 while(!eDiffCrossSection.eof())
192 {
193 G4double tDummy;
194 G4double eDummy;
195 eDiffCrossSection>>tDummy>>eDummy;
196 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
197
198 G4double tmp;
199 for (int j=0; j<6; j++)
200 {
201 eDiffCrossSection>> tmp;
202
203 eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
204
205 if (fasterCode)
206 {
207 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
208 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
209 }
210 else
211 {
212 // SI - only if eof is not reached !
213 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
214 eVecm[tDummy].push_back(eDummy);
215 }
216
217 }
218 }
219 //
220
221 // *** PROTON
222 tableFile[proton] = fileProton;
223 lowEnergyLimit[proton] = 50. * keV;
224 highEnergyLimit[proton] = 10. * GeV;
225
226 // Cross section
228 tableP->LoadData(fileProton);
229 tableData[proton] = tableP;
230
231 // Final state
232 std::ostringstream pFullFileName;
233
234 if (fasterCode) pFullFileName << path << "/microelec/sigmadiff_cumulated_inelastic_p_Si.dat";
235 else pFullFileName << path << "/microelec/sigmadiff_inelastic_p_Si.dat";
236
237 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
238
239 if (!pDiffCrossSection)
240 {
241 if (fasterCode) G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
242 FatalException,"Missing data file:/microelec/sigmadiff_cumulated_inelastic_p_Si.dat");
243
244 else G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
245 FatalException,"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
246 }
247
248 pTdummyVec.push_back(0.);
249 while(!pDiffCrossSection.eof())
250 {
251 G4double tDummy;
252 G4double eDummy;
253 pDiffCrossSection>>tDummy>>eDummy;
254 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
255 for (int j=0; j<6; j++)
256 {
257 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
258
259 if (fasterCode)
260 {
261 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
262 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
263 }
264 else
265 {
266 // SI - only if eof is not reached !
267 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
268 pVecm[tDummy].push_back(eDummy);
269 }
270 }
271 }
272
273 if (particle==electronDef)
274 {
275 SetLowEnergyLimit(lowEnergyLimit[electron]);
276 SetHighEnergyLimit(highEnergyLimit[electron]);
277 }
278
279 if (particle==protonDef)
280 {
281 SetLowEnergyLimit(lowEnergyLimit[proton]);
282 SetHighEnergyLimit(highEnergyLimit[proton]);
283 }
284
285 if( verboseLevel>0 )
286 {
287 G4cout << "MicroElec Inelastic model is initialized " << G4endl
288 << "Energy range: "
289 << LowEnergyLimit() / keV << " keV - "
290 << HighEnergyLimit() / MeV << " MeV for "
291 << particle->GetParticleName()
292 << " with mass (amu) " << particle->GetPDGMass()/proton_mass_c2
293 << " and charge " << particle->GetPDGCharge()
294 << G4endl << G4endl ;
295 }
296
297 //
298 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
299
300 if (isInitialised) { return; }
302 isInitialised = true;
303}
304
305//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
306
308 const G4ParticleDefinition* particleDefinition,
309 G4double ekin,
310 G4double,
311 G4double)
312{
313 if (verboseLevel > 3)
314 G4cout << "Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" << G4endl;
315
316 G4double density = material->GetTotNbOfAtomsPerVolume();
317
318 // Calculate total cross section for model
319 G4double lowLim = 0;
320 G4double highLim = 0;
321 G4double sigma=0;
322
323 const G4String& particleName = particleDefinition->GetParticleName();
324 G4String nameLocal = particleName ;
325
326 G4double Zeff2 = 1.0;
327 G4double Mion_c2 = particleDefinition->GetPDGMass();
328
329 if (Mion_c2 > proton_mass_c2)
330 {
331 G4ionEffectiveCharge EffCharge ;
332 G4double Zeff = EffCharge.EffectiveCharge(particleDefinition, material,ekin);
333 Zeff2 = Zeff*Zeff;
334
335 if (verboseLevel > 3)
336 G4cout << "Before scaling : " << G4endl
337 << "Particle : " << nameLocal << ", mass : " << Mion_c2/proton_mass_c2 << "*mp, charge " << Zeff
338 << ", Ekin (eV) = " << ekin/eV << G4endl ;
339
340 ekin *= proton_mass_c2/Mion_c2 ;
341 nameLocal = "proton" ;
342
343 if (verboseLevel > 3)
344 G4cout << "After scaling : " << G4endl
345 << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl ;
346 }
347
348 if (material == nistSi || material->GetBaseMaterial() == nistSi)
349 {
350 auto pos1 = lowEnergyLimit.find(nameLocal);
351 if (pos1 != lowEnergyLimit.end())
352 {
353 lowLim = pos1->second;
354 }
355
356 auto pos2 = highEnergyLimit.find(nameLocal);
357 if (pos2 != highEnergyLimit.end())
358 {
359 highLim = pos2->second;
360 }
361
362 if (ekin >= lowLim && ekin < highLim)
363 {
364 auto pos = tableData.find(nameLocal);
365 if (pos != tableData.end())
366 {
367 G4MicroElecCrossSectionDataSet* table = pos->second;
368 if (table != nullptr)
369 {
370 sigma = table->FindValue(ekin);
371 }
372 }
373 else
374 {
375 G4Exception("G4MicroElecInelasticModel::CrossSectionPerVolume","em0002",
376 FatalException,"Model not applicable to particle type.");
377 }
378 }
379 else
380 {
381 if (nameLocal!="e-")
382 {
383 // G4cout << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl;
384 // G4cout << "### Warning: particle energy out of bounds! ###" << G4endl;
385 }
386 }
387
388 if (verboseLevel > 3)
389 {
390 G4cout << "---> Kinetic energy (eV)=" << ekin/eV << G4endl;
391 G4cout << " - Cross section per Si atom (cm^2)=" << sigma*Zeff2/cm2 << G4endl;
392 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./cm) << G4endl;
393 }
394
395 } // if (SiMaterial)
396 return sigma*density*Zeff2;
397}
398
399//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400
401void G4MicroElecInelasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
402 const G4MaterialCutsCouple* couple,
403 const G4DynamicParticle* particle,
404 G4double,
405 G4double)
406{
407
408 if (verboseLevel > 3)
409 G4cout << "Calling SampleSecondaries() of G4MicroElecInelasticModel" << G4endl;
410
411 G4double lowLim = 0;
412 G4double highLim = 0;
413
414 G4double ekin = particle->GetKineticEnergy();
415 G4double k = ekin ;
416
417 G4ParticleDefinition* PartDef = particle->GetDefinition();
418 const G4String& particleName = PartDef->GetParticleName();
419 G4String nameLocal2 = particleName ;
420 G4double particleMass = particle->GetDefinition()->GetPDGMass();
421
422 if (particleMass > proton_mass_c2)
423 {
424 k *= proton_mass_c2/particleMass ;
425 PartDef = G4Proton::ProtonDefinition();
426 nameLocal2 = "proton" ;
427 }
428
429 auto pos1 = lowEnergyLimit.find(nameLocal2);
430 if (pos1 != lowEnergyLimit.end())
431 {
432 lowLim = pos1->second;
433 }
434
435 auto pos2 = highEnergyLimit.find(nameLocal2);
436 if (pos2 != highEnergyLimit.end())
437 {
438 highLim = pos2->second;
439 }
440
441 if (k >= lowLim && k < highLim)
442 {
443 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
444 G4double totalEnergy = ekin + particleMass;
445 G4double pSquare = ekin * (totalEnergy + particleMass);
446 G4double totalMomentum = std::sqrt(pSquare);
447
448 G4int Shell = 0;
449
450 /* if (!fasterCode)*/ Shell = RandomSelect(k,nameLocal2);
451
452 // SI: The following protection is necessary to avoid infinite loops :
453 // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
454 // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
455 // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
456 // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
457
458 G4double bindingEnergy = SiStructure.Energy(Shell);
459
460 if (verboseLevel > 3)
461 {
462 G4cout << "---> Kinetic energy (eV)=" << k/eV << G4endl ;
463 G4cout << "Shell: " << Shell << ", energy: " << bindingEnergy/eV << G4endl;
464 }
465
466 // sample deexcitation
467 std::size_t secNumberInit = 0; // need to know at a certain point the energy of secondaries
468 std::size_t secNumberFinal = 0; // So I'll make the difference and then sum the energies
469
470 //SI: additional protection if tcs interpolation method is modified
471 if (k<bindingEnergy) return;
472
473 G4int Z = 14;
474
475 if(fAtomDeexcitation && Shell > 2) {
476
478
479 if (Shell == 4)
480 {
482 }
483 else if (Shell == 3)
484 {
486 }
487
488 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
489 secNumberInit = fvect->size();
490 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
491 secNumberFinal = fvect->size();
492 }
493
494 G4double secondaryKinetic=-1000*eV;
495
496 if (!fasterCode)
497 {
498 secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef,k,Shell);
499 }
500 else
501 {
502 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(PartDef,k,Shell);
503 }
504
505 if (verboseLevel > 3)
506 {
507 G4cout << "Ionisation process" << G4endl;
508 G4cout << "Shell: " << Shell << " Kin. energy (eV)=" << k/eV
509 << " Sec. energy (eV)=" << secondaryKinetic/eV << G4endl;
510 }
511
512 G4ThreeVector deltaDirection =
513 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
514 Z, Shell,
515 couple->GetMaterial());
516
518 {
519 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
520 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
521 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
522 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
523 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
524 finalPx /= finalMomentum;
525 finalPy /= finalMomentum;
526 finalPz /= finalMomentum;
527
528 G4ThreeVector direction;
529 direction.set(finalPx,finalPy,finalPz);
530
532 }
533 else
535
536 // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
537 G4double deexSecEnergy = 0;
538 for (std::size_t j=secNumberInit; j < secNumberFinal; ++j) {
539 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
540
541 fParticleChangeForGamma->SetProposedKineticEnergy(ekin-bindingEnergy-secondaryKinetic);
542 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy-deexSecEnergy);
543
544 if (secondaryKinetic>0)
545 {
546 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
547 fvect->push_back(dp);
548 }
549 }
550}
551
552//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553
554G4double G4MicroElecInelasticModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
555 G4double k, G4int shell)
556{
557 if (particleDefinition == G4Electron::ElectronDefinition())
558 {
559 G4double maximumEnergyTransfer=0.;
560 if ((k+SiStructure.Energy(shell))/2. > k) maximumEnergyTransfer=k;
561 else maximumEnergyTransfer = (k+SiStructure.Energy(shell))/2.;
562
563 G4double crossSectionMaximum = 0.;
564
565 G4double minEnergy = SiStructure.Energy(shell);
566 G4double maxEnergy = maximumEnergyTransfer;
567 G4int nEnergySteps = 100;
568
569 G4double value(minEnergy);
570 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
571 G4int step(nEnergySteps);
572 while (step>0)
573 {
574 step--;
575 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
576 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
577 value*=stpEnergy;
578 }
579
580 G4double secondaryElectronKineticEnergy=0.;
581 do
582 {
583 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
584 } while(G4UniformRand()*crossSectionMaximum >
585 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
586
587 return secondaryElectronKineticEnergy;
588 }
589
590 if (particleDefinition == G4Proton::ProtonDefinition())
591 {
592 G4double maximumEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
593 G4double crossSectionMaximum = 0.;
594
595 G4double minEnergy = SiStructure.Energy(shell);
596 G4double maxEnergy = maximumEnergyTransfer;
597 G4int nEnergySteps = 100;
598
599 G4double value(minEnergy);
600 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
601 G4int step(nEnergySteps);
602 while (step>0)
603 {
604 step--;
605 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
606 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
607 value*=stpEnergy;
608 }
609
610 G4double secondaryElectronKineticEnergy = 0.;
611 do
612 {
613 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
614
615 } while(G4UniformRand()*crossSectionMaximum >
616 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
617 return secondaryElectronKineticEnergy;
618 }
619
620 return 0;
621}
622
623//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
624
625// The following section is not used anymore but is kept for memory
626// GetAngularDistribution()->SampleDirectionForShell is used instead
627
628/*void G4MicroElecInelasticModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
629 G4double k,
630 G4double secKinetic,
631 G4double & cosTheta,
632 G4double & phi )
633 {
634 if (particleDefinition == G4Electron::ElectronDefinition())
635 {
636 phi = twopi * G4UniformRand();
637 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
638 cosTheta = std::sqrt(1.-sin2O);
639 }
640
641 if (particleDefinition == G4Proton::ProtonDefinition())
642 {
643 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
644 phi = twopi * G4UniformRand();
645 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
646 }
647
648 else
649 {
650 G4double maxSecKinetic = 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
651 phi = twopi * G4UniformRand();
652 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
653 }
654 }
655 */
656
657//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
658
660 G4double k,
661 G4double energyTransfer,
662 G4int LevelIndex)
663{
664 G4double sigma = 0.;
665
666 if (energyTransfer >= SiStructure.Energy(LevelIndex))
667 {
668 G4double valueT1 = 0;
669 G4double valueT2 = 0;
670 G4double valueE21 = 0;
671 G4double valueE22 = 0;
672 G4double valueE12 = 0;
673 G4double valueE11 = 0;
674
675 G4double xs11 = 0;
676 G4double xs12 = 0;
677 G4double xs21 = 0;
678 G4double xs22 = 0;
679
680 if (particleDefinition == G4Electron::ElectronDefinition())
681 {
682 // k should be in eV and energy transfer eV also
683 auto t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
684 auto t1 = t2-1;
685 // The following condition avoids situations where energyTransfer >last vector element
686 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
687 {
688 auto e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
689 auto e11 = e12-1;
690 auto e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
691 auto e21 = e22-1;
692
693 valueT1 =*t1;
694 valueT2 =*t2;
695 valueE21 =*e21;
696 valueE22 =*e22;
697 valueE12 =*e12;
698 valueE11 =*e11;
699
700 xs11 = eDiffCrossSectionData[LevelIndex][valueT1][valueE11];
701 xs12 = eDiffCrossSectionData[LevelIndex][valueT1][valueE12];
702 xs21 = eDiffCrossSectionData[LevelIndex][valueT2][valueE21];
703 xs22 = eDiffCrossSectionData[LevelIndex][valueT2][valueE22];
704 }
705 }
706
707 if (particleDefinition == G4Proton::ProtonDefinition())
708 {
709 // k should be in eV and energy transfer eV also
710 auto t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
711 auto t1 = t2-1;
712 if (energyTransfer <= pVecm[(*t1)].back() && energyTransfer <= pVecm[(*t2)].back() )
713 {
714 auto e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
715 auto e11 = e12-1;
716
717 auto e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
718 auto e21 = e22-1;
719
720 valueT1 =*t1;
721 valueT2 =*t2;
722 valueE21 =*e21;
723 valueE22 =*e22;
724 valueE12 =*e12;
725 valueE11 =*e11;
726
727 xs11 = pDiffCrossSectionData[LevelIndex][valueT1][valueE11];
728 xs12 = pDiffCrossSectionData[LevelIndex][valueT1][valueE12];
729 xs21 = pDiffCrossSectionData[LevelIndex][valueT2][valueE21];
730 xs22 = pDiffCrossSectionData[LevelIndex][valueT2][valueE22];
731 }
732 }
733
734 sigma = QuadInterpolator( valueE11, valueE12,
735 valueE21, valueE22,
736 xs11, xs12,
737 xs21, xs22,
738 valueT1, valueT2,
739 k, energyTransfer);
740 }
741 return sigma;
742}
743
744//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
745
746G4double G4MicroElecInelasticModel::Interpolate(G4double e1,
747 G4double e2,
748 G4double e,
749 G4double xs1,
750 G4double xs2)
751{
752 G4double value = 0.;
753
754 // Log-log interpolation by default
755 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
756 && !fasterCode)
757 {
758 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
759 G4double b = std::log10(xs2) - a*std::log10(e2);
760 G4double sigma = a*std::log10(e) + b;
761 value = (std::pow(10.,sigma));
762
763 }
764
765 // Switch to log-lin interpolation for faster code
766 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
767 {
768 G4double d1 = std::log10(xs1);
769 G4double d2 = std::log10(xs2);
770 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
771 }
772
773 // Switch to lin-lin interpolation for faster code
774 // in case one of xs1 or xs2 (=cum proba) value is zero
775 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0)) // && fasterCode)
776 {
777 G4double d1 = xs1;
778 G4double d2 = xs2;
779 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
780 }
781
782 return value;
783}
784
785//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
786
787G4double G4MicroElecInelasticModel::QuadInterpolator(G4double e11, G4double e12,
788 G4double e21, G4double e22,
789 G4double xs11, G4double xs12,
790 G4double xs21, G4double xs22,
791 G4double t1, G4double t2,
792 G4double t, G4double e)
793{
794 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
795 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
796 G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
797 return value;
798}
799
800//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
801
802G4int G4MicroElecInelasticModel::RandomSelect(G4double k, const G4String& particle )
803{
804 G4int level = 0;
805
806 auto pos = tableData.find(particle);
807 if (pos != tableData.cend())
808 {
809 G4MicroElecCrossSectionDataSet* table = pos->second;
810
811 if (table != nullptr)
812 {
813 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
814 const G4int n = (G4int)table->NumberOfComponents();
815 G4int i(n);
816 G4double value = 0.;
817
818 while (i>0)
819 {
820 --i;
821 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
822 value += valuesBuffer[i];
823 }
824
825 value *= G4UniformRand();
826
827 i = n;
828
829 while (i > 0)
830 {
831 --i;
832
833 if (valuesBuffer[i] > value)
834 {
835 delete[] valuesBuffer;
836 return i;
837 }
838 value -= valuesBuffer[i];
839 }
840
841 if (valuesBuffer) delete[] valuesBuffer;
842
843 }
844 }
845 else
846 {
847 G4Exception("G4MicroElecInelasticModel::RandomSelect","em0002",FatalException,"Model not applicable to particle type.");
848 }
849
850 return level;
851}
852
853//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
854
855G4double G4MicroElecInelasticModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
856 G4double k,
857 G4int shell)
858{
859
860 G4double secondaryElectronKineticEnergy = 0.;
861 G4double random = G4UniformRand();
862 secondaryElectronKineticEnergy = TransferedEnergy(particleDefinition,
863 k / eV,
864 shell,
865 random) * eV
866 - SiStructure.Energy(shell);
867
868 if (secondaryElectronKineticEnergy < 0.)
869 return 0.;
870 //
871 return secondaryElectronKineticEnergy;
872}
873
874//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
875
877 G4double k,
878 G4int ionizationLevelIndex,
879 G4double random)
880{
881 G4double nrj = 0.;
882
883 G4double valueK1 = 0;
884 G4double valueK2 = 0;
885 G4double valuePROB21 = 0;
886 G4double valuePROB22 = 0;
887 G4double valuePROB12 = 0;
888 G4double valuePROB11 = 0;
889
890 G4double nrjTransf11 = 0;
891 G4double nrjTransf12 = 0;
892 G4double nrjTransf21 = 0;
893 G4double nrjTransf22 = 0;
894
895 G4double maximumEnergyTransfer1 = 0;
896 G4double maximumEnergyTransfer2 = 0;
897 G4double maximumEnergyTransferP = 4.* (electron_mass_c2 / proton_mass_c2) * k;
898 G4double bindingEnergy = SiStructure.Energy(ionizationLevelIndex)*1e6;
899
900 if (particleDefinition == G4Electron::ElectronDefinition())
901 {
902 // k should be in eV
903 auto k2 = std::upper_bound(eTdummyVec.begin(),
904 eTdummyVec.end(),
905 k);
906 auto k1 = k2 - 1;
907
908 /*
909 G4cout << "----> k=" << k
910 << " " << *k1
911 << " " << *k2
912 << " " << random
913 << " " << ionizationLevelIndex
914 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
915 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
916 << G4endl;
917 */
918
919 // SI : the following condition avoids situations where random >last vector element
920 if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
921 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
922 {
923 auto prob12 =
924 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
925 eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
926 random);
927 auto prob11 = prob12 - 1;
928 auto prob22 =
929 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
930 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
931 random);
932 auto prob21 = prob22 - 1;
933
934 valueK1 = *k1;
935 valueK2 = *k2;
936 valuePROB21 = *prob21;
937 valuePROB22 = *prob22;
938 valuePROB12 = *prob12;
939 valuePROB11 = *prob11;
940
941 // The following condition avoid getting transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
942 if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
943 else nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
944 if(valuePROB12 == 1)
945 {
946 if ((valueK1+bindingEnergy)/2. > valueK1) maximumEnergyTransfer1=valueK1;
947 else maximumEnergyTransfer1 = (valueK1+bindingEnergy)/2.;
948
949 nrjTransf12 = maximumEnergyTransfer1;
950 }
951 else nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
952
953 if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
954 else nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
955 if(valuePROB22 == 1)
956 {
957 if ((valueK2+bindingEnergy)/2. > valueK2) maximumEnergyTransfer2=valueK2;
958 else maximumEnergyTransfer2 = (valueK2+bindingEnergy)/2.;
959
960 nrjTransf22 = maximumEnergyTransfer2;
961 }
962 else
963 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
964 }
965 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
966 if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
967 {
968 auto prob22 =
969 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
970 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
971 random);
972
973 auto prob21 = prob22 - 1;
974
975 valueK1 = *k1;
976 valueK2 = *k2;
977 valuePROB21 = *prob21;
978 valuePROB22 = *prob22;
979
980 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
981 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
982
983 G4double interpolatedvalue2 = Interpolate(valuePROB21,
984 valuePROB22,
985 random,
986 nrjTransf21,
987 nrjTransf22);
988
989 // zeros are explicitly set
990 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
991 return value;
992 }
993 }
994 //
995 else if (particleDefinition == G4Proton::ProtonDefinition())
996 {
997 // k should be in eV
998 auto k2 = std::upper_bound(pTdummyVec.begin(),
999 pTdummyVec.end(),
1000 k);
1001
1002 auto k1 = k2 - 1;
1003
1004 /*
1005 G4cout << "----> k=" << k
1006 << " " << *k1
1007 << " " << *k2
1008 << " " << random
1009 << " " << ionizationLevelIndex
1010 << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1011 << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1012 << G4endl;
1013 */
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
1025 auto prob11 = prob12 - 1;
1026
1027 auto prob22 =
1028 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1029 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1030 random);
1031
1032 auto prob21 = prob22 - 1;
1033
1034 valueK1 = *k1;
1035 valueK2 = *k2;
1036 valuePROB21 = *prob21;
1037 valuePROB22 = *prob22;
1038 valuePROB12 = *prob12;
1039 valuePROB11 = *prob11;
1040
1041 // The following condition avoid getting transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
1042 if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
1043 else nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1044 if(valuePROB12 == 1) nrjTransf12 = maximumEnergyTransferP;
1045 else nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1046 if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
1047 else nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1048 if(valuePROB22 == 1) nrjTransf22 = maximumEnergyTransferP;
1049 else nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1050
1051 }
1052
1053 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1054 if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1055 {
1056 auto prob22 =
1057 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1058 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1059 random);
1060
1061 auto prob21 = prob22 - 1;
1062
1063 valueK1 = *k1;
1064 valueK2 = *k2;
1065 valuePROB21 = *prob21;
1066 valuePROB22 = *prob22;
1067
1068 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1069 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1070
1071 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1072 valuePROB22,
1073 random,
1074 nrjTransf21,
1075 nrjTransf22);
1076
1077 // zeros are explicitly set
1078 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1079
1080 return value;
1081 }
1082 }
1083 // End electron and proton cases
1084
1085 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1086 * 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
G4AtomicShellEnumerator
const char * G4FindDataDir(const char *)
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
const G4int Z[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
double z() const
Hep3Vector unit() const
double x() const
double y() const
void set(double x, double y, double z)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:88
static G4Electron * Electron()
Definition: G4Electron.cc:93
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:228
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:204
G4bool LoadData(const G4String &argFileName) override
G4double FindValue(G4double e, G4int componentId=0) const override
size_t NumberOfComponents(void) const override
const G4VEMDataSet * GetComponent(G4int componentId) const override
G4double TransferedEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4ParticleChangeForGamma * fParticleChangeForGamma
G4MicroElecInelasticModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="MicroElecInelasticModel")
G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax) override
G4double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
G4double Energy(G4int level)
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
static G4NistManager * Instance()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:87
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:746
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:600
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:124
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:753
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:802
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:607
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)