328{
330
332 const G4double minKinEnergy = std::min(cut, tmax);
333 const G4double maxKinEnergy = std::min(maxEnergy, tmax);
334 if(minKinEnergy >= maxKinEnergy) { return; }
335
336
337
338
339 G4double totEnergy = kineticEnergy + mass;
340 G4double etot2 = totEnergy*totEnergy;
341 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
342
346 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
347
350
351
352 do {
354 deltaKinEnergy = minKinEnergy*maxKinEnergy
355 /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
356
357 f = 1.0 - beta2*deltaKinEnergy/tmax;
358 if( 0.0 < spin ) {
359 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
360 f += f1;
361 }
362
363
364 } while( fmax*rndm[1] > f);
365
366
367
368
369 G4double x = formfact*deltaKinEnergy;
370 if(x > 1.e-6) {
371
374 if( 0.0 < spin ) {
375 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
376 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
377 }
378 if(grej > 1.1) {
379 G4cout <<
"### G4LindhardSorensenIonModel WARNING: grej= " << grej
381 << " Ekin(MeV)= " << kineticEnergy
382 << " delEkin(MeV)= " << deltaKinEnergy
384 }
385 if(rndmEngineMod->
flat() > grej) {
return; }
386 }
387
389
391
394
395 deltaDirection =
397
398 } else {
399
401 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
402 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
404 cost = std::min(cost, 1.0);
405 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
406
408
409 deltaDirection.
set(sint*std::cos(phi),sint*std::sin(phi), cost) ;
411 }
412
413
414
415
416
417
418
419
420
421
422
423
425
426 vdp->push_back(delta);
427
428
429 kineticEnergy -= deltaKinEnergy;
431 finalP = finalP.
unit();
432
435}
G4GLOB_DLL std::ostream G4cout
void set(double x, double y, double z)
Hep3Vector & rotateUz(const Hep3Vector &)
virtual void flatArray(const int size, double *vect)=0
const G4ThreeVector & GetMomentumDirection() const
G4ThreeVector GetMomentum() const
G4double GetTotalMomentum() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
const G4String & GetParticleName() const
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
G4int SelectRandomAtomNumber(const G4Material *) const