48 :XSTableElectron(0),XSTablePositron(0),
49 theDeltaTable(0),energyGrid(0)
58 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
60 theDeltaTable =
new std::map<const G4Material*,G4PhysicsFreeVector*>;
95 std::map<const G4Material*,G4PhysicsFreeVector*>::iterator k;
110 G4cout <<
"G4PenelopeIonisationXSHandler. Tables have been cleared"
125 G4Exception(
"G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
134 G4Exception(
"G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
136 "The Cross Section Table for e- was not initialized correctly!");
139 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
150 G4Exception(
"G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
152 "The Cross Section Table for e+ was not initialized correctly!");
155 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
172 G4Exception(
"G4PenelopeIonisationXSHandler::BuildXSTable()",
183 G4cout <<
"G4PenelopeIonisationXSHandler: going to build cross section table " <<
G4endl;
188 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
208 size_t numberOfOscillators = theTable->size();
213 ed <<
"Energy Grid looks not initialized" <<
G4endl;
215 G4Exception(
"G4PenelopeIonisationXSHandler::BuildXSTable()",
222 for (
size_t bin=0;bin<
nBins;bin++)
229 for (
size_t iosc=0;iosc<numberOfOscillators;iosc++)
243 ed <<
"Problem in calculating the shell XS for shell # "
245 G4Exception(
"G4PenelopeIonisationXSHandler::BuildXSTable()",
250 if (tempStorage->size() != 6)
253 ed <<
"Problem in calculating the shell XS " <<
G4endl;
254 ed <<
"Result has dimension " << tempStorage->size() <<
" instead of 6" <<
G4endl;
255 G4Exception(
"G4PenelopeIonisationXSHandler::BuildXSTable()",
260 XH0 += stre*(*tempStorage)[0];
261 XH1 += stre*(*tempStorage)[1];
262 XH2 += stre*(*tempStorage)[2];
263 XS0 += stre*(*tempStorage)[3];
264 XS1 += stre*(*tempStorage)[4];
265 XS2 += stre*(*tempStorage)[5];
298 G4Exception(
"G4PenelopeIonisationXSHandler::GetDensityCorrection()",
300 "Delta Table not initialized. Was Initialise() run?");
305 G4cout <<
"G4PenelopeIonisationXSHandler::GetDensityCorrection()" <<
G4endl;
314 result = vec->
Value(logene);
319 ed <<
"Unable to build table for " << mat->
GetName() <<
G4endl;
320 G4Exception(
"G4PenelopeIonisationXSHandler::GetDensityCorrection()",
334 size_t numberOfOscillators = theTable->size();
339 ed <<
"Energy Grid for Delta table looks not initialized" <<
G4endl;
341 G4Exception(
"G4PenelopeIonisationXSHandler::BuildDeltaTable()",
348 for (
size_t bin=0;bin<
nBins;bin++)
357 G4double TST = totalZ/(gamSq*plasmaSq);
362 for (
size_t i=0;i<numberOfOscillators;i++)
381 for (
size_t i=0;i<numberOfOscillators;i++)
397 wl2 = 0.5*(wl2l+wl2u);
399 for (
size_t i=0;i<numberOfOscillators;i++)
409 if ((wl2u-wl2l)>1e-12*wl2)
415 for (
size_t i=0;i<numberOfOscillators;i++)
420 std::log(1.0+(wl2/(wri*wri)));
422 delta = (delta/totalZ)-wl2/(gamSq*plasmaSq);
425 theVector->
PutValue(bin,std::log(energy),delta);
446 for (
size_t i=0;i<6;i++)
447 result->push_back(0.);
451 if (energy < ionEnergy)
458 G4double gamma = 1.0+energy/electron_mass_c2;
460 G4double beta = (gammaSq-1.0)/gammaSq;
461 G4double pielr2 =
pi*classic_electr_radius*classic_electr_radius;
462 G4double constant = pielr2*2.0*electron_mass_c2/beta;
463 G4double XHDT0 = std::log(gammaSq)-beta;
465 G4double cpSq = energy*(energy+2.0*electron_mass_c2);
467 G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
476 G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
481 if (resEne > 1e-6*energy)
482 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
485 QM = resEne*resEne/(beta*2.0*electron_mass_c2);
486 QM = QM*(1.0-0.5*QM/electron_mass_c2);
490 SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
517 if (wl < wu-(1e-5*
eV))
519 H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) +
520 (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
521 amol*(wu-wl)/(ee*ee);
522 H1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) +
523 (2.0-amol)*std::log((ee-wu)/(ee-wl)) +
524 amol*(wu*wu-wl*wl)/(2.0*ee*ee);
525 H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) -
526 (wl*(2.0*ee-wl)/(ee-wl)) +
527 (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) +
528 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
533 if (wl > wu-(1e-5*
eV))
535 (*result)[0] = constant*H0;
536 (*result)[1] = constant*H1;
537 (*result)[2] = constant*H2;
538 (*result)[3] = constant*S0;
539 (*result)[4] = constant*S1;
540 (*result)[5] = constant*S2;
544 S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) +
545 (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
546 amol*(wu-wl)/(ee*ee);
547 S1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) +
548 (2.0-amol)*std::log((ee-wu)/(ee-wl)) +
549 amol*(wu*wu-wl*wl)/(2.0*ee*ee);
550 S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) -
551 (wl*(2.0*ee-wl)/(ee-wl)) +
552 (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) +
553 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
555 (*result)[0] = constant*H0;
556 (*result)[1] = constant*H1;
557 (*result)[2] = constant*H2;
558 (*result)[3] = constant*S0;
559 (*result)[4] = constant*S1;
560 (*result)[5] = constant*S2;
578 for (
size_t i=0;i<6;i++)
579 result->push_back(0.);
583 if (energy < ionEnergy)
590 G4double gamma = 1.0+energy/electron_mass_c2;
592 G4double beta = (gammaSq-1.0)/gammaSq;
593 G4double pielr2 =
pi*classic_electr_radius*classic_electr_radius;
594 G4double constant = pielr2*2.0*electron_mass_c2/beta;
595 G4double XHDT0 = std::log(gammaSq)-beta;
597 G4double cpSq = energy*(energy+2.0*electron_mass_c2);
599 G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
600 G4double g12 = (gamma+1.0)*(gamma+1.0);
602 G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq-1.0);
604 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g12;
605 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/g12;
614 G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
619 if (resEne > 1e-6*energy)
620 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
623 QM = resEne*resEne/(beta*2.0*electron_mass_c2);
624 QM = QM*(1.0-0.5*QM/electron_mass_c2);
628 SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
656 if (wl < wu-(1e-5*
eV))
660 H0 += (1.0/wl) - (1.0/wu)- bha1*std::log(wu/wl)/energy
661 + bha2*(wu-wl)/energySq
662 - bha3*(wuSq-wlSq)/(2.0*energySq*
energy)
663 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
664 H1 += std::log(wu/wl) - bha1*(wu-wl)/energy
665 + bha2*(wuSq-wlSq)/(2.0*energySq)
666 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*
energy)
667 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
668 H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
669 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
670 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
671 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
676 if (wl > wu-(1e-5*
eV))
678 (*result)[0] = constant*H0;
679 (*result)[1] = constant*H1;
680 (*result)[2] = constant*H2;
681 (*result)[3] = constant*S0;
682 (*result)[4] = constant*S1;
683 (*result)[5] = constant*S2;
690 S0 += (1.0/wl) - (1.0/wu) - bha1*std::log(wu/wl)/energy
691 + bha2*(wu-wl)/energySq
692 - bha3*(wuSq-wlSq)/(2.0*energySq*
energy)
693 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
695 S1 += std::log(wu/wl) - bha1*(wu-wl)/energy
696 + bha2*(wuSq-wlSq)/(2.0*energySq)
697 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*
energy)
698 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
700 S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
701 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
702 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
703 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
705 (*result)[0] = constant*H0;
706 (*result)[1] = constant*H1;
707 (*result)[2] = constant*H2;
708 (*result)[3] = constant*S0;
709 (*result)[4] = constant*S1;
710 (*result)[5] = constant*S2;
std::map< std::pair< const G4Material *, G4double >, G4PenelopeCrossSection * > * XSTableElectron
G4PenelopeOscillatorTable * GetOscillatorTableIonisation(const G4Material *)
G4double GetIonisationEnergy()
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
std::ostringstream G4ExceptionDescription
G4double GetPlasmaEnergySquared(const G4Material *)
Returns the squared plasma energy.
G4DataVector * ComputeShellCrossSectionsElectron(G4PenelopeOscillator *, G4double energy, G4double cut, G4double delta)
std::map< std::pair< const G4Material *, G4double >, G4PenelopeCrossSection * > * XSTablePositron
const G4String & GetName() const
void AddShellCrossSectionPoint(size_t binNumber, size_t shellID, G4double energy, G4double xs)
G4PhysicsLogVector * energyGrid
const G4PenelopeCrossSection * GetCrossSectionTableForCouple(const G4ParticleDefinition *, const G4Material *, const G4double cut) const
Returns the table of cross sections for the given particle, given material and given cut as a G4Penel...
size_t GetVectorLength() const
G4double GetCutoffRecoilResonantEnergy()
G4double GetLowEdgeEnergy(size_t binNumber) const
const G4String & GetParticleName() const
G4double GetOscillatorStrength()
void AddCrossSectionPoint(size_t binNumber, G4double energy, G4double XH0, G4double XH1, G4double XH2, G4double XS0, G4double XS1, G4double XS2)
Public interface for the master thread.
G4DataVector * ComputeShellCrossSectionsPositron(G4PenelopeOscillator *, G4double energy, G4double cut, G4double delta)
std::map< const G4Material *, G4PhysicsFreeVector * > * theDeltaTable
G4GLOB_DLL std::ostream G4cout
void BuildDeltaTable(const G4Material *)
G4PenelopeIonisationXSHandler(size_t nBins=200)
Constructor. nBins is the number of intervals in the energy grid. By default the energy grid goes fro...
std::ostream & tab(std::ostream &)
static G4PenelopeOscillatorManager * GetOscillatorManager()
G4double GetResonanceEnergy() const
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4PenelopeOscillatorManager * oscManager
static G4Positron * Positron()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
G4double GetDensityCorrection(const G4Material *, const G4double energy) const
Returns the density coeection for the material at the given energy.
static G4Electron * Electron()
void NormalizeShellCrossSections()
void BuildXSTable(const G4Material *, G4double cut, const G4ParticleDefinition *, G4bool isMaster=true)
This can be inkoved only by the master.
G4double GetTotalZ(const G4Material *)
These are cumulative for the molecule Returns the total Z for the molecule.
virtual ~G4PenelopeIonisationXSHandler()
Destructor. Clean all tables.