83 MeanFreePath = (*theMeanFreePathTable)(aMaterial->
GetIndex())->Value(kineticEnergy);
94 if( path ==
nullptr ) {
98 "variable G4LEDATA not defined");
114 std::vector<G4Material*>::const_iterator matite;
115 for( matite = materialTable->begin(); matite != materialTable->end(); matite++ ) {
120 std::string dirName = std::string(path) +
"/lepts/";
121 std::string fnParam = dirName + mateName +
"." + aParticleName +
".param.dat";
122 std::string baseName = std::string(path) +
"/lepts/" + mateName +
"." + aParticleName;
124 if( !bData )
continue;
127 std::string fnIXS = baseName +
".IXS.dat";
129 std::map< G4int, std::vector<G4double> > integralXS =
ReadIXS(fnIXS, aMaterial);
132 if( integralXS.empty() ) {
133 G4cerr <<
" Integral cross sections will be set to 0. for material " << mateName <<
G4endl;
135 ptrVector->PutValue(0,
DBL_MAX);
136 ptrVector->PutValue(1,
DBL_MAX);
138 std::size_t matIdx = aMaterial->
GetIndex();
144 std::map< G4int, std::vector<G4double> >::const_iterator itei;
145 for( itei = integralXS.begin(); itei != integralXS.end(); itei++ ){
146 G4cout <<
GetName() <<
" : " << (*itei).first <<
" INTEGRALXS NDATA " << (*itei).second.size() <<
G4endl;
152 std::string fnDXS = baseName +
".DXS.dat";
153 std::string fnRMT = baseName +
".RMT.dat";
154 std::string fnEloss = baseName +
".Eloss.dat";
155 std::string fnEloss2 = baseName +
".Eloss2.dat";
158 if( !
theDiffXS[aMaterial]->IsFileFound() ) {
162 (
G4String(
"File not found :" + fnDXS).c_str()));
173 (
G4String(
"File not found :" + fnEloss).c_str()));
186 std::size_t matIdx = aMaterial->
GetIndex();
190 LowEdgeEnergy = ptrVector->Energy(ii);
202 G4double eVEnergy = LowEdgeEnergy/CLHEP::eV;
206 if(eVEnergy < integralXS[0][1] ) {
213 if( eVEnergy > integralXS[0][jj]) {
219 aa = integralXS[0][Bin];
220 bb = integralXS[0][Bin+1];
221 crossSection = (integralXS[
theXSType][Bin] + (integralXS[
theXSType][Bin+1]-integralXS[
theXSType][Bin])/(bb-aa)*(eVEnergy-aa) ) * 1.e-16*CLHEP::cm2;
223 if(
verboseLevel >= 3 )
G4cout <<
" crossSection " << crossSection <<
" " <<integralXS[
theXSType][Bin] <<
" + " << (integralXS[
theXSType][Bin+1]-integralXS[
theXSType][Bin]) <<
" / " << (bb-aa) <<
" *" << (eVEnergy-aa) <<
" * " << 1.e-16*CLHEP::cm2 <<
G4endl;;
226 SIGMA = NbOfMoleculesPerVolume * crossSection;
228 <<
" Bin " << Bin <<
" TOTAL " << aa <<
" " << bb
237 ptrVector->PutValue(ii, fValue);
251 x =
theDiffXS[aMaterial]->SampleAngleMT(e, el);
257 G4double Pi = std::sqrt( std::pow( (Ei/27.2/137),2) +2*Ei/27.2);
258 G4double Pd = std::sqrt( std::pow( (Ed/27.2/137),2) +2*Ed/27.2);
265 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*Pi*Pd);
266 if( co > 1. ) co = 1.;
278 G4double sinTeta = std::sqrt(1.0-cosTeta*cosTeta);
280 G4double dirx = sinTeta*std::cos(Phi) , diry = sinTeta*std::sin(Phi) , dirz = cosTeta ;
299 G4double sinTeta = std::sqrt(1.0-cosTeta*cosTeta);
301 G4double dirx = sinTeta*std::cos(Phi) , diry = sinTeta*std::sin(Phi) , dirz = cosTeta ;
314 el =
theElostDistr[aMaterial]->Sample(eMin/CLHEP::eV, eMax/CLHEP::eV)*CLHEP::eV;
317 if(
verboseLevel >= 2 )
G4cout << aMaterial->
GetName() <<
"SampleEnergyLoss/eV " << el/CLHEP::eV <<
" eMax/eV " << eMax/CLHEP::eV <<
" "
327 std::ifstream fin(fnParam);
328 if (!fin.is_open()) {
332 (
G4String(
"File not found: ")+ fnParam).c_str());
338 fin >> IonisPot >> IonisPotInt;
339 if(
verboseLevel >= 1 )
G4cout <<
"Read param (" << fnParam <<
")\t IonisPot: " << IonisPot
340 <<
" IonisPotInt: " << IonisPotInt <<
G4endl;
348 for(
G4int ii = 0; ii < nelem; ++ii ) {
349 MolecularMass += aMaterial->
GetElement(ii)->
GetA()*atomsV[ii]/CLHEP::g;
357 <<
" IonisPotInt: " << IonisPotInt/CLHEP::eV <<
" eV"
358 <<
" MolecularMass " << MolecularMass/(CLHEP::g/CLHEP::mole) <<
" g/mole" <<
G4endl;
366 std::map< G4int, std::vector<G4double> > integralXS;
369 std::ifstream fin(fnIXS);
370 if (!fin.is_open()) {
374 (
G4String(
"File not found: ")+ fnIXS).c_str());
378 G4int nXSdat, nXSsub;
379 fin >> nXSdat >> nXSsub;
381 <<
" nXSsub: " << nXSsub <<
G4endl;
386 for (
G4int ip=0; ip<=nXSsub; ip++) {
387 integralXS[ip].push_back(0.);
389 for (
G4int ie=1; ie<=nXSdat; ie++) {
390 for (
G4int ip=0; ip<=nXSsub; ip++) {
392 integralXS[ip].push_back(xsdat);
const char * G4FindDataDir(const char *)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::vector< G4Material * > G4MaterialTable
G4GLOB_DLL std::ostream G4cerr
G4GLOB_DLL std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
G4double GetDensity() const
const G4Element * GetElement(G4int iel) const
const G4int * GetAtomsVector() const
static std::size_t GetNumberOfMaterials()
std::size_t GetIndex() const
static G4MaterialTable * GetMaterialTable()
std::size_t GetNumberOfElements() const
const G4String & GetName() const
const G4String & GetParticleName() const
void insertAt(std::size_t, G4PhysicsVector *)
void SetHighEnergyLimit(G4double)
void SetLowEnergyLimit(G4double)
const G4String & GetName() const
G4bool ReadParam(G4String fileName, const G4Material *aMaterial)
std::map< const G4Material *, G4LEPTSDiffXS * > theDiffXS
G4ThreeVector SampleNewDirection(const G4Material *aMaterial, G4ThreeVector Dir, G4double e, G4double el)
G4double SampleEnergyLoss(const G4Material *aMaterial, G4double eMin, G4double eMax)
std::map< const G4Material *, G4double > theMolecularMass
~G4VLEPTSModel() override
G4double GetMeanFreePath(const G4Material *mate, const G4ParticleDefinition *aParticle, G4double kineticEnergy)
std::map< const G4Material *, G4int > theNXSsub
G4VLEPTSModel(const G4String &processName)
std::map< const G4Material *, G4LEPTSElossDistr * > theElostDistr
std::map< const G4Material *, G4double > theIonisPotInt
std::map< const G4Material *, G4double > theIonisPot
G4PhysicsTable * theMeanFreePathTable
G4double theHighestEnergyLimit
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
virtual std::map< G4int, std::vector< G4double > > ReadIXS(G4String fileName, const G4Material *aMaterial)
G4double theLowestEnergyLimit
std::map< const G4Material *, G4LEPTSDistribution * > theRMTDistr
void BuildMeanFreePathTable(const G4Material *aMaterial, std::map< G4int, std::vector< G4double > > &integralXS)
std::map< const G4Material *, G4int > theNXSdat
G4double SampleAngle(const G4Material *aMaterial, G4double e, G4double el)