G4SBBremTable Class Reference

#include <G4SBBremTable.hh>

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

Detailed Description

Definition at line 63 of file G4SBBremTable.hh.

Constructor & Destructor Documentation

G4SBBremTable()

G4SBBremTable::G4SBBremTable

(

)

Definition at line 64 of file G4SBBremTable.cc.

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

~G4SBBremTable()

G4SBBremTable::~G4SBBremTable

(

)

Definition at line 69 of file G4SBBremTable.cc.

{
  ClearSamplingTables();
}

Member Function Documentation

ClearSamplingTables()

void G4SBBremTable::ClearSamplingTables

(

)

Definition at line 452 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;
}

Referenced by ~G4SBBremTable().

Initialize()

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

Definition at line 74 of file G4SBBremTable.cc.

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

Referenced by G4SeltzerBergerModel::Initialise().

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 84 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 std::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;
}

Referenced by G4SeltzerBergerModel::SampleSecondaries().

---

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

• G4SBBremTable.hh
• G4SBBremTable.cc