Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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 63 of file G4SBBremTable.cc.

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

◆ ~G4SBBremTable()

G4SBBremTable::~G4SBBremTable ( )

Definition at line 68 of file G4SBBremTable.cc.

69{
71}
void ClearSamplingTables()

Member Function Documentation

◆ ClearSamplingTables()

void G4SBBremTable::ClearSamplingTables ( )

Definition at line 460 of file G4SBBremTable.cc.

460 {
461 for (G4int iz=0; iz<fMaxZet+1; ++iz) {
462 if (fSBSamplingTables[iz]) {
463 for (G4int iee=0; iee<fNumElEnergy; ++iee) {
464 if (fSBSamplingTables[iz]->fTablesPerEnergy[iee]) {
465 fSBSamplingTables[iz]->fTablesPerEnergy[iee]->fSTable.clear();
466 fSBSamplingTables[iz]->fTablesPerEnergy[iee]->fCumCutValues.clear();
467 }
468 }
469 fSBSamplingTables[iz]->fTablesPerEnergy.clear();
470 fSBSamplingTables[iz]->fGammaECuts.clear();
471 fSBSamplingTables[iz]->fLogGammaECuts.clear();
472 fSBSamplingTables[iz]->fMatCutIndxToGamCutIndx.clear();
473 //
474 delete fSBSamplingTables[iz];
475 fSBSamplingTables[iz] = nullptr;
476 }
477 }
478 fSBSamplingTables.clear();
479 fElEnergyVect.clear();
480 fLElEnergyVect.clear();
481 fKappaVect.clear();
482 fLKappaVect.clear();
483 fMaxZet = -1;
484}
int G4int
Definition: G4Types.hh:85

Referenced by ~G4SBBremTable().

◆ Initialize()

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

Definition at line 73 of file G4SBBremTable.cc.

74{
75 fUsedLowEenergy = lowe;
76 fUsedHighEenergy = highe;
77 BuildSamplingTables();
78 InitSamplingTables();
79// Dump();
80}

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 83 of file G4SBBremTable.cc.

90{
91 static const G4double kAlpha2Pi = CLHEP::twopi*CLHEP::fine_structure_const;
92 const G4double zet = (G4double)iZet;
93 const G4int izet = std::max(std::min(fMaxZet, iZet),1);
94 G4double eGamma = 0.;
95 // this should be checked in the caller
96 // if (eekin<=gcut) return kappa;
97 const G4double lElEnergy = leekin;
98 const SamplingTablePerZ* stZ = fSBSamplingTables[izet];
99 // get the gamma cut of this Z that corresponds to the current mat-cuts
100 const size_t gamCutIndx = stZ->fMatCutIndxToGamCutIndx[matCutIndx];
101 // gcut was not found: should never happen (only in verbose mode)
102 if (gamCutIndx >= stZ->fNumGammaCuts || stZ->fGammaECuts[gamCutIndx]!=gcut) {
103 G4String msg = " Gamma cut="+std::to_string(gcut) + " [MeV] was not found ";
104 msg += "in case of Z = " + std::to_string(iZet) + ". ";
105 G4Exception("G4SBBremTable::SampleEnergyTransfer()","em0X",FatalException,
106 msg.c_str());
107 }
108 const G4double lGCut = stZ->fLogGammaECuts[gamCutIndx];
109 // get the random engine
110 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
111 // find lower e- energy bin
112 G4bool isCorner = false; // indicate that the lower edge e- energy < gam-gut
113 G4bool isSimply = false; // simply sampling: isCorner+lower egde is selected
114 G4int elEnergyIndx = stZ->fMaxElEnergyIndx;
115 // only if e- ekin is below the maximum value(use table at maximum otherwise)
116 if (eekin < fElEnergyVect[elEnergyIndx]) {
117 const G4double val = (lElEnergy-fLogMinElEnergy)*fILDeltaElEnergy;
118 elEnergyIndx = (G4int)val;
119 G4double pIndxH = val-elEnergyIndx;
120 // check if we are at limiting case: lower edge e- energy < gam-gut
121 if (fElEnergyVect[elEnergyIndx]<=gcut) {
122 // recompute the probability of taking the higher e- energy bin()
123 pIndxH = (lElEnergy-lGCut)/(fLElEnergyVect[elEnergyIndx+1]-lGCut);
124 isCorner = true;
125 }
126 //
127 if (rndmEngine->flat()<pIndxH) {
128 ++elEnergyIndx; // take the table at the higher e- energy bin
129 } else if (isCorner) { // take the table at the lower e- energy bin
130 // special sampling need to be done if lower edge e- energy < gam-gut:
131 // actually, we "sample" from a table "built" at the gamm-cut i.e. delta
132 isSimply = true;
133 }
134 }
135 // should never happen under normal conditions but add protection
136 if (!stZ->fTablesPerEnergy[elEnergyIndx]) {
137 return 0.;
138 }
139 // Do the photon energy sampling:
140 const STable *st = stZ->fTablesPerEnergy[elEnergyIndx];
141 const std::vector<G4double>& cVect = st->fCumCutValues;
142 const std::vector<STPoint>& pVect = st->fSTable;
143 const G4double minVal = cVect[gamCutIndx];
144
145 // should never happen under normal conditions but add protection
146 if (minVal >= 1.) {
147 return 0.;
148 }
149 // some transfomrmtion variables used in the looop
150 const G4double lCurKappaC = lGCut - leekin;
151 const G4double lUsedKappaC = lGCut - fLElEnergyVect[elEnergyIndx];
152 // dielectric (always) and e+ correction suppressions (if the primary is e+)
153 G4double suppression = 1.;
154 G4double rndm[2];
155 // rejection loop starts here
156 do {
157 rndmEngine->flatArray(2, rndm);
158 G4double kappa = 1.0;
159 if (!isSimply) {
160 const G4double cumRV = rndm[0]*(1.-minVal)+minVal;
161 // find lower index of the values in the Cumulative Function: use linear
162 // instead of binary search because it's faster in our case
163 const G4int cumLIndx = LinSearch(pVect, fNumKappa, cumRV)-1;
164// const G4int cumLIndx = std::lower_bound( pVect.begin(), pVect.end(), cumRV,
165// [](const STPoint& p, const double& cumV) {
166// return p.fCum<cumV; } ) - pVect.begin() -1;
167 const STPoint& stPL = pVect[cumLIndx];
168 const G4double pA = stPL.fParA;
169 const G4double pB = stPL.fParB;
170 const G4double cumL = stPL.fCum;
171 const G4double cumH = pVect[cumLIndx+1].fCum;
172 const G4double lKL = fLKappaVect[cumLIndx];
173 const G4double lKH = fLKappaVect[cumLIndx+1];
174 const G4double dm1 = (cumRV-cumL)/(cumH-cumL);
175 const G4double dm2 = (1.+pA+pB)*dm1;
176 const G4double dm3 = 1.+dm1*(pA+pB*dm1);
177 // kappa sampled at E_i e- energy
178 const G4double lKappa = lKL+dm2/dm3*(lKH-lKL);
179 // transform lKappa to [log(gcut/eekin),0] form [log(gcut/E_i),0]
180 kappa = G4Exp(lKappa*lCurKappaC/lUsedKappaC);
181 } else {
182// const G4double upLimit = std::min(1.*CLHEP::eV,eekin-gcut);
183// kappa = (gcut+rndm[0]*upLimit)/eekin;
184 kappa = 1.-rndm[0]*(1.-gcut/eekin);
185 }
186 // compute the emitted photon energy: k
187 eGamma = kappa*eekin;
188 const G4double invEGamma = 1./eGamma;
189 // compute dielectric suppression: 1/(1+[gk_p/k]^2)
190 suppression = 1./(1.+dielSupConst*invEGamma*invEGamma);
191 // add positron correction if particle is e+
192 if (!isElectron) {
193 const G4double e1 = eekin - gcut;
194 const G4double iBeta1 = (e1 + CLHEP::electron_mass_c2)
195 / std::sqrt(e1*(e1 + 2.*CLHEP::electron_mass_c2));
196 const G4double e2 = eekin - eGamma;
197 const G4double iBeta2 = (e2 + CLHEP::electron_mass_c2)
198 / std::sqrt(e2*(e2 + 2.*CLHEP::electron_mass_c2));
199 const G4double dum = kAlpha2Pi*zet*(iBeta1 - iBeta2);
200 suppression = (dum > -12.) ? suppression*G4Exp(dum) : 0.;
201 }
202 } while (rndm[1] > suppression);
203 // end of rejection loop
204 // return the sampled photon energy value k
205 return eGamma;
206}
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
virtual double flat()=0
virtual void flatArray(const int size, double *vect)=0

Referenced by G4SeltzerBergerModel::SampleSecondaries().


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