G4DiffuseElastic Class Reference

#include <G4DiffuseElastic.hh>

Inheritance diagram for G4DiffuseElastic:

Public Member Functions
	G4DiffuseElastic ()
virtual	~G4DiffuseElastic ()
virtual G4bool	IsApplicable (const G4HadProjectile &, G4Nucleus &)
void	Initialise ()
void	InitialiseOnFly (G4double Z, G4double A)
void	BuildAngleTable ()
virtual G4double	SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4double	NeutronTuniform (G4int Z)
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	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 ()
Public Member Functions inherited from G4HadronElastic
	G4HadronElastic (const G4String &name="hElasticLHEP")
	~G4HadronElastic () override
G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) override
G4double	SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override
G4double	GetSlopeCof (const G4int pdg)
void	SetLowestEnergyLimit (G4double value)
G4double	LowestEnergyLimit () const
G4double	ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
void	ModelDescription (std::ostream &) const override
Public Member Functions inherited from G4HadronicInteraction
	G4HadronicInteraction (const G4String &modelName="HadronicModel")
virtual	~G4HadronicInteraction ()
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)
G4int	GetVerboseLevel () const
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
virtual const std::pair< G4double, G4double >	GetFatalEnergyCheckLevels () const
virtual std::pair< G4double, G4double >	GetEnergyMomentumCheckLevels () const
void	SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
virtual void	BuildPhysicsTable (const G4ParticleDefinition &)
virtual void	InitialiseModel ()
	G4HadronicInteraction (const G4HadronicInteraction &right)=delete
const G4HadronicInteraction &	operator= (const G4HadronicInteraction &right)=delete
G4bool	operator== (const G4HadronicInteraction &right) const =delete
G4bool	operator!= (const G4HadronicInteraction &right) const =delete

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

Detailed Description

Definition at line 58 of file G4DiffuseElastic.hh.

Constructor & Destructor Documentation

G4DiffuseElastic()

G4DiffuseElastic::G4DiffuseElastic

(

)

Definition at line 75 of file G4DiffuseElastic.cc.

  : G4HadronElastic("DiffuseElastic"), fParticle(0)
{
  SetMinEnergy( 0.01*MeV );
  SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
 
  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 = 300;  // Increased from the original 200 to have no wider log-energy-bins up to 10 PeV 
  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;
}

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

~G4DiffuseElastic()

G4DiffuseElastic::~G4DiffuseElastic ( )

virtual

Definition at line 117 of file G4DiffuseElastic.cc.

{
  if ( fEnergyVector ) 
  {
    delete fEnergyVector;
    fEnergyVector = 0;
  }
  for ( std::vector<G4PhysicsTable*>::iterator it = fAngleBank.begin();
        it != fAngleBank.end(); ++it ) 
  {
    if ( (*it) ) (*it)->clearAndDestroy();
 
    delete *it;
    *it = 0;
  }
  fAngleTable = 0;
}

Member Function Documentation

BesselJone()

G4double G4DiffuseElastic::BesselJone ( G4double z )

inline

Definition at line 329 of file G4DiffuseElastic.hh.

{
  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;
}

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

BesselJzero()

G4double G4DiffuseElastic::BesselJzero ( G4double z )

inline

Definition at line 277 of file G4DiffuseElastic.hh.

{
  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;
}

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

BesselOneByArg()

G4double G4DiffuseElastic::BesselOneByArg ( G4double z )

inline

Definition at line 404 of file G4DiffuseElastic.hh.

{
  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;
}

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

BuildAngleTable()

void G4DiffuseElastic::BuildAngleTable

(

)

Definition at line 1001 of file G4DiffuseElastic.cc.

{
  G4int i, j;
  G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
  G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
 
  G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
  
  fAngleTable = new G4PhysicsTable( fEnergyBin );
 
  for( i = 0; i < fEnergyBin; i++)
  {
    kinE        = fEnergyVector->GetLowEdgeEnergy(i);
    partMom     = std::sqrt( kinE*(kinE + 2*m1) );
 
    fWaveVector = partMom/hbarc;
 
    G4double kR     = fWaveVector*fNuclearRadius;
    G4double kR2    = kR*kR;
    G4double kRmax  = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
    G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1
    // G4double kRlim  = 1.2;
    // G4double kRlim2 = kRlim*kRlim/kR2;
 
    alphaMax = kRmax*kRmax/kR2;
 
 
    // if (alphaMax > 4.) alphaMax = 4.;  // vmg05-02-09: was pi2 
    // if ( alphaMax > 4. || alphaMax < 1. ) alphaMax = 15.;  // vmg27.11.14 
 
    // if ( alphaMax > 4. || alphaMax < 1. ) alphaMax = CLHEP::pi*CLHEP::pi;  // vmg06.01.15 
 
    // G4cout<<"alphaMax = "<<alphaMax<<", ";
 
    if ( alphaMax >= CLHEP::pi*CLHEP::pi ) alphaMax = CLHEP::pi*CLHEP::pi;   // vmg21.10.15 
 
    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 );
 
    G4double delth = alphaMax/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 = delth*(j-1);
      // if(alpha1 < kRlim2) alpha1 = kRlim2;
      alpha2 = alpha1 + delth;
 
      // if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
      if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false;
 
      delta = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
      // delta = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
 
      sum += delta;
      
      angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
      //      G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2 << " delta "<< delta <<"; sum = "<<sum<<G4endl;
    }
    fAngleTable->insertAt(i, angleVector);
 
    // delete[] angleVector; 
    // delete[] angleBins; 
  }
  return;
}

Referenced by Initialise(), and InitialiseOnFly().

CalculateAm()

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

inline

Definition at line 450 of file G4DiffuseElastic.hh.

{
  G4double k   = momentum/CLHEP::hbarc;
  G4double ch  = 1.13 + 3.76*n*n;
  G4double zn  = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
  G4double zn2 = zn*zn;
  fAm          = ch/zn2;
 
  return fAm;
}

Referenced by BuildAngleTable(), GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), GetDiffuseElasticSumXsc(), and TestAngleTable().

CalculateNuclearRad()

G4double G4DiffuseElastic::CalculateNuclearRad ( G4double A )

inline

Definition at line 465 of file G4DiffuseElastic.hh.

{
  G4double R, r0, a11, a12, a13, a2, a3;
 
  a11 = 1.26;  // 1.08, 1.16
  a12 = 1.;  // 1.08, 1.16
  a13 = 1.12;  // 1.08, 1.16
  a2 = 1.1;
  a3 = 1.;
 
  // Special rms radii for light nucleii
 
  if (A < 50.)
  {
    if     (std::abs(A-1.) < 0.5)                         return 0.89*CLHEP::fermi; // p
    else if(std::abs(A-2.) < 0.5)                         return 2.13*CLHEP::fermi; // d
    else if(  // std::abs(Z-1.) < 0.5 && 
std::abs(A-3.) < 0.5) return 1.80*CLHEP::fermi; // t
 
    // else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96CLHEP::fermi; // He3
    else if( // std::abs(Z-2.) < 0.5 && 
std::abs(A-4.) < 0.5) return 1.68*CLHEP::fermi; // He4
 
    else if(  // std::abs(Z-3.) < 0.5
        std::abs(A-7.) < 0.5   )                         return 2.40*CLHEP::fermi; // Li7
    else if(  // std::abs(Z-4.) < 0.5 
std::abs(A-9.) < 0.5)                         return 2.51*CLHEP::fermi; // Be9
 
    else if( 10.  < A && A <= 16. ) r0  = a11*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;   // 1.08CLHEP::fermi;
    else if( 15.  < A && A <= 20. ) r0  = a12*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
    else if( 20.  < A && A <= 30. ) r0  = a13*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
    else                            r0  = a2*CLHEP::fermi;
 
    R = r0*G4Pow::GetInstance()->A13(A);
  }
  else
  {
    r0 = a3*CLHEP::fermi;
 
    R  = r0*G4Pow::GetInstance()->powA(A, 0.27);
  }
  fNuclearRadius = R;
  return R;
  /*
  G4double r0;
  if( A < 50. )
  {
    if( A > 10. ) r0  = 1.16*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;   // 1.08CLHEP::fermi;
    else          r0  = 1.1*CLHEP::fermi;
    fNuclearRadius = r0*G4Pow::GetInstance()->A13(A);
  }
  else
  {
    r0 = 1.7*CLHEP::fermi;   // 1.7*CLHEP::fermi;
    fNuclearRadius = r0*G4Pow::GetInstance()->powA(A, 0.27); // 0.27);
  }
  return fNuclearRadius;
  */
}

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

CalculateParticleBeta()

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

inline

Definition at line 425 of file G4DiffuseElastic.hh.

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

Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), and GetDiffuseElasticSumXsc().

CalculateZommerfeld()

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

inline

Definition at line 439 of file G4DiffuseElastic.hh.

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

Referenced by BuildAngleTable(), GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), GetDiffuseElasticSumXsc(), and TestAngleTable().

DampFactor()

G4double G4DiffuseElastic::DampFactor ( G4double z )

inline

Definition at line 381 of file G4DiffuseElastic.hh.

{
  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;
}

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

GetCoulombElasticXsc()

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

inline

Definition at line 529 of file G4DiffuseElastic.hh.

{
  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;
}

Referenced by GetInvCoulombElasticXsc().

GetCoulombIntegralXsc()

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

inline

Definition at line 578 of file G4DiffuseElastic.hh.

{
  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);
           xsc          /= (1 - c1 + am)*(1 - c2 + am);
 
  return xsc;
}

GetCoulombTotalXsc()

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

inline

Definition at line 553 of file G4DiffuseElastic.hh.

{
  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;
}

GetDiffElasticProb()

G4double G4DiffuseElastic::GetDiffElasticProb

(

G4double

theta

)

Definition at line 380 of file G4DiffuseElastic.cc.

{
  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*G4Pow::GetInstance()->A13(A);
 
  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 if (fParticle == theNeutron)
  {
    diffuse =  0.63*fermi; // 1.63*fermi; //
    G4double k0 = 1*GeV/hbarc;
    diffuse *= k0/fWaveVector;
 
    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 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;  
 
  G4double lambda = 15.; // 15 ok
 
  //  G4double kgamma    = fWaveVector*gamma;   // wavek*delta;
 
  G4double kgamma    = lambda*(1.-G4Exp(-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.-G4Exp(-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;
}

Referenced by GetDiffuseElasticXsc().

GetDiffElasticSumProb()

G4double G4DiffuseElastic::GetDiffElasticSumProb

(

G4double

theta

)

Definition at line 468 of file G4DiffuseElastic.cc.

{
  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*G4Pow::GetInstance()->A13(A);
 
  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 if (fParticle == theNeutron)
  {
    diffuse = 0.63*fermi;
    // diffuse = 0.6*fermi;
    G4double k0 = 1*GeV/hbarc;
    diffuse *= k0/fWaveVector;
    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.-G4Exp(-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.-G4Exp(-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;
}

Referenced by GetDiffuseElasticSumXsc().

GetDiffElasticSumProbA()

G4double G4DiffuseElastic::GetDiffElasticSumProbA

(

G4double

alpha

)

Definition at line 574 of file G4DiffuseElastic.cc.

{
  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*G4Pow::GetInstance()->A13(A);
 
  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 if ( fParticle == theNeutron )
  {
    diffuse = 0.63*fermi;
    // diffuse = 0.6*fermi;
    // G4double k0 = 0.8*GeV/hbarc;
    // diffuse *= k0/fWaveVector;
    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.-G4Exp(-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. - G4Exp( -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;
}

Referenced by GetIntegrandFunction().

GetDiffuseElasticSumXsc()

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

Definition at line 244 of file G4DiffuseElastic.cc.

{
  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;
}

Referenced by GetInvElasticSumXsc().

GetDiffuseElasticXsc()

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

Definition at line 173 of file G4DiffuseElastic.cc.

{
  fParticle      = particle;
  fWaveVector    = momentum/hbarc;
  fAtomicWeight  = A;
  fAddCoulomb    = false;
  fNuclearRadius = CalculateNuclearRad(A);
 
  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
 
  return sigma;
}

Referenced by GetInvElasticXsc().

GetIntegrandFunction()

G4double G4DiffuseElastic::GetIntegrandFunction

(

G4double

theta

)

Definition at line 680 of file G4DiffuseElastic.cc.

{
  G4double result;
 
  result  = GetDiffElasticSumProbA(alpha);
 
  // result *= 2*pi*std::sin(theta);
 
  return result;
}

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

GetInvCoulombElasticXsc()

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

Definition at line 331 of file G4DiffuseElastic.cc.

{
  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;
}

GetInvElasticSumXsc()

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

Definition at line 279 of file G4DiffuseElastic.cc.

{
  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;
}

GetInvElasticXsc()

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

Definition at line 194 of file G4DiffuseElastic.cc.

{
  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;
}

GetNuclearRadius()

G4double G4DiffuseElastic::GetNuclearRadius ( )

inline

Definition at line 191 of file G4DiffuseElastic.hh.

{return fNuclearRadius;};

GetScatteringAngle()

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

Definition at line 1094 of file G4DiffuseElastic.cc.

{
 G4double x1, x2, y1, y2, randAngle;
 
  if( iAngle == 0 )
  {
    randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
    // iAngle++;
  }
  else
  {
    if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
    {
      iAngle = G4int((*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;
}

Referenced by SampleTableThetaCMS().

Initialise()

void G4DiffuseElastic::Initialise

(

)

Definition at line 139 of file G4DiffuseElastic.cc.

{
 
  // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
 
  const G4ElementTable* theElementTable = G4Element::GetElementTable();
 
  std::size_t jEl, numOfEl = G4Element::GetNumberOfElements();
 
  for( jEl = 0; jEl < numOfEl; ++jEl) // application element loop
  {
    fAtomicNumber = (*theElementTable)[jEl]->GetZ();     // atomic number
    fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( fAtomicNumber ) );
    fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
 
    if( verboseLevel > 0 ) 
    {   
      G4cout<<"G4DiffuseElastic::Initialise() the element: "
        <<(*theElementTable)[jEl]->GetName()<<G4endl;
    }
    fElementNumberVector.push_back(fAtomicNumber);
    fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
 
    BuildAngleTable();
    fAngleBank.push_back(fAngleTable);
  }  
  return;
}

InitialiseOnFly()

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

Definition at line 975 of file G4DiffuseElastic.cc.

{
  fAtomicNumber  = Z;     // atomic number
  fAtomicWeight  = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) );
 
  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
 
  if( verboseLevel > 0 )    
  {
    G4cout<<"G4DiffuseElastic::InitialiseOnFly() the element with Z = "
      <<Z<<"; and A = "<<A<<G4endl;
  }
  fElementNumberVector.push_back(fAtomicNumber);
 
  BuildAngleTable();
 
  fAngleBank.push_back(fAngleTable);
 
  return;
}

Referenced by SampleTableThetaCMS().

IntegralElasticProb()

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

Definition at line 696 of file G4DiffuseElastic.cc.

{
  G4double result;
  fParticle      = particle;
  fWaveVector    = momentum/hbarc;
  fAtomicWeight  = A;
 
  fNuclearRadius = CalculateNuclearRad(A);
 
 
  G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
 
  // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 
  result = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 
 
  return result;
}

IsApplicable()

G4bool G4DiffuseElastic::IsApplicable ( const G4HadProjectile & projectile, G4Nucleus & nucleus )

inlinevirtual

Reimplemented from G4HadronicInteraction.

Definition at line 232 of file G4DiffuseElastic.hh.

{
  if( ( projectile.GetDefinition() == G4Proton::Proton() ||
        projectile.GetDefinition() == G4Neutron::Neutron() ||
        projectile.GetDefinition() == G4PionPlus::PionPlus() ||
        projectile.GetDefinition() == G4PionMinus::PionMinus() ||
        projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
        projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
 
        nucleus.GetZ_asInt() >= 2 ) return true;
  else                              return false;
}

NeutronTuniform()

G4double G4DiffuseElastic::NeutronTuniform

(

G4int

Z

)

Definition at line 826 of file G4DiffuseElastic.cc.

{
  G4double elZ  = G4double(Z);
  elZ -= 1.;
  // G4double Tkin = 20.*G4Exp(-elZ/10.) + 1.;
  G4double Tkin = 12.*G4Exp(-elZ/10.) + 1.;
  return Tkin;
}

Referenced by SampleInvariantT().

SampleInvariantT()

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

virtual

Reimplemented from G4HadronicInteraction.

Definition at line 788 of file G4DiffuseElastic.cc.

{
  fParticle = aParticle;
  G4double m1 = fParticle->GetPDGMass(), t;
  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();
  
  if( aParticle == theNeutron)
  {
    G4double Tmax = NeutronTuniform( Z );
    G4double pCMS2 = momentumCMS*momentumCMS;
    G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1;
 
    if( Tkin <= Tmax )
    {
      t = 4.*pCMS2*G4UniformRand();
      // G4cout<<Tkin<<", "<<Tmax<<", "<<std::sqrt(t)<<";   ";
      return t;
    }
  }
  
  t = SampleTableT( aParticle,  momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
 
  return t;
}

SampleT()

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

Definition at line 721 of file G4DiffuseElastic.cc.

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

Referenced by SampleThetaLab().

SampleTableT()

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

Definition at line 840 of file G4DiffuseElastic.cc.

{
  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;
}

Referenced by SampleInvariantT().

SampleTableThetaCMS()

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

Definition at line 855 of file G4DiffuseElastic.cc.

{
  std::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 << "G4DiffuseElastic: 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++)
  {
    if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
  }
  if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;   // kinE is more then theMaxEnergy
  if ( iMomentum < 0 )           iMomentum = 0; // against negative index, kinE < theMinEnergy
 
  // G4cout<<"iMomentum = "<<iMomentum<<G4endl;
 
  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;
    // G4cout<<"randAngle = "<<randAngle<<G4endl;
  }
  // G4double angle = randAngle;
  // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
 
  if(randAngle < 0.) randAngle = 0.;
 
  return randAngle;
}

Referenced by SampleTableT().

SampleThetaCMS()

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

Definition at line 734 of file G4DiffuseElastic.cc.

{
  G4int i, iMax = 100;  
  G4double norm, theta1, theta2, thetaMax;
  G4double result = 0., sum = 0.;
 
  fParticle      = particle;
  fWaveVector    = momentum/hbarc;
  fAtomicWeight  = A;
 
  fNuclearRadius = CalculateNuclearRad(A);
  
  thetaMax = 10.174/fWaveVector/fNuclearRadius;
 
  if (thetaMax > pi) thetaMax = pi;
 
  G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
 
  // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 
  norm = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., thetaMax );
 
  norm *= G4UniformRand();
 
  for(i = 1; i <= iMax; i++)
  {
    theta1 = (i-1)*thetaMax/iMax; 
    theta2 = i*thetaMax/iMax;
    sum   += integral.Legendre10(this,&G4DiffuseElastic::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;
}

Referenced by SampleT().

SampleThetaLab()

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

Definition at line 1137 of file G4DiffuseElastic.cc.

{
  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 << "G4DiffuseElastic: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;
}

SetHEModelLowLimit()

void G4DiffuseElastic::SetHEModelLowLimit ( G4double value )

inline

Definition at line 256 of file G4DiffuseElastic.hh.

{
  lowEnergyLimitHE = value;
}

SetLowestEnergyLimit()

void G4DiffuseElastic::SetLowestEnergyLimit ( G4double value )

inline

Definition at line 266 of file G4DiffuseElastic.hh.

{
  lowestEnergyLimit = value;
}

SetPlabLowLimit()

void G4DiffuseElastic::SetPlabLowLimit ( G4double value )

inline

Definition at line 251 of file G4DiffuseElastic.hh.

{
  plabLowLimit = value;
}

SetQModelLowLimit()

void G4DiffuseElastic::SetQModelLowLimit ( G4double value )

inline

Definition at line 261 of file G4DiffuseElastic.hh.

{
  lowEnergyLimitQ = value;
}

SetRecoilKinEnergyLimit()

void G4DiffuseElastic::SetRecoilKinEnergyLimit ( G4double value )

inline

Definition at line 246 of file G4DiffuseElastic.hh.

{
  lowEnergyRecoilLimit = value;
}

TestAngleTable()

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

Definition at line 1344 of file G4DiffuseElastic.cc.

{
  fAtomicNumber  = Z;     // atomic number
  fAtomicWeight  = A;     // number of nucleons
  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
  
     
  
  G4cout<<"G4DiffuseElastic::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<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(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, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
    deltaL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
    deltaAG  = integral.AdaptiveGauss(this, &G4DiffuseElastic::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, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
  sumL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
  sumAG  = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction, 
                                       0., alpha2,epsilon);
  G4cout<<G4endl;
  G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
            <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
  */
  return;
}

ThetaCMStoThetaLab()

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

Definition at line 1225 of file G4DiffuseElastic.cc.

{
  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;
}

ThetaLabToThetaCMS()

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

Definition at line 1285 of file G4DiffuseElastic.cc.

{
  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:

• G4DiffuseElastic.hh
• G4DiffuseElastic.cc