165{
166
168
169
171
173
174
176
177
179 if(nShells > (
G4int)fProbabilities.size()) { fProbabilities.resize(nShells); }
182 for(i=0; i<nShells; ++i) {
183
185
186 fProbabilities[i] = totprob;
187 }
188
189
190 static const G4int nlooplim = 1000;
192
195
198
199 do {
200 ++nloop;
201
202
205
206
207 for(i=0; i<nShells; ++i) { if(xprob <= fProbabilities[i]) { break; } }
208
210 lv1.
set(0.0,0.0,energy,energy);
211
212
213
214
215
216
220
221
222 G4double eTotMomentum = sqrt(eKinEnergy*(eKinEnergy + 2*electron_mass_c2));
225 G4double sintet = sqrt((1 - costet)*(1 + costet));
226 lv2.
set(eTotMomentum*sintet*cos(phi),eTotMomentum*sintet*sin(phi),
227 eTotMomentum*costet,eKinEnergy + electron_mass_c2);
230
231 gamEnergy0 = lv1.
e();
232
233
234
235
236
237 G4double E0_m = gamEnergy0/electron_mass_c2;
238
239
240
241
242
243
245
249 G4double alpha2 = alpha1 + 0.5*(1 - epsilon0sq);
250
251 do {
252 ++nloop;
253
254 if(nloop > nlooplim) { return; }
255
256
258
259 if ( alpha1 > alpha2*rndm[0] ) {
262
263 } else {
264 epsilonsq = epsilon0sq + (1.- epsilon0sq)*rndm[1];
266 }
267
269 sint2 = onecost*(2.-onecost);
270 greject = 1. -
epsilon*sint2/(1.+ epsilonsq);
271
272
273 } while (greject < rndm[2]);
274 gamEnergy1 =
epsilon*gamEnergy0;
275
276
277 lv2.
set(0.0,0.0,0.0,electron_mass_c2);
278 lv2 += lv1;
279
280
281
282
283 if(sint2 < 0.0) { sint2 = 0.0; }
284 costet = 1. - onecost;
285 sintet = sqrt(sint2);
286 phi = twopi * rndmEngineMod->
flat();
287
288
289
290
294 lv1.
set(gamEnergy1*v.
x(),gamEnergy1*v.
y(),gamEnergy1*v.
z(),gamEnergy1);
295 lv2 -= lv1;
296
297
299 eKinEnergy = lv2.
e() - electron_mass_c2 - ePotEnergy;
300
301
302
303 } while ( eKinEnergy < 0.0 );
304
305
306
307
308
310 gamEnergy1 = lv1.
e();
315 } else {
317 gamEnergy1 = 0.0;
318 }
320
321
322
323
324
330 fvect->push_back(dp);
331 } else { eKinEnergy = 0.0; }
332
335
336
337
338 if(fAtomDeexcitation) {
344 G4int nbefore = fvect->size();
346 G4int nafter = fvect->size();
347
348 for (
G4int j=nbefore; j<nafter; ++j) {
349 G4double e = ((*fvect)[j])->GetKineticEnergy();
350 if(esec + e > edep) {
351
352 e = edep - esec;
353 ((*fvect)[j])->SetKineticEnergy(e);
354 esec += e;
355
356
357
358
359
360
361
362
363
364
365
366 for (
G4int jj=nafter-1; jj>j; --jj) {
367 delete (*fvect)[jj];
368 fvect->pop_back();
369 }
370 break;
371 }
372 esec += e;
373 }
374 edep -= esec;
375 }
376 }
377 if(std::abs(energy - gamEnergy1 - eKinEnergy - esec - edep) > eV) {
378 G4cout <<
"### G4KleinNishinaModel dE(eV)= "
379 << (
energy - gamEnergy1 - eKinEnergy - esec - edep)/eV
380 << " shell= " << i
381 << " E(keV)= " << energy/keV
382 << " Ebind(keV)= " << bindingEnergy/keV
383 << " Eg(keV)= " << gamEnergy1/keV
384 << " Ee(keV)= " << eKinEnergy/keV
385 << " Esec(keV)= " << esec/keV
386 << " Edep(keV)= " << edep/keV
388 }
389
390 if(edep > 0.0) {
392 }
393}
double epsilon(double density, double 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(G4double Px, G4double Py, G4double Pz)
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)