#include <G4GSMottCorrection.hh>

Data Structures
struct	DataPerDelta

struct	DataPerEkin

struct	DataPerMaterial

Public Member Functions
	G4GSMottCorrection (G4bool iselectron=true)

void	GetMottCorrectionFactors (G4double logekin, G4double beta2, G4int matindx, G4double &mcToScr, G4double &mcToQ1, G4double &mcToG2PerG1)

G4double	GetMottRejectionValue (G4double logekin, G4double G4beta2, G4double q1, G4double cost, G4int matindx, G4int &ekindx, G4int &deltindx)

void	Initialise ()

	~G4GSMottCorrection ()

Static Public Member Functions
static G4int	GetMaxZet ()

Private Member Functions
void	AllocateDataPerMaterial (DataPerMaterial *)

void	ClearMCDataPerElement ()

void	ClearMCDataPerMaterial ()

void	DeAllocateDataPerMaterial (DataPerMaterial *)

void	InitMCDataMaterial (const G4Material *)

void	InitMCDataPerElement ()

void	InitMCDataPerMaterials ()

void	LoadMCDataElement (const G4Element *)

void	ReadCompressedFile (std::string fname, std::istringstream &iss)

Private Attributes
G4double	fInvDelAngle

G4double	fInvDelBeta2

G4double	fInvDelDelta

G4double	fInvLogDelEkin

G4bool	fIsElectron

G4double	fLogMinEkin

G4double	fMaxEkin

std::vector< DataPerMaterial * >	fMCDataPerElement

std::vector< DataPerMaterial * >	fMCDataPerMaterial

G4double	fMinBeta2

Static Private Attributes
static const std::string	gElemSymbols []

static constexpr G4double	gMaxBeta2 = 0.9999

static constexpr G4double	gMaxDelta = 0.9

static constexpr G4int	gMaxZet = 98

static constexpr G4double	gMidEkin = 100.*CLHEP::keV

static constexpr G4double	gMinEkin = 1.*CLHEP::keV

static constexpr G4int	gNumAngle = 32

static constexpr G4int	gNumBeta2 = 16

static constexpr G4int	gNumDelta = 28

static constexpr G4int	gNumEkin = 31

Detailed Description

Definition at line 87 of file G4GSMottCorrection.hh.

Constructor & Destructor Documentation

◆ G4GSMottCorrection()

G4GSMottCorrection::G4GSMottCorrection ( G4bool iselectron = true )

Definition at line 71 of file G4GSMottCorrection.cc.

                                                        : fIsElectron(iselectron) {
  // init grids related data member values
  fMaxEkin        = CLHEP::electron_mass_c2*(1./std::sqrt(1.-gMaxBeta2)-1.);
  fLogMinEkin     = G4Log(gMinEkin);
  fInvLogDelEkin  = (gNumEkin-gNumBeta2)/G4Log(gMidEkin/gMinEkin);
  G4double pt2    = gMidEkin*(gMidEkin+2.0*CLHEP::electron_mass_c2);
  fMinBeta2       = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
  fInvDelBeta2    = (gNumBeta2-1.)/(gMaxBeta2-fMinBeta2);
  fInvDelDelta    = (gNumDelta-1.)/gMaxDelta;
  fInvDelAngle    = gNumAngle-1.;
}

References CLHEP::electron_mass_c2, fInvDelAngle, fInvDelBeta2, fInvDelDelta, fInvLogDelEkin, fLogMinEkin, fMaxEkin, fMinBeta2, G4Log(), gMaxBeta2, gMaxDelta, gMidEkin, gMinEkin, gNumAngle, gNumBeta2, gNumDelta, and gNumEkin.

◆ ~G4GSMottCorrection()

G4GSMottCorrection::~G4GSMottCorrection ( )

Definition at line 84 of file G4GSMottCorrection.cc.

                                        {
  ClearMCDataPerElement();
  ClearMCDataPerMaterial();
}

References ClearMCDataPerElement(), and ClearMCDataPerMaterial().

Member Function Documentation

◆ AllocateDataPerMaterial()

void G4GSMottCorrection::AllocateDataPerMaterial ( DataPerMaterial * data )

private

Definition at line 506 of file G4GSMottCorrection.cc.

                                                                      {
  data->fDataPerEkin = new DataPerEkin*[gNumEkin]();
  for (G4int iek=0; iek<gNumEkin; ++iek) {
    DataPerEkin *perEkin   = new DataPerEkin();
    perEkin->fDataPerDelta = new DataPerDelta*[gNumDelta]();
    for (G4int idel=0; idel<gNumDelta; ++idel) {
      DataPerDelta *perDelta       = new DataPerDelta();
      perDelta->fRejFuntion        = new double[gNumAngle]();
      perEkin->fDataPerDelta[idel] = perDelta;
    }
    data->fDataPerEkin[iek] = perEkin;
  }
}

References G4GSMottCorrection::DataPerEkin::fDataPerDelta, G4GSMottCorrection::DataPerMaterial::fDataPerEkin, G4GSMottCorrection::DataPerDelta::fRejFuntion, gNumAngle, gNumDelta, and gNumEkin.

Referenced by InitMCDataMaterial(), and LoadMCDataElement().

◆ ClearMCDataPerElement()

void G4GSMottCorrection::ClearMCDataPerElement ( )

private

Definition at line 535 of file G4GSMottCorrection.cc.

                                               {
  for (size_t i=0; i<fMCDataPerElement.size(); ++i) {
    if (fMCDataPerElement[i]) {
      DeAllocateDataPerMaterial(fMCDataPerElement[i]);
      delete fMCDataPerElement[i];
    }
  }
  fMCDataPerElement.clear();
}

References DeAllocateDataPerMaterial(), and fMCDataPerElement.

Referenced by ~G4GSMottCorrection().

◆ ClearMCDataPerMaterial()

void G4GSMottCorrection::ClearMCDataPerMaterial ( )

private

Definition at line 545 of file G4GSMottCorrection.cc.

                                                {
  for (size_t i=0; i<fMCDataPerMaterial.size(); ++i) {
    if (fMCDataPerMaterial[i]) {
      DeAllocateDataPerMaterial(fMCDataPerMaterial[i]);
      delete fMCDataPerMaterial[i];
    }
  }
  fMCDataPerMaterial.clear();
}

References DeAllocateDataPerMaterial(), and fMCDataPerMaterial.

Referenced by Initialise(), and ~G4GSMottCorrection().

◆ DeAllocateDataPerMaterial()

void G4GSMottCorrection::DeAllocateDataPerMaterial ( DataPerMaterial * data )

private

Definition at line 520 of file G4GSMottCorrection.cc.

                                                                        {
  for (G4int iek=0; iek<gNumEkin; ++iek) {
    DataPerEkin *perEkin = data->fDataPerEkin[iek]; //new DataPerEkin();
    for (G4int idel=0; idel<gNumDelta; ++idel) {
      DataPerDelta *perDelta = perEkin->fDataPerDelta[idel];
      delete [] perDelta->fRejFuntion;
      delete perDelta;
    }
    delete [] perEkin->fDataPerDelta;
    delete perEkin;
  }
  delete [] data->fDataPerEkin;
}

References G4GSMottCorrection::DataPerEkin::fDataPerDelta, G4GSMottCorrection::DataPerMaterial::fDataPerEkin, G4GSMottCorrection::DataPerDelta::fRejFuntion, gNumDelta, and gNumEkin.

Referenced by ClearMCDataPerElement(), and ClearMCDataPerMaterial().

◆ GetMaxZet()

static G4int G4GSMottCorrection::GetMaxZet ( )

inlinestatic

Definition at line 101 of file G4GSMottCorrection.hh.

101{ return gMaxZet; }

G4GSMottCorrection::gMaxZet

static constexpr G4int gMaxZet

Definition: G4GSMottCorrection.hh:164

References gMaxZet.

Referenced by G4GoudsmitSaundersonTable::InitMoliereMSCParams().

◆ GetMottCorrectionFactors()

void G4GSMottCorrection::GetMottCorrectionFactors	(	G4double	logekin,
		G4double	beta2,
		G4int	matindx,
		G4double &	mcToScr,
		G4double &	mcToQ1,
		G4double &	mcToG2PerG1
	)

Definition at line 90 of file G4GSMottCorrection.cc.

                                                                                           {
  G4int    ekinIndxLow = 0;
  G4double remRfaction = 0.;
  if (beta2>=gMaxBeta2) {
    ekinIndxLow = gNumEkin - 1;
    // remRfaction = -1.
  } else if (beta2>=fMinBeta2) {  // linear interpolation on \beta^2
    remRfaction   = (beta2 - fMinBeta2) * fInvDelBeta2;
    ekinIndxLow   = (G4int)remRfaction;
    remRfaction  -= ekinIndxLow;
    ekinIndxLow  += (gNumEkin - gNumBeta2);
  } else if (logekin>=fLogMinEkin) {
    remRfaction   = (logekin - fLogMinEkin) * fInvLogDelEkin;
    ekinIndxLow   = (G4int)remRfaction;
    remRfaction  -= ekinIndxLow;
  } // the defaults otherwise i.e. use the lowest energy values when ekin is smaller than the minum ekin
  //
  DataPerEkin *perEkinLow  = fMCDataPerMaterial[matindx]->fDataPerEkin[ekinIndxLow];
  mcToScr      = perEkinLow->fMCScreening;
  mcToQ1       = perEkinLow->fMCFirstMoment;
  mcToG2PerG1  = perEkinLow->fMCSecondMoment;
  if (remRfaction>0.) {
    DataPerEkin *perEkinHigh = fMCDataPerMaterial[matindx]->fDataPerEkin[ekinIndxLow+1];
    mcToScr      += remRfaction*(perEkinHigh->fMCScreening    - perEkinLow->fMCScreening);
    mcToQ1       += remRfaction*(perEkinHigh->fMCFirstMoment  - perEkinLow->fMCFirstMoment);
    mcToG2PerG1  += remRfaction*(perEkinHigh->fMCSecondMoment - perEkinLow->fMCSecondMoment);
  }
}

References fInvDelBeta2, fInvLogDelEkin, fLogMinEkin, fMCDataPerMaterial, G4GSMottCorrection::DataPerEkin::fMCFirstMoment, G4GSMottCorrection::DataPerEkin::fMCScreening, G4GSMottCorrection::DataPerEkin::fMCSecondMoment, fMinBeta2, gMaxBeta2, gNumBeta2, and gNumEkin.

◆ GetMottRejectionValue()

double G4GSMottCorrection::GetMottRejectionValue	(	G4double	logekin,
		G4double	G4beta2,
		G4double	q1,
		G4double	cost,
		G4int	matindx,
		G4int &	ekindx,
		G4int &	deltindx
	)

Definition at line 122 of file G4GSMottCorrection.cc.

                                                                                                {
  G4double val   = 1.0;
  G4double delta = q1/(0.5+q1);
  // check if converged to 1 for all angles => accept cost
  if (delta>=gMaxDelta) {
    return val;
  }
  //
  // check if kinetic energy index needs to be determined
  if (ekindx<0) {
    G4int    ekinIndxLow  = 0;
    G4double probIndxHigh = 0.;  // will be the prob. of taking the ekinIndxLow+1 bin
    if (beta2>gMaxBeta2) {
      ekinIndxLow = gNumEkin - 1;
      // probIndxHigh = -1.
    } else if (beta2>=fMinBeta2) {    // linear interpolation on \beta^2
      probIndxHigh  = (beta2 - fMinBeta2) * fInvDelBeta2;
      ekinIndxLow   = (G4int)probIndxHigh;
      probIndxHigh -= ekinIndxLow;
      ekinIndxLow  += (gNumEkin - gNumBeta2);
    } else if (logekin>fLogMinEkin) { // linear interpolation on \ln(E_{kin})
      probIndxHigh  = (logekin - fLogMinEkin) * fInvLogDelEkin;
      ekinIndxLow   = (G4int)probIndxHigh;
      probIndxHigh -= ekinIndxLow;
    } // the defaults otherwise i.e. use the lowest energy values when ekin is smaller than the minum ekin
    //
    // check if need to take the higher ekin index
    if (G4UniformRand()<probIndxHigh) {
      ++ekinIndxLow;
    }
    // set kinetic energy grid index
    ekindx = ekinIndxLow;
  }
  // check if delta value index needs to be determined (note: in case of single scattering deltindx will be set to 0 by
  // by the caller but the ekindx will be -1: kinetic energy index is not known but the delta index is known)
  if (deltindx<0) {
    // note: delta is for sure < gMaxDelta at this point ( and minimum delta value is 0)
    G4double probIndxHigh = delta*fInvDelDelta;  // will be the prob. of taking the deltIndxLow+1 bin
    G4int    deltIndxLow  = (G4int)probIndxHigh;
    probIndxHigh         -= deltIndxLow;
    // check if need to take the higher delta index
    if (G4UniformRand()<probIndxHigh) {
      ++deltIndxLow;
    }
    // set the delta value grid index
    deltindx = deltIndxLow;
  }
  //
  // get the corresponding distribution
  DataPerDelta *perDelta  = fMCDataPerMaterial[matindx]->fDataPerEkin[ekindx]->fDataPerDelta[deltindx];
  //
  // determine lower index of the angular bin
  G4double ang         = std::sqrt(0.5*(1.-cost)); // sin(0.5\theta) in [0,1]
  G4double remRfaction = ang*fInvDelAngle;
  G4int    angIndx     = (G4int)remRfaction;
  remRfaction         -= angIndx;
  if (angIndx<gNumAngle-2) { // normal case: linear interpolation
    val          = remRfaction*(perDelta->fRejFuntion[angIndx+1]-perDelta->fRejFuntion[angIndx]) + perDelta->fRejFuntion[angIndx];
  } else {   // last bin
    G4double dum = ang-1.+1./fInvDelAngle;
    val          = perDelta->fSA + dum*(perDelta->fSB + dum*(perDelta->fSC + dum*perDelta->fSD));
  }
  return val;
}

References fInvDelAngle, fInvDelBeta2, fInvDelDelta, fInvLogDelEkin, fLogMinEkin, fMCDataPerMaterial, fMinBeta2, G4GSMottCorrection::DataPerDelta::fRejFuntion, G4GSMottCorrection::DataPerDelta::fSA, G4GSMottCorrection::DataPerDelta::fSB, G4GSMottCorrection::DataPerDelta::fSC, G4GSMottCorrection::DataPerDelta::fSD, G4UniformRand, gMaxBeta2, gMaxDelta, gNumAngle, gNumBeta2, and gNumEkin.

◆ Initialise()

void G4GSMottCorrection::Initialise ( )

Definition at line 189 of file G4GSMottCorrection.cc.

                                    {
  // load Mott-correction data for each elements that belongs to materials that are used in the detector
  InitMCDataPerElement();
  // clrea Mott-correction data per material
  ClearMCDataPerMaterial();
  // initialise Mott-correction data for the materials that are used in the detector
  InitMCDataPerMaterials();
}

References ClearMCDataPerMaterial(), InitMCDataPerElement(), and InitMCDataPerMaterials().

◆ InitMCDataMaterial()

void G4GSMottCorrection::InitMCDataMaterial ( const G4Material * mat )

private

Definition at line 349 of file G4GSMottCorrection.cc.

                                                                 {
  constexpr G4double const1   = 7821.6;      // [cm2/g]
  constexpr G4double const2   = 0.1569;      // [cm2 MeV2 / g]
  constexpr G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
 
  G4double constFactor        = CLHEP::electron_mass_c2*CLHEP::fine_structure_const/0.88534;
  constFactor                *= constFactor;  // (mc^2)^2\alpha^2/( C_{TF}^2)
  // allocate memory
  DataPerMaterial *perMat     = new DataPerMaterial();
  AllocateDataPerMaterial(perMat);
  fMCDataPerMaterial[mat->GetIndex()] = perMat;
  //
  const G4ElementVector* elemVect           = mat->GetElementVector();
  const G4int            numElems           = mat->GetNumberOfElements();
  const G4double*        nbAtomsPerVolVect  = mat->GetVecNbOfAtomsPerVolume();
  G4double               totNbAtomsPerVol   = mat->GetTotNbOfAtomsPerVolume();
  //
  // 1. Compute material dependent part of Moliere's b_c \chi_c^2
  //    (with \xi=1 (i.e. total sub-threshold scattering power correction)
  G4double moliereBc  = 0.0;
  G4double moliereXc2 = 0.0;
  G4double zs         = 0.0;
  G4double ze         = 0.0;
  G4double zx         = 0.0;
  G4double sa         = 0.0;
  G4double xi         = 1.0;
  for (G4int ielem=0; ielem<numElems; ++ielem) {
    G4double zet = (*elemVect)[ielem]->GetZ();
    if (zet>gMaxZet) {
      zet = (G4double)gMaxZet;
    }
    G4double iwa  = (*elemVect)[ielem]->GetN();
    G4double ipz  = nbAtomsPerVolVect[ielem]/totNbAtomsPerVol;
    G4double dum  = ipz*zet*(zet+xi);
    zs           += dum;
    ze           += dum*(-2.0/3.0)*G4Log(zet);
    zx           += dum*G4Log(1.0+3.34*finstrc2*zet*zet);
    sa           += ipz*iwa;
  }
  G4double density = mat->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
  //
  moliereBc  = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs);  //[1/cm]
  moliereXc2 = const2*density*zs/sa;  // [MeV2/cm]
  // change to Geant4 internal units of 1/length and energ2/length
  moliereBc  *= 1.0/CLHEP::cm;
  moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
  //
  // 2. loop over the kinetic energy grid
  for (G4int iek=0; iek<gNumEkin; ++iek) {
    // 2./a. set current kinetic energy and pt2 value
      G4double ekin          = G4Exp(fLogMinEkin+iek/fInvLogDelEkin);
      G4double pt2  = ekin*(ekin+2.0*CLHEP::electron_mass_c2);
      if (ekin>gMidEkin) {
        G4double b2   = fMinBeta2+(iek-(gNumEkin-gNumBeta2))/fInvDelBeta2;
        ekin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-b2)-1.);
        pt2  = ekin*(ekin+2.0*CLHEP::electron_mass_c2);
      }
    // 2./b. loop over the elements at the current kinetic energy point
    for (G4int ielem=0; ielem<numElems; ++ielem) {
      const G4Element *elem = (*elemVect)[ielem];
      G4double zet  = elem->GetZ();
      if (zet>gMaxZet) {
        zet = (G4double)gMaxZet;
      }
      G4int izet         = G4lrint(zet);
      // xi should be one i.e. z(z+1) since total sub-threshold scattering power correction
      G4double nZZPlus1  = nbAtomsPerVolVect[ielem]*zet*(zet+1.0)/totNbAtomsPerVol;
      G4double Z23       = std::pow(zet,2./3.);
      //
      DataPerEkin *perElemPerEkin  = fMCDataPerElement[izet]->fDataPerEkin[iek];
      DataPerEkin *perMatPerEkin   = perMat->fDataPerEkin[iek];
      //
      // 2./b./(i) Add the 3 Mott-correction factors
      G4double mcScrCF = perElemPerEkin->fMCScreening;     // \kappa_i[1.13+3.76(\alpha Z_i)^2] with \kappa_i=scr_mc/scr_sr
      // compute the screening parameter correction factor (Z_i contribution to the material)
      // src_{mc} = C \exp\left[ \frac{ \sum_i n_i Z_i(Z_i+1)\ln[Z_{i}^{2/3}\kappa_i(1.13+3.76(\alpha Z_i)^2)] } {\sum_i n_i Z_i(Z_i+1)}
      // with C = \frac{(mc^2)^\alpha^2} {4(pc)^2 C_{TF}^2} = constFactor/(4*(pc)^2)
      // here we compute the \sum_i n_i Z_i(Z_i+1)\ln[Z_{i}^{2/3}\kappa_i(1.13+3.76(\alpha Z_i)^2)] part
      perMatPerEkin->fMCScreening += nZZPlus1*G4Log(Z23*mcScrCF);
      // compute the corrected screening parameter for the current Z_i and E_{kin}
      // src(Z_i)_{mc} = \frac{(mc^2)^\alpha^2 Z_i^{2/3}} {4(pc)^2 C_{TF}^2} \kappa_i[1.13+3.76(\alpha Z_i)^2]
      mcScrCF *= constFactor*Z23/(4.*pt2);
      // compute first moment correction factor
      // q1_{mc} = \frac{ \sum_i n_i Z_i(Z_i+1) A_i  B_i } {\sum_i n_i Z_i(Z_i+1)} \frac{1}{C}
      // where:
      // A_i(src(Z_i)_{mc}) = [\ln(1+1/src(Z_i)_{mc}) - 1/(1+src(Z_i)_{mc})]; where \sigma(Z_i)_{tr1}^(sr) = A_i(src(Z_i)_{mc}) [2\pi r_0 Z_i mc^2/(pc)\beta]^2
      // B_i = \beta_i \gamma_i with beta_i(Z_i) = \sigma(Z_i)_{tr1}^(PWA)/\sigma(Z_i,src(Z_i)_{mc})_{tr1}^(sr)
      // and \gamma_i = \sigma(Z_i)_{el}^(MC-DCS)/\sigma(Z_i,src(Z_i)_{mc})_{el}^(sr)
      // C(src_{mc}) = [\ln(1+1/src_{mc}) - 1/(1+src_{mc})]; where \sigma_{tr1}^(sr) = C(src_{mc}) [2\pi r_0 Z_i mc^2/(pc)\beta]^2
      // A_i x B_i is stored in file per e-/e+, E_{kin} and Z_i
      // here we compute the \sum_i n_i Z_i(Z_i+1) A_i  B_i part
      perMatPerEkin->fMCFirstMoment += nZZPlus1*(G4Log(1.+1./mcScrCF)-1./(1.+mcScrCF))*perElemPerEkin->fMCFirstMoment;
      // compute the second moment correction factor
      // [G2/G1]_{mc} = \frac{ \sum_i n_i Z_i(Z_i+1) A_i } {\sum_i n_i Z_i(Z_i+1)} \frac{1}{C}
      // with A_i(Z_i) = G2(Z_i)^{PWA}/G1(Z_i)^{PWA} and C=G2(Z_i,scr_{mc})^{sr}/G1(Z_i,scr_{mc})^{sr}}
      // here we compute the \sum_i n_i Z_i(Z_i+1) A_i part
      perMatPerEkin->fMCSecondMoment += nZZPlus1*perElemPerEkin->fMCSecondMoment;
      //
      // 2./b./(ii) Go for the rejection funtion part
      // I. loop over delta values
      for (G4int idel=0; idel<gNumDelta; ++idel) {
        DataPerDelta *perMatPerDelta  = perMatPerEkin->fDataPerDelta[idel];
        DataPerDelta *perElemPerDelta = perElemPerEkin->fDataPerDelta[idel];
        // I./a. loop over angles (i.e. the \sin(0.5\theta) values) and add the rejection function
        for (G4int iang=0; iang<gNumAngle; ++iang) {
          perMatPerDelta->fRejFuntion[iang] += nZZPlus1*perElemPerDelta->fRejFuntion[iang];
        }
        // I./b. get the last bin spline parameters and add them (a+bx+cx^2+dx^3)
        perMatPerDelta->fSA += nZZPlus1*perElemPerDelta->fSA;
        perMatPerDelta->fSB += nZZPlus1*perElemPerDelta->fSB;
        perMatPerDelta->fSC += nZZPlus1*perElemPerDelta->fSC;
        perMatPerDelta->fSD += nZZPlus1*perElemPerDelta->fSD;
      }
      //
      // 2./b./(iii) When the last element has been added:
      if (ielem==numElems-1) {
        //
        // 1. the remaining part of the sreening correction and divide the corrected screening par. with Moliere's one:
        //    (Moliere screening parameter = moliereXc2/(4(pc)^2 moliereBc) )
        G4double dumScr   = G4Exp(perMatPerEkin->fMCScreening/zs);
        perMatPerEkin->fMCScreening = constFactor*dumScr*moliereBc/moliereXc2;
        //
        // 2. the remaining part of the first moment correction and divide by the one computed by using the corrected
        //    screening parameter (= (mc^2)^\alpha^2/(4(pc)^2C_{TF}^2) dumScr
        G4double scrCorTed = constFactor*dumScr/(4.*pt2);
        G4double dum0      = G4Log(1.+1./scrCorTed);
        perMatPerEkin->fMCFirstMoment = perMatPerEkin->fMCFirstMoment/(zs*(dum0-1./(1.+scrCorTed)));
        //
        // 3. the remaining part of the second moment correction and divide by the one computed by using the corrected
        //    screening parameter
        G4double G2PerG1   =  3.*(1.+scrCorTed)*((1.+2.*scrCorTed)*dum0-2.)/((1.+scrCorTed)*dum0-1.);
        perMatPerEkin->fMCSecondMoment = perMatPerEkin->fMCSecondMoment/(zs*G2PerG1);
        //
        // 4. scale the maximum of the rejection function to unity and correct the last bin spline parameters as well
        // I. loop over delta values
        for (G4int idel=0; idel<gNumDelta; ++idel) {
          DataPerDelta *perMatPerDelta  = perMatPerEkin->fDataPerDelta[idel];
          G4double maxVal = -1.;
          // II. llop over angles
          for (G4int iang=0; iang<gNumAngle; ++iang) {
            if (perMatPerDelta->fRejFuntion[iang]>maxVal)
              maxVal = perMatPerDelta->fRejFuntion[iang];
          }
          for (G4int iang=0; iang<gNumAngle; ++iang) {
            perMatPerDelta->fRejFuntion[iang] /=maxVal;
          }
          perMatPerDelta->fSA /= maxVal;
          perMatPerDelta->fSB /= maxVal;
          perMatPerDelta->fSC /= maxVal;
          perMatPerDelta->fSD /= maxVal;
        }
      }
    }
  }
}

Referenced by InitMCDataPerMaterials().

◆ InitMCDataPerElement()

void G4GSMottCorrection::InitMCDataPerElement ( )

private

Definition at line 199 of file G4GSMottCorrection.cc.

                                              {
  // do it only once
  if (fMCDataPerElement.size()<gMaxZet+1) {
    fMCDataPerElement.resize(gMaxZet+1,nullptr);
  }
  // loop over all materials, for those that are used check the list of elements and load data from file if the
  // corresponding data has not been loaded yet
  G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
  size_t numMatCuts = thePCTable->GetTableSize();
  for (size_t imc=0; imc<numMatCuts; ++imc) {
    const G4MaterialCutsCouple *matCut =  thePCTable->GetMaterialCutsCouple(imc);
    if (!matCut->IsUsed()) {
      continue;
    }
    const G4Material      *mat      = matCut->GetMaterial();
    const G4ElementVector *elemVect = mat->GetElementVector();
    //
    size_t numElems = elemVect->size();
    for (size_t ielem=0; ielem<numElems; ++ielem) {
      const G4Element *elem = (*elemVect)[ielem];
      G4int izet = G4lrint(elem->GetZ());
      if (izet>gMaxZet) {
        izet = gMaxZet;
      }
      if (!fMCDataPerElement[izet]) {
        LoadMCDataElement(elem);
      }
    }
  }
}

References elem, fMCDataPerElement, G4lrint(), G4Material::GetElementVector(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), gMaxZet, G4MaterialCutsCouple::IsUsed(), and LoadMCDataElement().

Referenced by Initialise().

◆ InitMCDataPerMaterials()

void G4GSMottCorrection::InitMCDataPerMaterials ( )

private

Definition at line 231 of file G4GSMottCorrection.cc.

                                                {
  // prepare size of the container
  size_t numMaterials = G4Material::GetNumberOfMaterials();
  if (fMCDataPerMaterial.size()!=numMaterials) {
    fMCDataPerMaterial.resize(numMaterials);
  }
  // init. Mott-correction data for the Materials that are used in the geometry
  G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
  size_t numMatCuts = thePCTable->GetTableSize();
  for (size_t imc=0; imc<numMatCuts; ++imc) {
    const G4MaterialCutsCouple *matCut =  thePCTable->GetMaterialCutsCouple(imc);
    if (!matCut->IsUsed()) {
      continue;
    }
    const G4Material *mat = matCut->GetMaterial();
    if (!fMCDataPerMaterial[mat->GetIndex()]) {
      InitMCDataMaterial(mat);
    }
  }
}

References fMCDataPerMaterial, G4Material::GetIndex(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4Material::GetNumberOfMaterials(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), InitMCDataMaterial(), and G4MaterialCutsCouple::IsUsed().

Referenced by Initialise().

◆ LoadMCDataElement()

void G4GSMottCorrection::LoadMCDataElement ( const G4Element * elem )

private

Definition at line 254 of file G4GSMottCorrection.cc.

                                                                {
  // allocate memory
  G4int izet = elem->GetZasInt();
  if (izet>gMaxZet) {
    izet = gMaxZet;
  }
  DataPerMaterial *perElem = new DataPerMaterial();
  AllocateDataPerMaterial(perElem);
  fMCDataPerElement[izet]  = perElem;
  //
  // load data from file
  char* tmppath = std::getenv("G4LEDATA");
  if (!tmppath) {
    G4Exception("G4GSMottCorrection::LoadMCDataElement()","em0006",
        FatalException,
        "Environment variable G4LEDATA not defined");
    return;
  }
  std::string path(tmppath);
  if (fIsElectron) {
    path += "/msc_GS/MottCor/el/";
  } else {
    path += "/msc_GS/MottCor/pos/";
  }
  std::string fname = path+"rej_"+gElemSymbols[izet-1];
  std::istringstream infile(std::ios::in);
  ReadCompressedFile(fname, infile);
  // check if file is open !!!
  for (G4int iek=0; iek<gNumEkin; ++iek) {
    DataPerEkin *perEkin = perElem->fDataPerEkin[iek];
    // 1. get the 3 Mott-correction factors for the current kinetic energy
    infile >> perEkin->fMCScreening;
    infile >> perEkin->fMCFirstMoment;
    infile >> perEkin->fMCSecondMoment;
    // 2. load each data per delta:
    for (G4int idel=0; idel<gNumDelta; ++idel) {
      DataPerDelta *perDelta = perEkin->fDataPerDelta[idel];
      // 2./a.  : first the rejection function values
      for (G4int iang=0; iang<gNumAngle; ++iang) {
        infile >> perDelta->fRejFuntion[iang];
      }
      // 2./b. : then the 4 spline parameter for the last bin
      infile >> perDelta->fSA;
      infile >> perDelta->fSB;
      infile >> perDelta->fSC;
      infile >> perDelta->fSD;
    }
  }
}

Referenced by InitMCDataPerElement().

◆ ReadCompressedFile()

void G4GSMottCorrection::ReadCompressedFile	(	std::string	fname,
		std::istringstream &	iss
	)

private

Definition at line 305 of file G4GSMottCorrection.cc.

                                                                                  {
  std::string *dataString = nullptr;
  std::string 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)
     int 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 std::string((char*)uncompdata, (long)complen);
     // delete the uncompressed data buffer
     delete [] uncompdata;
  } else {
    std::string msg = "  Problem while trying to read " + compfilename + " data file.\n";
    G4Exception("G4GSMottCorrection::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 LoadMCDataElement().

Field Documentation

◆ fInvDelAngle

G4double G4GSMottCorrection::fInvDelAngle

private

Definition at line 176 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), and GetMottRejectionValue().

◆ fInvDelBeta2

G4double G4GSMottCorrection::fInvDelBeta2

private

Definition at line 174 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), GetMottCorrectionFactors(), GetMottRejectionValue(), and InitMCDataMaterial().

◆ fInvDelDelta

G4double G4GSMottCorrection::fInvDelDelta

private

Definition at line 175 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), and GetMottRejectionValue().

◆ fInvLogDelEkin

G4double G4GSMottCorrection::fInvLogDelEkin

private

Definition at line 172 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), GetMottCorrectionFactors(), GetMottRejectionValue(), and InitMCDataMaterial().

◆ fIsElectron

G4bool G4GSMottCorrection::fIsElectron

private

Definition at line 159 of file G4GSMottCorrection.hh.

Referenced by LoadMCDataElement().

◆ fLogMinEkin

G4double G4GSMottCorrection::fLogMinEkin

private

Definition at line 171 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), GetMottCorrectionFactors(), GetMottRejectionValue(), and InitMCDataMaterial().

◆ fMaxEkin

G4double G4GSMottCorrection::fMaxEkin

private

Definition at line 170 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection().

◆ fMCDataPerElement

std::vector<DataPerMaterial*> G4GSMottCorrection::fMCDataPerElement

private

Definition at line 180 of file G4GSMottCorrection.hh.

Referenced by ClearMCDataPerElement(), InitMCDataMaterial(), InitMCDataPerElement(), and LoadMCDataElement().

◆ fMCDataPerMaterial

std::vector<DataPerMaterial*> G4GSMottCorrection::fMCDataPerMaterial

private

Definition at line 181 of file G4GSMottCorrection.hh.

Referenced by ClearMCDataPerMaterial(), GetMottCorrectionFactors(), GetMottRejectionValue(), InitMCDataMaterial(), and InitMCDataPerMaterials().

◆ fMinBeta2

G4double G4GSMottCorrection::fMinBeta2

private

Definition at line 173 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), GetMottCorrectionFactors(), GetMottRejectionValue(), and InitMCDataMaterial().

◆ gElemSymbols

const std::string G4GSMottCorrection::gElemSymbols

staticprivate

Initial value:

= {"H","He","Li","Be","B" ,
 "C" ,"N" ,"O" ,"F" ,"Ne","Na","Mg","Al","Si","P" , "S","Cl","Ar","K" ,"Ca","Sc",
 "Ti","V" ,"Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb",
 "Sr","Y" ,"Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I" ,
 "Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm",
 "Yb","Lu","Hf","Ta","W" ,"Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At",
 "Rn","Fr","Ra","Ac","Th","Pa","U" ,"Np","Pu","Am","Cm","Bk","Cf"}

Definition at line 178 of file G4GSMottCorrection.hh.

Referenced by LoadMCDataElement().

◆ gMaxBeta2

constexpr G4double G4GSMottCorrection::gMaxBeta2 = 0.9999

staticconstexprprivate

Definition at line 167 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), GetMottCorrectionFactors(), and GetMottRejectionValue().

◆ gMaxDelta

constexpr G4double G4GSMottCorrection::gMaxDelta = 0.9

staticconstexprprivate

Definition at line 168 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), and GetMottRejectionValue().

◆ gMaxZet

constexpr G4int G4GSMottCorrection::gMaxZet = 98

staticconstexprprivate

Definition at line 164 of file G4GSMottCorrection.hh.

Referenced by GetMaxZet(), InitMCDataMaterial(), InitMCDataPerElement(), and LoadMCDataElement().

◆ gMidEkin

constexpr G4double G4GSMottCorrection::gMidEkin = 100.*CLHEP::keV

staticconstexprprivate

Definition at line 166 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), and InitMCDataMaterial().

◆ gMinEkin

constexpr G4double G4GSMottCorrection::gMinEkin = 1.*CLHEP::keV

staticconstexprprivate

Definition at line 165 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection().

◆ gNumAngle

constexpr G4int G4GSMottCorrection::gNumAngle = 32

staticconstexprprivate

Definition at line 163 of file G4GSMottCorrection.hh.

Referenced by AllocateDataPerMaterial(), G4GSMottCorrection(), GetMottRejectionValue(), InitMCDataMaterial(), and LoadMCDataElement().

◆ gNumBeta2

constexpr G4int G4GSMottCorrection::gNumBeta2 = 16

staticconstexprprivate

Definition at line 161 of file G4GSMottCorrection.hh.

Referenced by G4GSMottCorrection(), GetMottCorrectionFactors(), GetMottRejectionValue(), and InitMCDataMaterial().

◆ gNumDelta

constexpr G4int G4GSMottCorrection::gNumDelta = 28

staticconstexprprivate

Definition at line 162 of file G4GSMottCorrection.hh.

Referenced by AllocateDataPerMaterial(), DeAllocateDataPerMaterial(), G4GSMottCorrection(), InitMCDataMaterial(), and LoadMCDataElement().

◆ gNumEkin

constexpr G4int G4GSMottCorrection::gNumEkin = 31

staticconstexprprivate

Definition at line 160 of file G4GSMottCorrection.hh.

Referenced by AllocateDataPerMaterial(), DeAllocateDataPerMaterial(), G4GSMottCorrection(), GetMottCorrectionFactors(), GetMottRejectionValue(), InitMCDataMaterial(), and LoadMCDataElement().

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

source/processes/electromagnetic/standard/include/G4GSMottCorrection.hh
source/processes/electromagnetic/standard/src/G4GSMottCorrection.cc

Data Structures

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ G4GSMottCorrection()

◆ ~G4GSMottCorrection()

Member Function Documentation

◆ AllocateDataPerMaterial()

◆ ClearMCDataPerElement()

◆ ClearMCDataPerMaterial()

◆ DeAllocateDataPerMaterial()

◆ GetMaxZet()

◆ GetMottCorrectionFactors()

◆ GetMottRejectionValue()

◆ Initialise()

◆ InitMCDataMaterial()

◆ InitMCDataPerElement()

◆ InitMCDataPerMaterials()

◆ LoadMCDataElement()

◆ ReadCompressedFile()

Field Documentation

◆ fInvDelAngle

◆ fInvDelBeta2

◆ fInvDelDelta

◆ fInvLogDelEkin

◆ fIsElectron

◆ fLogMinEkin

◆ fMaxEkin

◆ fMCDataPerElement

◆ fMCDataPerMaterial

◆ fMinBeta2

◆ gElemSymbols

◆ gMaxBeta2

◆ gMaxDelta

◆ gMaxZet

◆ gMidEkin

◆ gMinEkin

◆ gNumAngle

◆ gNumBeta2

◆ gNumDelta

◆ gNumEkin