#include <G4GoudsmitSaundersonTable.hh>

Classes
struct	GSMSCAngularDtr

Public Member Functions
	G4GoudsmitSaundersonTable (G4bool iselectron)

	~G4GoudsmitSaundersonTable ()

void	Initialise (G4double lownergylimit, G4double highenergylimit)

void	LoadMSCData ()

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)

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)

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

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	SetOptionMottCorrection (G4bool val)

void	SetOptionPWACorrection (G4bool val)

G4double	ComputeScatteringPowerCorrection (const G4MaterialCutsCouple *matcut, G4double ekin)

void	InitSCPCorrection ()

Detailed Description

Definition at line 83 of file G4GoudsmitSaundersonTable.hh.

Constructor & Destructor Documentation

◆ G4GoudsmitSaundersonTable()

G4GoudsmitSaundersonTable::G4GoudsmitSaundersonTable ( G4bool iselectron )

Definition at line 112 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;
}

◆ ~G4GoudsmitSaundersonTable()

G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable ( )

Definition at line 132 of file G4GoudsmitSaundersonTable.cc.

                                                      {
  for (std::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 (std::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 (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
    if (fSCPCPerMatCuts[imc]) {
      fSCPCPerMatCuts[imc]->fVSCPC.clear();
      delete fSCPCPerMatCuts[imc];
    }
  }
  fSCPCPerMatCuts.clear();
  //
  gIsInitialised = false;
}

Member Function Documentation

◆ ComputeScatteringPowerCorrection()

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

Definition at line 604 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;
  std::size_t lindx   = (std::size_t)remaining;
  remaining          -= lindx;
  std::size_t imax    = fSCPCPerMatCuts[imc]->fVSCPC.size()-1;
  if (lindx>=imax) {
    corFactor = fSCPCPerMatCuts[imc]->fVSCPC[imax];
  } else {
    corFactor = fSCPCPerMatCuts[imc]->fVSCPC[lindx] + remaining*(fSCPCPerMatCuts[imc]->fVSCPC[lindx+1]-fSCPCPerMatCuts[imc]->fVSCPC[lindx]);
  }
  return corFactor;
}

◆ GetGSAngularDtr()

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

Definition at line 360 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;
}

Referenced by SampleCosTheta().

◆ GetMoliereBc()

G4double G4GoudsmitSaundersonTable::GetMoliereBc ( G4int matindx )

inline

Definition at line 122 of file G4GoudsmitSaundersonTable.hh.

122{ return gMoliereBc[matindx]; }

Referenced by InitSCPCorrection().

◆ GetMoliereXc2()

G4double G4GoudsmitSaundersonTable::GetMoliereXc2 ( G4int matindx )

inline

Definition at line 124 of file G4GoudsmitSaundersonTable.hh.

124{ return gMoliereXc2[matindx]; }

Referenced by InitSCPCorrection().

◆ GetMottCorrectionFactors()

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

Definition at line 536 of file G4GoudsmitSaundersonTable.cc.

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

◆ Initialise()

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

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

◆ InitSCPCorrection()

void G4GoudsmitSaundersonTable::InitSCPCorrection ( )

Definition at line 625 of file G4GoudsmitSaundersonTable.cc.

                                                  {
  // get the material-cuts table
  G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
  std::size_t numMatCuts            = thePCTable->GetTableSize();
  // clear container if any
  for (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
    if (fSCPCPerMatCuts[imc]) {
      fSCPCPerMatCuts[imc]->fVSCPC.clear();
      delete fSCPCPerMatCuts[imc];
      fSCPCPerMatCuts[imc] = nullptr;
    }
  }
  //
  // set size of the container and create the corresponding data structures
  fSCPCPerMatCuts.resize(numMatCuts,nullptr);
  // loop over the material-cuts and create scattering power correction data structure for each
  for (G4int imc=0; imc<(G4int)numMatCuts; ++imc) {
    const G4MaterialCutsCouple *matCut =  thePCTable->GetMaterialCutsCouple(imc);
    // 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 = (G4int)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;
    }
  }
}

Referenced by Initialise().

◆ LoadMSCData()

void G4GoudsmitSaundersonTable::LoadMSCData ( )

Definition at line 438 of file G4GoudsmitSaundersonTable.cc.

                                            {
  gGSMSCAngularDistributions1.resize(gLAMBNUM*gQNUM1,nullptr);
  const G4String str1 = G4EmParameters::Instance()->GetDirLEDATA() + "/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) {
      auto 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 = G4EmParameters::Instance()->GetDirLEDATA() + "/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) {
        auto 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();
  }
}

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

Referenced by Sampling().

◆ SampleGSSRCosTheta()

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

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

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

◆ SetOptionMottCorrection()

void G4GoudsmitSaundersonTable::SetOptionMottCorrection ( G4bool val )

inline

Definition at line 131 of file G4GoudsmitSaundersonTable.hh.

131{ fIsMottCorrection = val; }

◆ SetOptionPWACorrection()

void G4GoudsmitSaundersonTable::SetOptionPWACorrection ( G4bool val )

inline

Definition at line 133 of file G4GoudsmitSaundersonTable.hh.

133{ fIsPWACorrection = val; }

◆ SingleScattering()

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

Definition at line 507 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;
}

Referenced by Sampling().

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

Classes

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ G4GoudsmitSaundersonTable()

◆ ~G4GoudsmitSaundersonTable()

Member Function Documentation

◆ ComputeScatteringPowerCorrection()

◆ GetGSAngularDtr()

◆ GetMoliereBc()

◆ GetMoliereXc2()

◆ GetMottCorrectionFactors()

◆ Initialise()

◆ InitSCPCorrection()

◆ LoadMSCData()

◆ SampleCosTheta()

◆ SampleGSSRCosTheta()

◆ Sampling()

◆ SetOptionMottCorrection()

◆ SetOptionPWACorrection()

◆ SingleScattering()