G4PenelopeIonisationXSHandler Class Reference

#include <G4PenelopeIonisationXSHandler.hh>

Public Member Functions
void	BuildXSTable (const G4Material *, G4double cut, const G4ParticleDefinition *, G4bool isMaster=true)
	This can be inkoved only by the master. More...
	G4PenelopeIonisationXSHandler (const G4PenelopeIonisationXSHandler &)=delete
	G4PenelopeIonisationXSHandler (size_t nBins=200)
const G4PenelopeCrossSection *	GetCrossSectionTableForCouple (const G4ParticleDefinition *, const G4Material *, const G4double cut) const
G4double	GetDensityCorrection (const G4Material *, const G4double energy) const
	Returns the density coeection for the material at the given energy. More...
G4PenelopeIonisationXSHandler &	operator= (const G4PenelopeIonisationXSHandler &right)=delete
void	SetVerboseLevel (G4int vl)
	Setter for the verbosity level. More...
virtual	~G4PenelopeIonisationXSHandler ()
	Destructor. Clean all tables. More...

Private Member Functions
void	BuildDeltaTable (const G4Material *)
G4DataVector *	ComputeShellCrossSectionsElectron (G4PenelopeOscillator *, G4double energy, G4double cut, G4double delta)
G4DataVector *	ComputeShellCrossSectionsPositron (G4PenelopeOscillator *, G4double energy, G4double cut, G4double delta)

Private Attributes
std::map< const G4Material *, G4PhysicsFreeVector * > *	fDeltaTable
G4PhysicsLogVector *	fEnergyGrid
size_t	fNBins
G4PenelopeOscillatorManager *	fOscManager
G4int	fVerboseLevel
std::map< std::pair< const G4Material *, G4double >, G4PenelopeCrossSection * > *	fXSTableElectron
std::map< std::pair< const G4Material *, G4double >, G4PenelopeCrossSection * > *	fXSTablePositron

Detailed Description

Definition at line 59 of file G4PenelopeIonisationXSHandler.hh.

Constructor & Destructor Documentation

G4PenelopeIonisationXSHandler() [1/2]

G4PenelopeIonisationXSHandler::G4PenelopeIonisationXSHandler ( size_t  nBins = 200 )

explicit

Constructor. nBins is the number of intervals in the energy grid. By default the energy grid goes from 100 eV to 100 GeV.

Definition at line 46 of file G4PenelopeIonisationXSHandler.cc.

  :fXSTableElectron(nullptr),fXSTablePositron(nullptr),
   fDeltaTable(nullptr),fEnergyGrid(nullptr)
{
  fNBins = nb;
  G4double LowEnergyLimit = 100.0*eV;
  G4double HighEnergyLimit = 100.0*GeV;
  fOscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
  fXSTableElectron = new 
    std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
  fXSTablePositron = new 
    std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
 
  fDeltaTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
  fEnergyGrid = new G4PhysicsLogVector(LowEnergyLimit,
                      HighEnergyLimit, 
                      fNBins-1); //one hidden bin is added
  fVerboseLevel = 0;
}

References eV, fDeltaTable, fEnergyGrid, fNBins, fOscManager, fVerboseLevel, fXSTableElectron, fXSTablePositron, G4PenelopeOscillatorManager::GetOscillatorManager(), and GeV.

~G4PenelopeIonisationXSHandler()

G4PenelopeIonisationXSHandler::~G4PenelopeIonisationXSHandler ( )

virtual

Destructor. Clean all tables.

Definition at line 68 of file G4PenelopeIonisationXSHandler.cc.

{
  if (fXSTableElectron)
    {
      for (auto& item : (*fXSTableElectron))
    {     
      //G4PenelopeCrossSection* tab = i->second;
      delete item.second;
    }
      delete fXSTableElectron;
      fXSTableElectron = nullptr;
    }
 
  if (fXSTablePositron)
    {
      for (auto& item : (*fXSTablePositron))
    {
      //G4PenelopeCrossSection* tab = i->second;
      delete item.second;
    }
      delete fXSTablePositron;
      fXSTablePositron = nullptr;
    }
  if (fDeltaTable)
    {      
      for (auto& item : (*fDeltaTable))
    delete item.second;     
      delete fDeltaTable;
      fDeltaTable = nullptr;
    } 
  if (fEnergyGrid)
    delete fEnergyGrid;
  
  if (fVerboseLevel > 2)
    G4cout << "G4PenelopeIonisationXSHandler. Tables have been cleared" 
       << G4endl;
}

References fDeltaTable, fEnergyGrid, fVerboseLevel, fXSTableElectron, fXSTablePositron, G4cout, and G4endl.

G4PenelopeIonisationXSHandler() [2/2]

G4PenelopeIonisationXSHandler::G4PenelopeIonisationXSHandler ( const G4PenelopeIonisationXSHandler &  )

delete

Member Function Documentation

BuildDeltaTable()

void G4PenelopeIonisationXSHandler::BuildDeltaTable ( const G4Material *  mat )

private

Definition at line 321 of file G4PenelopeIonisationXSHandler.cc.

{
  G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(mat);
  G4double plasmaSq = fOscManager->GetPlasmaEnergySquared(mat);
  G4double totalZ = fOscManager->GetTotalZ(mat);
  size_t numberOfOscillators = theTable->size();
 
  if (fEnergyGrid->GetVectorLength() != fNBins) 
    {
      G4ExceptionDescription ed;
      ed << "Energy Grid for Delta table looks not initialized" << G4endl;
      ed << fNBins << " " << fEnergyGrid->GetVectorLength() << G4endl;
      G4Exception("G4PenelopeIonisationXSHandler::BuildDeltaTable()",
          "em2030",FatalException,ed);
    }
 
  G4PhysicsFreeVector* theVector = new G4PhysicsFreeVector(fNBins);
 
  //loop on the energy grid
  for (size_t bin=0;bin<fNBins;bin++)
    {
      G4double delta = 0.;
      G4double energy = fEnergyGrid->GetLowEdgeEnergy(bin);
 
      //Here calculate delta
      G4double gam = 1.0+(energy/electron_mass_c2);
      G4double gamSq = gam*gam;
 
      G4double TST = totalZ/(gamSq*plasmaSq);
      G4double wl2 = 0;
      G4double fdel = 0;
 
      //loop on oscillators
      for (size_t i=0;i<numberOfOscillators;i++)
    {
      G4PenelopeOscillator* theOsc = (*theTable)[i];
      G4double wri = theOsc->GetResonanceEnergy();
      fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2);
    }      
      if (fdel >= TST) //if fdel < TST, delta = 0
    {
      //get last oscillator
      G4PenelopeOscillator* theOsc = (*theTable)[numberOfOscillators-1]; 
      wl2 = theOsc->GetResonanceEnergy()*theOsc->GetResonanceEnergy();
      
      //First iteration
      G4bool loopAgain = false;
      do
        {
          loopAgain = false;
          wl2 += wl2;
          fdel = 0.;
          for (size_t i=0;i<numberOfOscillators;i++)
        {
          G4PenelopeOscillator* theOscLocal1 = (*theTable)[i];
          G4double wri = theOscLocal1->GetResonanceEnergy();
          fdel += theOscLocal1->GetOscillatorStrength()/(wri*wri+wl2);
        }
          if (fdel > TST)
        loopAgain = true;
        }while(loopAgain);
 
      G4double wl2l = 0;
      G4double wl2u = wl2;
      //second iteration
      do
        {        
          loopAgain = false;
          wl2 = 0.5*(wl2l+wl2u);
          fdel = 0;
          for (size_t i=0;i<numberOfOscillators;i++)
        {
          G4PenelopeOscillator* theOscLocal2 = (*theTable)[i];
          G4double wri = theOscLocal2->GetResonanceEnergy();
          fdel += theOscLocal2->GetOscillatorStrength()/(wri*wri+wl2);
        }
          if (fdel > TST)
        wl2l = wl2;
          else
        wl2u = wl2;
          if ((wl2u-wl2l)>1e-12*wl2)
        loopAgain = true;
        }while(loopAgain);
      
      //Eventually get density correction
      delta = 0.;
      for (size_t i=0;i<numberOfOscillators;i++)
        {
          G4PenelopeOscillator* theOscLocal3 = (*theTable)[i];
          G4double wri = theOscLocal3->GetResonanceEnergy();
          delta += theOscLocal3->GetOscillatorStrength()*
        G4Log(1.0+(wl2/(wri*wri)));          
        }
      delta = (delta/totalZ)-wl2/(gamSq*plasmaSq);
    }
      energy = std::max(1e-9*eV,energy); //prevents log(0)
      theVector->PutValue(bin,G4Log(energy),delta);
    }
  fDeltaTable->insert(std::make_pair(mat,theVector));
  return;
}

References source.hepunit::electron_mass_c2, G4INCL::KinematicsUtils::energy(), eV, FatalException, fDeltaTable, fEnergyGrid, fNBins, fOscManager, G4endl, G4Exception(), G4Log(), G4InuclParticleNames::gam, G4PhysicsVector::GetLowEdgeEnergy(), G4PenelopeOscillator::GetOscillatorStrength(), G4PenelopeOscillatorManager::GetOscillatorTableIonisation(), G4PenelopeOscillatorManager::GetPlasmaEnergySquared(), G4PenelopeOscillator::GetResonanceEnergy(), G4PenelopeOscillatorManager::GetTotalZ(), G4PhysicsVector::GetVectorLength(), G4INCL::Math::max(), and G4PhysicsFreeVector::PutValue().

Referenced by BuildXSTable().

BuildXSTable()

void G4PenelopeIonisationXSHandler::BuildXSTable	(	const G4Material *	mat,
		G4double	cut,
		const G4ParticleDefinition *	part,
		G4bool	isMaster = true
	)

This can be inkoved only by the master.

Definition at line 158 of file G4PenelopeIonisationXSHandler.cc.

{
  //Just to check
  if (!isMaster)
    G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
                "em0100",FatalException,"Worker thread in this method");
  
  //
  //This method fills the G4PenelopeCrossSection containers for electrons or positrons
  //and for the given material/cut couple. The calculation is done as sum over the 
  //individual shells.
  //Equivalent of subroutines EINaT and PINaT of Penelope
  //
  if (fVerboseLevel > 2)
    {
      G4cout << "G4PenelopeIonisationXSHandler: going to build cross section table " << G4endl;
      G4cout << "for " << part->GetParticleName() << " in " << mat->GetName() << G4endl;
      G4cout << "Cut= " << cut/keV << " keV" << G4endl;
    }
 
  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
  //Check if the table already exists
  if (part == G4Electron::Electron())
    {
      if (fXSTableElectron->count(theKey)) //table already built    
    return;
    }
  if (part == G4Positron::Positron())
    {
      if (fXSTablePositron->count(theKey)) //table already built    
    return;
    }
  
  //check if the material has been built
  if (!(fDeltaTable->count(mat)))
    BuildDeltaTable(mat);
 
 
  //Tables have been already created (checked by GetCrossSectionTableForCouple)
  G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(mat);
  size_t numberOfOscillators = theTable->size();
 
  if (fEnergyGrid->GetVectorLength() != fNBins) 
    {
      G4ExceptionDescription ed;
      ed << "Energy Grid looks not initialized" << G4endl;
      ed << fNBins << " " << fEnergyGrid->GetVectorLength() << G4endl;
      G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
          "em2030",FatalException,ed);
    }
 
  G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(fNBins,numberOfOscillators);
 
  //loop on the energy grid
  for (size_t bin=0;bin<fNBins;bin++)
    {
       G4double energy = fEnergyGrid->GetLowEdgeEnergy(bin);
       G4double XH0=0, XH1=0, XH2=0;
       G4double XS0=0, XS1=0, XS2=0;
   
       //oscillator loop
       for (size_t iosc=0;iosc<numberOfOscillators;iosc++)
     {
       G4DataVector* tempStorage = nullptr;
 
       G4PenelopeOscillator* theOsc = (*theTable)[iosc];
       G4double delta = GetDensityCorrection(mat,energy);
       if (part == G4Electron::Electron())       
         tempStorage = ComputeShellCrossSectionsElectron(theOsc,energy,cut,delta);
       else if (part == G4Positron::Positron())
         tempStorage = ComputeShellCrossSectionsPositron(theOsc,energy,cut,delta);
       //check results are all right
       if (!tempStorage)
         {
           G4ExceptionDescription ed;
           ed << "Problem in calculating the shell XS for shell # " 
          << iosc << G4endl;
           G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
               "em2031",FatalException,ed);            
           delete XSEntry;
           return;
         }
       if (tempStorage->size() != 6)
         {
           G4ExceptionDescription ed;
           ed << "Problem in calculating the shell XS " << G4endl;
           ed << "Result has dimension " << tempStorage->size() << " instead of 6" << G4endl;
           G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
               "em2031",FatalException,ed); 
         }
       G4double stre = theOsc->GetOscillatorStrength();
 
       XH0 += stre*(*tempStorage)[0];
       XH1 += stre*(*tempStorage)[1];
       XH2 += stre*(*tempStorage)[2];
       XS0 += stre*(*tempStorage)[3];
       XS1 += stre*(*tempStorage)[4];
       XS2 += stre*(*tempStorage)[5];
       XSEntry->AddShellCrossSectionPoint(bin,iosc,energy,stre*(*tempStorage)[0]);
       if (tempStorage)
         {
           delete tempStorage;
           tempStorage = nullptr;
         }
     }       
       XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
    }
  //Do (only once) the final normalization
  XSEntry->NormalizeShellCrossSections();
 
  //Insert in the appropriate table
  if (part == G4Electron::Electron())      
    fXSTableElectron->insert(std::make_pair(theKey,XSEntry));
  else if (part == G4Positron::Positron())
    fXSTablePositron->insert(std::make_pair(theKey,XSEntry));
  else
    delete XSEntry;
  
  return;
}

References G4PenelopeCrossSection::AddCrossSectionPoint(), G4PenelopeCrossSection::AddShellCrossSectionPoint(), BuildDeltaTable(), ComputeShellCrossSectionsElectron(), ComputeShellCrossSectionsPositron(), G4Electron::Electron(), G4INCL::KinematicsUtils::energy(), FatalException, fDeltaTable, fEnergyGrid, fNBins, fOscManager, fVerboseLevel, fXSTableElectron, fXSTablePositron, G4cout, G4endl, G4Exception(), GetDensityCorrection(), G4PhysicsVector::GetLowEdgeEnergy(), G4Material::GetName(), G4PenelopeOscillator::GetOscillatorStrength(), G4PenelopeOscillatorManager::GetOscillatorTableIonisation(), G4ParticleDefinition::GetParticleName(), G4PhysicsVector::GetVectorLength(), keV, G4PenelopeCrossSection::NormalizeShellCrossSections(), and G4Positron::Positron().

Referenced by G4PenelopeIonisationModel::ComputeDEDXPerVolume(), G4PenelopeIonisationCrossSection::CrossSection(), G4PenelopeIonisationModel::CrossSectionPerVolume(), and G4PenelopeIonisationModel::Initialise().

ComputeShellCrossSectionsElectron()

G4DataVector * G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsElectron ( G4PenelopeOscillator *  theOsc, G4double  energy, G4double  cut, G4double  delta  )

private

Definition at line 424 of file G4PenelopeIonisationXSHandler.cc.

{
  //
  //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for 
  //the given oscillator/cut and at the given energy.
  //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
  //Equivalent of subroutines EINaT1 of Penelope
  //
  // Results are _per target electron_
  //
  G4DataVector* result = new G4DataVector();
  for (size_t i=0;i<6;i++)
    result->push_back(0.);
  G4double ionEnergy = theOsc->GetIonisationEnergy();
  
  //return a set of zero's if the energy it too low to excite the current oscillator
  if (energy < ionEnergy)
    return result;
 
  G4double H0=0.,H1=0.,H2=0.;
  G4double S0=0.,S1=0.,S2=0.;
 
  //Define useful constants to be used in the calculation
  G4double gamma = 1.0+energy/electron_mass_c2;
  G4double gammaSq = gamma*gamma;
  G4double beta = (gammaSq-1.0)/gammaSq;
  G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2
  G4double constant = pielr2*2.0*electron_mass_c2/beta;
  G4double XHDT0 = G4Log(gammaSq)-beta;
 
  G4double cpSq = energy*(energy+2.0*electron_mass_c2);
  G4double cp = std::sqrt(cpSq);
  G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
 
  //
  // Distant interactions
  //
  G4double resEne = theOsc->GetResonanceEnergy();
  G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
  if (energy > resEne)
    {
      G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
      G4double cp1 = std::sqrt(cp1Sq);
      
      //Distant longitudinal interactions
      G4double QM = 0;
      if (resEne > 1e-6*energy)
    QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
      else
    {
      QM = resEne*resEne/(beta*2.0*electron_mass_c2);
      QM = QM*(1.0-0.5*QM/electron_mass_c2);
    }
      G4double SDL1 = 0;
      if (QM < cutoffEne)
    SDL1 = G4Log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
      
      //Distant transverse interactions
      if (SDL1)
    {
      G4double SDT1 = std::max(XHDT0-delta,0.0);
      G4double SD1 = SDL1+SDT1;
      if (cut > resEne)
        {
          S1 = SD1; //XS1
          S0 = SD1/resEne; //XS0
          S2 = SD1*resEne; //XS2
        }
      else
        {
          H1 = SD1; //XH1
          H0 = SD1/resEne; //XH0
          H2 = SD1*resEne; //XH2
        }
    }
    }
  //
  // Close collisions (Moller's cross section)
  //
  G4double wl = std::max(cut,cutoffEne);
  G4double ee = energy + ionEnergy;
  G4double wu = 0.5*ee;
  if (wl < wu-(1e-5*eV))
    {
      H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + 
    (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu))/ee + 
    amol*(wu-wl)/(ee*ee);
      H1 += G4Log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + 
    (2.0-amol)*G4Log((ee-wu)/(ee-wl)) + 
    amol*(wu*wu-wl*wl)/(2.0*ee*ee);
      H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - 
    (wl*(2.0*ee-wl)/(ee-wl)) + 
    (3.0-amol)*ee*G4Log((ee-wu)/(ee-wl)) +  
    amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
      wu = wl;
    }
  wl = cutoffEne;
  
  if (wl > wu-(1e-5*eV))
    {
      (*result)[0] = constant*H0;
      (*result)[1] = constant*H1;
      (*result)[2] = constant*H2;
      (*result)[3] = constant*S0;
      (*result)[4] = constant*S1;
      (*result)[5] = constant*S2;
      return result;
    }
 
  S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + 
    (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
    amol*(wu-wl)/(ee*ee);
  S1 += G4Log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + 
    (2.0-amol)*G4Log((ee-wu)/(ee-wl)) + 
    amol*(wu*wu-wl*wl)/(2.0*ee*ee);
  S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - 
    (wl*(2.0*ee-wl)/(ee-wl)) + 
    (3.0-amol)*ee*G4Log((ee-wu)/(ee-wl)) + 
    amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
 
  (*result)[0] = constant*H0;
  (*result)[1] = constant*H1;
  (*result)[2] = constant*H2;
  (*result)[3] = constant*S0;
  (*result)[4] = constant*S1;
  (*result)[5] = constant*S2;
  return result;
}

References anonymous_namespace{G4PionRadiativeDecayChannel.cc}::beta, source.hepunit::classic_electr_radius, cp, source.hepunit::electron_mass_c2, G4INCL::KinematicsUtils::energy(), eV, G4Log(), G4PenelopeOscillator::GetCutoffRecoilResonantEnergy(), G4PenelopeOscillator::GetIonisationEnergy(), G4PenelopeOscillator::GetResonanceEnergy(), G4INCL::Math::max(), and pi.

Referenced by BuildXSTable().

ComputeShellCrossSectionsPositron()

G4DataVector * G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsPositron ( G4PenelopeOscillator *  theOsc, G4double  energy, G4double  cut, G4double  delta  )

private

Definition at line 556 of file G4PenelopeIonisationXSHandler.cc.

{
  //
  //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for 
  //the given oscillator/cut and at the given energy.
  //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
  //Equivalent of subroutines PINaT1 of Penelope
  //
  // Results are _per target electron_
  //
  G4DataVector* result = new G4DataVector();
  for (size_t i=0;i<6;i++)
    result->push_back(0.);
  G4double ionEnergy = theOsc->GetIonisationEnergy();
  
  //return a set of zero's if the energy it too low to excite the current oscillator
  if (energy < ionEnergy)
    return result;
 
  G4double H0=0.,H1=0.,H2=0.;
  G4double S0=0.,S1=0.,S2=0.;
 
  //Define useful constants to be used in the calculation
  G4double gamma = 1.0+energy/electron_mass_c2;
  G4double gammaSq = gamma*gamma;
  G4double beta = (gammaSq-1.0)/gammaSq;
  G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2
  G4double constant = pielr2*2.0*electron_mass_c2/beta;
  G4double XHDT0 = G4Log(gammaSq)-beta;
 
  G4double cpSq = energy*(energy+2.0*electron_mass_c2);
  G4double cp = std::sqrt(cpSq);
  G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
  G4double g12 = (gamma+1.0)*(gamma+1.0);
  //Bhabha coefficients
  G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq-1.0);
  G4double bha2 = amol*(3.0+1.0/g12);
  G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g12;
  G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/g12;
 
  //
  // Distant interactions
  //
  G4double resEne = theOsc->GetResonanceEnergy();
  G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
  if (energy > resEne)
    {
      G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
      G4double cp1 = std::sqrt(cp1Sq);
      
      //Distant longitudinal interactions
      G4double QM = 0;
      if (resEne > 1e-6*energy)
    QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
      else
    {
      QM = resEne*resEne/(beta*2.0*electron_mass_c2);
      QM = QM*(1.0-0.5*QM/electron_mass_c2);
    }
      G4double SDL1 = 0;
      if (QM < cutoffEne)
    SDL1 = G4Log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
      
      //Distant transverse interactions
      if (SDL1)
    {
      G4double SDT1 = std::max(XHDT0-delta,0.0);
      G4double SD1 = SDL1+SDT1;
      if (cut > resEne)
        {
          S1 = SD1; //XS1
          S0 = SD1/resEne; //XS0
          S2 = SD1*resEne; //XS2
        }
      else
        {
          H1 = SD1; //XH1
          H0 = SD1/resEne; //XH0
          H2 = SD1*resEne; //XH2
        }
    }
    }
  //
  // Close collisions (Bhabha's cross section)
  //
  G4double wl = std::max(cut,cutoffEne);
  G4double wu = energy; 
  G4double energySq = energy*energy;
  if (wl < wu-(1e-5*eV))
    {
      G4double wlSq = wl*wl;
      G4double wuSq = wu*wu;
      H0 += (1.0/wl) - (1.0/wu)- bha1*G4Log(wu/wl)/energy  
    + bha2*(wu-wl)/energySq  
    - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
    + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
      H1 += G4Log(wu/wl) - bha1*(wu-wl)/energy
    + bha2*(wuSq-wlSq)/(2.0*energySq)
    - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy)
    + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
      H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
    + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
    - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
    + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
      wu = wl;
    }
  wl = cutoffEne;
  
  if (wl > wu-(1e-5*eV))
    {
      (*result)[0] = constant*H0;
      (*result)[1] = constant*H1;
      (*result)[2] = constant*H2;
      (*result)[3] = constant*S0;
      (*result)[4] = constant*S1;
      (*result)[5] = constant*S2;
      return result;
    }
 
  G4double wlSq = wl*wl;
  G4double wuSq = wu*wu;
 
  S0 += (1.0/wl) - (1.0/wu) - bha1*G4Log(wu/wl)/energy 
    + bha2*(wu-wl)/energySq  
    - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
    + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
 
  S1 += G4Log(wu/wl) - bha1*(wu-wl)/energy
    + bha2*(wuSq-wlSq)/(2.0*energySq)
    - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy)
    + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
 
  S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
    + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
    - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
    + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
 
 (*result)[0] = constant*H0;
 (*result)[1] = constant*H1;
 (*result)[2] = constant*H2;
 (*result)[3] = constant*S0;
 (*result)[4] = constant*S1;
 (*result)[5] = constant*S2;
 
 return result;
}

References anonymous_namespace{G4PionRadiativeDecayChannel.cc}::beta, source.hepunit::classic_electr_radius, cp, source.hepunit::electron_mass_c2, G4INCL::KinematicsUtils::energy(), eV, G4Log(), G4PenelopeOscillator::GetCutoffRecoilResonantEnergy(), G4PenelopeOscillator::GetIonisationEnergy(), G4PenelopeOscillator::GetResonanceEnergy(), G4INCL::Math::max(), and pi.

Referenced by BuildXSTable().

GetCrossSectionTableForCouple()

const G4PenelopeCrossSection * G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple	(	const G4ParticleDefinition *	part,
		const G4Material *	mat,
		const G4double	cut
	)		const

Returns the table of cross sections for the given particle, given material and given cut as a G4PenelopeCrossSection* pointer.

Definition at line 109 of file G4PenelopeIonisationXSHandler.cc.

{
  if (part != G4Electron::Electron() && part != G4Positron::Positron())
    {
      G4ExceptionDescription ed;
      ed << "Invalid particle: " << part->GetParticleName() << G4endl;
      G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
          "em0001",FatalException,ed);
      return nullptr;
    }
 
  if (part == G4Electron::Electron())
    {
      if (!fXSTableElectron)
    {
      G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
              "em0028",FatalException,  
              "The Cross Section Table for e- was not initialized correctly!");
      return nullptr;
    }
      std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
      if (fXSTableElectron->count(theKey)) //table already built    
    return fXSTableElectron->find(theKey)->second;
      else   
        return nullptr;  
    }
  
  if (part == G4Positron::Positron())
    {
      if (!fXSTablePositron)
    {
      G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
              "em0028",FatalException,  
              "The Cross Section Table for e+ was not initialized correctly!");
      return nullptr;
    }
      std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
      if (fXSTablePositron->count(theKey)) //table already built    
    return fXSTablePositron->find(theKey)->second;
      else
        return nullptr;  
   }
  return nullptr;
}

References G4Electron::Electron(), FatalException, fXSTableElectron, fXSTablePositron, G4endl, G4Exception(), G4ParticleDefinition::GetParticleName(), and G4Positron::Positron().

Referenced by G4PenelopeIonisationModel::ComputeDEDXPerVolume(), G4PenelopeIonisationCrossSection::CrossSection(), G4PenelopeIonisationModel::CrossSectionPerVolume(), G4PenelopeIonisationModel::SampleFinalStateElectron(), and G4PenelopeIonisationModel::SampleFinalStatePositron().

GetDensityCorrection()

G4double G4PenelopeIonisationXSHandler::GetDensityCorrection	(	const G4Material *	mat,
		const G4double	energy
	)		const

Returns the density coeection for the material at the given energy.

Definition at line 284 of file G4PenelopeIonisationXSHandler.cc.

{
  G4double result = 0;
  if (!fDeltaTable)
    {
      G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()",
          "em2032",FatalException,
          "Delta Table not initialized. Was Initialise() run?");
      return 0;
    }
  if (energy <= 0*eV)
    {
      G4cout << "G4PenelopeIonisationXSHandler::GetDensityCorrection()" << G4endl;
      G4cout << "Invalid energy " << energy/eV << " eV " << G4endl;
      return 0;
    }
  G4double logene = G4Log(energy);
 
  if (fDeltaTable->count(mat))
    {
      const G4PhysicsFreeVector* vec = fDeltaTable->find(mat)->second;
      result = vec->Value(logene); //the table has delta vs. ln(E)      
    }
  else
    {
      G4ExceptionDescription ed;      
      ed << "Unable to build table for " << mat->GetName() << G4endl;
      G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()",
          "em2033",FatalException,ed);
    }
 
  return result;
}

References G4INCL::KinematicsUtils::energy(), eV, FatalException, fDeltaTable, G4cout, G4endl, G4Exception(), G4Log(), G4Material::GetName(), and G4PhysicsVector::Value().

Referenced by BuildXSTable(), G4PenelopeIonisationModel::SampleFinalStateElectron(), and G4PenelopeIonisationModel::SampleFinalStatePositron().

operator=()

G4PenelopeIonisationXSHandler & G4PenelopeIonisationXSHandler::operator= ( const G4PenelopeIonisationXSHandler &  right )

delete

SetVerboseLevel()

void G4PenelopeIonisationXSHandler::SetVerboseLevel ( G4int  vl )

inline

Setter for the verbosity level.

Definition at line 77 of file G4PenelopeIonisationXSHandler.hh.

{fVerboseLevel = vl;};

References fVerboseLevel.

Referenced by G4PenelopeIonisationModel::Initialise().

Field Documentation

fDeltaTable

std::map<const G4Material*,G4PhysicsFreeVector*>* G4PenelopeIonisationXSHandler::fDeltaTable

private

Definition at line 105 of file G4PenelopeIonisationXSHandler.hh.

Referenced by BuildDeltaTable(), BuildXSTable(), G4PenelopeIonisationXSHandler(), GetDensityCorrection(), and ~G4PenelopeIonisationXSHandler().

fEnergyGrid

G4PhysicsLogVector* G4PenelopeIonisationXSHandler::fEnergyGrid

private

Definition at line 108 of file G4PenelopeIonisationXSHandler.hh.

Referenced by BuildDeltaTable(), BuildXSTable(), G4PenelopeIonisationXSHandler(), and ~G4PenelopeIonisationXSHandler().

fNBins

size_t G4PenelopeIonisationXSHandler::fNBins

private

Definition at line 111 of file G4PenelopeIonisationXSHandler.hh.

Referenced by BuildDeltaTable(), BuildXSTable(), and G4PenelopeIonisationXSHandler().

fOscManager

G4PenelopeOscillatorManager* G4PenelopeIonisationXSHandler::fOscManager

private

Definition at line 98 of file G4PenelopeIonisationXSHandler.hh.

Referenced by BuildDeltaTable(), BuildXSTable(), and G4PenelopeIonisationXSHandler().

fVerboseLevel

G4int G4PenelopeIonisationXSHandler::fVerboseLevel

private

Definition at line 110 of file G4PenelopeIonisationXSHandler.hh.

Referenced by BuildXSTable(), G4PenelopeIonisationXSHandler(), SetVerboseLevel(), and ~G4PenelopeIonisationXSHandler().

fXSTableElectron

std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>* G4PenelopeIonisationXSHandler::fXSTableElectron

private

Definition at line 101 of file G4PenelopeIonisationXSHandler.hh.

Referenced by BuildXSTable(), G4PenelopeIonisationXSHandler(), GetCrossSectionTableForCouple(), and ~G4PenelopeIonisationXSHandler().

fXSTablePositron

std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>* G4PenelopeIonisationXSHandler::fXSTablePositron

private

Definition at line 102 of file G4PenelopeIonisationXSHandler.hh.

Referenced by BuildXSTable(), G4PenelopeIonisationXSHandler(), GetCrossSectionTableForCouple(), and ~G4PenelopeIonisationXSHandler().

---

The documentation for this class was generated from the following files:

• source/processes/electromagnetic/lowenergy/include/G4PenelopeIonisationXSHandler.hh
• source/processes/electromagnetic/lowenergy/src/G4PenelopeIonisationXSHandler.cc