G4GoudsmitSaundersonMscModel Class Reference

#include <G4GoudsmitSaundersonMscModel.hh>

Inheritance diagram for G4GoudsmitSaundersonMscModel:

Public Member Functions
	G4GoudsmitSaundersonMscModel (const G4String &nam="GoudsmitSaunderson")
virtual	~G4GoudsmitSaundersonMscModel ()
virtual void	Initialise (const G4ParticleDefinition *, const G4DataVector &)
void	StartTracking (G4Track *)
virtual G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double, G4double, G4double)
virtual G4ThreeVector &	SampleScattering (const G4DynamicParticle *, G4double safety)
virtual G4double	ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep)
virtual G4double	ComputeGeomPathLength (G4double truePathLength)
virtual G4double	ComputeTrueStepLength (G4double geomStepLength)
Public Member Functions inherited from G4VMscModel
	G4VMscModel (const G4String &nam)
virtual	~G4VMscModel ()
virtual G4double	ComputeTruePathLengthLimit (const G4Track &track, G4double &stepLimit)
virtual G4double	ComputeGeomPathLength (G4double truePathLength)
virtual G4double	ComputeTrueStepLength (G4double geomPathLength)
virtual G4ThreeVector &	SampleScattering (const G4DynamicParticle *, G4double safety)
virtual void	SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double tmax)
void	SetStepLimitType (G4MscStepLimitType)
void	SetLateralDisplasmentFlag (G4bool val)
void	SetRangeFactor (G4double)
void	SetGeomFactor (G4double)
void	SetSkin (G4double)
void	SetSampleZ (G4bool)
G4VEnergyLossProcess *	GetIonisation () const
void	SetIonisation (G4VEnergyLossProcess *, const G4ParticleDefinition *part)
G4double	ComputeSafety (const G4ThreeVector &position, G4double limit)
G4double	ComputeGeomLimit (const G4Track &, G4double &presafety, G4double limit)
G4double	GetDEDX (const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
G4double	GetRange (const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
G4double	GetEnergy (const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
G4double	GetTransportMeanFreePath (const G4ParticleDefinition *part, G4double kinEnergy)
Public Member Functions inherited from G4VEmModel
	G4VEmModel (const G4String &nam)
virtual	~G4VEmModel ()
virtual void	Initialise (const G4ParticleDefinition *, const G4DataVector &)=0
virtual void	SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
virtual G4double	ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
virtual G4double	CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double	ChargeSquareRatio (const G4Track &)
virtual G4double	GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual G4double	GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual void	StartTracking (G4Track *)
virtual void	CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
virtual G4double	Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
virtual G4double	MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *)
virtual void	SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual void	DefineForRegion (const G4Region *)
void	InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
G4double	ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
G4double	CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double	ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4int	SelectIsotopeNumber (const G4Element *)
const G4Element *	SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element *	SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void	SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
void	SetCrossSectionTable (G4PhysicsTable *)
G4PhysicsTable *	GetCrossSectionTable ()
G4VEmFluctuationModel *	GetModelOfFluctuations ()
G4VEmAngularDistribution *	GetAngularDistribution ()
void	SetAngularDistribution (G4VEmAngularDistribution *)
G4double	HighEnergyLimit () const
G4double	LowEnergyLimit () const
G4double	HighEnergyActivationLimit () const
G4double	LowEnergyActivationLimit () const
G4double	PolarAngleLimit () const
G4double	SecondaryThreshold () const
G4bool	LPMFlag () const
G4bool	DeexcitationFlag () const
G4bool	ForceBuildTableFlag () const
void	SetHighEnergyLimit (G4double)
void	SetLowEnergyLimit (G4double)
void	SetActivationHighEnergyLimit (G4double)
void	SetActivationLowEnergyLimit (G4double)
G4bool	IsActive (G4double kinEnergy)
void	SetPolarAngleLimit (G4double)
void	SetSecondaryThreshold (G4double)
void	SetLPMFlag (G4bool val)
void	SetDeexcitationFlag (G4bool val)
void	ForceBuildTable (G4bool val)
G4double	MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
const G4String &	GetName () const
void	SetCurrentCouple (const G4MaterialCutsCouple *)
const G4Element *	GetCurrentElement () const

Additional Inherited Members
Protected Member Functions inherited from G4VMscModel
G4ParticleChangeForMSC *	GetParticleChangeForMSC (const G4ParticleDefinition *p=0)
G4double	ConvertTrueToGeom (G4double &tLength, G4double &gLength)
Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss *	GetParticleChangeForLoss ()
G4ParticleChangeForGamma *	GetParticleChangeForGamma ()
virtual G4double	MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
const G4MaterialCutsCouple *	CurrentCouple () const
void	SetCurrentElement (const G4Element *)
Protected Attributes inherited from G4VMscModel
G4double	facrange
G4double	facgeom
G4double	facsafety
G4double	skin
G4double	dtrl
G4double	lambdalimit
G4double	geomMin
G4double	geomMax
G4ThreeVector	fDisplacement
G4MscStepLimitType	steppingAlgorithm
G4bool	samplez
G4bool	latDisplasment
Protected Attributes inherited from G4VEmModel
G4VParticleChange *	pParticleChange
G4PhysicsTable *	xSectionTable
const std::vector< G4double > *	theDensityFactor
const std::vector< G4int > *	theDensityIdx

Detailed Description

Definition at line 73 of file G4GoudsmitSaundersonMscModel.hh.

Constructor & Destructor Documentation

G4GoudsmitSaundersonMscModel()

G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel

(

const G4String &

nam = "GoudsmitSaunderson"

)

Definition at line 97 of file G4GoudsmitSaundersonMscModel.cc.

  : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV)
{ 
  currentKinEnergy=currentRange=skindepth=par1=par2=par3
    =zPathLength=truePathLength
    =tausmall=taulim=tlimit=charge=lambdalimit=tPathLength=lambda0=lambda1
    =lambda11=mass=0.0;
  currentMaterialIndex = -1;
 
  fr=0.02,rangeinit=0.,masslimite=0.6*MeV,
  particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
  tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
  tlimitminfix=1.e-6*mm;stepmin=tlimitminfix;lambdalimit=1.*mm;smallstep=1.e10;
  theManager=G4LossTableManager::Instance();
  inside=false;insideskin=false;
  samplez=false;
  firstStep = true; 
 
  GSTable = new G4GoudsmitSaundersonTable();
 
  if(ener[0] < 0.0){ 
    G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
    LoadELSEPAXSections();
  }
}

~G4GoudsmitSaundersonMscModel()

G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel ( )

virtual

Definition at line 123 of file G4GoudsmitSaundersonMscModel.cc.

{
  delete GSTable;
}

Member Function Documentation

ComputeCrossSectionPerAtom()

G4double G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  particle, G4double  KineticEnergy, G4double  AtomicNumber, G4double  , G4double  , G4double    )

virtual

Reimplemented from G4VEmModel.

Definition at line 139 of file G4GoudsmitSaundersonMscModel.cc.

{  
  G4double kinEnergy = kineticEnergy;
  if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
  if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
 
  G4double cs(0.0), cs0(0.0);
  CalculateIntegrals(p,Z,kinEnergy,cs0,cs);
  
  return cs;
}  

ComputeGeomPathLength()

G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength ( G4double  truePathLength )

virtual

Reimplemented from G4VMscModel.

Definition at line 646 of file G4GoudsmitSaundersonMscModel.cc.

{
  firstStep = false; 
  par1 = -1. ;  
  par2 = par3 = 0. ;  
 
  //  do the true -> geom transformation
  zPathLength = tPathLength;
 
  // z = t for very small tPathLength
  if(tPathLength < tlimitminfix) { return zPathLength; }
 
  // this correction needed to run MSC with eIoni and eBrem inactivated
  // and makes no harm for a normal run
  if(tPathLength > currentRange)
    { tPathLength = currentRange; }
 
  G4double tau   = tPathLength/lambda1 ;
 
  if ((tau <= tausmall) || insideskin) {
    zPathLength  = tPathLength;
    if(zPathLength > lambda1) { zPathLength = lambda1; }
    return zPathLength;
  }
 
  G4double zmean = tPathLength;
  if (tPathLength < currentRange*dtrl) {
    if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
    else             zmean = lambda1*(1.-exp(-tau));
  } else if(currentKinEnergy < mass || tPathLength == currentRange) {
    par1 = 1./currentRange ;
    par2 = 1./(par1*lambda1) ;
    par3 = 1.+par2 ;
    if(tPathLength < currentRange)
      zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
    else
      zmean = 1./(par1*par3) ;
  } else {
    G4double T1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
 
    lambda11 = GetTransportMeanFreePath(particle,T1);
 
    par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
    par2 = 1./(par1*lambda1) ;
    par3 = 1.+par2 ;
    zmean = (1.-exp(par3*log(lambda11/lambda1)))/(par1*par3) ;
  }
 
  zPathLength = zmean ;
  //  sample z
  if(samplez) {
 
    const G4double  ztmax = 0.99;
    G4double zt = zmean/tPathLength ;
 
    if (tPathLength > stepmin && zt < ztmax) {
 
      G4double u,cz1;
      if(zt >= 0.333333333) {
 
        G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
        cz1 = 1.+cz ;
        G4double u0 = cz/cz1 ;
        G4double grej ;
        do {
      u = exp(log(G4UniformRand())/cz1) ;
      grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
    } while (grej < G4UniformRand()) ;
 
      } else {
        cz1 = 1./zt-1.;
        u = 1.-exp(log(G4UniformRand())/cz1) ;
      }
      zPathLength = tPathLength*u ;
    }
  }
  if(zPathLength > lambda1) zPathLength = lambda1;
  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda1 << G4endl;
 
  return zPathLength;
}

ComputeTruePathLengthLimit()

G4double G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit ( const G4Track &  track, G4double &  currentMinimalStep  )

virtual

Reimplemented from G4VMscModel.

Definition at line 453 of file G4GoudsmitSaundersonMscModel.cc.

{
  tPathLength = currentMinimalStep;
  const G4DynamicParticle* dp = track.GetDynamicParticle();
  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
  G4StepStatus stepStatus = sp->GetStepStatus();
  currentCouple = track.GetMaterialCutsCouple();
  SetCurrentCouple(currentCouple); 
  currentMaterialIndex = currentCouple->GetIndex();
  currentKinEnergy = dp->GetKineticEnergy();
  currentRange = GetRange(particle,currentKinEnergy,currentCouple);
 
  lambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
 
  // stop here if small range particle
  if(inside || tPathLength < tlimitminfix) {
    return ConvertTrueToGeom(tPathLength, currentMinimalStep);
  }  
  if(tPathLength > currentRange) tPathLength = currentRange;
 
  G4double presafety = sp->GetSafety();
 
  //G4cout << "G4GS::StepLimit tPathLength= " 
  //     <<tPathLength<<" safety= " << presafety
  //       << " range= " <<currentRange<< " lambda= "<<lambda1
  //     << " Alg: " << steppingAlgorithm <<G4endl;
 
  // far from geometry boundary
  if(currentRange < presafety)
    {
      inside = true;
      return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
    }
 
  // standard  version
  //
  if (steppingAlgorithm == fUseDistanceToBoundary)
    {
      //compute geomlimit and presafety 
      G4double geomlimit = ComputeGeomLimit(track, presafety, tPathLength);
   
      // is far from boundary
      if(currentRange <= presafety)
    {
      inside = true;
      return ConvertTrueToGeom(tPathLength, currentMinimalStep);   
    }
 
      smallstep += 1.;
      insideskin = false;
 
      if(firstStep || stepStatus == fGeomBoundary)
    {
          rangeinit = currentRange;
          if(firstStep) smallstep = 1.e10;
          else  smallstep = 1.;
 
          //define stepmin here (it depends on lambda!)
          //rough estimation of lambda_elastic/lambda_transport
          G4double rat = currentKinEnergy/MeV ;
          rat = 1.e-3/(rat*(10.+rat)) ;
          //stepmin ~ lambda_elastic
          stepmin = rat*lambda1;
          skindepth = skin*stepmin;
          //define tlimitmin
          tlimitmin = 10.*stepmin;
          if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
 
      //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
      //     << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
      // constraint from the geometry 
      if((geomlimit < geombig) && (geomlimit > geommin))
        {
          if(stepStatus == fGeomBoundary)  
            tgeom = geomlimit/facgeom;
          else
            tgeom = 2.*geomlimit/facgeom;
        }
            else
              tgeom = geombig;
 
        }
 
      //step limit 
      tlimit = facrange*rangeinit;              
      if(tlimit < facsafety*presafety)
        tlimit = facsafety*presafety; 
 
      //lower limit for tlimit
      if(tlimit < tlimitmin) tlimit = tlimitmin;
 
      if(tlimit > tgeom) tlimit = tgeom;
 
      //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit  
      //      << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
 
      // shortcut
      if((tPathLength < tlimit) && (tPathLength < presafety) &&
         (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
    return ConvertTrueToGeom(tPathLength, currentMinimalStep);   
 
      // step reduction near to boundary
      if(smallstep < skin)
    {
      tlimit = stepmin;
      insideskin = true;
    }
      else if(geomlimit < geombig)
    {
      if(geomlimit > skindepth)
        {
          if(tlimit > geomlimit-0.999*skindepth)
        tlimit = geomlimit-0.999*skindepth;
        }
      else
        {
          insideskin = true;
          if(tlimit > stepmin) tlimit = stepmin;
        }
    }
 
      if(tlimit < stepmin) tlimit = stepmin;
 
      if(tPathLength > tlimit) tPathLength = tlimit; 
 
    }
    // for 'normal' simulation with or without magnetic field 
    //  there no small step/single scattering at boundaries
  else if(steppingAlgorithm == fUseSafety)
    {
      // compute presafety again if presafety <= 0 and no boundary
      // i.e. when it is needed for optimization purposes
      if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 
        presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
 
      // is far from boundary
      if(currentRange < presafety)
        {
          inside = true;
          return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
        }
 
      if(firstStep || stepStatus == fGeomBoundary)
    {
      rangeinit = currentRange;
      fr = facrange;
      // 9.1 like stepping for e+/e- only (not for muons,hadrons)
      if(mass < masslimite) 
        {
          if(lambda1 > currentRange)
        rangeinit = lambda1;
          if(lambda1 > lambdalimit)
        fr *= 0.75+0.25*lambda1/lambdalimit;
        }
 
      //lower limit for tlimit
      G4double rat = currentKinEnergy/MeV ;
      rat = 1.e-3/(rat*(10.+rat)) ;
      tlimitmin = 10.*lambda1*rat;
      if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
    }
      //step limit
      tlimit = fr*rangeinit;               
 
      if(tlimit < facsafety*presafety)
        tlimit = facsafety*presafety;
 
      //lower limit for tlimit
      if(tlimit < tlimitmin) tlimit = tlimitmin;
 
      if(tPathLength > tlimit) tPathLength = tlimit;
    }
  
  // version similar to 7.1 (needed for some experiments)
  else
    {
      if (stepStatus == fGeomBoundary)
    {
      if (currentRange > lambda1) tlimit = facrange*currentRange;
      else                        tlimit = facrange*lambda1;
 
      if(tlimit < tlimitmin) tlimit = tlimitmin;
      if(tPathLength > tlimit) tPathLength = tlimit;
    }
    }
  //G4cout << "tPathLength= " << tPathLength  
  // << " currentMinimalStep= " << currentMinimalStep << G4endl;
  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
}

ComputeTrueStepLength()

G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength ( G4double  geomStepLength )

virtual

Reimplemented from G4VMscModel.

Definition at line 731 of file G4GoudsmitSaundersonMscModel.cc.

{
  // step defined other than transportation 
  if(geomStepLength == zPathLength && tPathLength <= currentRange)
    return tPathLength;
 
  // t = z for very small step
  zPathLength = geomStepLength;
  tPathLength = geomStepLength;
  if(geomStepLength < tlimitminfix) return tPathLength;
  
  // recalculation
  if((geomStepLength > lambda1*tausmall) && !insideskin)
    {
      if(par1 <  0.)
    tPathLength = -lambda1*log(1.-geomStepLength/lambda1) ;
      else 
    {
      if(par1*par3*geomStepLength < 1.)
        tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
      else 
        tPathLength = currentRange;
    }  
    }
  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
 
  return tPathLength;
}

Initialise()

void G4GoudsmitSaundersonMscModel::Initialise ( const G4ParticleDefinition *  p, const G4DataVector &    )

virtual

Implements G4VEmModel.

Definition at line 128 of file G4GoudsmitSaundersonMscModel.cc.

{ 
  skindepth=skin*stepmin;
  SetParticle(p);
  fParticleChange = GetParticleChangeForMSC(p);
}

SampleScattering()

G4ThreeVector & G4GoudsmitSaundersonMscModel::SampleScattering ( const G4DynamicParticle *  dynParticle, G4double  safety  )

virtual

Reimplemented from G4VMscModel.

Definition at line 154 of file G4GoudsmitSaundersonMscModel.cc.

{
  fDisplacement.set(0.0,0.0,0.0);
  G4double kineticEnergy = dynParticle->GetKineticEnergy();
  if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)||
     (tPathLength/tausmall < lambda1)) { return fDisplacement; }
 
  ///////////////////////////////////////////
  // Effective energy 
  G4double eloss  = 0.0;
  if (tPathLength > currentRange*dtrl) {
    eloss = kineticEnergy - 
      GetEnergy(particle,currentRange-tPathLength,currentCouple);
  } else {
    eloss = tPathLength*GetDEDX(particle,kineticEnergy,currentCouple);
  }
  /*
  G4double ttau      = kineticEnergy/electron_mass_c2;
  G4double ttau2     = ttau*ttau;
  G4double epsilonpp = eloss/kineticEnergy;
  G4double cst1  = epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
  kineticEnergy *= (1 - cst1);
  */
  kineticEnergy -= 0.5*eloss;
 
  ///////////////////////////////////////////
  // additivity rule for mixture and compound xsection's
  const G4Material* mat = currentCouple->GetMaterial();
  const G4ElementVector* theElementVector = mat->GetElementVector();
  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
  G4int nelm = mat->GetNumberOfElements();
  G4double s0(0.0), s1(0.0);
  lambda0 = 0.0;
  for(G4int i=0;i<nelm;i++)
    { 
      CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,s0,s1);
      lambda0 += (theAtomNumDensityVector[i]*s0);
    } 
  if(lambda0>0.0) lambda0 =1./lambda0;
 
  // Newton-Raphson root's finding method of scrA from: 
  // Sig1(PWA)/Sig0(PWA)=g1=2*scrA*((1+scrA)*log(1+1/scrA)-1)
  G4double g1=0.0;
  if(lambda1>0.0) { g1 = lambda0/lambda1; }
 
  G4double logx0,x1,delta;
  G4double x0=g1*0.5;
  // V.Ivanchenko added limit of the loop
  for(G4int i=0;i<1000;++i)
    {  
      logx0=std::log(1.+1./x0);
      x1 = x0-(x0*((1.+x0)*logx0-1.0)-g1*0.5)/( (1.+2.*x0)*logx0-2.0);
 
      // V.Ivanchenko cut step size of iterative procedure
      if(x1 < 0.0)         { x1 = 0.5*x0; }
      else if(x1 > 2*x0)   { x1 = 2*x0; }
      else if(x1 < 0.5*x0) { x1 = 0.5*x0; }
      delta = std::fabs( x1 - x0 );    
      x0 = x1;
      if(delta < 1.0e-3*x1) { break;}
    }
  G4double scrA = x1;
 
  G4double lambdan=0.;
 
  if(lambda0>0.0) { lambdan=tPathLength/lambda0; }
  if(lambdan<=1.0e-12) { return fDisplacement; }
 
  //G4cout << "E(eV)= " << kineticEnergy/eV << " L0= " << lambda0
  //    << " L1= " << lambda1 << G4endl;
 
  G4double Qn1 = lambdan *g1;//2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
  G4double Qn12 = 0.5*Qn1;
  
  G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
  G4double cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
  G4double us=0.0,vs=0.0,ws=1.0,wss=0.,x_coord=0.0,y_coord=0.0,z_coord=1.0;
 
  G4double epsilon1=G4UniformRand();
  G4double expn = std::exp(-lambdan);
 
  if(epsilon1<expn)// no scattering 
    { return fDisplacement; }
  else if((epsilon1<((1.+lambdan)*expn))||(lambdan<1.))//single or plural scattering (Rutherford DCS's)
    {
      G4double xi=G4UniformRand();
      xi= 2.*scrA*xi/(1.-xi + scrA);
      if(xi<0.)xi=0.;
      else if(xi>2.)xi=2.;      
      ws=(1. - xi);
      wss=std::sqrt(xi*(2.-xi));      
      G4double phi0=CLHEP::twopi*G4UniformRand(); 
      us=wss*cos(phi0);
      vs=wss*sin(phi0);
    }
  else // multiple scattering
    {
      // Ref.2 subsection 4.4 "The best solution found"
      // Sample first substep scattering angle
      SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
      G4double phi1  = CLHEP::twopi*G4UniformRand();
      cosPhi1 = cos(phi1);
      sinPhi1 = sin(phi1);
 
      // Sample second substep scattering angle
      SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
      G4double phi2  = CLHEP::twopi*G4UniformRand();
      cosPhi2 = cos(phi2);
      sinPhi2 = sin(phi2);
 
      // Overall scattering direction
      us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
      vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
      ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2; 
      G4double sqrtA=sqrt(scrA);
      if(acos(ws)<sqrtA)//small angle approximation for theta less than screening angle
      {
    G4int i=0;
    do{i++;
      ws=1.+Qn12*log(G4UniformRand());
    }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
    if(i>=19)ws=cos(sqrtA);
    wss=std::sqrt((1.-ws*ws));      
    us=wss*std::cos(phi1);
    vs=wss*std::sin(phi1);
      }
    }
    
  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
  G4ThreeVector newDirection(us,vs,ws);
  newDirection.rotateUz(oldDirection);
  fParticleChange->ProposeMomentumDirection(newDirection);
 
  // corresponding to error less than 1% in the exact formula of <z>
  if(Qn1<0.02) { z_coord = 1.0 - Qn1*(0.5 - Qn1/6.); }
  else         { z_coord = (1.-std::exp(-Qn1))/Qn1; }
  G4double rr = zPathLength*std::sqrt((1.- z_coord*z_coord)/(1.-ws*ws));
  x_coord  = rr*us;
  y_coord  = rr*vs;
 
  // displacement is computed relatively to the end point
  z_coord -= 1.0;
  z_coord *= zPathLength;
  /*
    G4cout << "G4GS::SampleSecondaries: e(MeV)= " << kineticEnergy
    << " sinTheta= " << sqrt(1.0 - ws*ws) 
    << " trueStep(mm)= " << tPathLength
    << " geomStep(mm)= " << zPathLength
    << G4endl;
  */
 
  fDisplacement.set(x_coord,y_coord,z_coord);
  fDisplacement.rotateUz(oldDirection);
 
  return fDisplacement;
}     

StartTracking()

void G4GoudsmitSaundersonMscModel::StartTracking ( G4Track *  track )

virtual

Reimplemented from G4VEmModel.

Definition at line 440 of file G4GoudsmitSaundersonMscModel.cc.

{
  SetParticle(track->GetDynamicParticle()->GetDefinition());
  firstStep = true; 
  inside = false;
  insideskin = false;
  tlimit = geombig;
}

---

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

• G4GoudsmitSaundersonMscModel.hh
• G4GoudsmitSaundersonMscModel.cc