#include <G4GoudsmitSaundersonTable.hh>

Data Structures
struct	GSMSCAngularDtr

struct	SCPCorrection

Public Member Functions
G4double	ComputeScatteringPowerCorrection (const G4MaterialCutsCouple *matcut, G4double ekin)

	G4GoudsmitSaundersonTable (G4bool iselectron)

GSMSCAngularDtr *	GetGSAngularDtr (G4double scra, G4double &lambdaval, G4double &qval, G4double &transfpar)

G4double	GetMoliereBc (G4int matindx)

G4double	GetMoliereXc2 (G4int matindx)

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

void	Initialise (G4double lownergylimit, G4double highenergylimit)

void	InitSCPCorrection ()

void	LoadMSCData ()

G4double	SampleCosTheta (G4double lambdaval, G4double qval, G4double scra, G4double lekin, G4double beta2, G4int matindx, GSMSCAngularDtr **gsDtr, G4int &mcekini, G4int &mcdelti, G4double &transfPar, G4bool isfirst)

G4double	SampleGSSRCosTheta (const GSMSCAngularDtr *gsDrt, G4double transfpar)

G4bool	Sampling (G4double lambdaval, G4double qval, G4double scra, G4double &cost, G4double &sint, G4double lekin, G4double beta2, G4int matindx, GSMSCAngularDtr **gsDtr, G4int &mcekini, G4int &mcdelti, G4double &transfPar, G4bool isfirst)

void	SetOptionMottCorrection (G4bool val)

void	SetOptionPWACorrection (G4bool val)

G4double	SingleScattering (G4double lambdaval, G4double scra, G4double lekin, G4double beta2, G4int matindx)

	~G4GoudsmitSaundersonTable ()

Private Member Functions
void	InitMoliereMSCParams ()

Private Attributes
G4double	fDeltaQ2

G4double	fHighEnergyLimit

G4double	fInvDeltaQ1

G4double	fInvDeltaQ2

G4double	fInvLogDeltaLambda

G4bool	fIsElectron

G4bool	fIsMottCorrection

G4bool	fIsPWACorrection

G4double	fLogDeltaLambda

G4double	fLogLambda0

G4double	fLowEnergyLimit

int	fNumSPCEbinPerDec

std::vector< SCPCorrection * >	fSCPCPerMatCuts

Static Private Attributes
static std::vector< GSMSCAngularDtr * >	gGSMSCAngularDistributions1

static std::vector< GSMSCAngularDtr * >	gGSMSCAngularDistributions2

static G4bool	gIsInitialised = false

static constexpr G4double	gLAMBMAX = 100000.0

static constexpr G4double	gLAMBMIN = 1.0

static constexpr G4int	gLAMBNUM = 64

static constexpr G4int	gNUMSCR1 = 201

static constexpr G4int	gNUMSCR2 = 51

static constexpr G4double	gQMAX1 = 0.99

static constexpr G4double	gQMAX2 = 7.99

static constexpr G4double	gQMIN1 = 0.001

static constexpr G4double	gQMIN2 = 0.99

static constexpr G4int	gQNUM1 = 15

static constexpr G4int	gQNUM2 = 32

Detailed Description

Definition at line 83 of file G4GoudsmitSaundersonTable.hh.

Constructor & Destructor Documentation

◆ G4GoudsmitSaundersonTable()

G4GoudsmitSaundersonTable::G4GoudsmitSaundersonTable ( G4bool iselectron )

Definition at line 111 of file G4GoudsmitSaundersonTable.cc.

                                                                      {
  fIsElectron         = iselectron;
  // set initial values: final values will be set in the Initialize method
  fLogLambda0         = 0.;              // will be set properly at init.
  fLogDeltaLambda     = 0.;              // will be set properly at init.
  fInvLogDeltaLambda  = 0.;              // will be set properly at init.
  fInvDeltaQ1         = 0.;              // will be set properly at init.
  fDeltaQ2            = 0.;              // will be set properly at init.
  fInvDeltaQ2         = 0.;              // will be set properly at init.
  //
  fLowEnergyLimit     =   0.1*CLHEP::keV; // will be set properly at init.
  fHighEnergyLimit    = 100.0*CLHEP::MeV; // will be set properly at init.
  //
  fIsMottCorrection   = false;            // will be set properly at init.
  fIsPWACorrection    = false;            // will be set properly at init.
  fMottCorrection     = nullptr;
  //
  fNumSPCEbinPerDec   = 3;
}

References fDeltaQ2, fHighEnergyLimit, fInvDeltaQ1, fInvDeltaQ2, fInvLogDeltaLambda, fIsElectron, fIsMottCorrection, fIsPWACorrection, fLogDeltaLambda, fLogLambda0, fLowEnergyLimit, fNumSPCEbinPerDec, CLHEP::keV, and CLHEP::MeV.

◆ ~G4GoudsmitSaundersonTable()

G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable ( )

Definition at line 131 of file G4GoudsmitSaundersonTable.cc.

                                                      {
  for (size_t i=0; i<gGSMSCAngularDistributions1.size(); ++i) {
    if (gGSMSCAngularDistributions1[i]) {
      delete [] gGSMSCAngularDistributions1[i]->fUValues;
      delete [] gGSMSCAngularDistributions1[i]->fParamA;
      delete [] gGSMSCAngularDistributions1[i]->fParamB;
      delete gGSMSCAngularDistributions1[i];
    }
  }
  gGSMSCAngularDistributions1.clear();
  for (size_t i=0; i<gGSMSCAngularDistributions2.size(); ++i) {
    if (gGSMSCAngularDistributions2[i]) {
      delete [] gGSMSCAngularDistributions2[i]->fUValues;
      delete [] gGSMSCAngularDistributions2[i]->fParamA;
      delete [] gGSMSCAngularDistributions2[i]->fParamB;
      delete gGSMSCAngularDistributions2[i];
    }
  }
  gGSMSCAngularDistributions2.clear();
  if (fMottCorrection) {
    delete fMottCorrection;
    fMottCorrection = nullptr;
  }
  // clear scp correction data
  for (size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
    if (fSCPCPerMatCuts[imc]) {
      fSCPCPerMatCuts[imc]->fVSCPC.clear();
      delete fSCPCPerMatCuts[imc];
    }
  }
  fSCPCPerMatCuts.clear();
  //
  gIsInitialised = false;
}

References fSCPCPerMatCuts, gGSMSCAngularDistributions1, gGSMSCAngularDistributions2, and gIsInitialised.

Member Function Documentation

◆ ComputeScatteringPowerCorrection()

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

Definition at line 611 of file G4GoudsmitSaundersonTable.cc.

                                                                                                                      {
  G4int    imc       = matcut->GetIndex();
  G4double corFactor = 1.0;
  if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin<=fSCPCPerMatCuts[imc]->fPrCut) {
    return corFactor;
  }
  // get the scattering power correction factor
  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 G4GoudsmitSaundersonMscModel::GetTransportMeanFreePath(), and G4GoudsmitSaundersonMscModel::GetTransportMeanFreePathOnly().

◆ GetGSAngularDtr()

G4GoudsmitSaundersonTable::GSMSCAngularDtr * G4GoudsmitSaundersonTable::GetGSAngularDtr	(	G4double	scra,
		G4double &	lambdaval,
		G4double &	qval,
		G4double &	transfpar
	)

Definition at line 359 of file G4GoudsmitSaundersonTable.cc.

                                                                                                            {
  GSMSCAngularDtr *dtr = nullptr;
  G4bool first         = false;
  // isotropic cost above gQMAX2 (i.e. dtr stays nullptr)
  if (qval<gQMAX2) {
    G4int    lamIndx  = -1; // lambda value index
    G4int    qIndx    = -1; // lambda value index
    // init to second grid Q values
    G4int    numQVal  = gQNUM2;
    G4double minQVal  = gQMIN2;
    G4double invDelQ  = fInvDeltaQ2;
    G4double pIndxH   = 0.; // probability of taking higher index
    // check if first or second grid needs to be used
    if (qval<gQMIN2) {  // first grid
      first = true;
      // protect against qval<gQMIN1
      if (qval<gQMIN1) {
        qval   = gQMIN1;
        qIndx  = 0;
        //pIndxH = 0.;
      }
      // set to first grid Q values
      numQVal  = gQNUM1;
      minQVal  = gQMIN1;
      invDelQ  = fInvDeltaQ1;
    }
    // make sure that lambda = s/lambda_el is in [gLAMBMIN,gLAMBMAX)
    // lambda<gLAMBMIN=1 is already handeled before so lambda>= gLAMBMIN for sure
    if (lambdaval>=gLAMBMAX) {
      lambdaval = gLAMBMAX-1.e-8;
      lamIndx   = gLAMBNUM-1;
    }
    G4double lLambda  = G4Log(lambdaval);
    //
    // determine lower lambda (=s/lambda_el) index: linear interp. on log(lambda) scale
    if (lamIndx<0) {
      pIndxH  = (lLambda-fLogLambda0)*fInvLogDeltaLambda;
      lamIndx = (G4int)(pIndxH);        // lower index of the lambda bin
      pIndxH  = pIndxH-lamIndx;       // probability of taking the higher index distribution
      if (G4UniformRand()<pIndxH) {
        ++lamIndx;
      }
    }
    //
    // determine lower Q (=s/lambda_el G1) index: linear interp. on Q
    if (qIndx<0) {
      pIndxH  = (qval-minQVal)*invDelQ;
      qIndx   = (G4int)(pIndxH);        // lower index of the Q bin
      pIndxH  = pIndxH-qIndx;
      if (G4UniformRand()<pIndxH) {
        ++qIndx;
      }
    }
    // set indx
    G4int indx = lamIndx*numQVal+qIndx;
    if (first) {
      dtr = gGSMSCAngularDistributions1[indx];
    } else {
      dtr = gGSMSCAngularDistributions2[indx];
    }
    // dtr might be nullptr that indicates isotropic cot distribution because:
    // - if the selected lamIndx, qIndx correspond to L(=s/lambda_el) and Q(=s/lambda_el G1) such that G1(=Q/L) > 1
    //   G1 should always be < 1 and if G1 is ~1 -> the dtr is isotropic (this can only happen in case of the 2. grid)
    //
    // compute the transformation parameter
    if (lambdaval>10.0) {
      transfpar = 0.5*(-2.77164+lLambda*( 2.94874-lLambda*(0.1535754-lLambda*0.00552888) ));
    } else {
      transfpar = 0.5*(1.347+lLambda*(0.209364-lLambda*(0.45525-lLambda*(0.50142-lLambda*0.081234))));
    }
    transfpar *= (lambdaval+4.0)*scra;
  }
  // return with the selected GS angular distribution that we need to sample cost from (if nullptr => isotropic cost)
  return dtr;
}

References fInvDeltaQ1, fInvDeltaQ2, fInvLogDeltaLambda, fLogLambda0, G4Log(), G4UniformRand, gGSMSCAngularDistributions1, gGSMSCAngularDistributions2, gLAMBMAX, gLAMBNUM, gQMAX2, gQMIN1, gQMIN2, gQNUM1, and gQNUM2.

Referenced by SampleCosTheta().

◆ GetMoliereBc()

G4double G4GoudsmitSaundersonTable::GetMoliereBc ( G4int matindx )

inline

Definition at line 122 of file G4GoudsmitSaundersonTable.hh.

122{ return gMoliereBc[matindx]; }

Referenced by G4GoudsmitSaundersonMscModel::CrossSectionPerVolume(), G4GoudsmitSaundersonMscModel::GetTransportMeanFreePath(), G4GoudsmitSaundersonMscModel::GetTransportMeanFreePathOnly(), and InitSCPCorrection().

◆ GetMoliereXc2()

G4double G4GoudsmitSaundersonTable::GetMoliereXc2 ( G4int matindx )

inline

Definition at line 124 of file G4GoudsmitSaundersonTable.hh.

124{ return gMoliereXc2[matindx]; }

Referenced by G4GoudsmitSaundersonMscModel::CrossSectionPerVolume(), G4GoudsmitSaundersonMscModel::GetTransportMeanFreePath(), G4GoudsmitSaundersonMscModel::GetTransportMeanFreePathOnly(), and InitSCPCorrection().

◆ GetMottCorrectionFactors()

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

Definition at line 543 of file G4GoudsmitSaundersonTable.cc.

                                                                                                   {
  if (fIsMottCorrection) {
    fMottCorrection->GetMottCorrectionFactors(logekin, beta2, matindx, mcToScr, mcToQ1, mcToG2PerG1);
  }
}

References fIsMottCorrection.

Referenced by G4GoudsmitSaundersonMscModel::CrossSectionPerVolume(), G4GoudsmitSaundersonMscModel::GetTransportMeanFreePath(), and G4GoudsmitSaundersonMscModel::GetTransportMeanFreePathOnly().

◆ Initialise()

void G4GoudsmitSaundersonTable::Initialise	(	G4double	lownergylimit,
		G4double	highenergylimit
	)

Definition at line 166 of file G4GoudsmitSaundersonTable.cc.

                                                                                           {
  fLowEnergyLimit     = lownergylimit;
  fHighEnergyLimit    = highenergylimit;
  G4double lLambdaMin = G4Log(gLAMBMIN);
  G4double lLambdaMax = G4Log(gLAMBMAX);
  fLogLambda0         = lLambdaMin;
  fLogDeltaLambda     = (lLambdaMax-lLambdaMin)/(gLAMBNUM-1.);
  fInvLogDeltaLambda  = 1./fLogDeltaLambda;
  fInvDeltaQ1         = 1./((gQMAX1-gQMIN1)/(gQNUM1-1.));
  fDeltaQ2            = (gQMAX2-gQMIN2)/(gQNUM2-1.);
  fInvDeltaQ2         = 1./fDeltaQ2;
  // load precomputed angular distributions and set up several values used during the sampling
  // these are particle independet => they go to static container: load them only onece
  if (!gIsInitialised) {
    // load pre-computed GS angular distributions (computed based on Screened-Rutherford DCS)
    LoadMSCData();
    gIsInitialised = true;
  }
  InitMoliereMSCParams();
  // Mott-correction: particle(e- or e+) dependet so init them
  if (fIsMottCorrection) {
    if (!fMottCorrection) {
      fMottCorrection = new G4GSMottCorrection(fIsElectron);
    }
    fMottCorrection->Initialise();
  }
  // init scattering power correction data; used only together with Mott-correction
  // (Moliere's parameters must be initialised before)
  if (fMottCorrection) {
    InitSCPCorrection();
  }
}

References fDeltaQ2, fHighEnergyLimit, fInvDeltaQ1, fInvDeltaQ2, fInvLogDeltaLambda, fIsElectron, fIsMottCorrection, fLogDeltaLambda, fLogLambda0, fLowEnergyLimit, G4Log(), gIsInitialised, gLAMBMAX, gLAMBMIN, gLAMBNUM, gQMAX1, gQMAX2, gQMIN1, gQMIN2, gQNUM1, gQNUM2, InitMoliereMSCParams(), InitSCPCorrection(), and LoadMSCData().

Referenced by G4GoudsmitSaundersonMscModel::Initialise().

◆ InitMoliereMSCParams()

void G4GoudsmitSaundersonTable::InitMoliereMSCParams ( )

private

Definition at line 553 of file G4GoudsmitSaundersonTable.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
 
   G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
   // get number of materials in the table
   size_t numMaterials = theMaterialTable->size();
   // make sure that we have long enough vectors
   if(gMoliereBc.size()<numMaterials) {
     gMoliereBc.resize(numMaterials);
     gMoliereXc2.resize(numMaterials);
   }
   G4double xi   = 1.0;
   G4int    maxZ = 200;
   if (fIsMottCorrection || fIsPWACorrection) {
     // xi   = 1.0;  <= always set to 1 from now on
     maxZ = G4GSMottCorrection::GetMaxZet();
   }
   //
   for (size_t imat=0; imat<numMaterials; ++imat) {
     const G4Material*      theMaterial     = (*theMaterialTable)[imat];
     const G4ElementVector* theElemVect     = theMaterial->GetElementVector();
     const G4int            numelems        = theMaterial->GetNumberOfElements();
     //
     const G4double*        theNbAtomsPerVolVect  = theMaterial->GetVecNbOfAtomsPerVolume();
     G4double               theTotNbAtomsPerVol   = theMaterial->GetTotNbOfAtomsPerVolume();
     //
     G4double zs = 0.0;
     G4double zx = 0.0;
     G4double ze = 0.0;
     G4double sa = 0.0;
     //
     for(G4int ielem = 0; ielem < numelems; ielem++) {
       G4double zet = (*theElemVect)[ielem]->GetZ();
       if (zet>maxZ) {
         zet = (G4double)maxZ;
       }
       G4double iwa  = (*theElemVect)[ielem]->GetN();
       G4double ipz  = theNbAtomsPerVolVect[ielem]/theTotNbAtomsPerVol;
       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 = theMaterial->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
     //
     gMoliereBc[theMaterial->GetIndex()]  = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs);  //[1/cm]
     gMoliereXc2[theMaterial->GetIndex()] = const2*density*zs/sa;  // [MeV2/cm]
     // change to Geant4 internal units of 1/length and energ2/length
     gMoliereBc[theMaterial->GetIndex()]  *= 1.0/CLHEP::cm;
     gMoliereXc2[theMaterial->GetIndex()] *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
   }
}

References CLHEP::cm, CLHEP::cm3, fIsMottCorrection, fIsPWACorrection, CLHEP::g, G4Exp(), G4Log(), G4Material::GetDensity(), G4Material::GetElementVector(), G4Material::GetIndex(), G4Material::GetMaterialTable(), G4GSMottCorrection::GetMaxZet(), G4Material::GetNumberOfElements(), G4Material::GetTotNbOfAtomsPerVolume(), G4Material::GetVecNbOfAtomsPerVolume(), maxZ, and CLHEP::MeV.

Referenced by Initialise().

◆ InitSCPCorrection()

void G4GoudsmitSaundersonTable::InitSCPCorrection ( )

Definition at line 632 of file G4GoudsmitSaundersonTable.cc.

                                                  {
  // get the material-cuts table
  G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
  size_t numMatCuts                 = thePCTable->GetTableSize();
  // clear container if any
  for (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 (size_t imc=0; imc<numMatCuts; ++imc) {
    const G4MaterialCutsCouple *matCut =  thePCTable->GetMaterialCutsCouple(imc);
    // get e- production cut in the current material-cuts in energy
    G4double limit;
    G4double ecut;
    if (fIsElectron) {
      ecut  = (*(thePCTable->GetEnergyCutsVector(idxG4ElectronCut)))[matCut->GetIndex()];
      limit = 2.*ecut;
    } else {
      ecut  = (*(thePCTable->GetEnergyCutsVector(idxG4PositronCut)))[matCut->GetIndex()];
      limit = ecut;
    }
    G4double min = std::max(limit,fLowEnergyLimit);
    G4double max = fHighEnergyLimit;
    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);
    G4double lmin        = G4Log(min);
    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;
    for (G4int ie=0; ie<numEbins; ++ie) {
      G4double ekin    = G4Exp(lmin+ie*ldel);
      G4double scpCorr = 1.0;
      // compute correction factor: I.Kawrakow NIMB 114(1996)307-326 (Eqs(32-37))
      if (ie>0) {
         G4double tau     = ekin/CLHEP::electron_mass_c2;
         G4double tauCut  = ecut/CLHEP::electron_mass_c2;
         // Moliere's screening parameter
         G4int    matindx = matCut->GetMaterial()->GetIndex();
         G4double A       = GetMoliereXc2(matindx)/(4.0*tau*(tau+2.)*GetMoliereBc(matindx));
         G4double gr      = (1.+2.*A)*G4Log(1.+1./A)-2.;
         G4double dum0    = (tau+2.)/(tau+1.);
         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.;
         }
         G4double z0 = matCut->GetMaterial()->GetIonisation()->GetZeffective();
         scpCorr     = 1.-gm*z0/(z0*(z0+1.));
      }
      fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCorr;
    }
  }
}

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

Referenced by Initialise().

◆ LoadMSCData()

void G4GoudsmitSaundersonTable::LoadMSCData ( )

Definition at line 437 of file G4GoudsmitSaundersonTable.cc.

                                            {
  char* path = std::getenv("G4LEDATA");
  if (!path) {
    G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
        FatalException,
        "Environment variable G4LEDATA not defined");
    return;
  }
  //
  gGSMSCAngularDistributions1.resize(gLAMBNUM*gQNUM1,nullptr);
  const G4String str1 = G4String(path) + "/msc_GS/GSGrid_1/gsDistr_";
  for (G4int il=0; il<gLAMBNUM; ++il) {
    G4String fname = str1 + std::to_string(il);
    std::ifstream infile(fname,std::ios::in);
    if (!infile.is_open()) {
      G4String msgc = "Cannot open file: " + fname;
      G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
          FatalException, msgc.c_str());
      return;
    }
    for (G4int iq=0; iq<gQNUM1; ++iq) {
      GSMSCAngularDtr *gsd = new GSMSCAngularDtr();
      infile >> gsd->fNumData;
      gsd->fUValues = new G4double[gsd->fNumData]();
      gsd->fParamA  = new G4double[gsd->fNumData]();
      gsd->fParamB  = new G4double[gsd->fNumData]();
      G4double ddummy;
      infile >> ddummy; infile >> ddummy;
      for (G4int i=0; i<gsd->fNumData; ++i) {
        infile >> gsd->fUValues[i];
        infile >> gsd->fParamA[i];
        infile >> gsd->fParamB[i];
      }
      gGSMSCAngularDistributions1[il*gQNUM1+iq] = gsd;
    }
    infile.close();
  }
  //
  // second grid
  gGSMSCAngularDistributions2.resize(gLAMBNUM*gQNUM2,nullptr);
  const G4String str2 = G4String(path) + "/msc_GS/GSGrid_2/gsDistr_";
  for (G4int il=0; il<gLAMBNUM; ++il) {
    G4String fname = str2 + std::to_string(il);
    std::ifstream infile(fname,std::ios::in);
    if (!infile.is_open()) {
      G4String msgc = "Cannot open file: " + fname;
      G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
          FatalException, msgc.c_str());
      return;
    }
    for (G4int iq=0; iq<gQNUM2; ++iq) {
      G4int numData;
      infile >> numData;
      if (numData>1) {
        GSMSCAngularDtr *gsd = new GSMSCAngularDtr();
        gsd->fNumData = numData;
        gsd->fUValues = new G4double[gsd->fNumData]();
        gsd->fParamA  = new G4double[gsd->fNumData]();
        gsd->fParamB  = new G4double[gsd->fNumData]();
        double ddummy;
        infile >> ddummy; infile >> ddummy;
        for (G4int i=0; i<gsd->fNumData; ++i) {
          infile >> gsd->fUValues[i];
          infile >> gsd->fParamA[i];
          infile >> gsd->fParamB[i];
        }
        gGSMSCAngularDistributions2[il*gQNUM2+iq] = gsd;
      } else {
        gGSMSCAngularDistributions2[il*gQNUM2+iq] = nullptr;
      }
    }
    infile.close();
  }
}

References FatalException, test::fname, G4GoudsmitSaundersonTable::GSMSCAngularDtr::fNumData, G4GoudsmitSaundersonTable::GSMSCAngularDtr::fParamA, G4GoudsmitSaundersonTable::GSMSCAngularDtr::fParamB, G4GoudsmitSaundersonTable::GSMSCAngularDtr::fUValues, G4Exception(), gGSMSCAngularDistributions1, gGSMSCAngularDistributions2, gLAMBNUM, gQNUM1, and gQNUM2.

Referenced by Initialise().

◆ SampleCosTheta()

G4double G4GoudsmitSaundersonTable::SampleCosTheta	(	G4double	lambdaval,
		G4double	qval,
		G4double	scra,
		G4double	lekin,
		G4double	beta2,
		G4int	matindx,
		GSMSCAngularDtr **	gsDtr,
		G4int &	mcekini,
		G4int &	mcdelti,
		G4double &	transfPar,
		G4bool	isfirst
	)

Definition at line 305 of file G4GoudsmitSaundersonTable.cc.

                                                                                        {
  G4double cost = 1.;
  // determine the base GS angular distribution if it is the first call (when sub-step sampling is used)
  if (isfirst) {
    *gsDtr = GetGSAngularDtr(scra, lambdaval, qval, transfPar);
  }
  // sample cost from the GS angular distribution (computed based on Screened-Rutherford DCS)
  cost = SampleGSSRCosTheta(*gsDtr, transfPar);
  // Mott-correction if it was requested by the user
  if (fIsMottCorrection && *gsDtr) {     // no Mott-correction in case of izotropic theta
    static const G4int nlooplim = 1000;
    G4int    nloop    =  0 ; // rejection loop counter
//    G4int    ekindx   = -1; // evaluate only in the first call
//    G4int    deltindx = -1 ; // evaluate only in the first call
    G4double val      = fMottCorrection->GetMottRejectionValue(lekin, beta2, qval, cost, matindx, mcekini, mcdelti);
    while (G4UniformRand()>val && ++nloop<nlooplim) {
      // sampling cos(theta)
      cost = SampleGSSRCosTheta(*gsDtr, transfPar);
      val  = fMottCorrection->GetMottRejectionValue(lekin, beta2, qval, cost, matindx, mcekini, mcdelti);
    };
  }
  return cost;
}

References fIsMottCorrection, G4UniformRand, GetGSAngularDtr(), and SampleGSSRCosTheta().

Referenced by Sampling().

◆ SampleGSSRCosTheta()

G4double G4GoudsmitSaundersonTable::SampleGSSRCosTheta	(	const GSMSCAngularDtr *	gsDrt,
		G4double	transfpar
	)

Definition at line 334 of file G4GoudsmitSaundersonTable.cc.

                                                                                                       {
  // check if isotropic theta (i.e. cost is uniform on [-1:1])
  if (!gsDtr) {
    return 1.-2.0*G4UniformRand();
  }
  //
  // sampling form the selected distribution
  G4double ndatm1 = gsDtr->fNumData-1.;
  G4double delta  = 1.0/ndatm1;
  // determine lower cumulative bin inidex
  G4double rndm   = G4UniformRand();
  G4int indxl     = rndm*ndatm1;
  G4double  aval  = rndm-indxl*delta;
  G4double  dum0  = delta*aval;
 
  G4double  dum1  = (1.0+gsDtr->fParamA[indxl]+gsDtr->fParamB[indxl])*dum0;
  G4double  dum2  = delta*delta + gsDtr->fParamA[indxl]*dum0 + gsDtr->fParamB[indxl]*aval*aval;
  G4double sample = gsDtr->fUValues[indxl] +  dum1/dum2 *(gsDtr->fUValues[indxl+1]-gsDtr->fUValues[indxl]);
  // transform back u to cos(theta) :
  // this is the sampled cos(theta) = (2.0*para*sample)/(1.0-sample+para)
  return 1.-(2.0*transfpar*sample)/(1.0-sample+transfpar);
}

References G4GoudsmitSaundersonTable::GSMSCAngularDtr::fNumData, G4GoudsmitSaundersonTable::GSMSCAngularDtr::fParamA, G4GoudsmitSaundersonTable::GSMSCAngularDtr::fParamB, G4GoudsmitSaundersonTable::GSMSCAngularDtr::fUValues, and G4UniformRand.

Referenced by SampleCosTheta().

◆ Sampling()

G4bool G4GoudsmitSaundersonTable::Sampling	(	G4double	lambdaval,
		G4double	qval,
		G4double	scra,
		G4double &	cost,
		G4double &	sint,
		G4double	lekin,
		G4double	beta2,
		G4int	matindx,
		GSMSCAngularDtr **	gsDtr,
		G4int &	mcekini,
		G4int &	mcdelti,
		G4double &	transfPar,
		G4bool	isfirst
	)

Definition at line 212 of file G4GoudsmitSaundersonTable.cc.

                                                                                {
  G4double rand0 = G4UniformRand();
  G4double expn  = G4Exp(-lambdaval);
  //
  // no scattering case
  if (rand0<expn) {
    cost = 1.0;
    sint = 0.0;
    return false;
  }
  //
  // single scattering case : sample from the single scattering PDF
  // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true)
  if (rand0<(1.+lambdaval)*expn) {
    // cost is sampled in SingleScattering()
    cost = SingleScattering(lambdaval, scra, lekin, beta2, matindx);
    // add protections
    if (cost<-1.0) cost = -1.0;
    if (cost>1.0)  cost =  1.0;
    // compute sin(theta) from the sampled cos(theta)
    G4double dum0 = 1.-cost;
    sint = std::sqrt(dum0*(2.0-dum0));
    return false;
  }
  //
  // handle this case:
  //      -lambdaval < 1 i.e. mean #elastic events along the step is < 1 but
  //       the currently sampled case is not 0 or 1 scattering. [Our minimal
  //       lambdaval (that we have precomputed, transformed angular distributions
  //       stored in a form of equally probabe intervalls together with rational
  //       interp. parameters) is 1.]
  //      -probability of having n elastic events follows Poisson stat. with
  //       lambdaval parameter.
  //      -the max. probability (when lambdaval=1) of having more than one
  //       elastic events is 0.2642411 and the prob of having 2,3,..,n elastic
  //       events decays rapidly with n. So set a max n to 10.
  //      -sampling of this cases is done in a one-by-one single elastic event way
  //       where the current #elastic event is sampled from the Poisson distr.
  if (lambdaval<1.0) {
    G4double prob, cumprob;
    prob = cumprob = expn;
    G4double curcost,cursint;
    // init cos(theta) and sin(theta) to the zero scattering values
    cost = 1.0;
    sint = 0.0;
    for (G4int iel=1; iel<10; ++iel) {
      // prob of having iel scattering from Poisson
      prob    *= lambdaval/(G4double)iel;
      cumprob += prob;
      //
      //sample cos(theta) from the singe scattering pdf:
      // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true)
      curcost       = SingleScattering(lambdaval, scra, lekin, beta2, matindx);
      G4double dum0 = 1.-curcost;
      cursint       = dum0*(2.0-dum0); // sin^2(theta)
      //
      // if we got current deflection that is not too small
      // then update cos(theta) sin(theta)
      if (cursint>1.0e-20) {
        cursint         = std::sqrt(cursint);
        G4double curphi = CLHEP::twopi*G4UniformRand();
        cost            = cost*curcost-sint*cursint*std::cos(curphi);
        sint            = std::sqrt(std::max(0.0, (1.0-cost)*(1.0+cost)));
      }
      //
      // check if we have done enough scattering i.e. sampling from the Poisson
      if (rand0<cumprob) {
        return false;
      }
    }
    // if reached the max iter i.e. 10
    return false;
  }
  //
  // multiple scattering case with lambdavalue >= 1:
  //   - use the precomputed and transformed Goudsmit-Saunderson angular
  //     distributions to sample cos(theta)
  //   - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true)
  cost = SampleCosTheta(lambdaval, qval, scra, lekin, beta2, matindx, gsDtr, mcekini, mcdelti, transfPar, isfirst);
  // add protections
  if (cost<-1.0)  cost = -1.0;
  if (cost> 1.0)  cost =  1.0;
  // compute cos(theta) and sin(theta) from the sampled 1-cos(theta)
  G4double dum0 = 1.0-cost;
  sint = std::sqrt(dum0*(2.0-dum0));
  // return true if it was msc
  return true;
}

References G4Exp(), G4UniformRand, G4INCL::Math::max(), SampleCosTheta(), SingleScattering(), and CLHEP::twopi.

Referenced by G4GoudsmitSaundersonMscModel::SampleMSC().

◆ SetOptionMottCorrection()

void G4GoudsmitSaundersonTable::SetOptionMottCorrection ( G4bool val )

inline

Definition at line 131 of file G4GoudsmitSaundersonTable.hh.

131{ fIsMottCorrection = val; }

References fIsMottCorrection.

Referenced by G4GoudsmitSaundersonMscModel::Initialise().

◆ SetOptionPWACorrection()

void G4GoudsmitSaundersonTable::SetOptionPWACorrection ( G4bool val )

inline

Definition at line 133 of file G4GoudsmitSaundersonTable.hh.

133{ fIsPWACorrection = val; }

References fIsPWACorrection.

Referenced by G4GoudsmitSaundersonMscModel::Initialise().

◆ SingleScattering()

G4double G4GoudsmitSaundersonTable::SingleScattering	(	G4double	lambdaval,
		G4double	scra,
		G4double	lekin,
		G4double	beta2,
		G4int	matindx
	)

Definition at line 514 of file G4GoudsmitSaundersonTable.cc.

                                                                    {
  G4double rand1 = G4UniformRand();
  // sample cost from the Screened-Rutherford DCS
  G4double cost  = 1.-2.0*scra*rand1/(1.0-rand1+scra);
  // Mott-correction if it was requested by the user
  if (fIsMottCorrection) {
    static const G4int nlooplim = 1000;  // rejection loop limit
    G4int    nloop    =  0 ; // loop counter
    G4int    ekindx   = -1 ; // evaluate only in the first call
    G4int    deltindx =  0 ; // single scattering case
    G4double q1       =  0.; // not used when deltindx = 0;
    // computing Mott rejection function value
    G4double val      = fMottCorrection->GetMottRejectionValue(lekin, beta2, q1, cost,
                                                               matindx, ekindx, deltindx);
    while (G4UniformRand()>val && ++nloop<nlooplim) {
      // sampling cos(theta) from the Screened-Rutherford DCS
      rand1 = G4UniformRand();
      cost  = 1.-2.0*scra*rand1/(1.0-rand1+scra);
      // computing Mott rejection function value
      val   = fMottCorrection->GetMottRejectionValue(lekin, beta2, q1, cost, matindx,
                                                     ekindx, deltindx);
    };
  }
  return cost;
}