G4NuclNuclDiffuseElastic Class Reference

#include <G4NuclNuclDiffuseElastic.hh>

Inheritance diagram for G4NuclNuclDiffuseElastic:

Public Member Functions
	G4NuclNuclDiffuseElastic ()
virtual	~G4NuclNuclDiffuseElastic ()
void	Initialise ()
void	InitialiseOnFly (G4double Z, G4double A)
void	BuildAngleTable ()
virtual G4double	SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
void	SetPlabLowLimit (G4double value)
void	SetHEModelLowLimit (G4double value)
void	SetQModelLowLimit (G4double value)
void	SetLowestEnergyLimit (G4double value)
void	SetRecoilKinEnergyLimit (G4double value)
G4double	SampleT (const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double	SampleTableT (const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double	SampleThetaCMS (const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double	SampleTableThetaCMS (const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double	GetScatteringAngle (G4int iMomentum, G4int iAngle, G4double position)
G4double	SampleThetaLab (const G4HadProjectile *aParticle, G4double tmass, G4double A)
G4double	GetDiffuseElasticXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double	GetInvElasticXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double	GetDiffuseElasticSumXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double	GetInvElasticSumXsc (const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double	IntegralElasticProb (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double	GetCoulombElasticXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
G4double	GetRutherfordXsc (G4double theta)
G4double	GetInvCoulombElasticXsc (const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double	GetCoulombTotalXsc (const G4ParticleDefinition *particle, G4double momentum, G4double Z)
G4double	GetCoulombIntegralXsc (const G4ParticleDefinition *particle, G4double momentum, G4double Z, G4double theta1, G4double theta2)
G4double	CalculateParticleBeta (const G4ParticleDefinition *particle, G4double momentum)
G4double	CalculateZommerfeld (G4double beta, G4double Z1, G4double Z2)
G4double	CalculateAm (G4double momentum, G4double n, G4double Z)
G4double	CalculateNuclearRad (G4double A)
G4double	ThetaCMStoThetaLab (const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
G4double	ThetaLabToThetaCMS (const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
void	TestAngleTable (const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
G4double	BesselJzero (G4double z)
G4double	BesselJone (G4double z)
G4double	DampFactor (G4double z)
G4double	BesselOneByArg (G4double z)
G4double	GetDiffElasticProb (G4double theta)
G4double	GetDiffElasticSumProb (G4double theta)
G4double	GetDiffElasticSumProbA (G4double alpha)
G4double	GetIntegrandFunction (G4double theta)
G4double	GetNuclearRadius ()
G4complex	GammaLogarithm (G4complex xx)
G4complex	GammaLogB2n (G4complex xx)
G4double	GetErf (G4double x)
G4double	GetCosHaPit2 (G4double t)
G4double	GetSinHaPit2 (G4double t)
G4double	GetCint (G4double x)
G4double	GetSint (G4double x)
G4complex	GetErfcComp (G4complex z, G4int nMax)
G4complex	GetErfcSer (G4complex z, G4int nMax)
G4complex	GetErfcInt (G4complex z)
G4complex	GetErfComp (G4complex z, G4int nMax)
G4complex	GetErfSer (G4complex z, G4int nMax)
G4double	GetExpCos (G4double x)
G4double	GetExpSin (G4double x)
G4complex	GetErfInt (G4complex z)
G4double	GetLegendrePol (G4int n, G4double x)
G4complex	TestErfcComp (G4complex z, G4int nMax)
G4complex	TestErfcSer (G4complex z, G4int nMax)
G4complex	TestErfcInt (G4complex z)
G4complex	CoulombAmplitude (G4double theta)
G4double	CoulombAmplitudeMod2 (G4double theta)
void	CalculateCoulombPhaseZero ()
G4double	CalculateCoulombPhase (G4int n)
void	CalculateRutherfordAnglePar ()
G4double	ProfileNear (G4double theta)
G4double	ProfileFar (G4double theta)
G4double	Profile (G4double theta)
G4complex	PhaseNear (G4double theta)
G4complex	PhaseFar (G4double theta)
G4complex	GammaLess (G4double theta)
G4complex	GammaMore (G4double theta)
G4complex	AmplitudeNear (G4double theta)
G4complex	AmplitudeFar (G4double theta)
G4complex	Amplitude (G4double theta)
G4double	AmplitudeMod2 (G4double theta)
G4complex	AmplitudeSim (G4double theta)
G4double	AmplitudeSimMod2 (G4double theta)
G4double	GetRatioSim (G4double theta)
G4double	GetRatioGen (G4double theta)
G4double	GetFresnelDiffuseXsc (G4double theta)
G4double	GetFresnelIntegrandXsc (G4double alpha)
G4complex	AmplitudeGla (G4double theta)
G4double	AmplitudeGlaMod2 (G4double theta)
G4complex	AmplitudeGG (G4double theta)
G4double	AmplitudeGGMod2 (G4double theta)
void	InitParameters (const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
void	InitDynParameters (const G4ParticleDefinition *theParticle, G4double partMom)
void	InitParametersGla (const G4DynamicParticle *aParticle, G4double partMom, G4double Z, G4double A)
G4double	GetHadronNucleonXscNS (G4ParticleDefinition *pParticle, G4double pTkin, G4ParticleDefinition *tParticle)
G4double	CalcMandelstamS (const G4double mp, const G4double mt, const G4double Plab)
G4double	GetProfileLambda ()
void	SetProfileLambda (G4double pl)
void	SetProfileDelta (G4double pd)
void	SetProfileAlpha (G4double pa)
void	SetCofLambda (G4double pa)
void	SetCofAlpha (G4double pa)
void	SetCofAlphaMax (G4double pa)
void	SetCofAlphaCoulomb (G4double pa)
void	SetCofDelta (G4double pa)
void	SetCofPhase (G4double pa)
void	SetCofFar (G4double pa)
void	SetEtaRatio (G4double pa)
void	SetMaxL (G4int l)
void	SetNuclearRadiusCof (G4double r)
G4double	GetCofAlphaMax ()
G4double	GetCofAlphaCoulomb ()
Public Member Functions inherited from G4HadronElastic
	G4HadronElastic (const G4String &name="hElasticLHEP")
virtual	~G4HadronElastic ()
virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void	SetLowestEnergyLimit (G4double value)
G4double	LowestEnergyLimit () const
G4double	ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
virtual void	Description () const
Public Member Functions inherited from G4HadronicInteraction
	G4HadronicInteraction (const G4String &modelName="HadronicModel")
virtual	~G4HadronicInteraction ()
virtual G4bool	IsApplicable (const G4HadProjectile &, G4Nucleus &)
G4double	GetMinEnergy () const
G4double	GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
void	SetMinEnergy (G4double anEnergy)
void	SetMinEnergy (G4double anEnergy, const G4Element *anElement)
void	SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
G4double	GetMaxEnergy () const
G4double	GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
void	SetMaxEnergy (const G4double anEnergy)
void	SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
void	SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
const G4HadronicInteraction *	GetMyPointer () const
virtual G4int	GetVerboseLevel () const
virtual void	SetVerboseLevel (G4int value)
const G4String &	GetModelName () const
void	DeActivateFor (const G4Material *aMaterial)
void	ActivateFor (const G4Material *aMaterial)
void	DeActivateFor (const G4Element *anElement)
void	ActivateFor (const G4Element *anElement)
G4bool	IsBlocked (const G4Material *aMaterial) const
G4bool	IsBlocked (const G4Element *anElement) const
void	SetRecoilEnergyThreshold (G4double val)
G4double	GetRecoilEnergyThreshold () const
G4bool	operator== (const G4HadronicInteraction &right) const
G4bool	operator!= (const G4HadronicInteraction &right) const
virtual const std::pair < G4double, G4double >	GetFatalEnergyCheckLevels () const
virtual std::pair< G4double, G4double >	GetEnergyMomentumCheckLevels () const
void	SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
virtual void	ModelDescription (std::ostream &outFile) const

Additional Inherited Members
Protected Member Functions inherited from G4HadronicInteraction
void	SetModelName (const G4String &nam)
G4bool	IsBlocked () const
void	Block ()
Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState	theParticleChange
G4int	verboseLevel
G4double	theMinEnergy
G4double	theMaxEnergy
G4bool	isBlocked

Detailed Description

Definition at line 59 of file G4NuclNuclDiffuseElastic.hh.

Constructor & Destructor Documentation

G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic

(

)

Definition at line 67 of file G4NuclNuclDiffuseElastic.cc.

References G4Alpha::Alpha(), G4Deuteron::Deuteron(), python.hepunit::GeV, python.hepunit::keV, python.hepunit::MeV, G4Neutron::Neutron(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4Proton::Proton(), G4HadronicInteraction::SetMaxEnergy(), G4HadronicInteraction::SetMinEnergy(), python.hepunit::TeV, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, and G4HadronicInteraction::verboseLevel.

  : G4HadronElastic("NNDiffuseElastic"), fParticle(0)
{
  SetMinEnergy( 50*MeV );
  SetMaxEnergy( 1.*TeV );
  verboseLevel = 0;
  lowEnergyRecoilLimit = 100.*keV;  
  lowEnergyLimitQ  = 0.0*GeV;  
  lowEnergyLimitHE = 0.0*GeV;  
  lowestEnergyLimit= 0.0*keV;  
  plabLowLimit     = 20.0*MeV;

  theProton   = G4Proton::Proton();
  theNeutron  = G4Neutron::Neutron();
  theDeuteron = G4Deuteron::Deuteron();
  theAlpha    = G4Alpha::Alpha();
  thePionPlus = G4PionPlus::PionPlus();
  thePionMinus= G4PionMinus::PionMinus();

  fEnergyBin = 200;
  fAngleBin  = 200;

  fEnergyVector =  new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
  fAngleTable = 0;

  fParticle = 0;
  fWaveVector = 0.;
  fAtomicWeight = 0.;
  fAtomicNumber = 0.;
  fNuclearRadius = 0.;
  fBeta = 0.;
  fZommerfeld = 0.;
  fAm = 0.;
  fAddCoulomb = false;
  // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle

  // Empirical parameters

  fCofAlphaMax     = 1.5;
  fCofAlphaCoulomb = 0.5;

  fProfileDelta  = 1.;
  fProfileAlpha  = 0.5;

  fCofLambda = 1.0;
  fCofDelta  = 0.04;
  fCofAlpha  = 0.095;

  fNuclearRadius1 = fNuclearRadius2 = fNuclearRadiusSquare = fNuclearRadiusCof 
    = fRutherfordRatio = fCoulombPhase0 = fHalfRutThetaTg = fHalfRutThetaTg2 
    = fRutherfordTheta = fProfileLambda = fCofPhase = fCofFar = fCofAlphaMax 
    = fCofAlphaCoulomb = fSumSigma = fEtaRatio = fReZ = 0.0;
  fMaxL = 0;

}

G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseElastic ( )

virtual

Definition at line 127 of file G4NuclNuclDiffuseElastic.cc.

References G4PhysicsTable::clearAndDestroy().

{
  if(fEnergyVector) delete fEnergyVector;

  if( fAngleTable )
  {
      fAngleTable->clearAndDestroy();
      delete fAngleTable ;
  }
}

Member Function Documentation

G4complex G4NuclNuclDiffuseElastic::Amplitude ( G4double  theta )

inline

Definition at line 973 of file G4NuclNuclDiffuseElastic.hh.

{
 
  G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
  // G4complex out = AmplitudeNear(theta);
  // G4complex out = AmplitudeFar(theta);
  return    out;
}

G4complex G4NuclNuclDiffuseElastic::AmplitudeFar ( G4double  theta )

inline

Definition at line 959 of file G4NuclNuclDiffuseElastic.hh.

{
  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
  G4complex out = G4complex(kappa/fWaveVector,0.);
  out *= ProfileFar(theta);
  out *= PhaseFar(theta);
  return out;
}

G4complex G4NuclNuclDiffuseElastic::AmplitudeGG

(

G4double

theta

)

Definition at line 1640 of file G4NuclNuclDiffuseElastic.cc.

References test::a, CoulombAmplitude(), G4cout, G4endl, and n.

{
  G4int n;
  G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
  G4double sinThetaH2 = sinThetaH*sinThetaH;
  G4complex out = G4complex(0.,0.); 
  G4complex im  = G4complex(0.,1.);

  a  = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
  b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;

  aTemp = a;

  for( n = 1; n < fMaxL; n++)
  {    
    T12b   = aTemp*std::exp(-b2/n)/n;         
    aTemp *= a;  
    out   += T12b; 
    G4cout<<"out = "<<out<<G4endl;  
  }
  out *= -4.*im*fWaveVector/CLHEP::pi;
  out += CoulombAmplitude(theta);
  return out;
}

G4double G4NuclNuclDiffuseElastic::AmplitudeGGMod2 ( G4double  theta )

inline

Definition at line 1064 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex out = AmplitudeGG(theta);
  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
  return   mod2;
}

G4complex G4NuclNuclDiffuseElastic::AmplitudeGla

(

G4double

theta

)

Definition at line 1613 of file G4NuclNuclDiffuseElastic.cc.

References test::b, CalculateCoulombPhase(), CoulombAmplitude(), GetLegendrePol(), and n.

{
  G4int n;
  G4double T12b, b, b2; // cosTheta = std::cos(theta);
  G4complex out = G4complex(0.,0.), shiftC, shiftN; 
  G4complex im  = G4complex(0.,1.);

  for( n = 0; n < fMaxL; n++)
  {
    shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
    // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
    b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
    b2 = b*b;
    T12b = fSumSigma*std::exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;         
    shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
    out +=  (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);   
  }
  out /= 2.*im*fWaveVector;
  out += CoulombAmplitude(theta);
  return out;
}

G4double G4NuclNuclDiffuseElastic::AmplitudeGlaMod2 ( G4double  theta )

inline

Definition at line 1053 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex out = AmplitudeGla(theta);
  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
  return   mod2;
}

G4double G4NuclNuclDiffuseElastic::AmplitudeMod2 ( G4double  theta )

inline

Definition at line 986 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex out = Amplitude(theta);
  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
  return   mod2;
}

G4complex G4NuclNuclDiffuseElastic::AmplitudeNear

(

G4double

theta

)

Definition at line 1557 of file G4NuclNuclDiffuseElastic.cc.

References CoulombAmplitude(), GammaLess(), GammaMore(), PhaseNear(), and ProfileNear().

{
  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
  G4complex out = G4complex(kappa/fWaveVector,0.);

  out *= PhaseNear(theta);

  if( theta <= fRutherfordTheta )
  {
    out *= GammaLess(theta) + ProfileNear(theta);
    // out *= GammaMore(theta) + ProfileNear(theta);
    out += CoulombAmplitude(theta);
  }
  else
  {
    out *= GammaMore(theta) + ProfileNear(theta);
    // out *= GammaLess(theta) + ProfileNear(theta);
  }
  return out;
}

G4complex G4NuclNuclDiffuseElastic::AmplitudeSim

(

G4double

theta

)

Definition at line 1582 of file G4NuclNuclDiffuseElastic.cc.

References CoulombAmplitude(), GetErfcInt(), and ProfileNear().

{
  G4double sinThetaR  = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
  G4double dTheta     = 0.5*(theta - fRutherfordTheta);
  G4double sindTheta  = std::sin(dTheta);
  G4double persqrt2   = std::sqrt(0.5);

  G4complex order     = G4complex(persqrt2,persqrt2);
  order              *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
  // order              *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;

  G4complex out;

  if ( theta <= fRutherfordTheta )
  {
    out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta);
  }
  else
  {
    out = 0.5*GetErfcInt(order)*ProfileNear(theta);
  }

  out *= CoulombAmplitude(theta);

  return out;
}

G4double G4NuclNuclDiffuseElastic::AmplitudeSimMod2 ( G4double  theta )

inline

Definition at line 1041 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex out = AmplitudeSim(theta);
  G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
  return   mod2;
}

G4double G4NuclNuclDiffuseElastic::BesselJone

(

G4double

z

)

Definition at line 2034 of file G4NuclNuclDiffuseElastic.cc.

Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().

{
  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;

  modvalue = std::fabs(value);

  if ( modvalue < 8.0 ) 
  {
    value2 = value*value;

    fact1  = value*(72362614232.0 + value2*(-7895059235.0 
                                  + value2*( 242396853.1
                                  + value2*(-2972611.439 
                                  + value2*( 15704.48260 
                                  + value2*(-30.16036606  ) ) ) ) ) );

    fact2  = 144725228442.0 + value2*(2300535178.0 
                            + value2*(18583304.74
                            + value2*(99447.43394 
                            + value2*(376.9991397 
                            + value2*1.0             ) ) ) );
    bessel = fact1/fact2;
  } 
  else 
  {
    arg    = 8.0/modvalue;

    value2 = arg*arg;

    shift  = modvalue - 2.356194491;

    fact1  = 1.0 + value2*( 0.183105e-2 
                 + value2*(-0.3516396496e-4
                 + value2*(0.2457520174e-5 
                 + value2*(-0.240337019e-6          ) ) ) );

    fact2 = 0.04687499995 + value2*(-0.2002690873e-3
                          + value2*( 0.8449199096e-5
                          + value2*(-0.88228987e-6
                          + value2*0.105787412e-6       ) ) );

    bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);

    if (value < 0.0) bessel = -bessel;
  }
  return bessel;
}

G4double G4NuclNuclDiffuseElastic::BesselJzero

(

G4double

z

)

Definition at line 1982 of file G4NuclNuclDiffuseElastic.cc.

Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().

{
  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;

  modvalue = std::fabs(value);

  if ( value < 8.0 && value > -8.0 )
  {
    value2 = value*value;

    fact1  = 57568490574.0 + value2*(-13362590354.0 
                           + value2*( 651619640.7 
                           + value2*(-11214424.18 
                           + value2*( 77392.33017 
                           + value2*(-184.9052456   ) ) ) ) );

    fact2  = 57568490411.0 + value2*( 1029532985.0 
                           + value2*( 9494680.718
                           + value2*(59272.64853
                           + value2*(267.8532712 
                           + value2*1.0               ) ) ) );

    bessel = fact1/fact2;
  } 
  else 
  {
    arg    = 8.0/modvalue;

    value2 = arg*arg;

    shift  = modvalue-0.785398164;

    fact1  = 1.0 + value2*(-0.1098628627e-2 
                 + value2*(0.2734510407e-4
                 + value2*(-0.2073370639e-5 
                 + value2*0.2093887211e-6    ) ) );

    fact2  = -0.1562499995e-1 + value2*(0.1430488765e-3
                              + value2*(-0.6911147651e-5 
                              + value2*(0.7621095161e-6
                              - value2*0.934945152e-7    ) ) );

    bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
  }
  return bessel;
}

G4double G4NuclNuclDiffuseElastic::BesselOneByArg ( G4double  z )

inline

Definition at line 414 of file G4NuclNuclDiffuseElastic.hh.

References test::x.

Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().

{
  G4double x2, result;
  
  if( std::fabs(x) < 0.01 )
  { 
   x     *= 0.5;
   x2     = x*x;
   result = 2. - x2 + x2*x2/6.;
  }
  else
  {
    result = BesselJone(x)/x; 
  }
  return result;
}

void G4NuclNuclDiffuseElastic::BuildAngleTable

(

)

Definition at line 958 of file G4NuclNuclDiffuseElastic.cc.

References GetFresnelIntegrandXsc(), G4PhysicsVector::GetLowEdgeEnergy(), G4ParticleDefinition::GetPDGMass(), InitDynParameters(), G4PhysicsTable::insertAt(), G4Integrator< T, F >::Legendre10(), python.hepunit::pi, and G4PhysicsFreeVector::PutValue().

Referenced by Initialise(), and InitialiseOnFly().

{
  G4int i, j;
  G4double partMom, kinE, m1 = fParticle->GetPDGMass();
  G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;

  // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl;

  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
  
  fAngleTable = new G4PhysicsTable(fEnergyBin);

  for( i = 0; i < fEnergyBin; i++)
  {
    kinE        = fEnergyVector->GetLowEdgeEnergy(i);

    // G4cout<<G4endl;
    // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl;

    partMom     = std::sqrt( kinE*(kinE + 2*m1) );

    InitDynParameters(fParticle, partMom);

    alphaMax = fRutherfordTheta*fCofAlphaMax;

    if(alphaMax > pi) alphaMax = pi;


    alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb;

    // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl;

    G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);

    // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );

    G4double delth = (alphaMax-alphaCoulomb)/fAngleBin;
        
    sum = 0.;

    // fAddCoulomb = false;
    fAddCoulomb = true;

    // for(j = 1; j < fAngleBin; j++)
    for(j = fAngleBin-1; j >= 1; j--)
    {
      // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
      // alpha2 = angleBins->GetLowEdgeEnergy(j);

      // alpha1 = alphaMax*(j-1)/fAngleBin;
      // alpha2 = alphaMax*( j )/fAngleBin;

      alpha1 = alphaCoulomb + delth*(j-1);
      // if(alpha1 < kRlim2) alpha1 = kRlim2;
      alpha2 = alpha1 + delth;

      delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc, alpha1, alpha2);
      // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);

      sum += delta;
      
      angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
      // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl;
    }
    fAngleTable->insertAt(i,angleVector);

    // delete[] angleVector; 
    // delete[] angleBins; 
  }
  return;
}

G4double G4NuclNuclDiffuseElastic::CalcMandelstamS ( const G4double  mp, const G4double  mt, const G4double  Plab  )

inline

Definition at line 1075 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetHadronNucleonXscNS().

{
  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
  G4double sMand  = mp*mp + mt*mt + 2*Elab*mt ;

  return sMand;
}

G4double G4NuclNuclDiffuseElastic::CalculateAm ( G4double  momentum, G4double  n, G4double  Z  )

inline

Definition at line 463 of file G4NuclNuclDiffuseElastic.hh.

References n.

Referenced by GetDiffuseElasticSumXsc(), InitDynParameters(), InitParameters(), InitParametersGla(), and TestAngleTable().

{
  G4double k   = momentum/CLHEP::hbarc;
  G4double ch  = 1.13 + 3.76*n*n;
  G4double zn  = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
  G4double zn2 = zn*zn;
  fAm          = ch/zn2;

  return fAm;
}

G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase ( G4int  n )

inline

Definition at line 839 of file G4NuclNuclDiffuseElastic.hh.

References z.

Referenced by AmplitudeGla().

{
  G4complex z        = G4complex(1. + n, fZommerfeld); 
  // G4complex gammalog = GammaLogarithm(z);
  G4complex gammalog = GammaLogB2n(z);
  return gammalog.imag();
}

void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero ( )

inline

Definition at line 826 of file G4NuclNuclDiffuseElastic.hh.

References z.

Referenced by InitDynParameters(), InitParameters(), and InitParametersGla().

{
  G4complex z        = G4complex(1,fZommerfeld); 
  // G4complex gammalog = GammaLogarithm(z);
  G4complex gammalog = GammaLogB2n(z);
  fCoulombPhase0     = gammalog.imag();
}

G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad ( G4double  A )

inline

Definition at line 478 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetDiffuseElasticSumXsc(), GetDiffuseElasticXsc(), Initialise(), InitialiseOnFly(), InitParameters(), InitParametersGla(), IntegralElasticProb(), SampleThetaCMS(), and TestAngleTable().

{
  G4double r0 = 1.*CLHEP::fermi, radius;
  // r0 *= 1.12;
  // r0 *= 1.44;
  r0 *= fNuclearRadiusCof;

  /*
  if( A < 50. )
  {
    if( A > 10. ) r0  = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi;   // 1.08*fermi;
    else          r0  = 1.1*CLHEP::fermi;

    radius = r0*std::pow(A, 1./3.);
  }
  else
  {
    r0 = 1.7*CLHEP::fermi;   // 1.7*fermi;

    radius = r0*std::pow(A, 0.27); // 0.27);
  }
  */
  radius = r0*std::pow(A, 1./3.);

  return radius;
}

G4double G4NuclNuclDiffuseElastic::CalculateParticleBeta ( const G4ParticleDefinition *  particle, G4double  momentum  )

inline

Definition at line 436 of file G4NuclNuclDiffuseElastic.hh.

References test::a, and G4ParticleDefinition::GetPDGMass().

Referenced by GetDiffuseElasticSumXsc().

{
  G4double mass = particle->GetPDGMass();
  G4double a    = momentum/mass;
  fBeta         = a/std::sqrt(1+a*a);

  return fBeta; 
}

void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar ( )

inline

Definition at line 853 of file G4NuclNuclDiffuseElastic.hh.

Referenced by InitDynParameters(), and InitParameters().

{
  fHalfRutThetaTg   = fZommerfeld/fProfileLambda;  // (fWaveVector*fNuclearRadius);
  fRutherfordTheta  = 2.*std::atan(fHalfRutThetaTg);
  fHalfRutThetaTg2  = fHalfRutThetaTg*fHalfRutThetaTg;
  // G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;

}

G4double G4NuclNuclDiffuseElastic::CalculateZommerfeld ( G4double  beta, G4double  Z1, G4double  Z2  )

inline

Definition at line 451 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetDiffuseElasticSumXsc(), InitDynParameters(), InitParameters(), InitParametersGla(), and TestAngleTable().

{
  fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;

  return fZommerfeld; 
}

G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude ( G4double  theta )

inline

Definition at line 792 of file G4NuclNuclDiffuseElastic.hh.

References z.

Referenced by AmplitudeGG(), AmplitudeGla(), AmplitudeNear(), and AmplitudeSim().

{
  G4complex ca;
  
  G4double sinHalfTheta  = std::sin(0.5*theta);
  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 
  sinHalfTheta2         += fAm;

  G4double order         = 2.*fCoulombPhase0 - fZommerfeld*std::log(sinHalfTheta2);
  G4complex z            = G4complex(0., order);
  ca                     = std::exp(z);

  ca                    *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);

  return ca; 
}

G4double G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2 ( G4double  theta )

inline

Definition at line 813 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex ca = CoulombAmplitude(theta);
  G4double out = ca.real()*ca.real() + ca.imag()*ca.imag();

  return  out;
}

G4double G4NuclNuclDiffuseElastic::DampFactor ( G4double  z )

inline

Definition at line 391 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().

{
  G4double df;
  G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials

  // x *= pi;

  if( std::fabs(x) < 0.01 )
  { 
    df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
  }
  else
  {
    df = x/std::sinh(x); 
  }
  return df;
}

G4complex G4NuclNuclDiffuseElastic::GammaLess

(

G4double

theta

)

Definition at line 1502 of file G4NuclNuclDiffuseElastic.cc.

References GetErfcInt().

Referenced by AmplitudeNear().

{
  G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);

  G4double u              = std::sqrt(0.5*fProfileLambda/sinThetaR);
  G4double kappa          = u/std::sqrt(CLHEP::pi);
  G4double dTheta         = theta - fRutherfordTheta;
  u                      *= dTheta;
  G4double u2             = u*u;
  G4double u2m2p3         = u2*2./3.;

  G4complex im            = G4complex(0.,1.);
  G4complex order         = G4complex(u,u);
  order                  /= std::sqrt(2.);

  G4complex gamma         = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi));
  G4complex a0            = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
  G4complex a1            = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
  G4complex out           = gamma*(1. - a1*dTheta) - a0;

  return out;
}

G4complex G4NuclNuclDiffuseElastic::GammaLogarithm

(

G4complex

xx

)

Definition at line 1958 of file G4NuclNuclDiffuseElastic.cc.

References z.

{
  const G4double cof[6] = { 76.18009172947146,     -86.50532032941677,
                             24.01409824083091,      -1.231739572450155,
                              0.1208650973866179e-2, -0.5395239384953e-5  } ;
  register G4int j;
  G4complex z = zz - 1.0;
  G4complex tmp = z + 5.5;
  tmp -= (z + 0.5) * std::log(tmp);
  G4complex ser = G4complex(1.000000000190015,0.);

  for ( j = 0; j <= 5; j++ )
  {
    z += 1.0;
    ser += cof[j]/z;
  }
  return -tmp + std::log(2.5066282746310005*ser);
}

G4complex G4NuclNuclDiffuseElastic::GammaLogB2n ( G4complex  xx )

inline

Definition at line 605 of file G4NuclNuclDiffuseElastic.hh.

References z.

{
  G4complex z1 = 12.*z;
  G4complex z2 = z*z;
  G4complex z3 = z2*z;
  G4complex z5 = z2*z3;
  G4complex z7 = z2*z5;

  z3 *= 360.;
  z5 *= 1260.;
  z7 *= 1680.;

  G4complex result  = (z-0.5)*std::log(z)-z+0.5*std::log(CLHEP::twopi);
            result += 1./z1 - 1./z3 +1./z5 -1./z7;
  return result;
}

G4complex G4NuclNuclDiffuseElastic::GammaMore

(

G4double

theta

)

Definition at line 1530 of file G4NuclNuclDiffuseElastic.cc.

References GetErfcInt().

Referenced by AmplitudeNear().

{
  G4double sinThetaR      = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);

  G4double u              = std::sqrt(0.5*fProfileLambda/sinThetaR);
  G4double kappa          = u/std::sqrt(CLHEP::pi);
  G4double dTheta         = theta - fRutherfordTheta;
  u                      *= dTheta;
  G4double u2             = u*u;
  G4double u2m2p3         = u2*2./3.;

  G4complex im            = G4complex(0.,1.);
  G4complex order         = G4complex(u,u);
  order                  /= std::sqrt(2.);
  G4complex gamma         = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi));
  G4complex a0            = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
  G4complex a1            = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
  G4complex out           = -gamma*(1. - a1*dTheta) - a0;

  return out;
}

G4double G4NuclNuclDiffuseElastic::GetCint ( G4double  x )

inline

Definition at line 762 of file G4NuclNuclDiffuseElastic.hh.

References GetCosHaPit2(), and G4Integrator< T, F >::Legendre96().

Referenced by GetRatioGen().

{
  G4double out;

  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;

  out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetCosHaPit2, 0., x );

  return out;
}

G4double G4NuclNuclDiffuseElastic::GetCofAlphaCoulomb ( )

inline

Definition at line 289 of file G4NuclNuclDiffuseElastic.hh.

{return fCofAlphaCoulomb;};

G4double G4NuclNuclDiffuseElastic::GetCofAlphaMax ( )

inline

Definition at line 288 of file G4NuclNuclDiffuseElastic.hh.

{return fCofAlphaMax;};

G4double G4NuclNuclDiffuseElastic::GetCosHaPit2 ( G4double  t )

inline

Definition at line 191 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetCint().

{return std::cos(CLHEP::halfpi*t*t);};

G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc ( const G4ParticleDefinition *  particle, G4double  theta, G4double  momentum, G4double  Z  )

inline

Definition at line 509 of file G4NuclNuclDiffuseElastic.hh.

References G4ParticleDefinition::GetPDGCharge(), n, and z.

Referenced by GetInvCoulombElasticXsc().

{
  G4double sinHalfTheta  = std::sin(0.5*theta);
  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
  G4double beta          = CalculateParticleBeta( particle, momentum);
  G4double z             = particle->GetPDGCharge();
  G4double n             = CalculateZommerfeld( beta, z, Z );
  G4double am            = CalculateAm( momentum, n, Z);
  G4double k             = momentum/CLHEP::hbarc;
  G4double ch            = 0.5*n/k;
  G4double ch2           = ch*ch;
  G4double xsc           = ch2/((sinHalfTheta2+am)*(sinHalfTheta2+am));

  return xsc;
}

G4double G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc ( const G4ParticleDefinition *  particle, G4double  momentum, G4double  Z, G4double  theta1, G4double  theta2  )

inline

Definition at line 574 of file G4NuclNuclDiffuseElastic.hh.

References plottest35::c1, G4ParticleDefinition::GetPDGCharge(), n, and z.

{
  G4double c1 = std::cos(theta1);
  //G4cout<<"c1 = "<<c1<<G4endl;
  G4double c2 = std::cos(theta2);
  // G4cout<<"c2 = "<<c2<<G4endl;
  G4double beta          = CalculateParticleBeta( particle, momentum);
  // G4cout<<"beta = "<<beta<<G4endl;
  G4double z             = particle->GetPDGCharge();
  G4double n             = CalculateZommerfeld( beta, z, Z );
  // G4cout<<"fZomerfeld = "<<n<<G4endl;
  G4double am            = CalculateAm( momentum, n, Z);
  // G4cout<<"cof Am = "<<am<<G4endl;
  G4double k             = momentum/CLHEP::hbarc;
  // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
  // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
  G4double ch            = n/k;
  G4double ch2           = ch*ch;
  am *= 2.;
  G4double xsc = ch2*CLHEP::twopi*(c1-c2)/((1 - c1 + am)*(1 - c2 + am));

  return xsc;
}

G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc ( const G4ParticleDefinition *  particle, G4double  momentum, G4double  Z  )

inline

Definition at line 548 of file G4NuclNuclDiffuseElastic.hh.

References G4cout, G4endl, G4ParticleDefinition::GetPDGCharge(), n, and z.

{
  G4double beta          = CalculateParticleBeta( particle, momentum);
  G4cout<<"beta = "<<beta<<G4endl;
  G4double z             = particle->GetPDGCharge();
  G4double n             = CalculateZommerfeld( beta, z, Z );
  G4cout<<"fZomerfeld = "<<n<<G4endl;
  G4double am            = CalculateAm( momentum, n, Z);
  G4cout<<"cof Am = "<<am<<G4endl;
  G4double k             = momentum/CLHEP::hbarc;
  G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
  G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
  G4double ch            = n/k;
  G4double ch2           = ch*ch;
  G4double xsc           = ch2*CLHEP::pi/(am +am*am);

  return xsc;
}

G4double G4NuclNuclDiffuseElastic::GetDiffElasticProb

(

G4double

theta

)

Definition at line 389 of file G4NuclNuclDiffuseElastic.cc.

References BesselJone(), BesselJzero(), BesselOneByArg(), DampFactor(), python.hepunit::fermi, G4InuclParticleNames::lambda, and python.hepunit::pi.

Referenced by GetDiffuseElasticXsc().

{
  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
  G4double delta, diffuse, gamma;
  G4double e1, e2, bone, bone2;

  // G4double wavek = momentum/hbarc;  // wave vector
  // G4double r0    = 1.08*fermi;
  // G4double rad   = r0*std::pow(A, 1./3.);

  G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
  G4double kr2   = kr*kr;
  G4double krt   = kr*theta;

  bzero      = BesselJzero(krt);
  bzero2     = bzero*bzero;    
  bone       = BesselJone(krt);
  bone2      = bone*bone;
  bonebyarg  = BesselOneByArg(krt);
  bonebyarg2 = bonebyarg*bonebyarg;  

  if (fParticle == theProton)
  {
    diffuse = 0.63*fermi;
    gamma   = 0.3*fermi;
    delta   = 0.1*fermi*fermi;
    e1      = 0.3*fermi;
    e2      = 0.35*fermi;
  }
  else // as proton, if were not defined 
  {
    diffuse = 0.63*fermi;
    gamma   = 0.3*fermi;
    delta   = 0.1*fermi*fermi;
    e1      = 0.3*fermi;
    e2      = 0.35*fermi;
  }
  G4double lambda = 15.; // 15 ok

  //  G4double kgamma    = fWaveVector*gamma;   // wavek*delta;

  G4double kgamma    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
  G4double kgamma2   = kgamma*kgamma;

  // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
  // G4double dk2t2 = dk2t*dk2t;
  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;

  G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;

  damp           = DampFactor(pikdt);
  damp2          = damp*damp;

  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
  G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;


  sigma  = kgamma2;
  // sigma  += dk2t2;
  sigma *= bzero2;
  sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
  sigma += kr2*bonebyarg2;
  sigma *= damp2;          // *rad*rad;

  return sigma;
}

G4double G4NuclNuclDiffuseElastic::GetDiffElasticSumProb

(

G4double

theta

)

Definition at line 466 of file G4NuclNuclDiffuseElastic.cc.

References BesselJone(), BesselJzero(), BesselOneByArg(), DampFactor(), python.hepunit::fermi, G4InuclParticleNames::lambda, and python.hepunit::pi.

Referenced by GetDiffuseElasticSumXsc().

{
  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
  G4double delta, diffuse, gamma;
  G4double e1, e2, bone, bone2;

  // G4double wavek = momentum/hbarc;  // wave vector
  // G4double r0    = 1.08*fermi;
  // G4double rad   = r0*std::pow(A, 1./3.);

  G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
  G4double kr2   = kr*kr;
  G4double krt   = kr*theta;

  bzero      = BesselJzero(krt);
  bzero2     = bzero*bzero;    
  bone       = BesselJone(krt);
  bone2      = bone*bone;
  bonebyarg  = BesselOneByArg(krt);
  bonebyarg2 = bonebyarg*bonebyarg;  

  if (fParticle == theProton)
  {
    diffuse = 0.63*fermi;
    // diffuse = 0.6*fermi;
    gamma   = 0.3*fermi;
    delta   = 0.1*fermi*fermi;
    e1      = 0.3*fermi;
    e2      = 0.35*fermi;
  }
  else // as proton, if were not defined 
  {
    diffuse = 0.63*fermi;
    gamma   = 0.3*fermi;
    delta   = 0.1*fermi*fermi;
    e1      = 0.3*fermi;
    e2      = 0.35*fermi;
  }
  G4double lambda = 15.; // 15 ok
  // G4double kgamma    = fWaveVector*gamma;   // wavek*delta;
  G4double kgamma    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;

  // G4cout<<"kgamma = "<<kgamma<<G4endl;

  if(fAddCoulomb)  // add Coulomb correction
  {
    G4double sinHalfTheta  = std::sin(0.5*theta);
    G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;

    kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
  }

  G4double kgamma2   = kgamma*kgamma;

  // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
  //   G4cout<<"dk2t = "<<dk2t<<G4endl;
  // G4double dk2t2 = dk2t*dk2t;
  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;

  G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;

  // G4cout<<"pikdt = "<<pikdt<<G4endl;

  damp           = DampFactor(pikdt);
  damp2          = damp*damp;

  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
  G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;

  sigma  = kgamma2;
  // sigma += dk2t2;
  sigma *= bzero2;
  sigma += mode2k2*bone2; 
  sigma += e2dk3t*bzero*bone;

  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
  sigma += kr2*bonebyarg2;  // correction at J1()/()

  sigma *= damp2;          // *rad*rad;

  return sigma;
}

G4double G4NuclNuclDiffuseElastic::GetDiffElasticSumProbA

(

G4double

alpha

)

Definition at line 561 of file G4NuclNuclDiffuseElastic.cc.

References BesselJone(), BesselJzero(), BesselOneByArg(), DampFactor(), python.hepunit::fermi, G4InuclParticleNames::lambda, and python.hepunit::pi.

Referenced by GetIntegrandFunction().

{
  G4double theta; 

  theta = std::sqrt(alpha);

  // theta = std::acos( 1 - alpha/2. );

  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
  G4double delta, diffuse, gamma;
  G4double e1, e2, bone, bone2;

  // G4double wavek = momentum/hbarc;  // wave vector
  // G4double r0    = 1.08*fermi;
  // G4double rad   = r0*std::pow(A, 1./3.);

  G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
  G4double kr2   = kr*kr;
  G4double krt   = kr*theta;

  bzero      = BesselJzero(krt);
  bzero2     = bzero*bzero;    
  bone       = BesselJone(krt);
  bone2      = bone*bone;
  bonebyarg  = BesselOneByArg(krt);
  bonebyarg2 = bonebyarg*bonebyarg;  

  if (fParticle == theProton)
  {
    diffuse = 0.63*fermi;
    // diffuse = 0.6*fermi;
    gamma   = 0.3*fermi;
    delta   = 0.1*fermi*fermi;
    e1      = 0.3*fermi;
    e2      = 0.35*fermi;
  }
  else // as proton, if were not defined 
  {
    diffuse = 0.63*fermi;
    gamma   = 0.3*fermi;
    delta   = 0.1*fermi*fermi;
    e1      = 0.3*fermi;
    e2      = 0.35*fermi;
  }
  G4double lambda = 15.; // 15 ok
  // G4double kgamma    = fWaveVector*gamma;   // wavek*delta;
  G4double kgamma    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;

  // G4cout<<"kgamma = "<<kgamma<<G4endl;

  if(fAddCoulomb)  // add Coulomb correction
  {
    G4double sinHalfTheta  = theta*0.5; // std::sin(0.5*theta);
    G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;

    kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
  }

  G4double kgamma2   = kgamma*kgamma;

  // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
  //   G4cout<<"dk2t = "<<dk2t<<G4endl;
  // G4double dk2t2 = dk2t*dk2t;
  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;

  G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;

  // G4cout<<"pikdt = "<<pikdt<<G4endl;

  damp           = DampFactor(pikdt);
  damp2          = damp*damp;

  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
  G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;

  sigma  = kgamma2;
  // sigma += dk2t2;
  sigma *= bzero2;
  sigma += mode2k2*bone2; 
  sigma += e2dk3t*bzero*bone;

  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
  sigma += kr2*bonebyarg2;  // correction at J1()/()

  sigma *= damp2;          // *rad*rad;

  return sigma;
}

G4double G4NuclNuclDiffuseElastic::GetDiffuseElasticSumXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 253 of file G4NuclNuclDiffuseElastic.cc.

References CalculateAm(), CalculateNuclearRad(), CalculateParticleBeta(), CalculateZommerfeld(), GetDiffElasticSumProb(), G4ParticleDefinition::GetPDGCharge(), python.hepunit::hbarc, and z.

Referenced by GetInvElasticSumXsc().

{
  fParticle      = particle;
  fWaveVector    = momentum/hbarc;
  fAtomicWeight  = A;
  fAtomicNumber  = Z;
  fNuclearRadius = CalculateNuclearRad(A);
  fAddCoulomb    = false;

  G4double z     = particle->GetPDGCharge();

  G4double kRt   = fWaveVector*fNuclearRadius*theta;
  G4double kRtC  = 1.9;

  if( z && (kRt > kRtC) )
  {
    fAddCoulomb = true;
    fBeta       = CalculateParticleBeta( particle, momentum);
    fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
    fAm         = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
  }
  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);

  return sigma;
}

G4double G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A
	)

Definition at line 182 of file G4NuclNuclDiffuseElastic.cc.

References CalculateNuclearRad(), GetDiffElasticProb(), and python.hepunit::hbarc.

Referenced by GetInvElasticXsc().

{
  fParticle      = particle;
  fWaveVector    = momentum/hbarc;
  fAtomicWeight  = A;
  fAddCoulomb    = false;
  fNuclearRadius = CalculateNuclearRad(A);

  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);

  return sigma;
}

G4double G4NuclNuclDiffuseElastic::GetErf ( G4double  x )

inline

Definition at line 626 of file G4NuclNuclDiffuseElastic.hh.

References z.

Referenced by GetErfComp(), and GetErfInt().

{
  G4double t, z, tmp, result;

  z   = std::fabs(x);
  t   = 1.0/(1.0+0.5*z);

  tmp = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
        t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
        t*(-0.82215223+t*0.17087277)))))))));

  if( x >= 0.) result = 1. - tmp;
  else         result = 1. + tmp; 
    
  return result;
}

G4complex G4NuclNuclDiffuseElastic::GetErfcComp ( G4complex  z, G4int  nMax  )

inline

Definition at line 647 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex erfcz = 1. - GetErfComp( z, nMax);
  return erfcz;
}

G4complex G4NuclNuclDiffuseElastic::GetErfcInt ( G4complex  z )

inline

Definition at line 667 of file G4NuclNuclDiffuseElastic.hh.

Referenced by AmplitudeSim(), GammaLess(), and GammaMore().

{
  G4complex erfcz = 1. - GetErfInt( z); // , nMax);
  return erfcz;
}

G4complex G4NuclNuclDiffuseElastic::GetErfComp	(	G4complex	z,
		G4int	nMax
	)

Definition at line 1415 of file G4NuclNuclDiffuseElastic.cc.

References GetErf(), n, and test::x.

{
  G4int n;
  G4double n2, cofn, shny, chny, fn, gn;

  G4double x = z.real();
  G4double y = z.imag();

  G4double outRe = 0., outIm = 0.;

  G4double twox  = 2.*x;
  G4double twoxy = twox*y;
  G4double twox2 = twox*twox;

  G4double cof1 = std::exp(-x*x)/CLHEP::pi;

  G4double cos2xy = std::cos(twoxy);
  G4double sin2xy = std::sin(twoxy);

  G4double twoxcos2xy = twox*cos2xy;
  G4double twoxsin2xy = twox*sin2xy;

  for( n = 1; n <= nMax; n++)
  {
    n2   = n*n;

    cofn = std::exp(-0.5*n2)/(n2+twox2);  // /(n2+0.5*twox2);

    chny = std::cosh(n*y);
    shny = std::sinh(n*y);

    fn  = twox - twoxcos2xy*chny + n*sin2xy*shny;
    gn  =        twoxsin2xy*chny + n*cos2xy*shny;

    fn *= cofn;
    gn *= cofn;

    outRe += fn;
    outIm += gn;
  }
  outRe *= 2*cof1;
  outIm *= 2*cof1;

  if(std::abs(x) < 0.0001)
  {
    outRe += GetErf(x);
    outIm += cof1*y;
  }
  else
  {
    outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
    outIm += cof1*sin2xy/twox;
  }
  return G4complex(outRe, outIm);
}

G4complex G4NuclNuclDiffuseElastic::GetErfcSer ( G4complex  z, G4int  nMax  )

inline

Definition at line 657 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex erfcz = 1. - GetErfSer( z, nMax);
  return erfcz;
}

G4complex G4NuclNuclDiffuseElastic::GetErfInt

(

G4complex

z

)

Definition at line 1476 of file G4NuclNuclDiffuseElastic.cc.

References GetErf(), GetExpCos(), GetExpSin(), G4Integrator< T, F >::Legendre96(), and test::x.

{
  G4double outRe, outIm;

  G4double x = z.real();
  G4double y = z.imag();
  fReZ       = x;

  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;

  outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
  outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );

  outRe *= 2./std::sqrt(CLHEP::pi);
  outIm *= 2./std::sqrt(CLHEP::pi);

  outRe += GetErf(x);

  return G4complex(outRe, outIm);
}

G4complex G4NuclNuclDiffuseElastic::GetErfSer ( G4complex  z, G4int  nMax  )

inline

Definition at line 714 of file G4NuclNuclDiffuseElastic.hh.

References test::a, test::b, n, and z.

{
  G4int n;
  G4double a =1., b = 1., tmp;
  G4complex sum = z, d = z;

  for( n = 1; n <= nMax; n++)
  {
    a *= 2.;
    b *= 2.*n +1.;
    d *= z*z;

    tmp = a/b;

    sum += tmp*d;
  }
  sum *= 2.*std::exp(-z*z)/std::sqrt(CLHEP::pi);

  return sum;
}

G4double G4NuclNuclDiffuseElastic::GetExpCos ( G4double  x )

inline

Definition at line 737 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetErfInt().

{
  G4double result;

  result = std::exp(x*x-fReZ*fReZ);
  result *= std::cos(2.*x*fReZ);
  return result;
}

G4double G4NuclNuclDiffuseElastic::GetExpSin ( G4double  x )

inline

Definition at line 748 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetErfInt().

{
  G4double result;

  result = std::exp(x*x-fReZ*fReZ);
  result *= std::sin(2.*x*fReZ);
  return result;
}

G4double G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc ( G4double  theta )

inline

Definition at line 1018 of file G4NuclNuclDiffuseElastic.hh.

{
  G4double ratio   = GetRatioGen(theta);
  G4double ruthXsc = GetRutherfordXsc(theta);
  G4double xsc     = ratio*ruthXsc;
  return xsc;
}

G4double G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc ( G4double  alpha )

inline

Definition at line 1030 of file G4NuclNuclDiffuseElastic.hh.

Referenced by BuildAngleTable().

{
  G4double theta = std::sqrt(alpha);
  G4double xsc     = GetFresnelDiffuseXsc(theta);
  return xsc;
}

G4double G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS	(	G4ParticleDefinition *	pParticle,
		G4double	pTkin,
		G4ParticleDefinition *	tParticle
	)

Definition at line 1811 of file G4NuclNuclDiffuseElastic.cc.

References CalcMandelstamS(), G4cout, G4endl, and G4ParticleDefinition::GetPDGMass().

Referenced by InitParametersGla().

{
  G4double xsection(0), /*Delta,*/ A0, B0;
  G4double hpXsc(0);
  G4double hnXsc(0);


  G4double targ_mass     = tParticle->GetPDGMass();
  G4double proj_mass     = pParticle->GetPDGMass(); 

  G4double proj_energy   = proj_mass + pTkin; 
  G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));

  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );

  sMand         /= CLHEP::GeV*CLHEP::GeV;  // in GeV for parametrisation
  proj_momentum /= CLHEP::GeV;
  proj_energy   /= CLHEP::GeV;
  proj_mass     /= CLHEP::GeV;
  G4double logS = std::log(sMand);

  // General PDG fit constants


  // fEtaRatio=Re[f(0)]/Im[f(0)]

  if( proj_momentum >= 1.2 )
  {
    fEtaRatio  = 0.13*(logS - 5.8579332)*std::pow(sMand,-0.18);
  }       
  else if( proj_momentum >= 0.6 )
  { 
    fEtaRatio = -75.5*(std::pow(proj_momentum,0.25)-0.95)/
      (std::pow(3*proj_momentum,2.2)+1);     
  }
  else 
  {
    fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
  }
  G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;

  // xsc
  
  if( proj_momentum >= 10. ) // high energy: pp = nn = np
    // if( proj_momentum >= 2.)
  {
    //Delta = 1.;

    //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;

    if( proj_momentum >= 10.)
    {
        B0 = 7.5;
        A0 = 100. - B0*std::log(3.0e7);

        xsection = A0 + B0*std::log(proj_energy) - 11
                  + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
                     0.93827*0.93827,-0.165);        //  mb
    }
  }
  else // low energy pp = nn != np
  {
      if(pParticle == tParticle) // pp or nn      // nn to be pp
      {
        if( proj_momentum < 0.73 )
        {
          hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
        }
        else if( proj_momentum < 1.05  )
        {
          hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
                         (std::log(proj_momentum/0.73));
        }
        else  // if( proj_momentum < 10.  )
        {
          hnXsc = 39.0 + 
              75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
        }
        xsection = hnXsc;
      }
      else  // pn to be np
      {
        if( proj_momentum < 0.8 )
        {
          hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
        }      
        else if( proj_momentum < 1.4 )
        {
          hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
        }
        else    // if( proj_momentum < 10.  )
        {
          hpXsc = 33.3+
              20.8*(std::pow(proj_momentum,2.0)-1.35)/
                 (std::pow(proj_momentum,2.50)+0.95);
        }
        xsection = hpXsc;
      }
  }
  xsection *= CLHEP::millibarn; // parametrised in mb
  G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
  return xsection;
}

G4double G4NuclNuclDiffuseElastic::GetIntegrandFunction

(

G4double

theta

)

Definition at line 657 of file G4NuclNuclDiffuseElastic.cc.

References GetDiffElasticSumProbA().

Referenced by IntegralElasticProb(), SampleThetaCMS(), and TestAngleTable().

{
  G4double result;

  result  = GetDiffElasticSumProbA(alpha);

  // result *= 2*pi*std::sin(theta);

  return result;
}

G4double G4NuclNuclDiffuseElastic::GetInvCoulombElasticXsc	(	const G4ParticleDefinition *	particle,
		G4double	tMand,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 340 of file G4NuclNuclDiffuseElastic.cc.

References CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), GetCoulombElasticXsc(), G4IonTable::GetIon(), G4ParticleTable::GetIonTable(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4He3::He3(), CLHEP::Hep3Vector::mag(), python.hepunit::pi, G4Triton::Triton(), and CLHEP::HepLorentzVector::vect().

{
  G4double m1 = particle->GetPDGMass();
  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));

  G4int iZ = static_cast<G4int>(Z+0.5);
  G4int iA = static_cast<G4int>(A+0.5);
  G4ParticleDefinition * theDef = 0;

  if      (iZ == 1 && iA == 1) theDef = theProton;
  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
  else if (iZ == 2 && iA == 4) theDef = theAlpha;
  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
 
  G4double tmass = theDef->GetPDGMass();

  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
  lv += lv1;

  G4ThreeVector bst = lv.boostVector();
  lv1.boost(-bst);

  G4ThreeVector p1 = lv1.vect();
  G4double ptot    = p1.mag();
  G4double ptot2 = ptot*ptot;
  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;

  if( cost >= 1.0 )      cost = 1.0;  
  else if( cost <= -1.0) cost = -1.0;
  
  G4double thetaCMS = std::acos(cost);

  G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );

  sigma *= pi/ptot2;

  return sigma;
}

G4double G4NuclNuclDiffuseElastic::GetInvElasticSumXsc	(	const G4ParticleDefinition *	particle,
		G4double	tMand,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 288 of file G4NuclNuclDiffuseElastic.cc.

References CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), GetDiffuseElasticSumXsc(), G4IonTable::GetIon(), G4ParticleTable::GetIonTable(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4He3::He3(), CLHEP::Hep3Vector::mag(), python.hepunit::pi, G4Triton::Triton(), and CLHEP::HepLorentzVector::vect().

{
  G4double m1 = particle->GetPDGMass();

  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));

  G4int iZ = static_cast<G4int>(Z+0.5);
  G4int iA = static_cast<G4int>(A+0.5);

  G4ParticleDefinition* theDef = 0;

  if      (iZ == 1 && iA == 1) theDef = theProton;
  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
  else if (iZ == 2 && iA == 4) theDef = theAlpha;
  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
 
  G4double tmass = theDef->GetPDGMass();

  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
  lv += lv1;

  G4ThreeVector bst = lv.boostVector();
  lv1.boost(-bst);

  G4ThreeVector p1 = lv1.vect();
  G4double ptot    = p1.mag();
  G4double ptot2   = ptot*ptot;
  G4double cost    = 1 - 0.5*std::fabs(tMand)/ptot2;

  if( cost >= 1.0 )      cost = 1.0;  
  else if( cost <= -1.0) cost = -1.0;
  
  G4double thetaCMS = std::acos(cost);

  G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );

  sigma *= pi/ptot2;

  return sigma;
}

G4double G4NuclNuclDiffuseElastic::GetInvElasticXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 203 of file G4NuclNuclDiffuseElastic.cc.

References CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), GetDiffuseElasticXsc(), G4IonTable::GetIon(), G4ParticleTable::GetIonTable(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4He3::He3(), CLHEP::Hep3Vector::mag(), python.hepunit::pi, G4Triton::Triton(), and CLHEP::HepLorentzVector::vect().

{
  G4double m1 = particle->GetPDGMass();
  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));

  G4int iZ = static_cast<G4int>(Z+0.5);
  G4int iA = static_cast<G4int>(A+0.5);
  G4ParticleDefinition * theDef = 0;

  if      (iZ == 1 && iA == 1) theDef = theProton;
  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
  else if (iZ == 2 && iA == 4) theDef = theAlpha;
  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
 
  G4double tmass = theDef->GetPDGMass();

  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
  lv += lv1;

  G4ThreeVector bst = lv.boostVector();
  lv1.boost(-bst);

  G4ThreeVector p1 = lv1.vect();
  G4double ptot    = p1.mag();
  G4double ptot2 = ptot*ptot;
  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;

  if( cost >= 1.0 )      cost = 1.0;  
  else if( cost <= -1.0) cost = -1.0;
  
  G4double thetaCMS = std::acos(cost);

  G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);

  sigma *= pi/ptot2;

  return sigma;
}

G4double G4NuclNuclDiffuseElastic::GetLegendrePol	(	G4int	n,
		G4double	x
	)

Definition at line 1389 of file G4NuclNuclDiffuseElastic.cc.

References test::x.

Referenced by AmplitudeGla().

{
  G4double legPol, epsilon = 1.e-6;
  G4double x = std::cos(theta);

  if     ( n  < 0 ) legPol = 0.;
  else if( n == 0 ) legPol = 1.;
  else if( n == 1 ) legPol = x;
  else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
  else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
  else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
  else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
  else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
  else           
  {
    // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;

    legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
  }
  return legPol; 
}

G4double G4NuclNuclDiffuseElastic::GetNuclearRadius ( )

inline

Definition at line 181 of file G4NuclNuclDiffuseElastic.hh.

{return fNuclearRadius;};

G4double G4NuclNuclDiffuseElastic::GetProfileLambda ( )

inline

Definition at line 270 of file G4NuclNuclDiffuseElastic.hh.

{return fProfileLambda;};

G4double G4NuclNuclDiffuseElastic::GetRatioGen

(

G4double

theta

)

Definition at line 1921 of file G4NuclNuclDiffuseElastic.cc.

References GetCint(), GetSint(), and Profile().

{
  G4double sinThetaR  = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
  G4double dTheta     = 0.5*(theta - fRutherfordTheta);
  G4double sindTheta  = std::sin(dTheta);

  G4double prof       = Profile(theta);
  G4double prof2      = prof*prof;
  // G4double profmod    = std::abs(prof);
  G4double order      = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;

  order = std::abs(order); // since sin changes sign!
  // G4cout<<"order = "<<order<<G4endl; 
 
  G4double cosFresnel = GetCint(order);  
  G4double sinFresnel = GetSint(order);  

  G4double out;

  if ( theta <= fRutherfordTheta )
  {
    out  = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2; 
    out += ( cosFresnel + sinFresnel - 1. )*prof;
  }
  else
  {
    out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
  }

  return out;
}

G4double G4NuclNuclDiffuseElastic::GetRatioSim ( G4double  theta )

inline

Definition at line 998 of file G4NuclNuclDiffuseElastic.hh.

{
  G4double sinThetaR  = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
  G4double dTheta     = 0.5*(theta - fRutherfordTheta);
  G4double sindTheta  = std::sin(dTheta);

  G4double order      = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
  // G4cout<<"order = "<<order<<G4endl;  
  G4double cosFresnel = 0.5 - GetCint(order);  
  G4double sinFresnel = 0.5 - GetSint(order);  

  G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );

  return out;
}

G4double G4NuclNuclDiffuseElastic::GetRutherfordXsc ( G4double  theta )

inline

Definition at line 532 of file G4NuclNuclDiffuseElastic.hh.

{
  G4double sinHalfTheta  = std::sin(0.5*theta);
  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;

  G4double ch2           = fRutherfordRatio*fRutherfordRatio;

  G4double xsc           = ch2/(sinHalfTheta2+fAm)/(sinHalfTheta2+fAm);

  return xsc;
}

G4double G4NuclNuclDiffuseElastic::GetScatteringAngle	(	G4int	iMomentum,
		G4int	iAngle,
		G4double	position
	)

Definition at line 1035 of file G4NuclNuclDiffuseElastic.cc.

References G4UniformRand.

Referenced by SampleTableThetaCMS().

{
 G4double x1, x2, y1, y2, randAngle;

  if( iAngle == 0 )
  {
    randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
    // iAngle++;
  }
  else
  {
    if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
    {
      iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
    }
    y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
    y2 = (*(*fAngleTable)(iMomentum))(iAngle);

    x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
    x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);

    if ( x1 == x2 )   randAngle = x2;
    else
    {
      if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
      else
      {
        randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
      }
    }
  }
  return randAngle;
}

G4double G4NuclNuclDiffuseElastic::GetSinHaPit2 ( G4double  t )

inline

Definition at line 192 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetSint().

{return std::sin(CLHEP::halfpi*t*t);};

G4double G4NuclNuclDiffuseElastic::GetSint ( G4double  x )

inline

Definition at line 777 of file G4NuclNuclDiffuseElastic.hh.

References GetSinHaPit2(), and G4Integrator< T, F >::Legendre96().

Referenced by GetRatioGen().

{
  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;

  G4double out = 
    integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetSinHaPit2, 0., x );

  return out;
}

void G4NuclNuclDiffuseElastic::InitDynParameters	(	const G4ParticleDefinition *	theParticle,
		G4double	partMom
	)

Definition at line 1716 of file G4NuclNuclDiffuseElastic.cc.

References test::a, CalculateAm(), CalculateCoulombPhaseZero(), CalculateRutherfordAnglePar(), CalculateZommerfeld(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4InuclParticleNames::lambda, and z.

Referenced by BuildAngleTable().

{
  G4double a = 0.;
  G4double z = theParticle->GetPDGCharge();
  G4double m1 = theParticle->GetPDGMass();

  fWaveVector = partMom/CLHEP::hbarc;

  G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;

  if( z )
  {
    a           = partMom/m1; // beta*gamma for m1
    fBeta       = a/std::sqrt(1+a*a);
    fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
    fRutherfordRatio = fZommerfeld/fWaveVector; 
    fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
  }
  fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
  fProfileDelta  = fCofDelta*fProfileLambda;
  fProfileAlpha  = fCofAlpha*fProfileLambda;

  CalculateCoulombPhaseZero();
  CalculateRutherfordAnglePar();

  return;
}

void G4NuclNuclDiffuseElastic::Initialise

(

)

Definition at line 142 of file G4NuclNuclDiffuseElastic.cc.

References BuildAngleTable(), CalculateNuclearRad(), G4cout, G4endl, G4ParticleDefinition::GetBaryonNumber(), G4Element::GetElementTable(), G4Element::GetNumberOfElements(), and G4HadronicInteraction::verboseLevel.

{

  // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );

  const G4ElementTable* theElementTable = G4Element::GetElementTable();
  size_t jEl, numOfEl = G4Element::GetNumberOfElements();

  // projectile radius

  G4double A1 = G4double( fParticle->GetBaryonNumber() );
  G4double R1 = CalculateNuclearRad(A1);

  for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
  {
    fAtomicNumber = (*theElementTable)[jEl]->GetZ();     // atomic number
    fAtomicWeight = (*theElementTable)[jEl]->GetN();     // number of nucleons

    fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
    fNuclearRadius += R1;

    if(verboseLevel > 0) 
    {   
      G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element: "
        <<(*theElementTable)[jEl]->GetName()<<G4endl;
    }
    fElementNumberVector.push_back(fAtomicNumber);
    fElementNameVector.push_back((*theElementTable)[jEl]->GetName());

    BuildAngleTable();
    fAngleBank.push_back(fAngleTable);
  }  
}

void G4NuclNuclDiffuseElastic::InitialiseOnFly	(	G4double	Z,
		G4double	A
	)

Definition at line 929 of file G4NuclNuclDiffuseElastic.cc.

References BuildAngleTable(), CalculateNuclearRad(), G4cout, G4endl, G4ParticleDefinition::GetBaryonNumber(), and G4HadronicInteraction::verboseLevel.

Referenced by SampleTableThetaCMS().

{
  fAtomicNumber  = Z;     // atomic number
  fAtomicWeight  = A;     // number of nucleons

  G4double A1 = G4double( fParticle->GetBaryonNumber() );
  G4double R1 = CalculateNuclearRad(A1);

  fNuclearRadius = CalculateNuclearRad(fAtomicWeight) + R1;
  
  if( verboseLevel > 0 )    
  {
    G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
      <<Z<<"; and A = "<<A<<G4endl;
  }
  fElementNumberVector.push_back(fAtomicNumber);

  BuildAngleTable();

  fAngleBank.push_back(fAngleTable);

  return;
}

void G4NuclNuclDiffuseElastic::InitParameters	(	const G4ParticleDefinition *	theParticle,
		G4double	partMom,
		G4double	Z,
		G4double	A
	)

Definition at line 1671 of file G4NuclNuclDiffuseElastic.cc.

References test::a, CalculateAm(), CalculateCoulombPhaseZero(), CalculateNuclearRad(), CalculateRutherfordAnglePar(), CalculateZommerfeld(), G4cout, G4endl, G4ParticleDefinition::GetBaryonNumber(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4InuclParticleNames::lambda, and z.

{
  fAtomicNumber  = Z;     // atomic number
  fAtomicWeight  = A;     // number of nucleons

  fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
  G4double A1     = G4double( theParticle->GetBaryonNumber() );   
  fNuclearRadius1 = CalculateNuclearRad(A1);
  // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
  fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;

  G4double a = 0.;
  G4double z = theParticle->GetPDGCharge();
  G4double m1 = theParticle->GetPDGMass();

  fWaveVector = partMom/CLHEP::hbarc;

  G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
  G4cout<<"kR = "<<lambda<<G4endl;

  if( z )
  {
    a           = partMom/m1; // beta*gamma for m1
    fBeta       = a/std::sqrt(1+a*a);
    fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
    fRutherfordRatio = fZommerfeld/fWaveVector; 
    fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
  }
  G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
  fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
  G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
  fProfileDelta  = fCofDelta*fProfileLambda;
  fProfileAlpha  = fCofAlpha*fProfileLambda;

  CalculateCoulombPhaseZero();
  CalculateRutherfordAnglePar();

  return;
}

void G4NuclNuclDiffuseElastic::InitParametersGla	(	const G4DynamicParticle *	aParticle,
		G4double	partMom,
		G4double	Z,
		G4double	A
	)

Definition at line 1751 of file G4NuclNuclDiffuseElastic.cc.

References test::a, CalculateAm(), CalculateCoulombPhaseZero(), CalculateNuclearRad(), CalculateZommerfeld(), G4cout, G4endl, G4ParticleDefinition::GetBaryonNumber(), G4DynamicParticle::GetDefinition(), GetHadronNucleonXscNS(), G4DynamicParticle::GetKineticEnergy(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and z.

{
  fAtomicNumber  = Z;     // target atomic number
  fAtomicWeight  = A;     // target number of nucleons

  fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
  G4double A1     = G4double( aParticle->GetDefinition()->GetBaryonNumber() );   
  fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
  fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
 

  G4double a  = 0., kR12;
  G4double z  = aParticle->GetDefinition()->GetPDGCharge();
  G4double m1 = aParticle->GetDefinition()->GetPDGMass();

  fWaveVector = partMom/CLHEP::hbarc;

  G4double pN = A1 - z;
  if( pN < 0. ) pN = 0.;

  G4double tN = A - Z;
  if( tN < 0. ) tN = 0.;

  G4double pTkin = aParticle->GetKineticEnergy();  
  pTkin /= A1;


  fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
              (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);

  G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl;
  G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiusSquare/CLHEP::millibarn<<" mb"<<G4endl;
  kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
  G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
  fMaxL = (G4int(kR12)+1)*4;
  G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;

  if( z )
  {
      a           = partMom/m1; // beta*gamma for m1
      fBeta       = a/std::sqrt(1+a*a);
      fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
      fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
  }

  CalculateCoulombPhaseZero();
 

  return;
}

G4double G4NuclNuclDiffuseElastic::IntegralElasticProb	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A
	)

Definition at line 673 of file G4NuclNuclDiffuseElastic.cc.

References CalculateNuclearRad(), GetIntegrandFunction(), python.hepunit::hbarc, and G4Integrator< T, F >::Legendre96().

{
  G4double result;
  fParticle      = particle;
  fWaveVector    = momentum/hbarc;
  fAtomicWeight  = A;

  fNuclearRadius = CalculateNuclearRad(A);


  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;

  // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 
  result = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 

  return result;
}

G4complex G4NuclNuclDiffuseElastic::PhaseFar ( G4double  theta )

inline

Definition at line 940 of file G4NuclNuclDiffuseElastic.hh.

References z.

{
  G4double twosigma = 2.*fCoulombPhase0; 
  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
  twosigma += fProfileLambda*theta - 0.25*CLHEP::pi;

  twosigma *= fCofPhase;

  G4complex z = G4complex(0., twosigma);

  return std::exp(z);
}

G4complex G4NuclNuclDiffuseElastic::PhaseNear ( G4double  theta )

inline

Definition at line 922 of file G4NuclNuclDiffuseElastic.hh.

References z.

Referenced by AmplitudeNear().

{
  G4double twosigma = 2.*fCoulombPhase0; 
  twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
  twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
  twosigma -= fProfileLambda*theta - 0.25*CLHEP::pi;

  twosigma *= fCofPhase;

  G4complex z = G4complex(0., twosigma);

  return std::exp(z);
}

G4double G4NuclNuclDiffuseElastic::Profile ( G4double  theta )

inline

Definition at line 903 of file G4NuclNuclDiffuseElastic.hh.

Referenced by GetRatioGen().

{
  G4double dTheta = fRutherfordTheta - theta;
  G4double result = 0., argument = 0.;

  if(std::abs(dTheta) < 0.001) result = 1.;
  else
  {
    argument = fProfileDelta*dTheta;
    result   = CLHEP::pi*argument;
    result  /= std::sinh(result);
  }
  return result;
}

G4double G4NuclNuclDiffuseElastic::ProfileFar ( G4double  theta )

inline

Definition at line 887 of file G4NuclNuclDiffuseElastic.hh.

{
  G4double dTheta   = fRutherfordTheta + theta;
  G4double argument = fProfileDelta*dTheta;

  G4double result   = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
  result           /= std::sinh(CLHEP::pi*argument);
  result           /= dTheta;

  return result;
}

G4double G4NuclNuclDiffuseElastic::ProfileNear ( G4double  theta )

inline

Definition at line 866 of file G4NuclNuclDiffuseElastic.hh.

Referenced by AmplitudeNear(), and AmplitudeSim().

{
  G4double dTheta = fRutherfordTheta - theta;
  G4double result = 0., argument = 0.;

  if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
  else
  {
    argument = fProfileDelta*dTheta;
    result   = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
    result  /= std::sinh(CLHEP::pi*argument);
    result  -= 1.;
    result  /= dTheta;
  }
  return result;
}

G4double G4NuclNuclDiffuseElastic::SampleInvariantT ( const G4ParticleDefinition *  p, G4double  plab, G4int  Z, G4int  A  )

virtual

Reimplemented from G4HadronElastic.

Definition at line 766 of file G4NuclNuclDiffuseElastic.cc.

References CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetPDGMass(), CLHEP::Hep3Vector::mag(), SampleTableT(), and CLHEP::HepLorentzVector::vect().

{
  fParticle = aParticle;
  G4double m1 = fParticle->GetPDGMass();
  G4double totElab = std::sqrt(m1*m1+p*p);
  G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
  G4LorentzVector lv1(p,0.0,0.0,totElab);
  G4LorentzVector  lv(0.0,0.0,0.0,mass2);   
  lv += lv1;

  G4ThreeVector bst = lv.boostVector();
  lv1.boost(-bst);

  G4ThreeVector p1 = lv1.vect();
  G4double momentumCMS = p1.mag();

  G4double t = SampleTableT( aParticle,  momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
  return t;
}

G4double G4NuclNuclDiffuseElastic::SampleT	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	A
	)

Definition at line 698 of file G4NuclNuclDiffuseElastic.cc.

References SampleThetaCMS().

Referenced by SampleThetaLab().

{
  G4double theta = SampleThetaCMS( aParticle,  p, A); // sample theta in cms
  G4double t     = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
  return t;
}

G4double G4NuclNuclDiffuseElastic::SampleTableT	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	Z,
		G4double	A
	)

Definition at line 791 of file G4NuclNuclDiffuseElastic.cc.

References SampleTableThetaCMS().

Referenced by SampleInvariantT().

{
  G4double alpha = SampleTableThetaCMS( aParticle,  p, Z, A); // sample theta2 in cms
  // G4double t     = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) );             // -t !!!
  G4double t     = p*p*alpha;             // -t !!!
  return t;
}

G4double G4NuclNuclDiffuseElastic::SampleTableThetaCMS	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	Z,
		G4double	A
	)

Definition at line 806 of file G4NuclNuclDiffuseElastic.cc.

References G4UniformRand, G4PhysicsVector::GetLowEdgeEnergy(), G4ParticleDefinition::GetPDGMass(), GetScatteringAngle(), InitialiseOnFly(), and position.

Referenced by SampleTableT().

{
  size_t iElement;
  G4int iMomentum, iAngle;  
  G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;  
  G4double m1 = particle->GetPDGMass();

  for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
  {
    if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
  }
  if ( iElement == fElementNumberVector.size() ) 
  {
    InitialiseOnFly(Z,A); // table preparation, if needed

    // iElement--;

    // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
    // << " is not found, return zero angle" << G4endl;
    // return 0.; // no table for this element
  }
  // G4cout<<"iElement = "<<iElement<<G4endl;

  fAngleTable = fAngleBank[iElement];

  G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;

  for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
  {
    // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl;
    
    if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
  }
  // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl;


  if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;   // kinE is more then theMaxEnergy
  if ( iMomentum < 0 )           iMomentum = 0; // against negative index, kinE < theMinEnergy


  if (iMomentum == fEnergyBin -1 || iMomentum == 0 )   // the table edges
  {
    position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();

    // G4cout<<"position = "<<position<<G4endl;

    for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
    {
      if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
    }
    if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;

    // G4cout<<"iAngle = "<<iAngle<<G4endl;

    randAngle = GetScatteringAngle(iMomentum, iAngle, position);

    // G4cout<<"randAngle = "<<randAngle<<G4endl;
  }
  else  // kinE inside between energy table edges
  {
    // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
    position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();

    // G4cout<<"position = "<<position<<G4endl;

    for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
    {
      // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
      if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
    }
    if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;

    // G4cout<<"iAngle = "<<iAngle<<G4endl;

    theta2  = GetScatteringAngle(iMomentum, iAngle, position);

    // G4cout<<"theta2 = "<<theta2<<G4endl;

    E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);

    // G4cout<<"E2 = "<<E2<<G4endl;
    
    iMomentum--;
    
    // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();

    // G4cout<<"position = "<<position<<G4endl;

    for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
    {
      // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
      if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
    }
    if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
    
    theta1  = GetScatteringAngle(iMomentum, iAngle, position);

    // G4cout<<"theta1 = "<<theta1<<G4endl;

    E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);

    // G4cout<<"E1 = "<<E1<<G4endl;

    W  = 1.0/(E2 - E1);
    W1 = (E2 - kinE)*W;
    W2 = (kinE - E1)*W;

    randAngle = W1*theta1 + W2*theta2;
    
    // randAngle = theta2;
  }
  // G4double angle = randAngle;
  // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
  // G4cout<<"randAngle = "<<randAngle/degree<<G4endl;

  return randAngle;
}

G4double G4NuclNuclDiffuseElastic::SampleThetaCMS	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	A
	)

Definition at line 712 of file G4NuclNuclDiffuseElastic.cc.

References CalculateNuclearRad(), G4UniformRand, GetIntegrandFunction(), python.hepunit::hbarc, G4Integrator< T, F >::Legendre10(), G4Integrator< T, F >::Legendre96(), python.hepunit::pi, and G4INCL::DeJongSpin::shoot().

Referenced by SampleT().

{
  G4int i, iMax = 100;  
  G4double norm, result, theta1, theta2, thetaMax, sum = 0.;

  fParticle      = particle;
  fWaveVector    = momentum/hbarc;
  fAtomicWeight  = A;

  fNuclearRadius = CalculateNuclearRad(A);
  
  thetaMax = 10.174/fWaveVector/fNuclearRadius;

  if (thetaMax > pi) thetaMax = pi;

  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;

  // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 
  norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );

  norm *= G4UniformRand();

  for(i = 1; i <= iMax; i++)
  {
    theta1 = (i-1)*thetaMax/iMax; 
    theta2 = i*thetaMax/iMax;
    sum   += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);

    if ( sum >= norm ) 
    {
      result = 0.5*(theta1 + theta2);
      break;
    }
  }
  if (i > iMax ) result = 0.5*(theta1 + theta2);

  G4double sigma = pi*thetaMax/iMax;

  result += G4RandGauss::shoot(0.,sigma);

  if(result < 0.) result = 0.;
  if(result > thetaMax) result = thetaMax;

  return result;
}

G4double G4NuclNuclDiffuseElastic::SampleThetaLab	(	const G4HadProjectile *	aParticle,
		G4double	tmass,
		G4double	A
	)

Definition at line 1078 of file G4NuclNuclDiffuseElastic.cc.

References CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), G4cout, G4endl, G4UniformRand, G4HadProjectile::Get4Momentum(), G4HadProjectile::GetDefinition(), G4ParticleDefinition::GetPDGMass(), G4HadProjectile::GetTotalMomentum(), python.hepunit::GeV, CLHEP::Hep3Vector::mag(), SampleT(), CLHEP::Hep3Vector::theta(), python.hepunit::twopi, CLHEP::HepLorentzVector::vect(), G4HadronicInteraction::verboseLevel, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

{
  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
  G4double m1 = theParticle->GetPDGMass();
  G4double plab = aParticle->GetTotalMomentum();
  G4LorentzVector lv1 = aParticle->Get4Momentum();
  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
  lv += lv1;

  G4ThreeVector bst = lv.boostVector();
  lv1.boost(-bst);

  G4ThreeVector p1 = lv1.vect();
  G4double ptot    = p1.mag();
  G4double tmax    = 4.0*ptot*ptot;
  G4double t       = 0.0;


  //
  // Sample t
  //
  
  t = SampleT( theParticle, ptot, A);

  // NaN finder
  if(!(t < 0.0 || t >= 0.0)) 
  {
    if (verboseLevel > 0) 
    {
      G4cout << "G4NuclNuclDiffuseElastic:WARNING: A = " << A 
         << " mom(GeV)= " << plab/GeV 
             << " S-wave will be sampled" 
         << G4endl; 
    }
    t = G4UniformRand()*tmax; 
  }
  if(verboseLevel>1)
  {
    G4cout <<" t= " << t << " tmax= " << tmax 
       << " ptot= " << ptot << G4endl;
  }
  // Sampling of angles in CM system

  G4double phi  = G4UniformRand()*twopi;
  G4double cost = 1. - 2.0*t/tmax;
  G4double sint;

  if( cost >= 1.0 ) 
  {
    cost = 1.0;
    sint = 0.0;
  }
  else if( cost <= -1.0) 
  {
    cost = -1.0;
    sint =  0.0;
  }
  else  
  {
    sint = std::sqrt((1.0-cost)*(1.0+cost));
  }    
  if (verboseLevel>1) 
  {
    G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
  }
  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
  v1 *= ptot;
  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));

  nlv1.boost(bst); 

  G4ThreeVector np1 = nlv1.vect();

    // G4double theta = std::acos( np1.z()/np1.mag() );  // degree;

  G4double theta = np1.theta();

  return theta;
}

void G4NuclNuclDiffuseElastic::SetCofAlpha ( G4double  pa )

inline

Definition at line 277 of file G4NuclNuclDiffuseElastic.hh.

{fCofAlpha = pa;};

void G4NuclNuclDiffuseElastic::SetCofAlphaCoulomb ( G4double  pa )

inline

Definition at line 279 of file G4NuclNuclDiffuseElastic.hh.

{fCofAlphaCoulomb = pa;};

void G4NuclNuclDiffuseElastic::SetCofAlphaMax ( G4double  pa )

inline

Definition at line 278 of file G4NuclNuclDiffuseElastic.hh.

{fCofAlphaMax = pa;};

void G4NuclNuclDiffuseElastic::SetCofDelta ( G4double  pa )

inline

Definition at line 281 of file G4NuclNuclDiffuseElastic.hh.

{fCofDelta = pa;};

void G4NuclNuclDiffuseElastic::SetCofFar ( G4double  pa )

inline

Definition at line 283 of file G4NuclNuclDiffuseElastic.hh.

{fCofFar = pa;};

void G4NuclNuclDiffuseElastic::SetCofLambda ( G4double  pa )

inline

Definition at line 275 of file G4NuclNuclDiffuseElastic.hh.

{fCofLambda = pa;};

void G4NuclNuclDiffuseElastic::SetCofPhase ( G4double  pa )

inline

Definition at line 282 of file G4NuclNuclDiffuseElastic.hh.

{fCofPhase = pa;};

void G4NuclNuclDiffuseElastic::SetEtaRatio ( G4double  pa )

inline

Definition at line 284 of file G4NuclNuclDiffuseElastic.hh.

{fEtaRatio = pa;};

void G4NuclNuclDiffuseElastic::SetHEModelLowLimit ( G4double  value )

inline

Definition at line 372 of file G4NuclNuclDiffuseElastic.hh.

{
  lowEnergyLimitHE = value;
}

void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit ( G4double  value )

inline

Definition at line 382 of file G4NuclNuclDiffuseElastic.hh.

{
  lowestEnergyLimit = value;
}

void G4NuclNuclDiffuseElastic::SetMaxL ( G4int  l )

inline

Definition at line 285 of file G4NuclNuclDiffuseElastic.hh.

{fMaxL = l;};

void G4NuclNuclDiffuseElastic::SetNuclearRadiusCof ( G4double  r )

inline

Definition at line 286 of file G4NuclNuclDiffuseElastic.hh.

{fNuclearRadiusCof = r;};

void G4NuclNuclDiffuseElastic::SetPlabLowLimit ( G4double  value )

inline

Definition at line 367 of file G4NuclNuclDiffuseElastic.hh.

{
  plabLowLimit = value;
}

void G4NuclNuclDiffuseElastic::SetProfileAlpha ( G4double  pa )

inline

Definition at line 274 of file G4NuclNuclDiffuseElastic.hh.

{fProfileAlpha = pa;};

void G4NuclNuclDiffuseElastic::SetProfileDelta ( G4double  pd )

inline

Definition at line 273 of file G4NuclNuclDiffuseElastic.hh.

{fProfileDelta = pd;};

void G4NuclNuclDiffuseElastic::SetProfileLambda ( G4double  pl )

inline

Definition at line 272 of file G4NuclNuclDiffuseElastic.hh.

References readPY::pl.

{fProfileLambda = pl;};

void G4NuclNuclDiffuseElastic::SetQModelLowLimit ( G4double  value )

inline

Definition at line 377 of file G4NuclNuclDiffuseElastic.hh.

{
  lowEnergyLimitQ = value;
}

void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit ( G4double  value )

inline

Definition at line 362 of file G4NuclNuclDiffuseElastic.hh.

{
  lowEnergyRecoilLimit = value;
}

void G4NuclNuclDiffuseElastic::TestAngleTable	(	const G4ParticleDefinition *	theParticle,
		G4double	partMom,
		G4double	Z,
		G4double	A
	)

Definition at line 1286 of file G4NuclNuclDiffuseElastic.cc.

References test::a, G4Integrator< T, F >::AdaptiveGauss(), CalculateAm(), CalculateNuclearRad(), CalculateZommerfeld(), python.hepunit::degree, G4cout, G4endl, GetIntegrandFunction(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), python.hepunit::hbarc, G4PhysicsTable::insertAt(), G4Integrator< T, F >::Legendre10(), G4Integrator< T, F >::Legendre96(), G4PhysicsFreeVector::PutValue(), and z.

{
  fAtomicNumber  = Z;     // atomic number
  fAtomicWeight  = A;     // number of nucleons
  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
  
     
  
  G4cout<<"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
      <<Z<<"; and A = "<<A<<G4endl;
 
  fElementNumberVector.push_back(fAtomicNumber);

 


  G4int i=0, j;
  G4double a = 0., z = theParticle->GetPDGCharge(),  m1 = fParticle->GetPDGMass();
  G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
  G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
  G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
  G4double epsilon = 0.001;

  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
  
  fAngleTable = new G4PhysicsTable(fEnergyBin);

  fWaveVector = partMom/hbarc;

  G4double kR     = fWaveVector*fNuclearRadius;
  G4double kR2    = kR*kR;
  G4double kRmax  = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
  G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1

  alphaMax = kRmax*kRmax/kR2;

  if (alphaMax > 4.) alphaMax = 4.;  // vmg05-02-09: was pi2 

  alphaCoulomb = kRcoul*kRcoul/kR2;

  if( z )
  {
      a           = partMom/m1; // beta*gamma for m1
      fBeta       = a/std::sqrt(1+a*a);
      fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
      fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
  }
  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);

  // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
    
  
  fAddCoulomb = false;

  for(j = 1; j < fAngleBin; j++)
  {
      // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
      // alpha2 = angleBins->GetLowEdgeEnergy(j);

    alpha1 = alphaMax*(j-1)/fAngleBin;
    alpha2 = alphaMax*( j )/fAngleBin;

    if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;

    deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
    deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
    deltaAG  = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 
                                       alpha1, alpha2,epsilon);

      // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
      //     <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;

    sumL10 += deltaL10;
    sumL96 += deltaL96;
    sumAG  += deltaAG;

    G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
            <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
      
    angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
  }
  fAngleTable->insertAt(i,angleVector);
  fAngleBank.push_back(fAngleTable);

  /*
  // Integral over all angle range - Bad accuracy !!!

  sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
  sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
  sumAG  = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 
                                       0., alpha2,epsilon);
  G4cout<<G4endl;
  G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
            <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
  */
  return;
}

G4complex G4NuclNuclDiffuseElastic::TestErfcComp ( G4complex  z, G4int  nMax  )

inline

Definition at line 678 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex miz = G4complex( z.imag(), -z.real() ); 
  G4complex erfcz = 1. - GetErfComp( miz, nMax);
  G4complex w = std::exp(-z*z)*erfcz;
  return w;
}

G4complex G4NuclNuclDiffuseElastic::TestErfcInt ( G4complex  z )

inline

Definition at line 702 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex miz = G4complex( z.imag(), -z.real() ); 
  G4complex erfcz = 1. - GetErfInt( miz); // , nMax);
  G4complex w = std::exp(-z*z)*erfcz;
  return w;
}

G4complex G4NuclNuclDiffuseElastic::TestErfcSer ( G4complex  z, G4int  nMax  )

inline

Definition at line 690 of file G4NuclNuclDiffuseElastic.hh.

{
  G4complex miz = G4complex( z.imag(), -z.real() ); 
  G4complex erfcz = 1. - GetErfSer( miz, nMax);
  G4complex w = std::exp(-z*z)*erfcz;
  return w;
}

G4double G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab	(	const G4DynamicParticle *	aParticle,
		G4double	tmass,
		G4double	thetaCMS
	)

Definition at line 1166 of file G4NuclNuclDiffuseElastic.cc.

References CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::Get4Momentum(), G4DynamicParticle::GetDefinition(), G4ParticleDefinition::GetPDGMass(), CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::theta(), python.hepunit::twopi, CLHEP::HepLorentzVector::vect(), G4HadronicInteraction::verboseLevel, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

{
  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
  G4double m1 = theParticle->GetPDGMass();
  // G4double plab = aParticle->GetTotalMomentum();
  G4LorentzVector lv1 = aParticle->Get4Momentum();
  G4LorentzVector lv(0.0,0.0,0.0,tmass);   

  lv += lv1;

  G4ThreeVector bst = lv.boostVector();

  lv1.boost(-bst);

  G4ThreeVector p1 = lv1.vect();
  G4double ptot    = p1.mag();

  G4double phi  = G4UniformRand()*twopi;
  G4double cost = std::cos(thetaCMS);
  G4double sint;

  if( cost >= 1.0 ) 
  {
    cost = 1.0;
    sint = 0.0;
  }
  else if( cost <= -1.0) 
  {
    cost = -1.0;
    sint =  0.0;
  }
  else  
  {
    sint = std::sqrt((1.0-cost)*(1.0+cost));
  }    
  if (verboseLevel>1) 
  {
    G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
  }
  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
  v1 *= ptot;
  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));

  nlv1.boost(bst); 

  G4ThreeVector np1 = nlv1.vect();


  G4double thetaLab = np1.theta();

  return thetaLab;
}

G4double G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS	(	const G4DynamicParticle *	aParticle,
		G4double	tmass,
		G4double	thetaLab
	)

Definition at line 1227 of file G4NuclNuclDiffuseElastic.cc.

References CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::Get4Momentum(), G4DynamicParticle::GetDefinition(), G4ParticleDefinition::GetPDGMass(), G4DynamicParticle::GetTotalMomentum(), CLHEP::Hep3Vector::theta(), python.hepunit::twopi, G4HadronicInteraction::verboseLevel, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

{
  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
  G4double m1 = theParticle->GetPDGMass();
  G4double plab = aParticle->GetTotalMomentum();
  G4LorentzVector lv1 = aParticle->Get4Momentum();
  G4LorentzVector lv(0.0,0.0,0.0,tmass);   

  lv += lv1;

  G4ThreeVector bst = lv.boostVector();

  // lv1.boost(-bst);

  // G4ThreeVector p1 = lv1.vect();
  // G4double ptot    = p1.mag();

  G4double phi  = G4UniformRand()*twopi;
  G4double cost = std::cos(thetaLab);
  G4double sint;

  if( cost >= 1.0 ) 
  {
    cost = 1.0;
    sint = 0.0;
  }
  else if( cost <= -1.0) 
  {
    cost = -1.0;
    sint =  0.0;
  }
  else  
  {
    sint = std::sqrt((1.0-cost)*(1.0+cost));
  }    
  if (verboseLevel>1) 
  {
    G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
  }
  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
  v1 *= plab;
  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));

  nlv1.boost(-bst); 

  G4ThreeVector np1 = nlv1.vect();


  G4double thetaCMS = np1.theta();

  return thetaCMS;
}

---

The documentation for this class was generated from the following files:

• G4NuclNuclDiffuseElastic.hh
• G4NuclNuclDiffuseElastic.cc