90{
91 static const G4double kAlpha2Pi = CLHEP::twopi*CLHEP::fine_structure_const;
93 const G4int izet = std::max(std::min(fMaxZet, iZet),1);
95
96
98 const SamplingTablePerZ* stZ = fSBSamplingTables[izet];
99
100 const size_t gamCutIndx = stZ->fMatCutIndxToGamCutIndx[matCutIndx];
101
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) + ". ";
106 msg.c_str());
107 }
108 const G4double lGCut = stZ->fLogGammaECuts[gamCutIndx];
109
111
114 G4int elEnergyIndx = stZ->fMaxElEnergyIndx;
115
116 if (eekin < fElEnergyVect[elEnergyIndx]) {
117 const G4double val = (lElEnergy-fLogMinElEnergy)*fILDeltaElEnergy;
118 elEnergyIndx = (
G4int)val;
120
121 if (fElEnergyVect[elEnergyIndx]<=gcut) {
122
123 pIndxH = (lElEnergy-lGCut)/(fLElEnergyVect[elEnergyIndx+1]-lGCut);
124 isCorner = true;
125 }
126
127 if (rndmEngine->
flat()<pIndxH) {
128 ++elEnergyIndx;
129 } else if (isCorner) {
130
131
132 isSimply = true;
133 }
134 }
135
136 if (!stZ->fTablesPerEnergy[elEnergyIndx]) {
137 return 0.;
138 }
139
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
146 if (minVal >= 1.) {
147 return 0.;
148 }
149
150 const G4double lCurKappaC = lGCut - leekin;
151 const G4double lUsedKappaC = lGCut - fLElEnergyVect[elEnergyIndx];
152
155
156 do {
159 if (!isSimply) {
160 const G4double cumRV = rndm[0]*(1.-minVal)+minVal;
161
162
163 const G4int cumLIndx = LinSearch(pVect, fNumKappa, cumRV)-1;
164
165
166
167 const STPoint& stPL = pVect[cumLIndx];
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
178 const G4double lKappa = lKL+dm2/dm3*(lKH-lKL);
179
180 kappa =
G4Exp(lKappa*lCurKappaC/lUsedKappaC);
181 } else {
182
183
184 kappa = 1.-rndm[0]*(1.-gcut/eekin);
185 }
186
187 eGamma = kappa*eekin;
188 const G4double invEGamma = 1./eGamma;
189
190 suppression = 1./(1.+dielSupConst*invEGamma*invEGamma);
191
192 if (!isElectron) {
194 const G4double iBeta1 = (e1 + CLHEP::electron_mass_c2)
195 / std::sqrt(e1*(e1 + 2.*CLHEP::electron_mass_c2));
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
204
205 return eGamma;
206}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
virtual void flatArray(const int size, double *vect)=0