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
273 G4cout <<
"Calling CrossSectionPerVolume() of G4MuElecInelasticModel" <<
G4endl;
306 G4cout <<
"Before scaling : " << G4endl
307 <<
"Particle : " << nameLocal <<
", mass : " << Mion_c2/
proton_mass_c2 <<
"*mp, charge " << Zeff
308 <<
", Ekin (eV) = " << ekin/
eV <<
G4endl ;
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 ;
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);
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;
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;
598 cosTheta = std::sqrt(1.-sin2O);
605 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
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;
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;
void set(double x, double y, double z)
std::vector< double > eTdummyVec
std::vector< ExP01TrackerHit * > a
static const G4double pos
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4VAtomDeexcitation * AtomDeexcitation()
static constexpr double MeV
void SetHighEnergyLimit(G4double)
G4int RandomSelect(G4double energy, const G4String &particle)
static constexpr double keV
const G4ThreeVector & GetMomentumDirection() const
const G4String & GetParticleName() const
G4double GetPDGCharge() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
TriDimensionMap pDiffCrossSectionData[7]
Hep3Vector & rotateUz(const Hep3Vector &)
virtual G4bool LoadData(const G4String &argFileName)
G4double GetPDGMass() const
G4double LowEnergyLimit() const
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
static constexpr double proton_mass_c2
G4MuElecSiStructure SiStructure
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)
virtual G4double FindValue(G4double e, G4int componentId=0) const
void SetDeexcitationFlag(G4bool val)
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
G4ParticleDefinition * GetDefinition() const
static constexpr double cm2
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
const XML_Char int const XML_Char * value
static G4Proton * ProtonDefinition()
TriDimensionMap eDiffCrossSectionData[7]
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
static constexpr double electron_mass_c2
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
G4MuElecInelasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MuElecInelasticModel")
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
static constexpr double twopi
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual ~G4MuElecInelasticModel()
static constexpr double eV
std::vector< double > pTdummyVec
static G4Electron * Electron()
G4double Energy(G4int level)
static G4Electron * ElectronDefinition()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void SetLowEnergyLimit(G4double)
void RandomizeEjectedElectronDirection(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4double outgoingParticleEnergy, G4double &cosTheta, G4double &phi)
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
G4double GetKineticEnergy() const
static constexpr double cm
G4GLOB_DLL std::ostream G4cout
static G4LossTableManager * Instance()
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double GetTotNbOfAtomsPerVolume() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4VAtomDeexcitation * fAtomDeexcitation
G4double bindingEnergy(G4int A, G4int Z)
G4double HighEnergyLimit() const
static constexpr double GeV
virtual size_t NumberOfComponents(void) const
const G4Material * GetBaseMaterial() const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
static G4NistManager * Instance()
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)