62 :
G4VEmModel(nam),fAtomDeexcitation(0),isInitialised(false)
66 G4cout <<
"*******************************************************************************" <<
G4endl;
67 G4cout <<
"*******************************************************************************" <<
G4endl;
68 G4cout <<
" The name of the class G4MuElecInelasticModel is changed to G4MicroElecInelasticModel. " <<
G4endl;
69 G4cout <<
" The obsolete class will be REMOVED with the next release of Geant4. " <<
G4endl;
70 G4cout <<
"*******************************************************************************" <<
G4endl;
71 G4cout <<
"*******************************************************************************" <<
G4endl;
86 G4cout <<
"MuElec inelastic model is constructed " <<
G4endl;
100 std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
121 G4cout <<
"Calling G4MuElecInelasticModel::Initialise()" <<
G4endl;
125 G4String fileElectron(
"microelec/sigma_inelastic_e_Si");
126 G4String fileProton(
"microelec/sigma_inelastic_p_Si");
136 char *path = getenv(
"G4LEDATA");
155 std::ostringstream eFullFileName;
156 eFullFileName << path <<
"/microelec/sigmadiff_inelastic_e_Si.dat";
157 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
159 if (!eDiffCrossSection)
161 G4Exception(
"G4MuElecInelasticModel::Initialise",
"em0003",
FatalException,
"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
165 while(!eDiffCrossSection.eof())
169 eDiffCrossSection>>tDummy>>eDummy;
171 for (
int j=0; j<6; j++)
176 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
178 eVecm[tDummy].push_back(eDummy);
202 std::ostringstream pFullFileName;
203 pFullFileName << path <<
"/microelec/sigmadiff_inelastic_p_Si.dat";
204 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
206 if (!pDiffCrossSection)
208 G4Exception(
"G4MuElecInelasticModel::Initialise",
"em0003",
FatalException,
"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
212 while(!pDiffCrossSection.eof())
216 pDiffCrossSection>>tDummy>>eDummy;
218 for (
int j=0; j<6; j++)
223 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
225 pVecm[tDummy].push_back(eDummy);
230 if (particle==electronDef)
236 if (particle==protonDef)
244 G4cout <<
"MuElec Inelastic model is initialized " << G4endl
249 <<
" with mass (amu) " << particle->
GetPDGMass()/proton_mass_c2
273 G4cout <<
"Calling CrossSectionPerVolume() of G4MuElecInelasticModel" <<
G4endl;
299 if (Mion_c2 > proton_mass_c2)
306 G4cout <<
"Before scaling : " << G4endl
307 <<
"Particle : " << nameLocal <<
", mass : " << Mion_c2/proton_mass_c2 <<
"*mp, charge " << Zeff
308 <<
", Ekin (eV) = " << ekin/
eV <<
G4endl ;
310 ekin *= proton_mass_c2/Mion_c2 ;
311 nameLocal =
"proton" ;
314 G4cout <<
"After scaling : " << G4endl
315 <<
"Particle : " << nameLocal <<
", Ekin (eV) = " << ekin/
eV <<
G4endl ;
321 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
325 lowLim = pos1->second;
328 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
332 highLim = pos2->second;
335 if (ekin >= lowLim && ekin < highLim)
337 std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
350 G4Exception(
"G4MuElecInelasticModel::CrossSectionPerVolume",
"em0002",
FatalException,
"Model not applicable to particle type.");
365 G4cout <<
" - Cross section per Si atom (cm^2)=" << sigma*Zeff2/
cm2 <<
G4endl;
366 G4cout <<
" - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./
cm) << G4endl;
370 return sigma*density*Zeff2;
385 G4cout <<
"Calling SampleSecondaries() of G4MuElecInelasticModel" <<
G4endl;
395 G4String nameLocal2 = particleName ;
398 if (particleMass > proton_mass_c2)
400 k *= proton_mass_c2/particleMass ;
402 nameLocal2 =
"proton" ;
405 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
410 lowLim = pos1->second;
413 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
418 highLim = pos2->second;
421 if (k >= lowLim && k < highLim)
424 G4double totalEnergy = ekin + particleMass;
425 G4double pSquare = ekin * (totalEnergy + particleMass);
426 G4double totalMomentum = std::sqrt(pSquare);
433 G4cout <<
"Shell: " << Shell <<
", energy: " << bindingEnergy/
eV <<
G4endl;
438 G4int secNumberInit = 0;
439 G4int secNumberFinal = 0;
455 secNumberInit = fvect->size();
457 secNumberFinal = fvect->size();
465 G4cout <<
"Shell: " << Shell <<
" Kin. energy (eV)=" << k/
eV
466 <<
" Sec. energy (eV)=" << secondaryKinetic/
eV <<
G4endl;
473 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
474 G4double dirX = sinTheta*std::cos(phi);
475 G4double dirY = sinTheta*std::sin(phi);
478 deltaDirection.rotateUz(primaryDirection);
482 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
484 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
485 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
486 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
487 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
488 finalPx /= finalMomentum;
489 finalPy /= finalMomentum;
490 finalPz /= finalMomentum;
493 direction.set(finalPx,finalPy,finalPz);
501 for (
G4int j=secNumberInit; j < secNumberFinal; j++) {
502 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
508 fvect->push_back(dp);
528 G4double maxEnergy = maximumEnergyTransfer;
529 G4int nEnergySteps = 100;
532 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
533 G4int step(nEnergySteps);
538 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
543 G4double secondaryElectronKineticEnergy=0.;
550 return secondaryElectronKineticEnergy;
556 G4double maximumEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
560 G4double maxEnergy = maximumEnergyTransfer;
561 G4int nEnergySteps = 100;
564 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
565 G4int step(nEnergySteps);
570 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
573 G4double secondaryElectronKineticEnergy = 0.;
580 return secondaryElectronKineticEnergy;
597 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
598 cosTheta = std::sqrt(1.-sin2O);
603 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
605 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
610 G4double maxSecKinetic = 4.* (electron_mass_c2 / particleDefinition->
GetPDGMass()) * k;
612 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
644 std::vector<double>::iterator t1 = t2-1;
646 if (energyTransfer <=
eVecm[(*t1)].back() && energyTransfer <=
eVecm[(*t2)].back() )
648 std::vector<double>::iterator e12 = std::upper_bound(
eVecm[(*t1)].begin(),
eVecm[(*t1)].end(), energyTransfer);
649 std::vector<double>::iterator e11 = e12-1;
651 std::vector<double>::iterator e22 = std::upper_bound(
eVecm[(*t2)].begin(),
eVecm[(*t2)].end(), energyTransfer);
652 std::vector<double>::iterator e21 = e22-1;
673 std::vector<double>::iterator t1 = t2-1;
674 if (energyTransfer <=
pVecm[(*t1)].back() && energyTransfer <=
pVecm[(*t2)].back() )
676 std::vector<double>::iterator e12 = std::upper_bound(
pVecm[(*t1)].begin(),
pVecm[(*t1)].end(), energyTransfer);
677 std::vector<double>::iterator e11 = e12-1;
679 std::vector<double>::iterator e22 = std::upper_bound(
pVecm[(*t2)].begin(),
pVecm[(*t2)].end(), energyTransfer);
680 std::vector<double>::iterator e21 = e22-1;
696 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
720 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
721 G4double b = std::log10(xs2) - a*std::log10(e2);
722 G4double sigma = a*std::log10(e) + b;
723 G4double value = (std::pow(10.,sigma));
748 std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
766 value += valuesBuffer[i];
777 if (valuesBuffer[i] > value)
779 delete[] valuesBuffer;
782 value -= valuesBuffer[i];
785 if (valuesBuffer)
delete[] valuesBuffer;
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4Electron * ElectronDefinition()
G4double LowEnergyLimit() const
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
virtual G4double FindValue(G4double e, G4int componentId=0) const
static G4LossTableManager * Instance()
virtual size_t NumberOfComponents(void) const
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static G4Proton * ProtonDefinition()
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
std::vector< double > pTdummyVec
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
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)
static G4NistManager * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
TriDimensionMap eDiffCrossSectionData[7]
G4GLOB_DLL std::ostream G4cout
TriDimensionMap pDiffCrossSectionData[7]
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
G4MuElecSiStructure SiStructure
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
G4double GetTotNbOfAtomsPerVolume() const
G4double GetPDGMass() const
void RandomizeEjectedElectronDirection(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4double outgoingParticleEnergy, G4double &cosTheta, G4double &phi)
G4VAtomDeexcitation * fAtomDeexcitation
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4MuElecInelasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MuElecInelasticModel")
const G4Material * GetBaseMaterial() const
std::vector< double > eTdummyVec
virtual ~G4MuElecInelasticModel()
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double Energy(G4int level)
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)
static const G4double pos
G4int RandomSelect(G4double energy, const G4String &particle)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit