G4WentzelVIRelXSection Class Reference

#include <G4WentzelVIRelXSection.hh>

Public Member Functions
	G4WentzelVIRelXSection ()
virtual	~G4WentzelVIRelXSection ()
void	Initialise (const G4ParticleDefinition *, G4double CosThetaLim)
void	SetupParticle (const G4ParticleDefinition *)
G4double	SetupTarget (G4int Z, G4double cut=DBL_MAX)
G4double	ComputeTransportCrossSectionPerAtom (G4double CosThetaMax)
G4ThreeVector	SampleSingleScattering (G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio=0.0)
G4double	ComputeNuclearCrossSection (G4double CosThetaMin, G4double CosThetaMax)
G4double	ComputeElectronCrossSection (G4double CosThetaMin, G4double CosThetaMax)
G4double	SetupKinematic (G4double kinEnergy, const G4Material *mat)
void	SetTargetMass (G4double value)
G4double	GetMomentumSquare () const
G4double	GetCosThetaNuc () const
G4double	GetCosThetaElec () const

Detailed Description

Definition at line 71 of file G4WentzelVIRelXSection.hh.

Constructor & Destructor Documentation

G4WentzelVIRelXSection::G4WentzelVIRelXSection

(

)

Definition at line 66 of file G4WentzelVIRelXSection.cc.

References python.hepunit::classic_electr_radius, G4Electron::Electron(), python.hepunit::electron_mass_c2, python.hepunit::eV, python.hepunit::fine_structure_const, G4NistManager::GetA27(), G4Pow::GetInstance(), Initialise(), G4NistManager::Instance(), python.hepunit::MeV, G4Positron::Positron(), G4Proton::Proton(), python.hepunit::proton_mass_c2, python.hepunit::twopi, test::x, and G4Pow::Z13().

                                               :
  numlimit(0.1),
  nwarnings(0),
  nwarnlimit(50),
  alpha2(fine_structure_const*fine_structure_const)
{
  fNistManager = G4NistManager::Instance();
  fG4pow = G4Pow::GetInstance();
  theElectron = G4Electron::Electron();
  thePositron = G4Positron::Positron();
  theProton   = G4Proton::Proton();
  lowEnergyLimit = 1.0*eV;
  G4double p0 = electron_mass_c2*classic_electr_radius;
  coeff = twopi*p0*p0;
  particle = 0;

  // Thomas-Fermi screening radii
  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282

  if(0.0 == ScreenRSquare[0]) {
    G4double a0 = electron_mass_c2/0.88534; 
    G4double constn = 6.937e-6/(MeV*MeV);

    ScreenRSquare[0] = alpha2*a0*a0;
    for(G4int j=1; j<100; ++j) {
      G4double x = a0*fG4pow->Z13(j);
      //ScreenRSquare[j] = 0.5*(1 + exp(-j*j*0.001))*alpha2*x*x;
      ScreenRSquare[j] = 0.5*alpha2*x*x;
      x = fNistManager->GetA27(j);
      FormFactor[j] = constn*x*x;
    } 
  }
  currentMaterial = 0;
  elecXSRatio = factB = factD = formfactA = screenZ = 0.0;
  cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;

  factB1= 0.5*CLHEP::pi*fine_structure_const;

  Initialise(theElectron, 1.0);
  targetMass = proton_mass_c2;
}

G4WentzelVIRelXSection::~G4WentzelVIRelXSection ( )

virtual

Definition at line 110 of file G4WentzelVIRelXSection.cc.

{}

Member Function Documentation

G4double G4WentzelVIRelXSection::ComputeElectronCrossSection ( G4double  CosThetaMin, G4double  CosThetaMax  )

inline

Definition at line 236 of file G4WentzelVIRelXSection.hh.

References G4INCL::Math::max().

Referenced by G4hCoulombScatteringModel::ComputeCrossSectionPerAtom().

{
  G4double xsec = 0.0;
  G4double cost1 = std::max(cosTMin,cosTetMaxElec);
  G4double cost2 = std::max(cosTMax,cosTetMaxElec);
  if(cost1 > cost2) {
    xsec = kinFactor*(cost1 - cost2)/((1.0 - cost1 + screenZ)*(1.0 - cost2 + screenZ));
  }
  return xsec;
}

G4double G4WentzelVIRelXSection::ComputeNuclearCrossSection ( G4double  CosThetaMin, G4double  CosThetaMax  )

inline

Definition at line 222 of file G4WentzelVIRelXSection.hh.

Referenced by G4hCoulombScatteringModel::ComputeCrossSectionPerAtom().

{
  G4double xsec = 0.0;
  if(cosTMax < cosTMin) {
    xsec = targetZ*kinFactor*(cosTMin - cosTMax)/
      ((1.0 - cosTMin + screenZ)*(1.0 - cosTMax + screenZ));
  }
  return xsec;
}

G4double G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom

(

G4double

CosThetaMax

)

Definition at line 193 of file G4WentzelVIRelXSection.cc.

References G4cout, G4endl, G4Log(), G4ParticleDefinition::GetParticleName(), G4INCL::Math::max(), and test::x.

Referenced by G4WentzelVIRelModel::ComputeCrossSectionPerAtom().

{
  G4double xsec = 0.0;
  if(cosTMax >= 1.0) { return xsec; }
 
  G4double xSection = 0.0;
  G4double x = 0; 
  G4double y = 0;
  G4double x1= 0;
  G4double x2= 0;
  G4double xlog = 0.0;

  G4double costm = std::max(cosTMax,cosTetMaxElec); 
  G4double fb = screenZ*factB;

  // scattering off electrons
  if(costm < 1.0) {
    x = (1.0 - costm)/screenZ;
    if(x < numlimit) { 
      x2 = 0.5*x*x;
      y  = x2*(1.0 - 1.3333333*x + 3*x2);
      if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
    } else { 
      x1= x/(1 + x);
      xlog = G4Log(1.0 + x);  
      y = xlog - x1; 
      if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
    }

    if(y < 0.0) {
      ++nwarnings;
      if(nwarnings < nwarnlimit) {
    G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
           << G4endl;
    G4cout << "y= " << y 
           << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 
           << " Z= " << targetZ << "  " 
           << particle->GetParticleName() << G4endl;
    G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ 
           << " x= " << x << G4endl;
      }
      y = 0.0;
    }
    xSection = y;
  }
  /* 
  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ 
     << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
     << " cut(MeV)= " << ecut/MeV  
     << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 
     << " zmaxN= " << (1.0 - cosThetaMax)/screenZ 
         << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
  */
  // scattering off nucleus
  if(cosTMax < 1.0) {
    x = (1.0 - cosTMax)/screenZ;
    if(x < numlimit) { 
      x2 = 0.5*x*x;
      y  = x2*(1.0 - 1.3333333*x + 3*x2); 
      if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
    } else { 
      x1= x/(1 + x);
      xlog = G4Log(1.0 + x);  
      y = xlog - x1; 
      if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
    }

    if(y < 0.0) {
      ++nwarnings;
      if(nwarnings < nwarnlimit) {
    G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
           << G4endl;
    G4cout << "y= " << y 
           << " e(MeV)= " << tkin << " Z= " << targetZ << "  " 
           << particle->GetParticleName() << G4endl;
    G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ 
           << " x= " << " x1= " << x1 <<G4endl;
      }
      y = 0.0;
    }
    xSection += y*targetZ; 
  }
  xSection *= kinFactor;
  /*
  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 
     << " screenZ= " << screenZ << " formF= " << formfactA 
     << " for " << particle->GetParticleName() 
     << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
     << " x= " << x 
     << G4endl;
  */
  return xSection; 
}

G4double G4WentzelVIRelXSection::GetCosThetaElec ( ) const

inline

Definition at line 214 of file G4WentzelVIRelXSection.hh.

{
  return cosTetMaxElec;
}

G4double G4WentzelVIRelXSection::GetCosThetaNuc ( ) const

inline

Definition at line 207 of file G4WentzelVIRelXSection.hh.

{
  return cosTetMaxNuc;
}

G4double G4WentzelVIRelXSection::GetMomentumSquare ( ) const

inline

Definition at line 200 of file G4WentzelVIRelXSection.hh.

Referenced by G4hCoulombScatteringModel::SampleSecondaries().

{
  return mom2;
}

void G4WentzelVIRelXSection::Initialise	(	const G4ParticleDefinition *	p,
		G4double	CosThetaLim
	)

Definition at line 115 of file G4WentzelVIRelXSection.cc.

References test::a, DBL_MAX, G4LossTableManager::FactorForAngleLimit(), G4LossTableManager::Instance(), and SetupParticle().

Referenced by G4WentzelVIRelXSection(), G4hCoulombScatteringModel::Initialise(), and G4WentzelVIRelModel::Initialise().

{
  SetupParticle(p);
  tkin = mom2 = momCM2 = 0.0;
  ecut = etag = DBL_MAX;
  targetZ = 0;
  cosThetaMax = CosThetaLim;
  G4double a = 
    G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
  factorA2 = 0.5*a*a;
  currentMaterial = 0;
}

G4ThreeVector G4WentzelVIRelXSection::SampleSingleScattering	(	G4double	CosThetaMin,
		G4double	CosThetaMax,
		G4double	elecRatio = 0.0
	)

Definition at line 290 of file G4WentzelVIRelXSection.cc.

References fm, G4UniformRand, G4INCL::Math::max(), CLHEP::Hep3Vector::set(), python.hepunit::twopi, and test::v.

Referenced by G4WentzelVIRelModel::SampleScattering(), and G4hCoulombScatteringModel::SampleSecondaries().

{
  G4ThreeVector v(0.0,0.0,1.0);
 
  G4double formf = formfactA;
  G4double cost1 = cosTMin;
  G4double cost2 = cosTMax;
  if(elecRatio > 0.0) {
    if(G4UniformRand() <= elecRatio) {
      formf = 0.0;
      cost1 = std::max(cost1,cosTetMaxElec);
      cost2 = std::max(cost2,cosTetMaxElec);
    }
  }
  if(cost1 < cost2) { return v; }

  G4double w1 = 1. - cost1 + screenZ;
  G4double w2 = 1. - cost2 + screenZ;
  G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;

  //G4double fm =  1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass);
  G4double fm =  1.0 + formf*z1;
  //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm );
  G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm );
  // "false" scattering
  if( G4UniformRand() > grej ) { return v; }
    // } 
  G4double cost = 1.0 - z1;

  if(cost > 1.0)       { cost = 1.0; }
  else if(cost < -1.0) { cost =-1.0; }
  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
  //G4cout << "sint= " << sint << G4endl;
  G4double phi  = twopi*G4UniformRand();
  G4double vx1 = sint*cos(phi);
  G4double vy1 = sint*sin(phi);

  // only direction is changed
  v.set(vx1,vy1,cost);
  return v;
}

void G4WentzelVIRelXSection::SetTargetMass ( G4double  value )

inline

Definition at line 192 of file G4WentzelVIRelXSection.hh.

Referenced by G4hCoulombScatteringModel::SampleSecondaries(), and SetupTarget().

{
  targetMass = value;
  factD = std::sqrt(mom2)/value;
}

G4double G4WentzelVIRelXSection::SetupKinematic ( G4double  kinEnergy, const G4Material *  mat  )

inline

Definition at line 176 of file G4WentzelVIRelXSection.hh.

References G4IonisParamMat::GetInvA23(), G4Material::GetIonisation(), and G4INCL::Math::max().

Referenced by G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4WentzelVIRelModel::ComputeCrossSectionPerAtom(), G4WentzelVIRelModel::ComputeGeomPathLength(), G4WentzelVIRelModel::ComputeTruePathLengthLimit(), and G4WentzelVIRelModel::ComputeTrueStepLength().

{
  if(ekin != tkin || mat != currentMaterial) { 
    currentMaterial = mat;
    tkin  = ekin;
    mom2  = tkin*(tkin + 2.0*mass);
    invbeta2 = 1.0 +  mass*mass/mom2;
    factB = spin/invbeta2;
    cosTetMaxNuc = 
    std::max(cosThetaMax,1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2);
  } 
  return cosTetMaxNuc;
}

void G4WentzelVIRelXSection::SetupParticle

(

const G4ParticleDefinition *

p

)

Definition at line 131 of file G4WentzelVIRelXSection.cc.

References python.hepunit::eplus, G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and G4ParticleDefinition::GetPDGSpin().

Referenced by Initialise(), and G4hCoulombScatteringModel::SetupParticle().

{
  particle = p;
  mass = particle->GetPDGMass();
  spin = particle->GetPDGSpin();
  if(0.0 != spin) { spin = 0.5; }
  G4double q = std::fabs(particle->GetPDGCharge()/eplus);
  chargeSquare = q*q;
  charge3 = chargeSquare*q;
  tkin = 0.0;
  currentMaterial = 0;
  targetZ = 0;
}

G4double G4WentzelVIRelXSection::SetupTarget	(	G4int	Z,
		G4double	cut = DBL_MAX
	)

Definition at line 148 of file G4WentzelVIRelXSection.cc.

References G4NistManager::GetAtomicMassAmu(), G4INCL::Math::min(), and SetTargetMass().

Referenced by G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4WentzelVIRelModel::ComputeCrossSectionPerAtom(), and G4WentzelVIRelModel::SampleScattering().

{
  G4double cosTetMaxNuc2 = cosTetMaxNuc;
  if(Z != targetZ || tkin != etag) {
    etag    = tkin; 
    targetZ = Z;
    if(targetZ > 99) { targetZ = 99; }
    SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2);
    //G4double tmass2 = targetMass*targetMass;
    //G4double etot = tkin + mass;
    //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass;
    //momCM2 = mom2*tmass2/invmass2;
    //gam0pcmp = (etot + targetMass)*targetMass/invmass2;
    //pcmp2    = tmass2/invmass2;

    kinFactor = coeff*targetZ*chargeSquare*invbeta2/mom2;

    screenZ = ScreenRSquare[targetZ]/mom2;
    if(Z > 1) {
      screenZ *= std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare);
      /*
      if(mass > MeV) {
    screenZ *= std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare);
      } else {
    G4double tau = tkin/mass;
    //  screenZ *= std::min(Z*invbeta2,
    screenZ *= std::min(Z*1.13,
        (1.13 +3.76*Z*Z*invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(Z)))));
      }
      */
    }
    if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
      cosTetMaxNuc2 = 0.0;
    }
    formfactA = FormFactor[targetZ]*mom2;

    cosTetMaxElec = 1.0;
    ComputeMaxElectronScattering(cut); 
  }
  return cosTetMaxNuc2;
} 

---

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

• G4WentzelVIRelXSection.hh
• G4WentzelVIRelXSection.cc