G4eDPWAElasticDCS Class Reference

#include <G4eDPWAElasticDCS.hh>

Data Structures
struct	OneSamplingTable
struct	SCPCorrection

Public Member Functions
void	ComputeCSPerAtom (G4int iz, G4double ekin, G4double &elcs, G4double &tr1cs, G4double &tr2cs, G4double mumin=0.0, G4double mumax=1.0)
G4double	ComputeScatteringPowerCorrection (const G4MaterialCutsCouple *matcut, G4double ekin)
	G4eDPWAElasticDCS (G4bool iselectron=true, G4bool isrestricted=false)
void	InitialiseForZ (std::size_t iz)
void	InitSCPCorrection (G4double lowEnergyLimit, G4double highEnergyLimit)
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)
	~G4eDPWAElasticDCS ()

Private Member Functions
void	BuildSmplingTableForZ (G4int iz)
void	ComputeMParams (const G4Material *mat, G4double &theBc, G4double &theXc2)
G4double	FindCumValue (G4double u, const OneSamplingTable &stable, const std::vector< G4double > &uvect)
const G4String &	FindDirectoryPath ()
void	LoadDCSForZ (G4int iz)
void	LoadGrid ()
void	ReadCompressedFile (G4String fname, std::istringstream &iss)
G4double	SampleMu (std::size_t izet, std::size_t ie, G4double r1, G4double muMin, G4double muMax)
G4double	SampleMu (std::size_t izet, std::size_t ie, G4double r1, G4double r2)

Private Attributes
std::vector< G4Physics2DVector * >	fDCS
std::vector< G4Physics2DVector * >	fDCSLow
G4bool	fIsElectron
G4bool	fIsRestrictedSamplingRequired
const G4int	fNumSPCEbinPerDec = 3
std::vector< std::vector< OneSamplingTable > * >	fSamplingTables
std::vector< SCPCorrection * >	fSCPCPerMatCuts

Static Private Attributes
static G4String	gDataDirectory = ""
static std::size_t	gIndxEnergyLim = 35
static G4double	gInvDelLogEkin = 1.0
static G4bool	gIsGridLoaded = false
static G4double	gLogMinEkin = 1.0
static constexpr std::size_t	gMaxZ = 103
static std::size_t	gNumEnergies = 106
static std::size_t	gNumThetas1 = 247
static std::size_t	gNumThetas2 = 128
static std::vector< G4double >	gTheEnergies
static std::vector< G4double >	gTheMus1
static std::vector< G4double >	gTheMus2
static std::vector< G4double >	gTheU1
static std::vector< G4double >	gTheU2
static const G4double	gWGL [8]
static const G4double	gXGL [8]

Detailed Description

Definition at line 107 of file G4eDPWAElasticDCS.hh.

Constructor & Destructor Documentation

G4eDPWAElasticDCS()

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

Definition at line 80 of file G4eDPWAElasticDCS.cc.

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

References fDCS, fDCSLow, fSamplingTables, and gMaxZ.

~G4eDPWAElasticDCS()

G4eDPWAElasticDCS::~G4eDPWAElasticDCS

(

)

Definition at line 89 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();
}

References fDCS, fDCSLow, fSamplingTables, and fSCPCPerMatCuts.

Member Function Documentation

BuildSmplingTableForZ()

void G4eDPWAElasticDCS::BuildSmplingTableForZ ( G4int  iz )

private

Definition at line 351 of file G4eDPWAElasticDCS.cc.

                                                      {
  // Check if it has already been done:
  if (fSamplingTables[iz]) return;
  // Do it otherwise:
  // allocate space
  std::vector<OneSamplingTable>* sTables = new std::vector<OneSamplingTable>(gNumEnergies);
  // read compressed sampling table data
  std::ostringstream oss;
  const G4String fname = fIsElectron ? "stables/el/" : "stables/pos/";
  oss << FindDirectoryPath() << fname << "stable_" << iz;
  std::istringstream fin(std::ios::in);
  ReadCompressedFile(oss.str(), fin);
  std::size_t ndata = 0;
  for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
    OneSamplingTable& aTable = (*sTables)[ie];
    // #data in this table
    fin >> ndata;
    aTable.SetSize(ndata, !fIsRestrictedSamplingRequired);
    // the A screening parameter value used for transformation of mu to u
    fin >> aTable.fScreenParA;
    // load data: Alias(W,I) + RatIn(Cum, A, B)
    if (!fIsRestrictedSamplingRequired)  {
      for (std::size_t id=0; id<ndata; ++id) {
        fin >> aTable.fW[id];
      }
      for (std::size_t id=0; id<ndata; ++id) {
        fin >> aTable.fI[id];
      }
    }
    for (std::size_t id=0; id<ndata; ++id) {
      fin >> aTable.fCum[id];
    }
    for (std::size_t id=0; id<ndata; ++id) {
      fin >> aTable.fA[id];
    }
    for (std::size_t id=0; id<ndata; ++id) {
      fin >> aTable.fB[id];
    }
  }
  fSamplingTables[iz] = sTables;
}

References G4eDPWAElasticDCS::OneSamplingTable::fA, G4eDPWAElasticDCS::OneSamplingTable::fB, G4eDPWAElasticDCS::OneSamplingTable::fCum, G4eDPWAElasticDCS::OneSamplingTable::fI, FindDirectoryPath(), fIsElectron, fIsRestrictedSamplingRequired, test::fname, fSamplingTables, G4eDPWAElasticDCS::OneSamplingTable::fScreenParA, G4eDPWAElasticDCS::OneSamplingTable::fW, gNumEnergies, ReadCompressedFile(), and G4eDPWAElasticDCS::OneSamplingTable::SetSize().

Referenced by InitialiseForZ().

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

References fDCS, fDCSLow, fIsElectron, G4Exp(), G4Log(), gIndxEnergyLim, gNumEnergies, gTheEnergies, gTheMus1, gTheMus2, gWGL, gXGL, G4INCL::Math::max(), G4INCL::Math::min(), CLHEP::twopi, and G4Physics2DVector::Value().

Referenced by G4eDPWACoulombScatteringModel::ComputeCrossSectionPerAtom().

ComputeMParams()

void G4eDPWAElasticDCS::ComputeMParams ( const G4Material *  mat, G4double &  theBc, G4double &  theXc2  )

private

Definition at line 670 of file G4eDPWAElasticDCS.cc.

                                                         {
   const G4double const1   = 7821.6;      // [cm2/g]
   const G4double const2   = 0.1569;      // [cm2 MeV2 / g]
   const G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
   //   G4double xi   = 1.0;
   const G4ElementVector* theElemVect     = mat->GetElementVector();
   const G4int            numelems        = mat->GetNumberOfElements();
   //
   const G4double*  theNbAtomsPerVolVect  = mat->GetVecNbOfAtomsPerVolume();
   G4double         theTotNbAtomsPerVol   = mat->GetTotNbOfAtomsPerVolume();
   //
   G4double zs = 0.0;
   G4double zx = 0.0;
   G4double ze = 0.0;
   G4double sa = 0.0;
   //
   for(G4int ielem = 0; ielem < numelems; ielem++) {
     const G4double zet  = (*theElemVect)[ielem]->GetZ();
     const G4double iwa  = (*theElemVect)[ielem]->GetN();
     const G4double ipz  = theNbAtomsPerVolVect[ielem]/theTotNbAtomsPerVol;
     const G4double dum  = ipz*zet*(zet+1.0);
     zs           += dum;
     ze           += dum*(-2.0/3.0)*G4Log(zet);
     zx           += dum*G4Log(1.0+3.34*finstrc2*zet*zet);
     sa           += ipz*iwa;
   }
   const G4double density = mat->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
   //
   theBc  = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs);  //[1/cm]
   theXc2 = const2*density*zs/sa;  // [MeV2/cm]
   // change to Geant4 internal units of 1/length and energ2/length
   theBc  *= 1.0/CLHEP::cm;
   theXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
}

References CLHEP::cm, CLHEP::cm3, CLHEP::g, G4Exp(), G4Log(), G4Material::GetDensity(), G4Material::GetElementVector(), G4Material::GetNumberOfElements(), G4Material::GetTotNbOfAtomsPerVolume(), G4Material::GetVecNbOfAtomsPerVolume(), and CLHEP::MeV.

Referenced by InitSCPCorrection().

ComputeScatteringPowerCorrection()

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

Definition at line 571 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;
  G4int    lindx       = (G4int)remaining;
  remaining           -= lindx;
  G4int    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;
}

References fSCPCPerMatCuts, G4Log(), G4MaterialCutsCouple::GetIndex(), and imax.

Referenced by G4eDPWACoulombScatteringModel::ComputeCrossSectionPerAtom().

FindCumValue()

G4double G4eDPWAElasticDCS::FindCumValue ( G4double  u, const OneSamplingTable &  stable, const std::vector< G4double > &  uvect  )

private

Definition at line 466 of file G4eDPWAElasticDCS.cc.

                                                                  {
  const std::size_t iLow = std::distance( uvect.begin(), std::upper_bound(uvect.begin(), uvect.end(), u) )-1;
  const G4double    tau  = (u-uvect[iLow])/(uvect[iLow+1]-uvect[iLow]); // Note: I could store 1/(fX[iLow+1]-fX[iLow])
  const G4double    dum0 = (1.0+stable.fA[iLow]*(1.0-tau)+stable.fB[iLow]);
  const G4double    dum1 = 2.0*stable.fB[iLow]*tau;
  const G4double    dum2 = 1.0 - std::sqrt(std::max(0.0, 1.0-2.0*dum1*tau/(dum0*dum0)));
  return std::min(stable.fCum[iLow+1], std::max(stable.fCum[iLow], stable.fCum[iLow]+dum0*dum2*(stable.fCum[iLow+1]-stable.fCum[iLow])/dum1 ));
}

References G4eDPWAElasticDCS::OneSamplingTable::fA, G4eDPWAElasticDCS::OneSamplingTable::fB, G4eDPWAElasticDCS::OneSamplingTable::fCum, G4INCL::Math::max(), and G4INCL::Math::min().

Referenced by SampleMu().

FindDirectoryPath()

const G4String & G4eDPWAElasticDCS::FindDirectoryPath ( )

private

Definition at line 500 of file G4eDPWAElasticDCS.cc.

                                                     {
  // check environment variable
  if (gDataDirectory.empty()) {
    const char* path = std::getenv("G4LEDATA");
    if (path) {
      std::ostringstream ost;
      ost << path << "/dpwa/";
      gDataDirectory = ost.str();
    } else {
      G4Exception("G4eDPWAElasticDCS::FindDirectoryPath()","em0006",
                  FatalException,
                  "Environment variable G4LEDATA not defined");
    }
  }
  return gDataDirectory;
}

References FatalException, G4Exception(), and gDataDirectory.

Referenced by BuildSmplingTableForZ(), LoadDCSForZ(), and LoadGrid().

InitialiseForZ()

void G4eDPWAElasticDCS::InitialiseForZ

(

std::size_t

iz

)

Definition at line 112 of file G4eDPWAElasticDCS.cc.

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

References BuildSmplingTableForZ(), gIsGridLoaded, LoadDCSForZ(), and LoadGrid().

Referenced by G4eDPWACoulombScatteringModel::Initialise().

InitSCPCorrection()

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

Definition at line 593 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 (std::size_t imc=0; imc<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;
    }
  }
}

References A, ComputeMParams(), CLHEP::electron_mass_c2, fIsElectron, fNumSPCEbinPerDec, fSCPCPerMatCuts, G4Exp(), G4Log(), G4lrint(), G4ProductionCutsTable::GetEnergyCutsVector(), G4MaterialCutsCouple::GetIndex(), G4Material::GetIonisation(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4IonisParamMat::GetZeffective(), idxG4ElectronCut, G4INCL::Math::max(), G4INCL::Math::min(), and G4InuclParticleNames::z0.

Referenced by G4eDPWACoulombScatteringModel::Initialise().

LoadDCSForZ()

void G4eDPWAElasticDCS::LoadDCSForZ ( G4int  iz )

private

Definition at line 175 of file G4eDPWAElasticDCS.cc.

                                            {
  // Check if it has already been done:
  if (fDCS[iz]) return;
  // Do it otherwise
  if (fIsElectron) {
    // e-
    // load the high energy part firt:
    // - with gNumThetas2 theta and gNumEnergies-gIndxEnergyLim energy values
    const std::size_t hNumEnergries = gNumEnergies-gIndxEnergyLim;
    G4Physics2DVector* v2DHigh = new G4Physics2DVector(gNumThetas2, hNumEnergries);
    v2DHigh->SetBicubicInterpolation(true);
    for (std::size_t it=0; it<gNumThetas2; ++it) {
      v2DHigh->PutX(it, gTheMus2[it]);
    }
    for (std::size_t ie=0; ie<hNumEnergries; ++ie) {
      v2DHigh->PutY(ie, gTheEnergies[gIndxEnergyLim+ie]);
    }
    std::ostringstream ossh;
    ossh << FindDirectoryPath() << "dcss/el/dcs_"<< iz<<"_h";
    std::istringstream finh(std::ios::in);
    ReadCompressedFile(ossh.str(), finh);
    G4double dum = 0.0;
    for (std::size_t it=0; it<gNumThetas2; ++it) {
      finh >> dum;
      for (std::size_t ie=0; ie<hNumEnergries; ++ie) {
        finh >> dum;
        v2DHigh->PutValue(it, ie, G4Log(dum*CLHEP::cm2/CLHEP::sr));
      }
    }
    // load the low energy part:
    // - with gNumThetas1 theta and gIndxEnergyLim+1 energy values (the +1 is
    //   for including the firts DCS from the higher part above for being
    //   able to perform interpolation between the high and low energy DCS set)
    G4Physics2DVector* v2DLow = new G4Physics2DVector(gNumThetas1, gIndxEnergyLim+1);
    v2DLow->SetBicubicInterpolation(true);
    for (std::size_t it=0; it<gNumThetas1; ++it) {
      v2DLow->PutX(it, gTheMus1[it]);
    }
    for (std::size_t ie=0; ie<gIndxEnergyLim+1; ++ie) {
      v2DLow->PutY(ie, gTheEnergies[ie]);
    }
    std::ostringstream ossl;
    ossl << FindDirectoryPath() << "dcss/el/dcs_"<< iz<<"_l";
    std::istringstream finl(std::ios::in);
    ReadCompressedFile(ossl.str(), finl);
    for (std::size_t it=0; it<gNumThetas1; ++it) {
      finl >> dum;
      for (std::size_t ie=0; ie<gIndxEnergyLim; ++ie) {
        finl >> dum;
        v2DLow->PutValue(it, ie, G4Log(dum*CLHEP::cm2/CLHEP::sr));
      }
    }
    // add the +1 part: interpolate the firts DCS from the high energy
    std::size_t ix = 0;
    std::size_t iy = 0;
    for (std::size_t it=0; it<gNumThetas1; ++it) {
      const G4double val = v2DHigh->Value(gTheMus1[it], gTheEnergies[gIndxEnergyLim], ix, iy);
      v2DLow->PutValue(it, gIndxEnergyLim, val);
    }
    // store
    fDCSLow[iz] = v2DLow;
    fDCS[iz]    = v2DHigh;
  } else {
    // e+
    G4Physics2DVector* v2D= new G4Physics2DVector(gNumThetas2, gNumEnergies);
    v2D->SetBicubicInterpolation(true);
    for (std::size_t it=0; it<gNumThetas2; ++it) {
      v2D->PutX(it, gTheMus2[it]);
    }
    for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
      v2D->PutY(ie, gTheEnergies[ie]);
    }
    std::ostringstream oss;
    oss << FindDirectoryPath() << "dcss/pos/dcs_"<< iz;
    std::istringstream fin(std::ios::in);
    ReadCompressedFile(oss.str(), fin);
    G4double dum = 0.0;
    for (std::size_t it=0; it<gNumThetas2; ++it) {
      fin >> dum;
      for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
        fin >> dum;
        v2D->PutValue(it, ie, G4Log(dum*CLHEP::cm2/CLHEP::sr));
      }
    }
    fDCS[iz]= v2D;
  }
}

References CLHEP::cm2, fDCS, fDCSLow, FindDirectoryPath(), fIsElectron, G4Log(), gIndxEnergyLim, gNumEnergies, gNumThetas1, gNumThetas2, gTheEnergies, gTheMus1, gTheMus2, G4Physics2DVector::PutValue(), G4Physics2DVector::PutX(), G4Physics2DVector::PutY(), ReadCompressedFile(), G4Physics2DVector::SetBicubicInterpolation(), CLHEP::sr, and G4Physics2DVector::Value().

Referenced by InitialiseForZ().

LoadGrid()

void G4eDPWAElasticDCS::LoadGrid ( )

private

Definition at line 123 of file G4eDPWAElasticDCS.cc.

                                 {
  G4String fname = FindDirectoryPath() + "grid.dat";
  std::ifstream infile(fname.c_str());
  if (!infile.is_open()) {
    G4String msg  =
         "    Problem while trying to read " + fname + " file.\n"+
         "    G4LEDATA version should be G4EMLOW7.12 or later.\n";
    G4Exception("G4eDPWAElasticDCS::ReadCompressedFile","em0006",
                FatalException,msg.c_str());
    return;
  }
  // read size
  infile >> gNumEnergies;
  infile >> gNumThetas1;
  infile >> gNumThetas2;
  // read the grids
  // - energy in [MeV]
  G4double dum = 0.0;
  gTheEnergies.resize(gNumEnergies);
  for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
    infile >> dum;
    gTheEnergies[ie] = G4Log(dum*CLHEP::MeV);
    if (gTheEnergies[ie]<G4Log(2.0*CLHEP::keV)) gIndxEnergyLim = ie;  // only for e-
  }
  ++gIndxEnergyLim;
  // store/set usefull logarithms of the kinetic energy grid
  gLogMinEkin    = gTheEnergies[0];
  gInvDelLogEkin = (gNumEnergies-1)/(gTheEnergies[gNumEnergies-1]-gTheEnergies[0]);
  // - theta1 in [deg.] (247): we store mu(theta) = 0.5[1-cos(theta)]
  gTheMus1.resize(gNumThetas1);
  gTheU1.resize(gNumThetas1);
  const double theA = 0.01;
  for (std::size_t it=0; it<gNumThetas1; ++it) {
    infile >> dum;
    gTheMus1[it] = 0.5*(1.0-std::cos(dum*CLHEP::degree));
    gTheU1[it]   = (theA+1.0)*gTheMus1[it]/(theA+gTheMus1[it]);
  }
  // - theta2 in [deg.] (128): we store mu(theta) = 0.5[1-cos(theta)]
  gTheMus2.resize(gNumThetas2);
  gTheU2.resize(gNumThetas2);
  for (std::size_t it=0; it<gNumThetas2; ++it) {
    infile >> dum;
    gTheMus2[it] = 0.5*(1.0-std::cos(dum*CLHEP::degree));
    gTheU2[it]   = (theA+1.0)*gTheMus2[it]/(theA+gTheMus2[it]);
 
  }
  infile.close();
  gIsGridLoaded = true;
}

References CLHEP::degree, FatalException, FindDirectoryPath(), test::fname, G4Exception(), G4Log(), gIndxEnergyLim, gInvDelLogEkin, gIsGridLoaded, gLogMinEkin, gNumEnergies, gNumThetas1, gNumThetas2, gTheEnergies, gTheMus1, gTheMus2, gTheU1, gTheU2, CLHEP::keV, and CLHEP::MeV.

Referenced by InitialiseForZ().

ReadCompressedFile()

void G4eDPWAElasticDCS::ReadCompressedFile ( G4String  fname, std::istringstream &  iss  )

private

Definition at line 521 of file G4eDPWAElasticDCS.cc.

                                                                           {
  G4String *dataString = nullptr;
  G4String compfilename(fname+".z");
  // create input stream with binary mode operation and positioning at the end of the file
  std::ifstream in(compfilename, std::ios::binary | std::ios::ate);
  if (in.good()) {
     // get current position in the stream (was set to the end)
     G4int fileSize = in.tellg();
     // set current position being the beginning of the stream
     in.seekg(0,std::ios::beg);
     // create (zlib) byte buffer for the data
     Bytef *compdata = new Bytef[fileSize];
     while(in) {
        in.read((char*)compdata, fileSize);
     }
     // create (zlib) byte buffer for the uncompressed data
     uLongf complen    = (uLongf)(fileSize*4);
     Bytef *uncompdata = new Bytef[complen];
     while (Z_OK!=uncompress(uncompdata, &complen, compdata, fileSize)) {
        // increase uncompressed byte buffer
        delete[] uncompdata;
        complen   *= 2;
        uncompdata = new Bytef[complen];
     }
     // delete the compressed data buffer
     delete [] compdata;
     // create a string from the uncompressed data (will be deallocated by the caller)
     dataString = new G4String((char*)uncompdata, (long)complen);
     // delete the uncompressed data buffer
     delete [] uncompdata;
  } else {
    G4String msg  =
         "    Problem while trying to read " + fname + " data file.\n"+
         "    G4LEDATA version should be G4EMLOW7.12 or later.\n";
    G4Exception("G4eDPWAElasticDCS::ReadCompressedFile","em0006",
                FatalException,msg.c_str());
    return;
  }
  // create the input string stream from the data string
  if (dataString) {
    iss.str(*dataString);
    in.close();
    delete dataString;
  }
}

References FatalException, test::fname, G4Exception(), uncompress(), and Z_OK.

Referenced by BuildSmplingTableForZ(), and LoadDCSForZ().

SampleCosineTheta()

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

Definition at line 403 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 std::size_t     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));
}

References gInvDelLogEkin, gLogMinEkin, gNumEnergies, gTheEnergies, G4INCL::Math::max(), G4INCL::Math::min(), and SampleMu().

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 426 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 std::size_t     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));
}

References gInvDelLogEkin, gLogMinEkin, gNumEnergies, gTheEnergies, G4INCL::Math::max(), G4INCL::Math::min(), and SampleMu().

Referenced by G4eDPWACoulombScatteringModel::SampleSecondaries().

SampleMu() [1/2]

G4double G4eDPWAElasticDCS::SampleMu ( std::size_t  izet, std::size_t  ie, G4double  r1, G4double  muMin, G4double  muMax  )

private

Definition at line 477 of file G4eDPWAElasticDCS.cc.

                                                                     {
  const OneSamplingTable& rtn = (*fSamplingTables[izet])[ie];
  const G4double theA    = rtn.fScreenParA;
  //
  const std::vector<G4double>& theUVect = (fIsElectron && ie<gIndxEnergyLim) ? gTheU1 : gTheU2;
  const G4double xiMin   = (muMin > 0.0) ? FindCumValue((theA+1.0)*muMin/(theA+muMin), rtn, theUVect) : 0.0;
  const G4double xiMax   = (muMax < 1.0) ? FindCumValue((theA+1.0)*muMax/(theA+muMax), rtn, theUVect) : 1.0;
  //
  const G4double    xi   = xiMin+r1*(xiMax-xiMin); // a smaple within the range
  const std::size_t iLow = std::distance( rtn.fCum.begin(), std::upper_bound(rtn.fCum.begin(), rtn.fCum.end(), xi) )-1;
  const G4double   delta = rtn.fCum[iLow + 1] - rtn.fCum[iLow];
  const G4double   aval  = xi - rtn.fCum[iLow];
 
  const G4double dum1    = (1.0 + rtn.fA[iLow] + rtn.fB[iLow]) * delta * aval;
  const G4double dum2    = delta * delta + rtn.fA[iLow] * delta * aval + rtn.fB[iLow] * aval * aval;
  const G4double    u    = theUVect[iLow] + dum1 / dum2 * (theUVect[iLow + 1] - theUVect[iLow]);
  return theA*u/(theA+1.0-u);
}

References G4eDPWAElasticDCS::OneSamplingTable::fA, G4eDPWAElasticDCS::OneSamplingTable::fB, G4eDPWAElasticDCS::OneSamplingTable::fCum, FindCumValue(), fIsElectron, fSamplingTables, G4eDPWAElasticDCS::OneSamplingTable::fScreenParA, gIndxEnergyLim, gTheU1, and gTheU2.

SampleMu() [2/2]

G4double G4eDPWAElasticDCS::SampleMu ( std::size_t  izet, std::size_t  ie, G4double  r1, G4double  r2  )

private

Definition at line 444 of file G4eDPWAElasticDCS.cc.

                                                                                  {
  OneSamplingTable& rtn = (*fSamplingTables[izet])[ie];
  // get the lower index of the bin by using the alias part
  const G4double rest = r1 * (rtn.fN - 1);
  std::size_t indxl   = (std::size_t)rest;
  const G4double dum0 = rest - indxl;
  if (rtn.fW[indxl] < dum0) indxl = rtn.fI[indxl];
  // sample value within the selected bin by using ratin based numerical inversion
  const G4double delta = rtn.fCum[indxl + 1] - rtn.fCum[indxl];
  const G4double aval  = r2 * delta;
 
  const G4double dum1  = (1.0 + rtn.fA[indxl] + rtn.fB[indxl]) * delta * aval;
  const G4double dum2  = delta * delta + rtn.fA[indxl] * delta * aval + rtn.fB[indxl] * aval * aval;
  const std::vector<G4double>& theUVect = (fIsElectron && ie<gIndxEnergyLim) ? gTheU1 : gTheU2;
  const G4double    u  = theUVect[indxl] + dum1 / dum2 * (theUVect[indxl + 1] - theUVect[indxl]);
  // transform back u to mu
  return rtn.fScreenParA*u/(rtn.fScreenParA+1.0-u);
}

References G4eDPWAElasticDCS::OneSamplingTable::fA, G4eDPWAElasticDCS::OneSamplingTable::fB, G4eDPWAElasticDCS::OneSamplingTable::fCum, G4eDPWAElasticDCS::OneSamplingTable::fI, fIsElectron, G4eDPWAElasticDCS::OneSamplingTable::fN, fSamplingTables, G4eDPWAElasticDCS::OneSamplingTable::fScreenParA, G4eDPWAElasticDCS::OneSamplingTable::fW, gIndxEnergyLim, gTheU1, and gTheU2.

Referenced by SampleCosineTheta(), and SampleCosineThetaRestricted().

Field Documentation

fDCS

std::vector<G4Physics2DVector*> G4eDPWAElasticDCS::fDCS

private

Definition at line 252 of file G4eDPWAElasticDCS.hh.

Referenced by ComputeCSPerAtom(), G4eDPWAElasticDCS(), LoadDCSForZ(), and ~G4eDPWAElasticDCS().

fDCSLow

std::vector<G4Physics2DVector*> G4eDPWAElasticDCS::fDCSLow

private

Definition at line 253 of file G4eDPWAElasticDCS.hh.

Referenced by ComputeCSPerAtom(), G4eDPWAElasticDCS(), LoadDCSForZ(), and ~G4eDPWAElasticDCS().

fIsElectron

G4bool G4eDPWAElasticDCS::fIsElectron

private

Definition at line 228 of file G4eDPWAElasticDCS.hh.

Referenced by BuildSmplingTableForZ(), ComputeCSPerAtom(), InitSCPCorrection(), LoadDCSForZ(), and SampleMu().

fIsRestrictedSamplingRequired

G4bool G4eDPWAElasticDCS::fIsRestrictedSamplingRequired

private

Definition at line 226 of file G4eDPWAElasticDCS.hh.

Referenced by BuildSmplingTableForZ().

fNumSPCEbinPerDec

const G4int G4eDPWAElasticDCS::fNumSPCEbinPerDec = 3

private

Definition at line 258 of file G4eDPWAElasticDCS.hh.

Referenced by InitSCPCorrection().

fSamplingTables

std::vector< std::vector<OneSamplingTable>* > G4eDPWAElasticDCS::fSamplingTables

private

Definition at line 255 of file G4eDPWAElasticDCS.hh.

Referenced by BuildSmplingTableForZ(), G4eDPWAElasticDCS(), SampleMu(), and ~G4eDPWAElasticDCS().

fSCPCPerMatCuts

std::vector<SCPCorrection*> G4eDPWAElasticDCS::fSCPCPerMatCuts

private

Definition at line 266 of file G4eDPWAElasticDCS.hh.

Referenced by ComputeScatteringPowerCorrection(), InitSCPCorrection(), and ~G4eDPWAElasticDCS().

gDataDirectory

G4String G4eDPWAElasticDCS::gDataDirectory = ""

staticprivate

Definition at line 232 of file G4eDPWAElasticDCS.hh.

Referenced by FindDirectoryPath().

gIndxEnergyLim

std::size_t G4eDPWAElasticDCS::gIndxEnergyLim = 35

staticprivate

Definition at line 237 of file G4eDPWAElasticDCS.hh.

Referenced by ComputeCSPerAtom(), LoadDCSForZ(), LoadGrid(), and SampleMu().

gInvDelLogEkin

G4double G4eDPWAElasticDCS::gInvDelLogEkin = 1.0

staticprivate

Definition at line 246 of file G4eDPWAElasticDCS.hh.

Referenced by LoadGrid(), SampleCosineTheta(), and SampleCosineThetaRestricted().

gIsGridLoaded

G4bool G4eDPWAElasticDCS::gIsGridLoaded = false

staticprivate

Definition at line 230 of file G4eDPWAElasticDCS.hh.

Referenced by InitialiseForZ(), and LoadGrid().

gLogMinEkin

G4double G4eDPWAElasticDCS::gLogMinEkin = 1.0

staticprivate

Definition at line 245 of file G4eDPWAElasticDCS.hh.

Referenced by LoadGrid(), SampleCosineTheta(), and SampleCosineThetaRestricted().

gMaxZ

constexpr std::size_t G4eDPWAElasticDCS::gMaxZ = 103

staticconstexprprivate

Definition at line 234 of file G4eDPWAElasticDCS.hh.

Referenced by G4eDPWAElasticDCS().

gNumEnergies

std::size_t G4eDPWAElasticDCS::gNumEnergies = 106

staticprivate

Definition at line 236 of file G4eDPWAElasticDCS.hh.

Referenced by BuildSmplingTableForZ(), ComputeCSPerAtom(), LoadDCSForZ(), LoadGrid(), SampleCosineTheta(), and SampleCosineThetaRestricted().

gNumThetas1

std::size_t G4eDPWAElasticDCS::gNumThetas1 = 247

staticprivate

Definition at line 238 of file G4eDPWAElasticDCS.hh.

Referenced by LoadDCSForZ(), and LoadGrid().

gNumThetas2

std::size_t G4eDPWAElasticDCS::gNumThetas2 = 128

staticprivate

Definition at line 239 of file G4eDPWAElasticDCS.hh.

Referenced by LoadDCSForZ(), and LoadGrid().

gTheEnergies

std::vector< G4double > G4eDPWAElasticDCS::gTheEnergies

staticprivate

Definition at line 240 of file G4eDPWAElasticDCS.hh.

Referenced by ComputeCSPerAtom(), LoadDCSForZ(), LoadGrid(), SampleCosineTheta(), and SampleCosineThetaRestricted().

gTheMus1

std::vector< G4double > G4eDPWAElasticDCS::gTheMus1

staticprivate

Definition at line 241 of file G4eDPWAElasticDCS.hh.

Referenced by ComputeCSPerAtom(), LoadDCSForZ(), and LoadGrid().

gTheMus2

std::vector< G4double > G4eDPWAElasticDCS::gTheMus2

staticprivate

Definition at line 242 of file G4eDPWAElasticDCS.hh.

Referenced by ComputeCSPerAtom(), LoadDCSForZ(), and LoadGrid().

gTheU1

std::vector< G4double > G4eDPWAElasticDCS::gTheU1

staticprivate

Definition at line 243 of file G4eDPWAElasticDCS.hh.

Referenced by LoadGrid(), and SampleMu().

gTheU2

std::vector< G4double > G4eDPWAElasticDCS::gTheU2

staticprivate

Definition at line 244 of file G4eDPWAElasticDCS.hh.

Referenced by LoadGrid(), and SampleMu().

gWGL

const G4double G4eDPWAElasticDCS::gWGL

staticprivate

Initial value:

= {
  5.06142681E-02, 1.11190517E-01, 1.56853323E-01, 1.81341892E-01,
  1.81341892E-01, 1.56853323E-01, 1.11190517E-01, 5.06142681E-02
}

Definition at line 250 of file G4eDPWAElasticDCS.hh.

Referenced by ComputeCSPerAtom().

gXGL

const G4double G4eDPWAElasticDCS::gXGL

staticprivate

Initial value:

= {
  1.98550718E-02, 1.01666761E-01, 2.37233795E-01, 4.08282679E-01,
  5.91717321E-01, 7.62766205E-01, 8.98333239E-01, 9.80144928E-01
}

Definition at line 249 of file G4eDPWAElasticDCS.hh.

Referenced by ComputeCSPerAtom().

---

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

• source/processes/electromagnetic/standard/include/G4eDPWAElasticDCS.hh
• source/processes/electromagnetic/standard/src/G4eDPWAElasticDCS.cc