G4SBBremTable Class Reference

#include <G4SBBremTable.hh>

Data Structures
struct	SamplingTablePerZ
struct	STable
struct	STPoint

Public Member Functions
void	ClearSamplingTables ()
	G4SBBremTable ()
void	Initialize (const G4double lowe, const G4double highe)
double	SampleEnergyTransfer (const G4double eekin, const G4double leekin, const G4double gcut, const G4double dielSupConst, const G4int izet, const G4int matCutIndx, const bool iselectron)
	~G4SBBremTable ()

Private Member Functions
void	BuildSamplingTables ()
void	InitSamplingTables ()
G4int	LinSearch (const std::vector< STPoint > &vect, const G4int size, const G4double val)
void	LoadSamplingTables (G4int iz)
void	LoadSTGrid ()
void	ReadCompressedFile (const G4String &fname, std::istringstream &iss)

Private Attributes
std::vector< G4double >	fElEnergyVect
G4double	fILDeltaElEnergy
std::vector< G4double >	fKappaVect
std::vector< G4double >	fLElEnergyVect
std::vector< G4double >	fLKappaVect
G4double	fLogMinElEnergy
G4int	fMaxZet
G4int	fNumElEnergy
G4int	fNumKappa
std::vector< SamplingTablePerZ * >	fSBSamplingTables
G4double	fUsedHighEenergy
G4double	fUsedLowEenergy

Detailed Description

Definition at line 63 of file G4SBBremTable.hh.

Constructor & Destructor Documentation

G4SBBremTable()

G4SBBremTable::G4SBBremTable

(

)

Definition at line 63 of file G4SBBremTable.cc.

 : fMaxZet(-1), fNumElEnergy(-1), fNumKappa(-1), fUsedLowEenergy(-1.),
   fUsedHighEenergy(-1.), fLogMinElEnergy(-1.), fILDeltaElEnergy(-1.)
{}

~G4SBBremTable()

G4SBBremTable::~G4SBBremTable

(

)

Definition at line 68 of file G4SBBremTable.cc.

{
  ClearSamplingTables();
}

References ClearSamplingTables().

Member Function Documentation

BuildSamplingTables()

void G4SBBremTable::BuildSamplingTables ( )

private

Definition at line 209 of file G4SBBremTable.cc.

                                        {
  // claer
  ClearSamplingTables();
  LoadSTGrid();
  // First elements and gamma cuts data structures need to be built:
  // loop over all material-cuts and add gamma cut to the list of elements
  const G4ProductionCutsTable
  *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
  // a temporary vector to store one element
  std::vector<size_t> vtmp(1,0);
  size_t numMatCuts = thePCTable->GetTableSize();
  for (size_t imc=0; imc<numMatCuts; ++imc) {
    const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc);
    if (!matCut->IsUsed()) {
      continue;
    }
    const G4Material*           mat = matCut->GetMaterial();
    const G4ElementVector* elemVect = mat->GetElementVector();
    const G4int              indxMC = matCut->GetIndex();
    const G4double gamCut = (*(thePCTable->GetEnergyCutsVector(0)))[indxMC];
    const size_t numElems = elemVect->size();
    for (size_t ielem=0; ielem<numElems; ++ielem) {
      const G4Element *elem = (*elemVect)[ielem];
      const G4int izet = std::max(std::min(fMaxZet, elem->GetZasInt()),1);
      if (!fSBSamplingTables[izet]) {
        // create data structure but do not load sampling tables yet: will be
        // loaded after we know what are the min/max e- energies where sampling
        // will be needed during the simulation for this Z
        // LoadSamplingTables(izet);
        fSBSamplingTables[izet] = new SamplingTablePerZ();
      }
      // add current gamma cut to the list of this element data (only if this
      // cut value is still not tehre)
      const std::vector<double> &cVect = fSBSamplingTables[izet]->fGammaECuts;
      size_t indx = std::find(cVect.begin(), cVect.end(), gamCut)-cVect.begin();
      if (indx==cVect.size()) {
        vtmp[0] = imc;
        fSBSamplingTables[izet]->fGamCutIndxToMatCutIndx.push_back(vtmp);
        fSBSamplingTables[izet]->fGammaECuts.push_back(gamCut);
        fSBSamplingTables[izet]->fLogGammaECuts.push_back(G4Log(gamCut));
        ++fSBSamplingTables[izet]->fNumGammaCuts;
      } else {
        fSBSamplingTables[izet]->fGamCutIndxToMatCutIndx[indx].push_back(imc);
      }
    }
  }
}

References ClearSamplingTables(), elem, fMaxZet, fSBSamplingTables, G4Log(), G4Material::GetElementVector(), G4ProductionCutsTable::GetEnergyCutsVector(), G4MaterialCutsCouple::GetIndex(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4MaterialCutsCouple::IsUsed(), LoadSTGrid(), G4INCL::Math::max(), and G4INCL::Math::min().

Referenced by Initialize().

ClearSamplingTables()

void G4SBBremTable::ClearSamplingTables

(

)

Definition at line 460 of file G4SBBremTable.cc.

                                        {
  for (G4int iz=0; iz<fMaxZet+1; ++iz) {
    if (fSBSamplingTables[iz]) {
      for (G4int iee=0; iee<fNumElEnergy; ++iee) {
        if (fSBSamplingTables[iz]->fTablesPerEnergy[iee]) {
          fSBSamplingTables[iz]->fTablesPerEnergy[iee]->fSTable.clear();
          fSBSamplingTables[iz]->fTablesPerEnergy[iee]->fCumCutValues.clear();
        }
      }
      fSBSamplingTables[iz]->fTablesPerEnergy.clear();
      fSBSamplingTables[iz]->fGammaECuts.clear();
      fSBSamplingTables[iz]->fLogGammaECuts.clear();
      fSBSamplingTables[iz]->fMatCutIndxToGamCutIndx.clear();
      //
      delete fSBSamplingTables[iz];
      fSBSamplingTables[iz] = nullptr;
    }
  }
  fSBSamplingTables.clear();
  fElEnergyVect.clear();
  fLElEnergyVect.clear();
  fKappaVect.clear();
  fLKappaVect.clear();
  fMaxZet = -1;
}

References fElEnergyVect, fKappaVect, fLElEnergyVect, fLKappaVect, fMaxZet, fNumElEnergy, and fSBSamplingTables.

Referenced by BuildSamplingTables(), and ~G4SBBremTable().

Initialize()

void G4SBBremTable::Initialize	(	const G4double	lowe,
		const G4double	highe
	)

Definition at line 73 of file G4SBBremTable.cc.

{
  fUsedLowEenergy  = lowe;
  fUsedHighEenergy = highe;
  BuildSamplingTables();
  InitSamplingTables();
//  Dump();
}

References BuildSamplingTables(), fUsedHighEenergy, fUsedLowEenergy, and InitSamplingTables().

Referenced by G4SeltzerBergerModel::Initialise().

InitSamplingTables()

void G4SBBremTable::InitSamplingTables ( )

private

Definition at line 257 of file G4SBBremTable.cc.

                                       {
  const size_t numMatCuts = G4ProductionCutsTable::GetProductionCutsTable()
                            ->GetTableSize();
  for (G4int iz=1; iz<fMaxZet+1; ++iz) {
    SamplingTablePerZ* stZ = fSBSamplingTables[iz];
    if (!stZ) continue;
    // Load-in sampling table data:
    LoadSamplingTables(iz);
    // init data
    for (G4int iee=0; iee<fNumElEnergy; ++iee) {
      if (!stZ->fTablesPerEnergy[iee])
        continue;
      const G4double elEnergy = fElEnergyVect[iee];
      // 1 indicates that gamma production is not possible at this e- energy
      stZ->fTablesPerEnergy[iee]->fCumCutValues.resize(stZ->fNumGammaCuts,1.);
      // sort gamma cuts and other members accordingly
      for (size_t i=0; i<stZ->fNumGammaCuts-1; ++i) {
        for (size_t j=i+1; j<stZ->fNumGammaCuts; ++j) {
          if (stZ->fGammaECuts[j]<stZ->fGammaECuts[i]) {
            G4double dum0                   = stZ->fGammaECuts[i];
            G4double dum1                   = stZ->fLogGammaECuts[i];
            std::vector<size_t>   dumv      = stZ->fGamCutIndxToMatCutIndx[i];
            stZ->fGammaECuts[i]             = stZ->fGammaECuts[j];
            stZ->fLogGammaECuts[i]          = stZ->fLogGammaECuts[j];
            stZ->fGamCutIndxToMatCutIndx[i] = stZ->fGamCutIndxToMatCutIndx[j];
            stZ->fGammaECuts[j]             = dum0;
            stZ->fLogGammaECuts[j]          = dum1;
            stZ->fGamCutIndxToMatCutIndx[j] = dumv;
          }
        }
      }
      // set couple indices to store the corresponding gamma cut index
      stZ->fMatCutIndxToGamCutIndx.resize(numMatCuts,-1);
      for (size_t i=0; i<stZ->fGamCutIndxToMatCutIndx.size(); ++i) {
        for (size_t j=0; j<stZ->fGamCutIndxToMatCutIndx[i].size(); ++j) {
          stZ->fMatCutIndxToGamCutIndx[stZ->fGamCutIndxToMatCutIndx[i][j]] = i;
        }
      }
      // clear temporary vector
      for (size_t i=0; i<stZ->fGamCutIndxToMatCutIndx.size(); ++i) {
        stZ->fGamCutIndxToMatCutIndx[i].clear();
      }
      stZ->fGamCutIndxToMatCutIndx.clear();
      //  init for each gamma cut that are below the e- energy
      for (size_t ic=0; ic<stZ->fNumGammaCuts; ++ic) {
        const G4double gamCut = stZ->fGammaECuts[ic];
        if (elEnergy>gamCut) {
          // find lower kappa index; compute the 'xi' i.e. cummulative value for
          // gamCut/elEnergy
          const G4double cutKappa = std::max(1.e-12, gamCut/elEnergy);
          const G4int       iKLow = (cutKappa>1.e-12)
          ? std::lower_bound(fKappaVect.begin(), fKappaVect.end(), cutKappa)
            - fKappaVect.begin() -1
          : 0;
          const STPoint* stpL = &(stZ->fTablesPerEnergy[iee]->fSTable[iKLow]);
          const STPoint* stpH = &(stZ->fTablesPerEnergy[iee]->fSTable[iKLow+1]);
          const G4double pA   = stpL->fParA;
          const G4double pB   = stpL->fParB;
          const G4double etaL = stpL->fCum;
          const G4double etaH = stpH->fCum;
          const G4double alph = G4Log(cutKappa/fKappaVect[iKLow])
                               /G4Log(fKappaVect[iKLow+1]/fKappaVect[iKLow]);
          const G4double dum  = pA*(alph-1.)-1.-pB;
          G4double val = etaL;
          if (alph==0.) {
            stZ->fTablesPerEnergy[iee]->fCumCutValues[ic] = val;
          } else {
            val = -(dum+std::sqrt(dum*dum-4.*pB*alph*alph))/(2.*pB*alph);
            val = val*(etaH-etaL)+etaL;
            stZ->fTablesPerEnergy[iee]->fCumCutValues[ic] = val;
          }
        }
      }
    }
  }
}

References G4SBBremTable::STPoint::fCum, fElEnergyVect, G4SBBremTable::SamplingTablePerZ::fGamCutIndxToMatCutIndx, G4SBBremTable::SamplingTablePerZ::fGammaECuts, fKappaVect, G4SBBremTable::SamplingTablePerZ::fLogGammaECuts, G4SBBremTable::SamplingTablePerZ::fMatCutIndxToGamCutIndx, fMaxZet, fNumElEnergy, G4SBBremTable::SamplingTablePerZ::fNumGammaCuts, G4SBBremTable::STPoint::fParA, G4SBBremTable::STPoint::fParB, fSBSamplingTables, G4SBBremTable::SamplingTablePerZ::fTablesPerEnergy, G4Log(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), LoadSamplingTables(), and G4INCL::Math::max().

Referenced by Initialize().

LinSearch()

G4int G4SBBremTable::LinSearch ( const std::vector< STPoint > &  vect, const G4int  size, const G4double  val  )

private

Definition at line 501 of file G4SBBremTable.cc.

                                                   {
  G4int i= 0;
  while (i + 3 < size) {
    if (vect [i + 0].fCum > val) return i + 0;
    if (vect [i + 1].fCum > val) return i + 1;
    if (vect [i + 2].fCum > val) return i + 2;
    if (vect [i + 3].fCum > val) return i + 3;
    i += 4;
  }
  while (i < size) {
    if (vect [i].fCum > val)
      break;
    ++i;
  }
  return i;
}

Referenced by SampleEnergyTransfer().

LoadSamplingTables()

void G4SBBremTable::LoadSamplingTables ( G4int  iz )

private

Definition at line 387 of file G4SBBremTable.cc.

                                               {
  // check if grid needs to be loaded first
  if (fMaxZet<0) {
    LoadSTGrid();
  }
  // load data for a given Z only once
  iz = std::max(std::min(fMaxZet, iz),1);
  char* path = std::getenv("G4LEDATA");
  if (!path) {
    G4Exception("G4SBBremTable::LoadSamplingTables()","em0006",
                FatalException, "Environment variable G4LEDATA not defined");
    return;
  }
  const G4String fname =  G4String(path) + "/brem_SB/SBTables/sTableSB_"
                        + std::to_string(iz);
  std::istringstream infile(std::ios::in);
  // read the compressed data file into the stream
  ReadCompressedFile(fname, infile);
  // the SamplingTablePerZ object was already created, set size of containers
  // then load sampling table data for each electron energies
  SamplingTablePerZ* zTable = fSBSamplingTables[iz];
  //
  // Determine min/max elektron kinetic energies and indices
  const G4double minGammaCut = zTable->fGammaECuts[ std::min_element(
                 std::begin(zTable->fGammaECuts),std::end(zTable->fGammaECuts))
                -std::begin(zTable->fGammaECuts)];
  const G4double elEmin = std::max(fUsedLowEenergy, minGammaCut);
  const G4double elEmax = fUsedHighEenergy;
  // find low/high elecrton energy indices where tables will be needed
  // low:
  zTable->fMinElEnergyIndx = 0;
  if (elEmin>=fElEnergyVect[fNumElEnergy-1]) {
    zTable->fMinElEnergyIndx = fNumElEnergy-1;
  } else {
    zTable->fMinElEnergyIndx = std::lower_bound(fElEnergyVect.begin(),
                        fElEnergyVect.end(), elEmin) - fElEnergyVect.begin() -1;
  }
  // high:
  zTable->fMaxElEnergyIndx = 0;
  if (elEmax>=fElEnergyVect[fNumElEnergy-1]) {
    zTable->fMaxElEnergyIndx = fNumElEnergy-1;
  } else {
    // lower + 1
    zTable->fMaxElEnergyIndx = std::lower_bound(fElEnergyVect.begin(),
                           fElEnergyVect.end(), elEmax) - fElEnergyVect.begin();
  }
  // protect
  if (zTable->fMaxElEnergyIndx<=zTable->fMinElEnergyIndx) {
    return;
  }
  // load sampling tables that are needed: file is already in the stream
  zTable->fTablesPerEnergy.resize(fNumElEnergy, nullptr);
  for (G4int iee=0; iee<fNumElEnergy; ++iee) {
    // go over data that are not needed
    if (iee<zTable->fMinElEnergyIndx || iee>zTable->fMaxElEnergyIndx) {
      for (G4int ik=0; ik<fNumKappa; ++ik) {
        G4double dum;
        infile >> dum; infile >> dum; infile >> dum;
      }
    } else { // load data that are needed
      zTable->fTablesPerEnergy[iee] = new STable();
      zTable->fTablesPerEnergy[iee]->fSTable.resize(fNumKappa);
      for (G4int ik=0; ik<fNumKappa; ++ik) {
        STPoint &stP = zTable->fTablesPerEnergy[iee]->fSTable[ik];
        infile >> stP.fCum;
        infile >> stP.fParA;
        infile >> stP.fParB;
      }
    }
  }
}

References FatalException, G4SBBremTable::STPoint::fCum, fElEnergyVect, G4SBBremTable::SamplingTablePerZ::fGammaECuts, G4SBBremTable::SamplingTablePerZ::fMaxElEnergyIndx, fMaxZet, G4SBBremTable::SamplingTablePerZ::fMinElEnergyIndx, test::fname, fNumElEnergy, fNumKappa, G4SBBremTable::STPoint::fParA, G4SBBremTable::STPoint::fParB, fSBSamplingTables, G4SBBremTable::SamplingTablePerZ::fTablesPerEnergy, fUsedHighEenergy, fUsedLowEenergy, G4Exception(), LoadSTGrid(), G4INCL::Math::max(), G4INCL::Math::min(), and ReadCompressedFile().

Referenced by InitSamplingTables().

LoadSTGrid()

void G4SBBremTable::LoadSTGrid ( )

private

Definition at line 335 of file G4SBBremTable.cc.

                               {
  char* path = std::getenv("G4LEDATA");
  if (!path) {
    G4Exception("G4SBBremTable::LoadSTGrid()","em0006",
                FatalException, "Environment variable G4LEDATA not defined");
    return;
  }
  const G4String fname =  G4String(path) + "/brem_SB/SBTables/grid";
  std::ifstream infile(fname,std::ios::in);
  if (!infile.is_open()) {
    G4String msgc = "Cannot open file: " + fname;
    G4Exception("G4SBBremTable::LoadSTGrid()","em0006",
                FatalException, msgc.c_str());
    return;
  }
  // get max Z, # electron energies and # kappa values
  infile >> fMaxZet;
  infile >> fNumElEnergy;
  infile >> fNumKappa;
  // allocate space for the data and get them in:
  // (1.) first eletron energy grid
  fElEnergyVect.resize(fNumElEnergy);
  fLElEnergyVect.resize(fNumElEnergy);
  for (G4int iee=0; iee<fNumElEnergy; ++iee) {
    G4double  dum;
    infile >> dum;
    fElEnergyVect[iee]  = dum*CLHEP::MeV;
    fLElEnergyVect[iee] = G4Log(fElEnergyVect[iee]);
  }
  // (2.) then the kappa grid
  fKappaVect.resize(fNumKappa);
  fLKappaVect.resize(fNumKappa);
  for (G4int ik=0; ik<fNumKappa; ++ik) {
    infile >> fKappaVect[ik];
    fLKappaVect[ik] = G4Log(fKappaVect[ik]);
  }
  // (3.) set size of the main container for sampling tables
  fSBSamplingTables.resize(fMaxZet+1,nullptr);
  // init electron energy grid related variables: use accurate values !!!
//  fLogMinElEnergy   = G4Log(fElEnergyVect[0]);
//  fILDeltaElEnergy  = 1./G4Log(fElEnergyVect[1]/fElEnergyVect[0]);
  const G4double elEmin  = 100.0*CLHEP::eV; //fElEnergyVect[0];
  const G4double elEmax  =  10.0*CLHEP::GeV;//fElEnergyVect[fNumElEnergy-1];
  fLogMinElEnergy  = G4Log(elEmin);
  fILDeltaElEnergy = 1./(G4Log(elEmax/elEmin)/(fNumElEnergy-1.0));
  // reset min/max energies if needed
  fUsedLowEenergy  = std::max(fUsedLowEenergy ,elEmin);
  fUsedHighEenergy = std::min(fUsedHighEenergy,elEmax);
  //
  infile.close();
}

References CLHEP::eV, FatalException, fElEnergyVect, fILDeltaElEnergy, fKappaVect, fLElEnergyVect, fLKappaVect, fLogMinElEnergy, fMaxZet, test::fname, fNumElEnergy, fNumKappa, fSBSamplingTables, fUsedHighEenergy, fUsedLowEenergy, G4Exception(), G4Log(), CLHEP::GeV, G4INCL::Math::max(), CLHEP::MeV, and G4INCL::Math::min().

Referenced by BuildSamplingTables(), and LoadSamplingTables().

ReadCompressedFile()

void G4SBBremTable::ReadCompressedFile ( const G4String &  fname, std::istringstream &  iss  )

private

Definition at line 521 of file G4SBBremTable.cc.

                                                              {
  std::string *dataString = nullptr;
  std::string compfilename(fname+".z");
  // create input stream with binary mode operation and positioning at the end
  // of the file
  std::ifstream in(compfilename, std::ios::binary | std::ios::ate);
  if (in.good()) {
     // get current position in the stream (was set to the end)
     int fileSize = in.tellg();
     // set current position being the beginning of the stream
     in.seekg(0,std::ios::beg);
     // create (zlib) byte buffer for the data
     Bytef *compdata = new Bytef[fileSize];
     while(in) {
        in.read((char*)compdata, fileSize);
     }
     // create (zlib) byte buffer for the uncompressed data
     uLongf complen    = (uLongf)(fileSize*4);
     Bytef *uncompdata = new Bytef[complen];
     while (Z_OK!=uncompress(uncompdata, &complen, compdata, fileSize)) {
        // increase uncompressed byte buffer
        delete[] uncompdata;
        complen   *= 2;
        uncompdata = new Bytef[complen];
     }
     // delete the compressed data buffer
     delete [] compdata;
     // create a string from the uncompressed data (will be deleted by the caller)
     dataString = new std::string((char*)uncompdata, (long)complen);
     // delete the uncompressed data buffer
     delete [] uncompdata;
  } else {
    std::string msg = "  Problem while trying to read "
                      + compfilename + " data file.\n";
    G4Exception("G4SBBremTable::ReadCompressedFile","em0006",
                FatalException,msg.c_str());
    return;
  }
  // create the input string stream from the data string
  if (dataString) {
    iss.str(*dataString);
    in.close();
    delete dataString;
  }
}

References FatalException, test::fname, G4Exception(), uncompress(), and Z_OK.

Referenced by LoadSamplingTables().

SampleEnergyTransfer()

double G4SBBremTable::SampleEnergyTransfer	(	const G4double	eekin,
		const G4double	leekin,
		const G4double	gcut,
		const G4double	dielSupConst,
		const G4int	izet,
		const G4int	matCutIndx,
		const bool	iselectron
	)

Definition at line 83 of file G4SBBremTable.cc.

{
  static const G4double kAlpha2Pi = CLHEP::twopi*CLHEP::fine_structure_const;
  const G4double zet = (G4double)iZet;
  const G4int   izet = std::max(std::min(fMaxZet, iZet),1);
  G4double eGamma    = 0.;
  // this should be checked in the caller
  // if (eekin<=gcut) return kappa;
  const G4double lElEnergy     = leekin;
  const SamplingTablePerZ* stZ = fSBSamplingTables[izet];
  // get the gamma cut of this Z that corresponds to the current mat-cuts
  const size_t gamCutIndx = stZ->fMatCutIndxToGamCutIndx[matCutIndx];
  // gcut was not found: should never happen (only in verbose mode)
  if (gamCutIndx >= stZ->fNumGammaCuts || stZ->fGammaECuts[gamCutIndx]!=gcut) {
    G4String msg = " Gamma cut="+std::to_string(gcut) + " [MeV] was not found ";
    msg += "in case of Z = " + std::to_string(iZet) + ". ";
    G4Exception("G4SBBremTable::SampleEnergyTransfer()","em0X",FatalException,
                msg.c_str());
  }
  const G4double lGCut = stZ->fLogGammaECuts[gamCutIndx];
  // get the random engine
  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
  // find lower e- energy bin
  G4bool isCorner = false; // indicate that the lower edge e- energy < gam-gut
  G4bool isSimply = false; // simply sampling: isCorner+lower egde is selected
  G4int elEnergyIndx   = stZ->fMaxElEnergyIndx;
  // only if e- ekin is below the maximum value(use table at maximum otherwise)
  if (eekin < fElEnergyVect[elEnergyIndx]) {
    const G4double val = (lElEnergy-fLogMinElEnergy)*fILDeltaElEnergy;
    elEnergyIndx       = (G4int)val;
    G4double pIndxH    = val-elEnergyIndx;
    // check if we are at limiting case: lower edge e- energy < gam-gut
    if (fElEnergyVect[elEnergyIndx]<=gcut) {
      // recompute the probability of taking the higher e- energy bin()
      pIndxH   = (lElEnergy-lGCut)/(fLElEnergyVect[elEnergyIndx+1]-lGCut);
      isCorner = true;
    }
    //
    if (rndmEngine->flat()<pIndxH) {
      ++elEnergyIndx;      // take the table at the higher e- energy bin
    } else if (isCorner) { // take the table at the lower  e- energy bin
      // special sampling need to be done if lower edge e- energy < gam-gut:
      // actually, we "sample" from a table "built" at the gamm-cut i.e. delta
      isSimply = true;
    }
  }
  // should never happen under normal conditions but add protection
  if (!stZ->fTablesPerEnergy[elEnergyIndx]) {
    return 0.;
  }
  // Do the photon energy sampling:
  const STable *st =  stZ->fTablesPerEnergy[elEnergyIndx];
  const std::vector<G4double>& cVect = st->fCumCutValues;
  const std::vector<STPoint>&  pVect = st->fSTable;
  const G4double minVal = cVect[gamCutIndx];
 
  // should never happen under normal conditions but add protection
  if (minVal >= 1.) {
    return 0.;
  }
  // some transfomrmtion variables used in the looop
  const G4double lCurKappaC  = lGCut - leekin;
  const G4double lUsedKappaC = lGCut - fLElEnergyVect[elEnergyIndx];
  // dielectric (always) and e+ correction suppressions (if the primary is e+)
  G4double suppression = 1.;
  G4double rndm[2];
  // rejection loop starts here
  do {
    rndmEngine->flatArray(2, rndm);
    G4double kappa = 1.0;
    if (!isSimply) {
      const G4double cumRV = rndm[0]*(1.-minVal)+minVal;
      // find lower index of the values in the Cumulative Function: use linear
      // instead of binary search because it's faster in our case
      const G4int cumLIndx = LinSearch(pVect, fNumKappa, cumRV)-1;
//      const G4int cumLIndx = std::lower_bound( pVect.begin(), pVect.end(), cumRV,
//                                    [](const STPoint& p, const double& cumV) {
//                                    return p.fCum<cumV; } ) - pVect.begin() -1;
      const STPoint& stPL  = pVect[cumLIndx];
      const G4double pA    = stPL.fParA;
      const G4double pB    = stPL.fParB;
      const G4double cumL  = stPL.fCum;
      const G4double cumH  = pVect[cumLIndx+1].fCum;
      const G4double lKL   = fLKappaVect[cumLIndx];
      const G4double lKH   = fLKappaVect[cumLIndx+1];
      const G4double dm1   = (cumRV-cumL)/(cumH-cumL);
      const G4double dm2   = (1.+pA+pB)*dm1;
      const G4double dm3   = 1.+dm1*(pA+pB*dm1);
      // kappa sampled at E_i e- energy
      const G4double lKappa = lKL+dm2/dm3*(lKH-lKL);
      // transform lKappa to [log(gcut/eekin),0] form [log(gcut/E_i),0]
      kappa  = G4Exp(lKappa*lCurKappaC/lUsedKappaC);
    } else {
//      const G4double upLimit = std::min(1.*CLHEP::eV,eekin-gcut);
//      kappa = (gcut+rndm[0]*upLimit)/eekin;
      kappa = 1.-rndm[0]*(1.-gcut/eekin);
    }
    // compute the emitted photon energy: k
    eGamma = kappa*eekin;
    const G4double invEGamma = 1./eGamma;
    // compute dielectric suppression: 1/(1+[gk_p/k]^2)
    suppression = 1./(1.+dielSupConst*invEGamma*invEGamma);
    // add positron correction if particle is e+
    if (!isElectron) {
      const G4double e1     = eekin - gcut;
      const G4double iBeta1 =  (e1 + CLHEP::electron_mass_c2)
                              / std::sqrt(e1*(e1 + 2.*CLHEP::electron_mass_c2));
      const G4double e2     = eekin - eGamma;
      const G4double iBeta2 =  (e2 + CLHEP::electron_mass_c2)
                              / std::sqrt(e2*(e2 + 2.*CLHEP::electron_mass_c2));
      const G4double dum    = kAlpha2Pi*zet*(iBeta1 - iBeta2);
      suppression = (dum > -12.) ? suppression*G4Exp(dum) : 0.;
    }
  } while (rndm[1] > suppression);
  // end of rejection loop
  // return the sampled photon energy value k
  return eGamma;
}

References e1, e2, CLHEP::electron_mass_c2, FatalException, G4SBBremTable::STPoint::fCum, G4SBBremTable::STable::fCumCutValues, fElEnergyVect, G4SBBremTable::SamplingTablePerZ::fGammaECuts, fILDeltaElEnergy, CLHEP::fine_structure_const, CLHEP::HepRandomEngine::flat(), CLHEP::HepRandomEngine::flatArray(), fLElEnergyVect, fLKappaVect, G4SBBremTable::SamplingTablePerZ::fLogGammaECuts, fLogMinElEnergy, G4SBBremTable::SamplingTablePerZ::fMatCutIndxToGamCutIndx, G4SBBremTable::SamplingTablePerZ::fMaxElEnergyIndx, fMaxZet, G4SBBremTable::SamplingTablePerZ::fNumGammaCuts, fNumKappa, G4SBBremTable::STPoint::fParA, G4SBBremTable::STPoint::fParB, fSBSamplingTables, G4SBBremTable::STable::fSTable, G4SBBremTable::SamplingTablePerZ::fTablesPerEnergy, G4Exception(), G4Exp(), G4InuclParticleNames::isElectron(), LinSearch(), G4INCL::Math::max(), G4INCL::Math::min(), and CLHEP::twopi.

Referenced by G4SeltzerBergerModel::SampleSecondaries().

Field Documentation

fElEnergyVect

std::vector<G4double> G4SBBremTable::fElEnergyVect

private

Definition at line 152 of file G4SBBremTable.hh.

Referenced by ClearSamplingTables(), InitSamplingTables(), LoadSamplingTables(), LoadSTGrid(), and SampleEnergyTransfer().

fILDeltaElEnergy

G4double G4SBBremTable::fILDeltaElEnergy

private

Definition at line 149 of file G4SBBremTable.hh.

Referenced by LoadSTGrid(), and SampleEnergyTransfer().

fKappaVect

std::vector<G4double> G4SBBremTable::fKappaVect

private

Definition at line 154 of file G4SBBremTable.hh.

Referenced by ClearSamplingTables(), InitSamplingTables(), and LoadSTGrid().

fLElEnergyVect

std::vector<G4double> G4SBBremTable::fLElEnergyVect

private

Definition at line 153 of file G4SBBremTable.hh.

Referenced by ClearSamplingTables(), LoadSTGrid(), and SampleEnergyTransfer().

fLKappaVect

std::vector<G4double> G4SBBremTable::fLKappaVect

private

Definition at line 155 of file G4SBBremTable.hh.

Referenced by ClearSamplingTables(), LoadSTGrid(), and SampleEnergyTransfer().

fLogMinElEnergy

G4double G4SBBremTable::fLogMinElEnergy

private

Definition at line 148 of file G4SBBremTable.hh.

Referenced by LoadSTGrid(), and SampleEnergyTransfer().

fMaxZet

G4int G4SBBremTable::fMaxZet

private

Definition at line 141 of file G4SBBremTable.hh.

Referenced by BuildSamplingTables(), ClearSamplingTables(), InitSamplingTables(), LoadSamplingTables(), LoadSTGrid(), and SampleEnergyTransfer().

fNumElEnergy

G4int G4SBBremTable::fNumElEnergy

private

Definition at line 142 of file G4SBBremTable.hh.

Referenced by ClearSamplingTables(), InitSamplingTables(), LoadSamplingTables(), and LoadSTGrid().

fNumKappa

G4int G4SBBremTable::fNumKappa

private

Definition at line 143 of file G4SBBremTable.hh.

Referenced by LoadSamplingTables(), LoadSTGrid(), and SampleEnergyTransfer().

fSBSamplingTables

std::vector<SamplingTablePerZ*> G4SBBremTable::fSBSamplingTables

private

Definition at line 158 of file G4SBBremTable.hh.

Referenced by BuildSamplingTables(), ClearSamplingTables(), InitSamplingTables(), LoadSamplingTables(), LoadSTGrid(), and SampleEnergyTransfer().

fUsedHighEenergy

G4double G4SBBremTable::fUsedHighEenergy

private

Definition at line 147 of file G4SBBremTable.hh.

Referenced by Initialize(), LoadSamplingTables(), and LoadSTGrid().

fUsedLowEenergy

G4double G4SBBremTable::fUsedLowEenergy

private

Definition at line 146 of file G4SBBremTable.hh.

Referenced by Initialize(), LoadSamplingTables(), and LoadSTGrid().

---

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

• source/processes/electromagnetic/standard/include/G4SBBremTable.hh
• source/processes/electromagnetic/standard/src/G4SBBremTable.cc