G4DiffuseElastic Class Reference

#include <G4DiffuseElastic.hh>

Inheritance diagram for G4DiffuseElastic:

Public Member Functions
	G4DiffuseElastic ()
virtual	~G4DiffuseElastic ()
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	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")
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 58 of file G4DiffuseElastic.hh.

Constructor & Destructor Documentation

G4DiffuseElastic::G4DiffuseElastic

(

)

Definition at line 67 of file G4DiffuseElastic.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("DiffuseElastic"), fParticle(0)
{
  SetMinEnergy( 0.01*GeV );
  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;
}

G4DiffuseElastic::~G4DiffuseElastic ( )

virtual

Definition at line 107 of file G4DiffuseElastic.cc.

References G4PhysicsTable::clearAndDestroy().

{
  if(fEnergyVector) delete fEnergyVector;

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

Member Function Documentation

G4double G4DiffuseElastic::BesselJone ( G4double  z )

inline

Definition at line 311 of file G4DiffuseElastic.hh.

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

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

  modvalue = 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 = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);

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

G4double G4DiffuseElastic::BesselJzero ( G4double  z )

inline

Definition at line 259 of file G4DiffuseElastic.hh.

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

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

  modvalue = 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 = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
  }
  return bessel;
}

G4double G4DiffuseElastic::BesselOneByArg ( G4double  z )

inline

Definition at line 386 of file G4DiffuseElastic.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 G4DiffuseElastic::BuildAngleTable

(

)

Definition at line 922 of file G4DiffuseElastic.cc.

References test::a, CalculateAm(), CalculateZommerfeld(), GetIntegrandFunction(), G4PhysicsVector::GetLowEdgeEnergy(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), python.hepunit::hbarc, G4PhysicsTable::insertAt(), G4Integrator< T, F >::Legendre10(), G4PhysicsFreeVector::PutValue(), and z.

Referenced by Initialise(), and InitialiseOnFly().

{
  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 

    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<<"; sum = "<<sum<<G4endl;
    }
    fAngleTable->insertAt(i,angleVector);

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

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

inline

Definition at line 432 of file G4DiffuseElastic.hh.

References n.

Referenced by BuildAngleTable(), GetDiffuseElasticSumXsc(), 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 G4DiffuseElastic::CalculateNuclearRad ( G4double  A )

inline

Definition at line 447 of file G4DiffuseElastic.hh.

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

{
  G4double r0;

  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;

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

    fNuclearRadius = r0*std::pow(A, 0.27); // 0.27);
  }
  return fNuclearRadius;
}

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

inline

Definition at line 407 of file G4DiffuseElastic.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; 
}

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

inline

Definition at line 421 of file G4DiffuseElastic.hh.

Referenced by BuildAngleTable(), GetDiffuseElasticSumXsc(), and TestAngleTable().

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

  return fZommerfeld; 
}

G4double G4DiffuseElastic::DampFactor ( G4double  z )

inline

Definition at line 363 of file G4DiffuseElastic.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;
}

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

inline

Definition at line 471 of file G4DiffuseElastic.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 G4DiffuseElastic::GetCoulombIntegralXsc ( const G4ParticleDefinition *  particle, G4double  momentum, G4double  Z, G4double  theta1, G4double  theta2  )

inline

Definition at line 520 of file G4DiffuseElastic.hh.

References plottest35::c1, G4cout, G4endl, 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);
           xsc          /= (1 - c1 + am)*(1 - c2 + am);

  return xsc;
}

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

inline

Definition at line 495 of file G4DiffuseElastic.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 G4DiffuseElastic::GetDiffElasticProb

(

G4double

theta

)

Definition at line 363 of file G4DiffuseElastic.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 G4DiffuseElastic::GetDiffElasticSumProb

(

G4double

theta

)

Definition at line 440 of file G4DiffuseElastic.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 G4DiffuseElastic::GetDiffElasticSumProbA

(

G4double

alpha

)

Definition at line 535 of file G4DiffuseElastic.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 G4DiffuseElastic::GetDiffuseElasticSumXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 227 of file G4DiffuseElastic.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 G4DiffuseElastic::GetDiffuseElasticXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A
	)

Definition at line 156 of file G4DiffuseElastic.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 G4DiffuseElastic::GetIntegrandFunction

(

G4double

theta

)

Definition at line 631 of file G4DiffuseElastic.cc.

References GetDiffElasticSumProbA().

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

{
  G4double result;

  result  = GetDiffElasticSumProbA(alpha);

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

  return result;
}

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

Definition at line 314 of file G4DiffuseElastic.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 G4DiffuseElastic::GetInvElasticSumXsc	(	const G4ParticleDefinition *	particle,
		G4double	tMand,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 262 of file G4DiffuseElastic.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 G4DiffuseElastic::GetInvElasticXsc	(	const G4ParticleDefinition *	particle,
		G4double	theta,
		G4double	momentum,
		G4double	A,
		G4double	Z
	)

Definition at line 177 of file G4DiffuseElastic.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 G4DiffuseElastic::GetNuclearRadius ( )

inline

Definition at line 186 of file G4DiffuseElastic.hh.

{return fNuclearRadius;};

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

Definition at line 1007 of file G4DiffuseElastic.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;
}

void G4DiffuseElastic::Initialise

(

)

Definition at line 122 of file G4DiffuseElastic.cc.

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

{

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

  const G4ElementTable* theElementTable = G4Element::GetElementTable();

  size_t jEl, numOfEl = G4Element::GetNumberOfElements();

  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);

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

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

Definition at line 896 of file G4DiffuseElastic.cc.

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

Referenced by SampleTableThetaCMS().

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

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

  BuildAngleTable();

  fAngleBank.push_back(fAngleTable);

  return;
}

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

Definition at line 647 of file G4DiffuseElastic.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<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;

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

  return result;
}

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

virtual

Reimplemented from G4HadronElastic.

Definition at line 738 of file G4DiffuseElastic.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 G4DiffuseElastic::SampleT	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	A
	)

Definition at line 672 of file G4DiffuseElastic.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 G4DiffuseElastic::SampleTableT	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	Z,
		G4double	A
	)

Definition at line 763 of file G4DiffuseElastic.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 G4DiffuseElastic::SampleTableThetaCMS	(	const G4ParticleDefinition *	aParticle,
		G4double	p,
		G4double	Z,
		G4double	A
	)

Definition at line 778 of file G4DiffuseElastic.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 << "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);

  return randAngle;
}

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

Definition at line 685 of file G4DiffuseElastic.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<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;
}

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

Definition at line 1050 of file G4DiffuseElastic.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 << "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;
}

void G4DiffuseElastic::SetHEModelLowLimit ( G4double  value )

inline

Definition at line 238 of file G4DiffuseElastic.hh.

{
  lowEnergyLimitHE = value;
}

void G4DiffuseElastic::SetLowestEnergyLimit ( G4double  value )

inline

Definition at line 248 of file G4DiffuseElastic.hh.

{
  lowestEnergyLimit = value;
}

void G4DiffuseElastic::SetPlabLowLimit ( G4double  value )

inline

Definition at line 233 of file G4DiffuseElastic.hh.

{
  plabLowLimit = value;
}

void G4DiffuseElastic::SetQModelLowLimit ( G4double  value )

inline

Definition at line 243 of file G4DiffuseElastic.hh.

{
  lowEnergyLimitQ = value;
}

void G4DiffuseElastic::SetRecoilKinEnergyLimit ( G4double  value )

inline

Definition at line 228 of file G4DiffuseElastic.hh.

{
  lowEnergyRecoilLimit = value;
}

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

Definition at line 1257 of file G4DiffuseElastic.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<<"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;
}

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

Definition at line 1138 of file G4DiffuseElastic.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 G4DiffuseElastic::ThetaLabToThetaCMS	(	const G4DynamicParticle *	aParticle,
		G4double	tmass,
		G4double	thetaLab
	)

Definition at line 1198 of file G4DiffuseElastic.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:

• G4DiffuseElastic.hh
• G4DiffuseElastic.cc