43 #ifndef G4DiffuseElastic_h
44 #define G4DiffuseElastic_h 1
46 #include <CLHEP/Units/PhysicalConstants.h>
280 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
282 modvalue = std::fabs(value);
284 if ( value < 8.0 && value > -8.0 )
286 value2 = value*value;
288 fact1 = 57568490574.0 + value2*(-13362590354.0
289 + value2*( 651619640.7
290 + value2*(-11214424.18
291 + value2*( 77392.33017
292 + value2*(-184.9052456 ) ) ) ) );
294 fact2 = 57568490411.0 + value2*( 1029532985.0
295 + value2*( 9494680.718
296 + value2*(59272.64853
297 + value2*(267.8532712
298 + value2*1.0 ) ) ) );
300 bessel = fact1/fact2;
308 shift = modvalue-0.785398164;
310 fact1 = 1.0 + value2*(-0.1098628627e-2
311 + value2*(0.2734510407e-4
312 + value2*(-0.2073370639e-5
313 + value2*0.2093887211e-6 ) ) );
315 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
316 + value2*(-0.6911147651e-5
317 + value2*(0.7621095161e-6
318 - value2*0.934945152e-7 ) ) );
320 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
332 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
334 modvalue = std::fabs(value);
336 if ( modvalue < 8.0 )
338 value2 = value*value;
340 fact1 = value*(72362614232.0 + value2*(-7895059235.0
341 + value2*( 242396853.1
342 + value2*(-2972611.439
343 + value2*( 15704.48260
344 + value2*(-30.16036606 ) ) ) ) ) );
346 fact2 = 144725228442.0 + value2*(2300535178.0
347 + value2*(18583304.74
348 + value2*(99447.43394
349 + value2*(376.9991397
350 + value2*1.0 ) ) ) );
351 bessel = fact1/fact2;
359 shift = modvalue - 2.356194491;
361 fact1 = 1.0 + value2*( 0.183105e-2
362 + value2*(-0.3516396496e-4
363 + value2*(0.2457520174e-5
364 + value2*(-0.240337019e-6 ) ) ) );
366 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
367 + value2*( 0.8449199096e-5
368 + value2*(-0.88228987e-6
369 + value2*0.105787412e-6 ) ) );
371 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
373 if (value < 0.0) bessel = -bessel;
389 if( std::fabs(x) < 0.01 )
391 df = 1./(1. + x/f2 + x*x/
f3 + x*x*x/
f4);
409 if( std::fabs(x) < 0.01 )
413 result = 2. - x2 + x2*x2/6.;
431 fBeta = a/std::sqrt(1+a*a);
442 fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
535 G4double sinHalfTheta = std::sin(0.5*theta);
536 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
544 G4double xsc = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
566 G4cout<<
"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<
G4endl;
601 xsc /= (1 - c1 + am)*(1 - c2 + am);
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
static c2_factory< G4double > c2
G4double IntegralElasticProb(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double CalculateNuclearRad(G4double A)
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
G4double GetCoulombTotalXsc(const G4ParticleDefinition *particle, G4double momentum, G4double Z)
G4double lowEnergyLimitHE
G4double BesselJzero(G4double z)
void InitialiseOnFly(G4double Z, G4double A)
G4double GetInvElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
std::vector< G4String > fElementNameVector
G4double NeutronTuniform(G4int Z)
G4ParticleDefinition * theProton
G4double GetDiffElasticSumProb(G4double theta)
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
G4double BesselJone(G4double z)
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double SampleThetaLab(const G4HadProjectile *aParticle, G4double tmass, G4double A)
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
static G4KaonMinus * KaonMinus()
void SetLowestEnergyLimit(G4double value)
const G4ParticleDefinition * thePionPlus
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
void SetRecoilKinEnergyLimit(G4double value)
const G4ParticleDefinition * thePionMinus
G4double A23(G4double A) const
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
const G4ParticleDefinition * GetDefinition() const
virtual ~G4DiffuseElastic()
G4double GetDiffElasticProb(G4double theta)
std::vector< G4PhysicsTable * > fAngleBank
void SetQModelLowLimit(G4double value)
static const double twopi
static G4Proton * Proton()
G4double GetNuclearRadius()
static G4PionPlus * PionPlus()
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double lowestEnergyLimit
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
static G4Neutron * Neutron()
G4ParticleDefinition * theAlpha
G4double DampFactor(G4double z)
G4PhysicsLogVector * fEnergyVector
G4double GetPDGMass() const
G4double A13(G4double A) const
G4double lowEnergyRecoilLimit
static G4PionMinus * PionMinus()
const G4double x[NPOINTSGL]
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double GetDiffElasticSumProbA(G4double alpha)
G4double GetCoulombIntegralXsc(const G4ParticleDefinition *particle, G4double momentum, G4double Z, G4double theta1, G4double theta2)
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
G4ParticleDefinition * theDeuteron
G4double BesselOneByArg(G4double z)
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
const G4ParticleDefinition * fParticle
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
void SetHEModelLowLimit(G4double value)
static G4KaonPlus * KaonPlus()
G4double GetPDGCharge() const
static const G4double alpha
std::vector< G4double > fElementNumberVector
void SetPlabLowLimit(G4double value)
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4PhysicsTable * fAngleTable
virtual G4bool IsApplicable(const G4HadProjectile &, G4Nucleus &)
static const double fermi
G4double GetIntegrandFunction(G4double theta)
G4ParticleDefinition * theNeutron