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;
101 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
120 if (verboseLevel > 3)
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");
164 eTdummyVec.push_back(0.);
165 while(!eDiffCrossSection.eof())
169 eDiffCrossSection>>tDummy>>eDummy;
170 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
171 for (
int j=0; j<6; j++)
173 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
176 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
178 eVecm[tDummy].push_back(eDummy);
188 tableFile[
proton] = fileProton;
198 tableData[
proton] = tableP;
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");
211 pTdummyVec.push_back(0.);
212 while(!pDiffCrossSection.eof())
216 pDiffCrossSection>>tDummy>>eDummy;
217 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
218 for (
int j=0; j<6; j++)
220 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
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
258 if (isInitialised) {
return; }
260 isInitialised =
true;
272 if (verboseLevel > 3)
273 G4cout <<
"Calling CrossSectionPerVolume() of G4MuElecInelasticModel" <<
G4endl;
305 if (verboseLevel > 3)
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" ;
313 if (verboseLevel > 3)
314 G4cout <<
"After scaling : " << G4endl
315 <<
"Particle : " << nameLocal <<
", Ekin (eV) = " << ekin/
eV <<
G4endl ;
321 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
322 pos1 = lowEnergyLimit.find(nameLocal);
323 if (pos1 != lowEnergyLimit.end())
325 lowLim = pos1->second;
328 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
329 pos2 = highEnergyLimit.find(nameLocal);
330 if (pos2 != highEnergyLimit.end())
332 highLim = pos2->second;
335 if (ekin >= lowLim && ekin < highLim)
337 std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
338 pos = tableData.find(nameLocal);
340 if (pos != tableData.end())
350 G4Exception(
"G4MuElecInelasticModel::CrossSectionPerVolume",
"em0002",
FatalException,
"Model not applicable to particle type.");
362 if (verboseLevel > 3)
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;
384 if (verboseLevel > 3)
385 G4cout <<
"Calling SampleSecondaries() of G4MuElecInelasticModel" <<
G4endl;
395 G4String nameLocal2 = particleName ;
402 nameLocal2 =
"proton" ;
405 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
406 pos1 = lowEnergyLimit.find(nameLocal2);
408 if (pos1 != lowEnergyLimit.end())
410 lowLim = pos1->second;
413 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
414 pos2 = highEnergyLimit.find(nameLocal2);
416 if (pos2 != highEnergyLimit.end())
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);
428 G4int Shell = RandomSelect(k,nameLocal2);
430 if (verboseLevel > 3)
433 G4cout <<
"Shell: " << Shell <<
", energy: " << bindingEnergy/
eV <<
G4endl;
438 G4int secNumberInit = 0;
439 G4int secNumberFinal = 0;
441 if(fAtomDeexcitation && Shell > 2) {
455 secNumberInit = fvect->size();
457 secNumberFinal = fvect->size();
460 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef,k,Shell);
462 if (verboseLevel > 3)
465 G4cout <<
"Shell: " << Shell <<
" Kin. energy (eV)=" << k/
eV
466 <<
" Sec. energy (eV)=" << secondaryKinetic/
eV <<
G4endl;
471 RandomizeEjectedElectronDirection(PartDef, k, secondaryKinetic, cosTheta, phi);
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);
522 if ((k+SiStructure.
Energy(shell))/2. > k) maximumEnergyTransfer=k;
523 else maximumEnergyTransfer = (k+SiStructure.
Energy(shell))/2.;
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.;
546 secondaryElectronKineticEnergy =
G4UniformRand() * (maximumEnergyTransfer-SiStructure.
Energy(shell));
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.;
576 secondaryElectronKineticEnergy =
G4UniformRand() * (maximumEnergyTransfer-SiStructure.
Energy(shell));
580 return secondaryElectronKineticEnergy;
588 void G4MuElecInelasticModel::RandomizeEjectedElectronDirection(
G4ParticleDefinition* particleDefinition,
598 cosTheta = std::sqrt(1.-sin2O);
605 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
612 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
625 if (energyTransfer >= SiStructure.
Energy(LevelIndex))
643 std::vector<double>::iterator
t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
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;
661 xs11 = eDiffCrossSectionData[LevelIndex][valueT1][valueE11];
662 xs12 = eDiffCrossSectionData[LevelIndex][valueT1][valueE12];
663 xs21 = eDiffCrossSectionData[LevelIndex][valueT2][valueE21];
664 xs22 = eDiffCrossSectionData[LevelIndex][valueT2][valueE22];
672 std::vector<double>::iterator
t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
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;
689 xs11 = pDiffCrossSectionData[LevelIndex][valueT1][valueE11];
690 xs12 = pDiffCrossSectionData[LevelIndex][valueT1][valueE12];
691 xs21 = pDiffCrossSectionData[LevelIndex][valueT2][valueE21];
692 xs22 = pDiffCrossSectionData[LevelIndex][valueT2][valueE22];
696 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
699 sigma = QuadInterpolator( valueE11, valueE12,
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);
736 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
737 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
738 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
748 std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
749 pos = tableData.find(particle);
751 if (pos != tableData.end())
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 &)
void set(double x, double y, double z)
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 GetKineticEnergy() const
G4double HighEnergyLimit() const
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static G4Proton * ProtonDefinition()
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
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
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
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
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
virtual ~G4MuElecInelasticModel()
const XML_Char int const XML_Char * value
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
G4double bindingEnergy(G4int A, G4int Z)
G4ParticleChangeForGamma * GetParticleChangeForGamma()