164{
165
167
168
170
172
173
175
176
178 if(nShells > (
G4int)fProbabilities.size()) { fProbabilities.resize(nShells); }
181 for(i=0; i<nShells; ++i) {
182
184
185 fProbabilities[i] = totprob;
186 }
187
188
189 static const G4int nlooplim = 1000;
191
194
197
198 do {
199 ++nloop;
200
201
204
205
206 for(i=0; i<nShells; ++i) { if(xprob <= fProbabilities[i]) { break; } }
207
209 lv1.
set(0.0,0.0,energy,energy);
210
211
212
213
214
215
219
220
221 G4double eTotMomentum = sqrt(eKinEnergy*(eKinEnergy + 2*electron_mass_c2));
224 G4double sintet = sqrt((1 - costet)*(1 + costet));
225 lv2.
set(eTotMomentum*sintet*cos(phi),eTotMomentum*sintet*sin(phi),
226 eTotMomentum*costet,eKinEnergy + electron_mass_c2);
229
230 gamEnergy0 = lv1.
e();
231
232
233
234
235
236 G4double E0_m = gamEnergy0/electron_mass_c2;
237
238
239
240
241
242
244
249
250 do {
251 ++nloop;
252
253 if(nloop > nlooplim) { return; }
254
255
257
258 if ( alpha1 >
alpha2*rndm[0] ) {
261
262 } else {
263 epsilonsq = epsilon0sq + (1.- epsilon0sq)*rndm[1];
265 }
266
268 sint2 = onecost*(2.-onecost);
269 greject = 1. -
epsilon*sint2/(1.+ epsilonsq);
270
271
272 } while (greject < rndm[2]);
273 gamEnergy1 =
epsilon*gamEnergy0;
274
275
276 lv2.
set(0.0,0.0,0.0,electron_mass_c2);
277 lv2 += lv1;
278
279
280
281
282 if(sint2 < 0.0) { sint2 = 0.0; }
283 costet = 1. - onecost;
284 sintet = sqrt(sint2);
285 phi = twopi * rndmEngineMod->
flat();
286
287
288
289
293 lv1.
set(gamEnergy1*v.
x(),gamEnergy1*v.
y(),gamEnergy1*v.
z(),gamEnergy1);
294 lv2 -= lv1;
295
296
298 eKinEnergy = lv2.
e() - electron_mass_c2 - ePotEnergy;
299
300
301
302 } while ( eKinEnergy < 0.0 );
303
304
305
306
307
309 gamEnergy1 = lv1.
e();
314 } else {
316 gamEnergy1 = 0.0;
317 }
319
320
321
322
323
328 fvect->push_back(dp);
329 } else { eKinEnergy = 0.0; }
330
333
334
335
336 if(nullptr != fAtomDeexcitation) {
345
346 for (
G4int j=nbefore; j<nafter; ++j) {
347 G4double e = ((*fvect)[j])->GetKineticEnergy();
348 if(esec + e > edep) {
349
350 e = edep - esec;
351 ((*fvect)[j])->SetKineticEnergy(e);
352 esec += e;
353
354
355
356
357
358
359
360
361
362
363
364 for (
G4int jj=nafter-1; jj>j; --jj) {
365 delete (*fvect)[jj];
366 fvect->pop_back();
367 }
368 break;
369 }
370 esec += e;
371 }
372 edep -= esec;
373 }
374 }
375 if(std::abs(energy - gamEnergy1 - eKinEnergy - esec - edep) > eV) {
376 G4cout <<
"### G4KleinNishinaModel dE(eV)= "
377 << (
energy - gamEnergy1 - eKinEnergy - esec - edep)/eV
378 << " shell= " << i
379 << " E(keV)= " << energy/keV
380 << " Ebind(keV)= " << bindingEnergy/keV
381 << " Eg(keV)= " << gamEnergy1/keV
382 << " Ee(keV)= " << eKinEnergy/keV
383 << " Esec(keV)= " << esec/keV
384 << " Edep(keV)= " << edep/keV
386 }
387
388 if(edep > 0.0) {
390 }
391}
G4double epsilon(G4double density, G4double temperature)
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
Hep3Vector boostVector() const
HepLorentzVector & boost(double, double, double)
void set(double x, double y, double z, double t)
virtual void flatArray(const int size, double *vect)=0
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
G4int GetNbOfAtomicShells() const
G4int GetNbOfShellElectrons(G4int index) const
G4double GetAtomicShell(G4int index) const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double energy(const ThreeVector &p, const G4double m)
G4double bindingEnergy(G4int A, G4int Z)