G4eDPWAElasticDCS Class Reference

#include <G4eDPWAElasticDCS.hh>

Public Member Functions
	G4eDPWAElasticDCS (G4bool iselectron=true, G4bool isrestricted=false)
	~G4eDPWAElasticDCS ()
void	InitialiseForZ (std::size_t iz)
void	ComputeCSPerAtom (G4int iz, G4double ekin, G4double &elcs, G4double &tr1cs, G4double &tr2cs, G4double mumin=0.0, G4double mumax=1.0)
G4double	SampleCosineTheta (std::size_t iz, G4double lekin, G4double r1, G4double r2, G4double r3)
G4double	SampleCosineThetaRestricted (std::size_t iz, G4double lekin, G4double r1, G4double r2, G4double costMax, G4double costMin)
G4double	ComputeScatteringPowerCorrection (const G4MaterialCutsCouple *matcut, G4double ekin)
void	InitSCPCorrection (G4double lowEnergyLimit, G4double highEnergyLimit)

Detailed Description

Definition at line 106 of file G4eDPWAElasticDCS.hh.

Constructor & Destructor Documentation

G4eDPWAElasticDCS()

G4eDPWAElasticDCS::G4eDPWAElasticDCS	(	G4bool	iselectron = true,
		G4bool	isrestricted = false
	)

Definition at line 82 of file G4eDPWAElasticDCS.cc.

: fIsRestrictedSamplingRequired(isrestricted), fIsElectron(iselectron) {
  fDCS.resize(gMaxZ+1, nullptr);
  fDCSLow.resize(gMaxZ+1, nullptr);
  fSamplingTables.resize(gMaxZ+1, nullptr);
}

~G4eDPWAElasticDCS()

G4eDPWAElasticDCS::~G4eDPWAElasticDCS

(

)

Definition at line 91 of file G4eDPWAElasticDCS.cc.

                                      {
  for (std::size_t i=0; i<fDCS.size(); ++i) {
    if (fDCS[i]) delete fDCS[i];
  }
  for (std::size_t i=0; i<fDCSLow.size(); ++i) {
    if (fDCSLow[i]) delete fDCSLow[i];
  }
  for (std::size_t i=0; i<fSamplingTables.size(); ++i) {
    if (fSamplingTables[i]) delete fSamplingTables[i];
  }
  // clear scp correction data
  for (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
    if (fSCPCPerMatCuts[imc]) {
      fSCPCPerMatCuts[imc]->fVSCPC.clear();
      delete fSCPCPerMatCuts[imc];
    }
  }
  fSCPCPerMatCuts.clear();
}

Member Function Documentation

ComputeCSPerAtom()

void G4eDPWAElasticDCS::ComputeCSPerAtom	(	G4int	iz,
		G4double	ekin,
		G4double &	elcs,
		G4double &	tr1cs,
		G4double &	tr2cs,
		G4double	mumin = 0.0,
		G4double	mumax = 1.0
	)

Definition at line 270 of file G4eDPWAElasticDCS.cc.

                                                                         {
  // init all cross section values to zero;
  elcs  = 0.0;
  tr1cs = 0.0;
  tr2cs = 0.0;
  // make sure that mu(theta) = 0.5[1-cos(theta)] limits have appropriate vals
  mumin = std::max(0.0, std::min(1.0, mumin));
  mumax = std::max(0.0, std::min(1.0, mumax));
  if (mumin>=mumax) return;
  // make sure that kin. energy is within the available range (10 eV-100MeV)  
  const G4double lekin = std::max(gTheEnergies[0], std::min(gTheEnergies[gNumEnergies-1], G4Log(ekin)));
  // if the lower, denser in theta, DCS set should be used
  const G4bool isLowerGrid = (fIsElectron && lekin<gTheEnergies[gIndxEnergyLim]);
  const std::vector<G4double>& theMuVector = (isLowerGrid) ? gTheMus1    : gTheMus2;
  const G4Physics2DVector*     the2DDCS    = (isLowerGrid) ? fDCSLow[iz] : fDCS[iz];
  // find lower/upper mu bin of integration:
  // 0.0 <= mumin < 1.0 for sure here
  const std::size_t iMuStart = (mumin == 0.0) ? 0 : std::distance( theMuVector.begin(), std::upper_bound(theMuVector.begin(), theMuVector.end(), mumin) )-1 ;
  // 0.0 < mumax <= 1.0 for sure here
  const std::size_t iMuEnd   = (mumax == 1.0) ? theMuVector.size()-2 : std::distance( theMuVector.begin(), std::upper_bound(theMuVector.begin(), theMuVector.end(), mumax) )-1 ;
  // perform numerical integration of the DCS over the given [mumin, mumax]
  // interval (where mu(theta) = 0.5[1-cos(theta)]) to get the elastic, first
  std::size_t ix = 0;
  std::size_t iy = 0;
  for (std::size_t imu=iMuStart; imu<=iMuEnd; ++imu) {
    G4double elcsPar   = 0.0;
    G4double tr1csPar  = 0.0;
    G4double tr2csPar  = 0.0;
    const G4double low = (imu==iMuStart) ? mumin     : theMuVector[imu];
    const G4double del = (imu==iMuEnd)   ? mumax-low : theMuVector[imu+1]-low;
    ix = imu;
    for (std::size_t igl=0; igl<8; ++igl) {
      const double mu  = low + del*gXGL[igl];
      const double dcs = G4Exp(the2DDCS->Value(mu, lekin, ix, iy));
      elcsPar  += gWGL[igl]*dcs;             // elastic
      tr1csPar += gWGL[igl]*dcs*mu;          // first transport
      tr2csPar += gWGL[igl]*dcs*mu*(1.0-mu); // second transport
    }
    elcs  += del*elcsPar;
    tr1cs += del*tr1csPar;
    tr2cs += del*tr2csPar;
  }
  elcs  *=  2.0*CLHEP::twopi;
  tr1cs *=  4.0*CLHEP::twopi;
  tr2cs *= 12.0*CLHEP::twopi;
}

Referenced by G4eDPWACoulombScatteringModel::ComputeCrossSectionPerAtom().

ComputeScatteringPowerCorrection()

G4double G4eDPWAElasticDCS::ComputeScatteringPowerCorrection	(	const G4MaterialCutsCouple *	matcut,
		G4double	ekin
	)

Definition at line 573 of file G4eDPWAElasticDCS.cc.

                                                                   {
  const G4int imc  = matcut->GetIndex();
  G4double corFactor = 1.0;
  if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin<=fSCPCPerMatCuts[imc]->fPrCut) {
    return corFactor;
  }
  // get the scattering power correction factor
  const G4double lekin = G4Log(ekin);
  G4double remaining   = (lekin-fSCPCPerMatCuts[imc]->fLEmin)*fSCPCPerMatCuts[imc]->fILDel;
  std::size_t lindx    = (G4int)remaining;
  remaining           -= lindx;
  std::size_t imax     = fSCPCPerMatCuts[imc]->fVSCPC.size()-1;
  if (lindx>=imax) {
    corFactor = fSCPCPerMatCuts[imc]->fVSCPC[imax];
  } else {
    corFactor = fSCPCPerMatCuts[imc]->fVSCPC[lindx] + remaining*(fSCPCPerMatCuts[imc]->fVSCPC[lindx+1]-fSCPCPerMatCuts[imc]->fVSCPC[lindx]);
  }
  return corFactor;
}

Referenced by G4eDPWACoulombScatteringModel::ComputeCrossSectionPerAtom().

InitialiseForZ()

void G4eDPWAElasticDCS::InitialiseForZ

(

std::size_t

iz

)

Definition at line 114 of file G4eDPWAElasticDCS.cc.

                                                   {
  if (!gIsGridLoaded) {
    LoadGrid();
  }
  LoadDCSForZ((G4int)iz);
  BuildSmplingTableForZ((G4int)iz);
}

Referenced by G4eDPWACoulombScatteringModel::Initialise().

InitSCPCorrection()

void G4eDPWAElasticDCS::InitSCPCorrection	(	G4double	lowEnergyLimit,
		G4double	highEnergyLimit
	)

Definition at line 595 of file G4eDPWAElasticDCS.cc.

                                                                    {
  // get the material-cuts table
  G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
  std::size_t numMatCuts            = thePCTable->GetTableSize();
  // clear container if any
  for (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
    if (fSCPCPerMatCuts[imc]) {
      fSCPCPerMatCuts[imc]->fVSCPC.clear();
      delete fSCPCPerMatCuts[imc];
      fSCPCPerMatCuts[imc] = nullptr;
    }
  }
  //
  // set size of the container and create the corresponding data structures
  fSCPCPerMatCuts.resize(numMatCuts,nullptr);
  // loop over the material-cuts and create scattering power correction data structure for each
  for (G4int imc=0; imc<(G4int)numMatCuts; ++imc) {
    const G4MaterialCutsCouple *matCut =  thePCTable->GetMaterialCutsCouple(imc);
    const G4Material* mat  = matCut->GetMaterial();
    // get e- production cut in the current material-cuts in energy
    const G4double ecut    = (*(thePCTable->GetEnergyCutsVector(idxG4ElectronCut)))[matCut->GetIndex()];
    const G4double limit   = fIsElectron ? 2.0*ecut : ecut;
    const G4double min     = std::max(limit,lowEnergyLimit);
    const G4double max     = highEnergyLimit;
    if (min>=max) {
      fSCPCPerMatCuts[imc] = new SCPCorrection();
      fSCPCPerMatCuts[imc]->fIsUse = false;
      fSCPCPerMatCuts[imc]->fPrCut = min;
      continue;
    }
    G4int numEbins         = fNumSPCEbinPerDec*G4lrint(std::log10(max/min));
    numEbins               = std::max(numEbins,3);
    const G4double lmin    = G4Log(min);
    const G4double ldel    = G4Log(max/min)/(numEbins-1.0);
    fSCPCPerMatCuts[imc]   = new SCPCorrection();
    fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbins,1.0);
    fSCPCPerMatCuts[imc]->fIsUse = true;
    fSCPCPerMatCuts[imc]->fPrCut = min;
    fSCPCPerMatCuts[imc]->fLEmin = lmin;
    fSCPCPerMatCuts[imc]->fILDel = 1./ldel;
    // compute Moliere material dependet parameetrs
    G4double moliereBc     = 0.0;
    G4double moliereXc2    = 0.0;
    ComputeMParams(mat, moliereBc, moliereXc2);
    // compute scattering power correction over the enrgy grid
    for (G4int ie=0; ie<numEbins; ++ie) {
      const G4double ekin = G4Exp(lmin+ie*ldel);
      G4double scpCorr    = 1.0;
      // compute correction factor: I.Kawrakow, Med.Phys.24,505-517(1997)(Eqs(32-37)
      if (ie>0) {
         const G4double tau     = ekin/CLHEP::electron_mass_c2;
         const G4double tauCut  = ecut/CLHEP::electron_mass_c2;
         // Moliere's screening parameter
         const G4double A       = moliereXc2/(4.0*tau*(tau+2.)*moliereBc);
         const G4double gr      = (1.+2.*A)*G4Log(1.+1./A)-2.;
         const G4double dum0    = (tau+2.)/(tau+1.);
         const G4double dum1    = tau+1.;
         G4double gm  = G4Log(0.5*tau/tauCut) + (1.+dum0*dum0)*G4Log(2.*(tau-tauCut+2.)/(tau+4.))
                        - 0.25*(tau+2.)*( tau+2.+2.*(2.*tau+1.)/(dum1*dum1))*
                          G4Log((tau+4.)*(tau-tauCut)/tau/(tau-tauCut+2.))
                        + 0.5*(tau-2*tauCut)*(tau+2.)*(1./(tau-tauCut)-1./(dum1*dum1));
         if (gm<gr) {
           gm = gm/gr;
         } else {
           gm = 1.;
         }
         const G4double z0 = matCut->GetMaterial()->GetIonisation()->GetZeffective();
         scpCorr = 1.-gm*z0/(z0*(z0+1.));
      }
      fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCorr;
    }
  }
}

Referenced by G4eDPWACoulombScatteringModel::Initialise().

SampleCosineTheta()

G4double G4eDPWAElasticDCS::SampleCosineTheta	(	std::size_t	iz,
		G4double	lekin,
		G4double	r1,
		G4double	r2,
		G4double	r3
	)

Definition at line 405 of file G4eDPWAElasticDCS.cc.

                                                               {
  lekin = std::max(gTheEnergies[0], std::min(gTheEnergies[gNumEnergies-1], lekin));
  // determine the discrete ekin sampling table to be used:
  //  - statistical interpolation (i.e. linear) on log energy scale
  const G4double      rem = (lekin-gLogMinEkin)*gInvDelLogEkin;
  const auto            k = (std::size_t)rem;
  const std::size_t iekin = (r1 < rem-k) ? k+1 : k;
  // sample the mu(t)=0.5(1-cos(t))
  const double mu    = SampleMu(iz, iekin, r2, r3);
  return std::max(-1.0, std::min(1.0, 1.0-2.0*mu));
}

Referenced by G4eDPWACoulombScatteringModel::SampleSecondaries().

SampleCosineThetaRestricted()

G4double G4eDPWAElasticDCS::SampleCosineThetaRestricted	(	std::size_t	iz,
		G4double	lekin,
		G4double	r1,
		G4double	r2,
		G4double	costMax,
		G4double	costMin
	)

Definition at line 428 of file G4eDPWAElasticDCS.cc.

                                                                                   {
  // costMin corresponds to mu-max while costMax to mu-min: mu(t)=0.5[1-cos(t)]
  lekin = std::max(gTheEnergies[0], std::min(gTheEnergies[gNumEnergies-1], lekin));
  // determine the discrete ekin sampling table to be used:
  //  - statistical interpolation (i.e. linear) on log energy scale
  const G4double      rem = (lekin-gLogMinEkin)*gInvDelLogEkin;
  const auto            k = (size_t)rem;
  const std::size_t iekin = (r1 < rem-k) ? k : k+1;
  // sample the mu(t)=0.5(1-cos(t))
  const G4double mu  = SampleMu(iz, iekin, r2, 0.5*(1.0-costMax), 0.5*(1.0-costMin));
  return std::max(-1.0, std::min(1.0, 1.0-2.0*mu));
}

Referenced by G4eDPWACoulombScatteringModel::SampleSecondaries().

---

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

• G4eDPWAElasticDCS.hh
• G4eDPWAElasticDCS.cc