64 :
G4VEmModel(nam),fAtomDeexcitation(0),isInitialised(false)
78 G4cout <<
"MicroElec inelastic model is constructed " <<
G4endl;
95 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
116 G4cout <<
"Calling G4MicroElecInelasticModel::Initialise()" <<
G4endl;
120 G4String fileElectron(
"microelec/sigma_inelastic_e_Si");
121 G4String fileProton(
"microelec/sigma_inelastic_p_Si");
131 char *path = getenv(
"G4LEDATA");
150 std::ostringstream eFullFileName;
151 eFullFileName << path <<
"/microelec/sigmadiff_inelastic_e_Si.dat";
152 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
154 if (!eDiffCrossSection)
156 G4Exception(
"G4MicroElecInelasticModel::Initialise",
"em0003",
FatalException,
"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
177 while(!eDiffCrossSection.eof())
181 eDiffCrossSection>>tDummy>>eDummy;
183 for (
int j=0; j<6; j++)
188 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
190 eVecm[tDummy].push_back(eDummy);
214 std::ostringstream pFullFileName;
215 pFullFileName << path <<
"/microelec/sigmadiff_inelastic_p_Si.dat";
216 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
218 if (!pDiffCrossSection)
220 G4Exception(
"G4MicroElecInelasticModel::Initialise",
"em0003",
FatalException,
"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
224 while(!pDiffCrossSection.eof())
228 pDiffCrossSection>>tDummy>>eDummy;
230 for (
int j=0; j<6; j++)
235 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
237 pVecm[tDummy].push_back(eDummy);
242 if (particle==electronDef)
248 if (particle==protonDef)
256 G4cout <<
"MicroElec Inelastic model is initialized " << G4endl
261 <<
" with mass (amu) " << particle->
GetPDGMass()/proton_mass_c2
285 G4cout <<
"Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" <<
G4endl;
311 if (Mion_c2 > proton_mass_c2)
318 G4cout <<
"Before scaling : " << G4endl
319 <<
"Particle : " << nameLocal <<
", mass : " << Mion_c2/proton_mass_c2 <<
"*mp, charge " << Zeff
320 <<
", Ekin (eV) = " << ekin/
eV <<
G4endl ;
322 ekin *= proton_mass_c2/Mion_c2 ;
323 nameLocal =
"proton" ;
326 G4cout <<
"After scaling : " << G4endl
327 <<
"Particle : " << nameLocal <<
", Ekin (eV) = " << ekin/
eV <<
G4endl ;
333 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
337 lowLim = pos1->second;
340 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
344 highLim = pos2->second;
347 if (ekin >= lowLim && ekin < highLim)
349 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
362 G4Exception(
"G4MicroElecInelasticModel::CrossSectionPerVolume",
"em0002",
FatalException,
"Model not applicable to particle type.");
377 G4cout <<
" - Cross section per Si atom (cm^2)=" << sigma*Zeff2/
cm2 <<
G4endl;
378 G4cout <<
" - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./
cm) << G4endl;
382 return sigma*density*Zeff2;
397 G4cout <<
"Calling SampleSecondaries() of G4MicroElecInelasticModel" <<
G4endl;
407 G4String nameLocal2 = particleName ;
410 if (particleMass > proton_mass_c2)
412 k *= proton_mass_c2/particleMass ;
414 nameLocal2 =
"proton" ;
417 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
422 lowLim = pos1->second;
425 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
430 highLim = pos2->second;
433 if (k >= lowLim && k < highLim)
436 G4double totalEnergy = ekin + particleMass;
437 G4double pSquare = ekin * (totalEnergy + particleMass);
438 G4double totalMomentum = std::sqrt(pSquare);
446 G4cout <<
"Shell: " << Shell <<
", energy: " << bindingEnergy/
eV <<
G4endl;
451 G4int secNumberInit = 0;
452 G4int secNumberFinal = 0;
470 secNumberInit = fvect->size();
472 secNumberFinal = fvect->size();
480 G4cout <<
"Shell: " << Shell <<
" Kin. energy (eV)=" << k/
eV
481 <<
" Sec. energy (eV)=" << secondaryKinetic/
eV <<
G4endl;
491 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
493 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
494 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
495 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
496 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
497 finalPx /= finalMomentum;
498 finalPy /= finalMomentum;
499 finalPz /= finalMomentum;
502 direction.set(finalPx,finalPy,finalPz);
510 for (
G4int j=secNumberInit; j < secNumberFinal; j++) {
511 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
517 fvect->push_back(dp);
537 G4double maxEnergy = maximumEnergyTransfer;
538 G4int nEnergySteps = 100;
541 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
542 G4int step(nEnergySteps);
547 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
552 G4double secondaryElectronKineticEnergy=0.;
559 return secondaryElectronKineticEnergy;
565 G4double maximumEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
569 G4double maxEnergy = maximumEnergyTransfer;
570 G4int nEnergySteps = 100;
573 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
574 G4int step(nEnergySteps);
579 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
582 G4double secondaryElectronKineticEnergy = 0.;
589 return secondaryElectronKineticEnergy;
657 std::vector<double>::iterator t1 = t2-1;
659 if (energyTransfer <=
eVecm[(*t1)].back() && energyTransfer <=
eVecm[(*t2)].back() )
661 std::vector<double>::iterator e12 = std::upper_bound(
eVecm[(*t1)].begin(),
eVecm[(*t1)].end(), energyTransfer);
662 std::vector<double>::iterator e11 = e12-1;
664 std::vector<double>::iterator e22 = std::upper_bound(
eVecm[(*t2)].begin(),
eVecm[(*t2)].end(), energyTransfer);
665 std::vector<double>::iterator e21 = e22-1;
686 std::vector<double>::iterator t1 = t2-1;
687 if (energyTransfer <=
pVecm[(*t1)].back() && energyTransfer <=
pVecm[(*t2)].back() )
689 std::vector<double>::iterator e12 = std::upper_bound(
pVecm[(*t1)].begin(),
pVecm[(*t1)].end(), energyTransfer);
690 std::vector<double>::iterator e11 = e12-1;
692 std::vector<double>::iterator e22 = std::upper_bound(
pVecm[(*t2)].begin(),
pVecm[(*t2)].end(), energyTransfer);
693 std::vector<double>::iterator e21 = e22-1;
709 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
733 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
734 G4double b = std::log10(xs2) - a*std::log10(e2);
735 G4double sigma = a*std::log10(e) + b;
736 G4double value = (std::pow(10.,sigma));
761 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
779 value += valuesBuffer[i];
790 if (valuesBuffer[i] > value)
792 delete[] valuesBuffer;
795 value -= valuesBuffer[i];
798 if (valuesBuffer)
delete[] valuesBuffer;
static G4Electron * ElectronDefinition()
std::vector< double > pTdummyVec
G4double LowEnergyLimit() const
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
TriDimensionMap eDiffCrossSectionData[7]
G4double HighEnergyLimit() const
G4MicroElecSiStructure SiStructure
G4MicroElecInelasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MicroElecInelasticModel")
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4Proton * ProtonDefinition()
G4VEmAngularDistribution * GetAngularDistribution()
virtual ~G4MicroElecInelasticModel()
G4ParticleDefinition * GetDefinition() const
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
static G4NistManager * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double Energy(G4int level)
void SetHighEnergyLimit(G4double)
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
std::vector< double > eTdummyVec
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
G4GLOB_DLL std::ostream G4cout
virtual size_t NumberOfComponents(void) const
const G4ThreeVector & GetMomentumDirection() const
G4int RandomSelect(G4double energy, const G4String &particle)
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual G4bool LoadData(const G4String &argFileName)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
TriDimensionMap pDiffCrossSectionData[7]
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
G4double GetTotNbOfAtomsPerVolume() const
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
G4double GetPDGMass() const
void SetAngularDistribution(G4VEmAngularDistribution *)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
const G4Material * GetBaseMaterial() const
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4VAtomDeexcitation * AtomDeexcitation()
void SetLowEnergyLimit(G4double)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
G4double GetPDGCharge() const
G4ThreeVector G4ParticleMomentum
G4double bindingEnergy(G4int A, G4int Z)
G4VAtomDeexcitation * fAtomDeexcitation
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
static const G4double pos
G4ParticleChangeForGamma * GetParticleChangeForGamma()
const G4Material * GetMaterial() const
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)