Geant4  10.02.p03
G4GoudsmitSaundersonMscModel Class Reference

#include <G4GoudsmitSaundersonMscModel.hh>

Inheritance diagram for G4GoudsmitSaundersonMscModel:
Collaboration diagram for G4GoudsmitSaundersonMscModel:

Public Member Functions

 G4GoudsmitSaundersonMscModel (const G4String &nam="GoudsmitSaunderson")
 
virtual ~G4GoudsmitSaundersonMscModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
void StartTracking (G4Track *)
 
G4double GetTransportMeanFreePath (const G4ParticleDefinition *, G4double)
 
void SingleScattering (G4double &cost, G4double &sint)
 
void SampleMSC ()
 
virtual G4ThreeVectorSampleScattering (const G4ThreeVector &, G4double safety)
 
virtual G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep)
 
virtual G4double ComputeGeomPathLength (G4double truePathLength)
 
virtual G4double ComputeTrueStepLength (G4double geomStepLength)
 
void SetOptionPWAScreening (G4bool opt)
 
- Public Member Functions inherited from G4VMscModel
 G4VMscModel (const G4String &nam)
 
virtual ~G4VMscModel ()
 
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)
 
G4VEnergyLossProcessGetIonisation () const
 
void SetIonisation (G4VEnergyLossProcess *, const G4ParticleDefinition *part)
 
G4double ComputeSafety (const G4ThreeVector &position, G4double limit=DBL_MAX)
 
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 InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
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 GetPartialCrossSection (const G4Material *, G4int, const G4ParticleDefinition *, G4double)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, 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 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 *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector *> *)
 
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 G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
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
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
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 SetForceBuildTable (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
const G4IsotopeGetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Private Member Functions

void SetParticle (const G4ParticleDefinition *p)
 
G4double GetLambda (G4double)
 
G4GoudsmitSaundersonMscModeloperator= (const G4GoudsmitSaundersonMscModel &right)
 
 G4GoudsmitSaundersonMscModel (const G4GoudsmitSaundersonMscModel &)
 
G4double GetTransportMeanFreePathOnly (const G4ParticleDefinition *, G4double)
 
G4double Randomizetlimit ()
 

Private Attributes

CLHEP::HepRandomEnginerndmEngineMod
 
G4double lowKEnergy
 
G4double highKEnergy
 
G4double currentKinEnergy
 
G4double currentRange
 
G4double fr
 
G4double rangeinit
 
G4double geombig
 
G4double geomlimit
 
G4double lambdalimit
 
G4double tlimit
 
G4double tgeom
 
G4int charge
 
G4int currentMaterialIndex
 
G4bool firstStep
 
G4double par1
 
G4double par2
 
G4double par3
 
G4double tlimitminfix2
 
G4double tausmall
 
G4double mass
 
G4double taulim
 
G4LossTableManagertheManager
 
const G4ParticleDefinitionparticle
 
G4ParticleChangeForMSC * fParticleChange
 
const G4MaterialCutsCouplecurrentCouple
 
G4bool fIsUsePWATotalXsecData
 
G4double presafety
 
G4double fZeff
 
G4double fLambda0
 
G4double fLambda1
 
G4double fScrA
 
G4double fG1
 
G4double fTheTrueStepLenght
 
G4double fTheTransportDistance
 
G4double fTheZPathLenght
 
G4ThreeVector fTheDisplacementVector
 
G4ThreeVector fTheNewDirection
 
G4bool fIsEndedUpOnBoundary
 
G4bool fIsMultipleSacettring
 
G4bool fIsSingleScattering
 
G4bool fIsEverythingWasDone
 
G4bool fIsNoScatteringInMSC
 
G4bool fIsNoDisplace
 
G4bool fIsInsideSkin
 
G4bool fIsWasOnBoundary
 
G4bool fIsFirstRealStep
 

Static Private Attributes

static G4GoudsmitSaundersonTablefgGSTable = 0
 
static G4PWATotalXsecTablefgPWAXsecTable = 0
 
static G4bool fgIsUseAccurate = true
 
static G4bool fgIsOptimizationOn = true
 

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 G4MaterialCutsCoupleCurrentCouple () 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
G4ElementDatafElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 111 of file G4GoudsmitSaundersonMscModel.hh.

Constructor & Destructor Documentation

◆ G4GoudsmitSaundersonMscModel() [1/2]

G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel ( const G4String nam = "GoudsmitSaunderson")

Definition at line 125 of file G4GoudsmitSaundersonMscModel.cc.

126  : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(1.*GeV) {
127  charge = 0;
129 
130  lambdalimit = 1.*mm ;
131  fr = 0.02 ;
132  rangeinit = 0.;
133  geombig = 1.e50*mm;
134  geomlimit = geombig;
135  tgeom = geombig;
136  tlimit = 1.e+10*mm;
137  presafety = 0.*mm;
138 
139  particle = 0;
141  firstStep = true;
142  currentKinEnergy = 0.0;
143  currentRange = 0.0;
144 
145  tlimitminfix2 = 1.*nm;
146  par1=par2=par3= 0.0;
147  tausmall = 1.e-16;
149  taulim = 1.e-6;
150 
151  currentCouple = 0;
152  fParticleChange = 0;
153 
154  // by default Moliere screeing parameter will be used but it can be set to the
155  // PWA screening before calling the Initialise method
156  fIsUsePWATotalXsecData = false;
157 
158  //*****************
159  fLambda0 = 0.0; // elastic mean free path
160  fLambda1 = 0.0; // first transport mean free path
161  fScrA = 0.0; // screening parameter
162  fG1 = 0.0; // first transport coef.
163  //******************
164  fTheTrueStepLenght = 0.;
166  fTheZPathLenght = 0.;
167  fTheDisplacementVector.set(0.,0.,0.);
168  fTheNewDirection.set(0.,0.,1.);
169 
170  fIsEverythingWasDone = false;
171  fIsMultipleSacettring = false;
172  fIsSingleScattering = false;
173  fIsEndedUpOnBoundary = false;
174  fIsNoScatteringInMSC = false;
175  fIsNoDisplace = false;
176  fIsInsideSkin = false;
177  fIsWasOnBoundary = false;
178  fIsFirstRealStep = false;
179 
180  fZeff = 1.;
181  //+++++++++++++
182 
183  rndmEngineMod = G4Random::getTheEngine();
184 
185  //=================== base
186  facsafety = 0.6;
187 }
void set(double x, double y, double z)
static G4LossTableManager * Instance()
const G4ParticleDefinition * particle
static const double nm
Definition: G4SIunits.hh:111
static const double GeV
Definition: G4SIunits.hh:214
float electron_mass_c2
Definition: hepunit.py:274
G4double facsafety
Definition: G4VMscModel.hh:179
static const double keV
Definition: G4SIunits.hh:213
G4VMscModel(const G4String &nam)
Definition: G4VMscModel.cc:58
static const double mm
Definition: G4SIunits.hh:114
const G4MaterialCutsCouple * currentCouple
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ~G4GoudsmitSaundersonMscModel()

G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel ( )
virtual

Definition at line 190 of file G4GoudsmitSaundersonMscModel.cc.

190  {
191  if(IsMaster()){
192  if(fgGSTable){
193  delete fgGSTable;
194  fgGSTable = 0;
195  }
196  if(fgPWAXsecTable){
197  delete fgPWAXsecTable;
198  fgPWAXsecTable = 0;
199  }
200  }
201 }
static G4PWATotalXsecTable * fgPWAXsecTable
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
static G4GoudsmitSaundersonTable * fgGSTable
Here is the call graph for this function:

◆ G4GoudsmitSaundersonMscModel() [2/2]

G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel ( const G4GoudsmitSaundersonMscModel )
private

Member Function Documentation

◆ ComputeGeomPathLength()

G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength ( G4double  truePathLength)
virtual

Reimplemented from G4VMscModel.

Definition at line 702 of file G4GoudsmitSaundersonMscModel.cc.

703 {
704  //
705  // convert true ->geom
706  // It is called from the step limitation ComputeTruePathLengthLimit if
707  // !fIsEverythingWasDone but protect:
708  //
709  par1 = -1. ;
710  par2 = par3 = 0. ;
711 
712  // if fIsEverythingWasDone = TRUE => fTheZPathLenght is already set
713  // so return with the already known value
714  // Otherwise:
715  if(!fIsEverythingWasDone) {
716  // this correction needed to run MSC with eIoni and eBrem inactivated
717  // and makes no harm for a normal run
719 
720  // do the true -> geom transformation
722 
723  // z = t for very small true-path-length
725 
727 
728  if (tau <= tausmall) {
730  } else if (fTheTrueStepLenght < currentRange*dtrl) {
731  if(tau < taulim) fTheZPathLenght = fTheTrueStepLenght*(1.-0.5*tau) ;
732  else fTheZPathLenght = fLambda1*(1.-G4Exp(-tau));
734  par1 = 1./currentRange ; // alpha =1/range_init for Ekin<mass
735  par2 = 1./(par1*fLambda1) ; // 1/(alphaxlambda01)
736  par3 = 1.+par2 ; // 1+1/
738  fTheZPathLenght = 1./(par1*par3) * (1.-std::pow(1.-par1*fTheTrueStepLenght,par3));
739  } else {
740  fTheZPathLenght = 1./(par1*par3);
741  }
742 
743  } else {
747 
748  par1 = (fLambda1-lambda1)/(fLambda1*fTheTrueStepLenght); // alpha
749  par2 = 1./(par1*fLambda1);
750  par3 = 1.+par2 ;
751  G4Pow *g4pow = G4Pow::GetInstance();
752  fTheZPathLenght = 1./(par1*par3) * (1.-g4pow->powA(1.-par1*fTheTrueStepLenght,par3));
753  }
754  }
756 
757  return fTheZPathLenght;
758 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double dtrl
Definition: G4VMscModel.hh:181
Definition: G4Pow.hh:56
G4double GetTransportMeanFreePathOnly(const G4ParticleDefinition *, G4double)
const G4ParticleDefinition * particle
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:315
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
double G4double
Definition: G4Types.hh:76
const G4MaterialCutsCouple * currentCouple
Here is the call graph for this function:

◆ ComputeTruePathLengthLimit()

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

Reimplemented from G4VMscModel.

Definition at line 374 of file G4GoudsmitSaundersonMscModel.cc.

376  {
377 
378  G4double skindepth = 0.;
379 
380  const G4DynamicParticle* dp = track.GetDynamicParticle();
381 
382  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
383  G4StepStatus stepStatus = sp->GetStepStatus();
384  currentCouple = track.GetMaterialCutsCouple();
389 
390 
391 
392  // *** Elastic and first transport mfp, fScrA and fG1 are also set in this !!!
394 
395  //
396  // Set initial values:
397  // : lengths are initialised to currentMinimalStep which is the true, minimum
398  // step length from all other physics
399  fTheTrueStepLenght = currentMinimalStep;
400  fTheTransportDistance = currentMinimalStep;
401  fTheZPathLenght = currentMinimalStep; // will need to be converted
402  fTheDisplacementVector.set(0.,0.,0.);
403  fTheNewDirection.set(0.,0.,1.);
404 
405  // Can everything be done in the step limit phase ?
406  fIsEverythingWasDone = false;
407  // Multiple scattering needs to be sample ?
408  fIsMultipleSacettring = false;
409  // Single scattering needs to be sample ?
410  fIsSingleScattering = false;
411  // Was zero deflection in multiple scattering sampling ?
412  fIsNoScatteringInMSC = false;
413  // Do not care about displacement in MSC sampling
414  // ( used only in the case of fgIsOptimizationOn = true)
415  fIsNoDisplace = false;
416 
417  // get pre-step point safety
418  presafety = sp->GetSafety();
419 
420  // Zeff = TotNbOfElectPerVolume/TotNbOfAtomsPerVolume
422 
423  // distance will take into account max-fluct.
424  G4double distance = currentRange;
425  distance *= (1.20-fZeff*(1.62e-2-9.22e-5*fZeff));
426 
427 
428 
429  //
430  // Possible optimization : if the distance is samller than the safety -> the
431  // particle will never leave this volume -> dispalcement
432  // as the effect of multiple elastic scattering can be skipped
433  // Important : this optimization can cause problems if one does scoring
434  // in a bigger volume since MSC won't be done deep inside the volume when
435  // distance < safety so don't use optimized-mode in such case.
436  if( fgIsOptimizationOn && (distance < presafety) ) {
437  // Indicate that we need to do MSC after transportation and no dispalcement.
438  fIsMultipleSacettring = true;
439  fIsNoDisplace = true;
441  //Compute geomlimit (and presafety) :
442  // - geomlimit will be:
443  // == the straight line distance to the boundary if currentRange is
444  // longer than that
445  // == a big value [geombig = 1.e50*mm] if currentRange is shorter than
446  // the straight line distance to the boundary
447  // - presafety will be updated as well
448  // So the particle can travell 'gemlimit' distance (along a straight
449  // line!) in its current direction:
450  // (1) before reaching a boundary (geomlimit < geombig) OR
451  // (2) before reaching its current range (geomlimit == geombig)
453 
454  // Record that the particle is on a boundary
455  if( (stepStatus==fGeomBoundary) || (stepStatus==fUndefined && presafety==0.0))
456  fIsWasOnBoundary = true;
457 
458  // Set skin depth = skin x elastic_mean_free_path
459  skindepth = skin*fLambda0;
460  // Init the flag that indicates that the particle are within a skindepth
461  // distance from a boundary
462  fIsInsideSkin = false;
463  // Check if we can try Single Scattering because we are within skindepth
464  // distance from/to a boundary OR the current minimum true-step-length is
465  // shorter than skindepth. NOTICE: the latest has only efficieny reasons
466  // because the MSC angular sampling is fine for any short steps but much
467  // faster to try single scattering in case of short steps.
468  if( (stepStatus == fGeomBoundary) || (presafety < skindepth) || (fTheTrueStepLenght < skindepth)){
469  // check if we are within skindepth distance from a boundary
470  if( (stepStatus == fGeomBoundary) || (presafety < skindepth))
471  fIsInsideSkin = true;
472  //Try single scattering:
473  // - sample distance to next single scattering interaction (sslimit)
474  // - compare to current minimum length
475  // == if sslimit is the shorter:
476  // - set the step length to sslimit
477  // - indicate that single scattering needs to be done
478  // == else : nothing to do
479  //- in both cases, the step length was very short so geometrical and
480  // true path length are the same
481  G4double sslimit = -1.*fLambda0*G4Log(rndmEngineMod->flat());
482  // compare to current minimum step length
483  if(sslimit < fTheTrueStepLenght) {
484  fTheTrueStepLenght = sslimit;
485  fIsSingleScattering = true;
486  }
487 
488  // short step -> true step length equal to geometrical path length
490  // Set taht everything is done in step-limit phase so no MSC call
491  // We will check if we need to perform the single-scattering angular
492  // sampling i.e. if single elastic scattering was the winer!
493  fIsEverythingWasDone = true;
494  } else {
495  // After checking we know that we cannot try single scattering so we will
496  // need to make an MSC step
497  // Indicate that we need to make and MSC step. We do not check if we can
498  // do it now i.e. if presafety>final_true_step_length so we let the
499  // fIsEverythingWasDone = false which indicates that we will perform
500  // MSC after transportation.
501  fIsMultipleSacettring = true;
502 
503  // Init the first-real-step falg: it will indicate if we do the first
504  // non-single scattering step in this volume with this particle
505  fIsFirstRealStep = false;
506  // If previously the partcile was on boundary it was within skin as
507  // well. When it is not within skin anymore it has just left the skin
508  // so we make the first real MSC step with the particle.
510  // reset the 'was on boundary' indicator flag
511  fIsWasOnBoundary = false;
512  fIsFirstRealStep = true;
513  }
514 
515  // If this is the first-real msc step (the partcile has just left the
516  // skin) or this is the first step with the particle (was born or
517  // primary):
518  // - set the initial range that will be used later to limit its step
519  // (only in this volume, because after boundary crossing at the
520  // first-real MSC step we will reset)
521  // - don't let the partcile to cross the volume just in one step
524  // If geomlimit < geombig than the particle might reach the boundary
525  // along its initial direction before losing its energy (in this step)
526  // Otherwise we can be sure that the particle will lose it energy
527  // before reaching the boundary along a starigth line so there is no
528  // geometrical limit appalied. [However, tgeom is set only in the
529  // first or the first-real MSC step. After the first or first real
530  // MSC step the direction will change tgeom won't guaranty anything!
531  // But we will try to end up within skindepth from the boundary using
532  // the actual value of geomlimit(See later at step reduction close to
533  // boundary).]
534  if(geomlimit < geombig){
535  // transfrom straight line distance to the boundary to real step
536  // length based on the mean values (using the prestep point
537  // first-transport mean free path i.e. no energy loss correction)
538  if((1.-geomlimit/fLambda1) > 0.)
540  // the 2-different case that could lead us here
541  if(firstStep)
542  tgeom = 2.*geomlimit/facgeom;
543  else
545  } else {
546  tgeom = geombig;
547  }
548  }
549 
550  // True step length limit from range factor. Noteice, that the initial
551  // range is used that was set at the first step or first-real MSC step
552  // in this volume with this particle.
554  // Take the minimum of the true step length limits coming from
555  // geometrical constraint or range-factor limitation
557 
558  // Step reduction close to boundary: we try to end up within skindepth
559  // from the boundary ( Notice: in case of mag. field it might not work
560  // because geomlimit is the straigth line distance to the boundary in
561  // the currect direction (if geomlimit<geombig) and mag. field can
562  // change the initial direction. So te particle might hit some boundary
563  // before in a different direction. However, here we restrict the true
564  // path length to this (straight line) lenght so the corresponding
565  // transport distance (straight line) will be even shorter than
566  // geomlimit-0.999*skindepth after the change of true->geom.
567  if(geomlimit < geombig)
568  tlimit = std::min(tlimit, geomlimit-0.999*skindepth);
569 
570  // randomize 1st step or 1st 'normal' step in volume
571  if(firstStep || fIsFirstRealStep) {
573  } else {
575  }
576  }
577  } else if (steppingAlgorithm == fUseSafety) { // THE ERROR_FREE stepping alg.
578  presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
580 
581  // Set skin depth = skin x elastic_mean_free_path
582  skindepth = skin*fLambda0;
583  // Check if we can try Single Scattering because we are within skindepth
584  // distance from/to a boundary OR the current minimum true-step-length is
585  // shorter than skindepth. NOTICE: the latest has only efficieny reasons
586  // because the MSC angular sampling is fine for any short steps but much
587  // faster to try single scattering in case of short steps.
588  if( (stepStatus == fGeomBoundary) || (presafety < skindepth) || (fTheTrueStepLenght < skindepth)){
589  //Try single scattering:
590  // - sample distance to next single scattering interaction (sslimit)
591  // - compare to current minimum length
592  // == if sslimit is the shorter:
593  // - set the step length to sslimit
594  // - indicate that single scattering needs to be done
595  // == else : nothing to do
596  //- in both cases, the step length was very short so geometrical and
597  // true path length are the same
598  G4double sslimit = -1.*fLambda0*G4Log(rndmEngineMod->flat());
599  // compare to current minimum step length
600  if(sslimit < fTheTrueStepLenght) {
601  fTheTrueStepLenght = sslimit;
602  fIsSingleScattering = true;
603  }
604 
605  // short step -> true step length equal to geometrical path length
607  // Set taht everything is done in step-limit phase so no MSC call
608  // We will check if we need to perform the single-scattering angular
609  // sampling i.e. if single elastic scattering was the winer!
610  fIsEverythingWasDone = true;
611  } else {
612  // After checking we know that we cannot try single scattering so we will
613  // need to make an MSC step
614  // Indicate that we need to make and MSC step.
615  fIsMultipleSacettring = true;
616  fIsEverythingWasDone = true;
617 
618  // limit from range factor
620 
621  // never let the particle go further than the safety if we are out of the skin
622  // if we are here we are out of the skin, presafety > 0.
625 
626  // make sure that we are still within the aplicability of condensed histry model
627  // i.e. true step length is not longer than first transport mean free path.
628  // We schould take into account energy loss along 0.5x lambda_transport1
629  // step length as well. So let it 0.4 x lambda_transport1
631  }
632  } else {
633 
634  // This is the default stepping algorithm: the fastest but the least
635  // accurate that corresponds to fUseSafety in Urban model. Note, that GS
636  // model can handle any short steps so we do not need the minimum limits
637 
638  // NO single scattering in case of skin or short steps (by defult the MSC
639  // model will be single or even no scattering in case of short steps
640  // compared to the elastic mean free path.)
641 
642  // indicate that MSC needs to be done (always and always after transportation)
643  fIsMultipleSacettring = true;
644 
645  if( stepStatus != fGeomBoundary ) {
646  presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
647  }
648 
649  // Far from boundary-> in optimized mode do not sample dispalcement.
650  if( (distance < presafety) && (fgIsOptimizationOn)) {
651  fIsNoDisplace = true;
652  } else {
653  // Urban like
654  if(firstStep || (stepStatus == fGeomBoundary)) {
655  rangeinit = currentRange;
656  fr = facrange;
657 // We don't use this: we won't converge to the single scattering results with
658 // decreasing range-factor.
659 // rangeinit = std::max(rangeinit, fLambda1);
660 // if(fLambda1 > lambdalimit) {
661 // fr *= (0.75+0.25*fLambda1/lambdalimit);
662 // }
663 
664  }
665 
666  //step limit
667  tlimit = std::max(fr*rangeinit, facsafety*presafety);
668 
669  // first step randomization
670  if(firstStep || stepStatus == fGeomBoundary) {
672  } else {
674  }
675  }
676  }
677 
678  // unset first-step
679  firstStep =false;
680 
681  // performe single scattering, multiple scattering if this later can be done safely here
684  // sample single scattering
685  G4double cost,sint,cosPhi,sinPhi;
686  SingleScattering(cost, sint);
688  sinPhi = std::sin(phi);
689  cosPhi = std::cos(phi);
690  fTheNewDirection.set(sint*cosPhi,sint*sinPhi,cost);
691  } else if(fIsMultipleSacettring){
692  // sample multiple scattering
693  SampleMSC(); // fTheZPathLenght, fTheDisplacementVector and fTheNewDirection will be set
694  } // and if single scattering but it was longer => nothing to do
695  } //else { do nothing here but after transportation
696 
697  return ConvertTrueToGeom(fTheTrueStepLenght,currentMinimalStep);
698 }
void set(double x, double y, double z)
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
G4double facgeom
Definition: G4VMscModel.hh:178
const G4Material * GetMaterial() const
G4double facrange
Definition: G4VMscModel.hh:177
G4double GetZeffective() const
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:247
G4double skin
Definition: G4VMscModel.hh:180
const G4ParticleDefinition * particle
virtual double flat()=0
G4double ComputeSafety(const G4ThreeVector &position, G4double limit=DBL_MAX)
Definition: G4VMscModel.hh:239
G4double GetKineticEnergy() const
void SingleScattering(G4double &cost, G4double &sint)
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:295
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:257
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double facsafety
Definition: G4VMscModel.hh:179
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:445
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:187
G4double GetTransportMeanFreePath(const G4ParticleDefinition *, G4double)
double G4double
Definition: G4Types.hh:76
static const double twopi
Definition: SystemOfUnits.h:54
const G4MaterialCutsCouple * currentCouple
Here is the call graph for this function:

◆ ComputeTrueStepLength()

G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength ( G4double  geomStepLength)
virtual

Reimplemented from G4VMscModel.

Definition at line 761 of file G4GoudsmitSaundersonMscModel.cc.

762 {
763  // init
764  fIsEndedUpOnBoundary = false;
765 
766  // step defined other than transportation
767  if(geomStepLength == fTheZPathLenght) {
768  return fTheTrueStepLenght;
769  }
770 
771  // else ::
772  // - set the flag that transportation was the winer so DoNothin in DOIT !!
773  // - convert geom -> true by using the mean value
774  fIsEndedUpOnBoundary = true; // OR LAST STEP
775 
776  fTheZPathLenght = geomStepLength;
777 
778  // was a short single scattering step
780  fTheTrueStepLenght = geomStepLength;
781  return fTheTrueStepLenght;
782  }
783 
784  // t = z for very small step
785  if(geomStepLength < tlimitminfix2) {
786  fTheTrueStepLenght = geomStepLength;
787  // recalculation
788  } else {
789  G4double tlength = geomStepLength;
790  if(geomStepLength > fLambda1*tausmall ) {
791  if(par1 < 0.) {
792  tlength = -fLambda1*G4Log(1.-geomStepLength/fLambda1) ;
793  } else {
794  if(par1*par3*geomStepLength < 1.) {
795  G4Pow *g4pow = G4Pow::GetInstance();
796  tlength = (1.-g4pow->powA( 1.-par1*par3*geomStepLength,1./par3))/par1;
797  } else {
798  tlength = currentRange;
799  }
800  }
801  if(tlength < geomStepLength) { tlength = geomStepLength; }
802  else if(tlength > fTheTrueStepLenght) { tlength = geomStepLength; }
803  }
804  fTheTrueStepLenght = tlength;
805  }
806  return fTheTrueStepLenght;
807 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
Definition: G4Pow.hh:56
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ GetLambda()

G4double G4GoudsmitSaundersonMscModel::GetLambda ( G4double  )
inlineprivate
Here is the caller graph for this function:

◆ GetTransportMeanFreePath()

G4double G4GoudsmitSaundersonMscModel::GetTransportMeanFreePath ( const G4ParticleDefinition ,
G4double  kineticEnergy 
)

Definition at line 246 of file G4GoudsmitSaundersonMscModel.cc.

247  {
248  // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
249  G4double efEnergy = kineticEnergy;
250  if(efEnergy<lowKEnergy) efEnergy=lowKEnergy;
251  if(efEnergy>highKEnergy)efEnergy=highKEnergy;
252 
255 
256  fLambda0 = 0.0; // elastic mean free path
257  fLambda1 = 0.0; // first transport mean free path
258  fScrA = 0.0; // screening parameter
259  fG1 = 0.0; // first transport coef.
260 
262  // use PWA total xsec data to determine screening and lambda0 lambda1
263  const G4ElementVector* theElementVector = mat->GetElementVector();
264  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
265  G4int nelm = mat->GetNumberOfElements();
266 
267  int elowindx = fgPWAXsecTable->GetPWATotalXsecForZet( G4lrint((*theElementVector)[0]->GetZ()) )
268  ->GetPWATotalXsecEnergyBinIndex(efEnergy);
269  for(G4int i=0;i<nelm;i++){
270  int zet = G4lrint((*theElementVector)[i]->GetZ());
271  // total elastic scattering cross section in Geant4 internal length2 units
272  int indx = (G4int)(1.5 + charge*1.5); // specify that want to get the total elastic scattering xsection
273  double elxsec = (fgPWAXsecTable->GetPWATotalXsecForZet(zet))->GetInterpXsec(efEnergy,elowindx,indx);
274  // first transport cross section in Geant4 internal length2 units
275  indx = (G4int)(2.5 + charge*1.5); // specify that want to get the first tarnsport xsection
276  double t1xsec = (fgPWAXsecTable->GetPWATotalXsecForZet(zet))->GetInterpXsec(efEnergy,elowindx,indx);
277  fLambda0 += (theAtomNumDensityVector[i]*elxsec);
278  fLambda1 += (theAtomNumDensityVector[i]*t1xsec);
279  }
280  if(fLambda0>0.0) { fLambda0 =1./fLambda0; }
281  if(fLambda1>0.0) { fLambda1 =1./fLambda1; }
282 
283  // determine screening parameter such that it gives back PWA G1
284  fG1=0.0;
285  if(fLambda1>0.0) { fG1 = fLambda0/fLambda1; }
286 
288  } else {
289  // Get SCREENING FROM MOLIER
290 
291  // below 1 keV it can give bananas so prevet (1 keV)
292  if(efEnergy<0.001) efEnergy = 0.001;
293 
294  // total mometum square in Geant4 internal energy2 units which is MeV2
295  G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2);
296  G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2);
297  G4int matindx = mat->GetIndex();
298  G4double bc = fgGSTable->GetMoliereBc(matindx);
299  fScrA = fgGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc);
300  // total elastic mean free path in Geant4 internal lenght units
301  fLambda0 = beta2/bc;//*(1.+fScrA);
302  fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/fScrA+1.0)-1.0);
304  }
305 
306  return fLambda1;
307 }
std::vector< G4Element * > G4ElementVector
const G4Material * GetMaterial() const
static G4PWATotalXsecTable * fgPWAXsecTable
size_t GetIndex() const
Definition: G4Material.hh:262
int G4int
Definition: G4Types.hh:78
Float_t mat
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
float electron_mass_c2
Definition: hepunit.py:274
static G4GoudsmitSaundersonTable * fgGSTable
const G4PWATotalXsecZ * GetPWATotalXsecForZet(G4int Z) const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4int GetPWATotalXsecEnergyBinIndex(G4double energy) const
G4double GetMoliereXc2(G4int matindx)
int G4lrint(double ad)
Definition: templates.hh:163
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4double GetMoliereBc(G4int matindx)
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
double G4double
Definition: G4Types.hh:76
const G4MaterialCutsCouple * currentCouple
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetTransportMeanFreePathOnly()

G4double G4GoudsmitSaundersonMscModel::GetTransportMeanFreePathOnly ( const G4ParticleDefinition ,
G4double  kineticEnergy 
)
private

Definition at line 311 of file G4GoudsmitSaundersonMscModel.cc.

312  {
313  // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
314  G4double efEnergy = kineticEnergy;
315  if(efEnergy<lowKEnergy) efEnergy=lowKEnergy;
316  if(efEnergy>highKEnergy)efEnergy=highKEnergy;
317 
320 
321  G4double lambda0 = 0.0; // elastc mean free path
322  G4double lambda1 = 0.0; // first transport mean free path
323  G4double scrA = 0.0; // screening parametr
324  G4double g1 = 0.0; // first transport mean free path
325 
327  // use PWA total xsec data to determine screening and lambda0 lambda1
328  const G4ElementVector* theElementVector = mat->GetElementVector();
329  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
330  G4int nelm = mat->GetNumberOfElements();
331 
332  int elowindx = fgPWAXsecTable->GetPWATotalXsecForZet( G4lrint((*theElementVector)[0]->GetZ()) )
333  ->GetPWATotalXsecEnergyBinIndex(efEnergy);
334  for(G4int i=0;i<nelm;i++){
335  int zet = G4lrint((*theElementVector)[i]->GetZ());
336  // first transport cross section in Geant4 internal length2 units
337  int indx = (G4int)(2.5 + charge*1.5); // specify that want to get the first tarnsport xsection
338  double t1xsec = (fgPWAXsecTable->GetPWATotalXsecForZet(zet))->GetInterpXsec(efEnergy,elowindx,indx);
339  lambda1 += (theAtomNumDensityVector[i]*t1xsec);
340  }
341  if(lambda1>0.0) { lambda1 =1./lambda1; }
342 
343  } else {
344  // Get SCREENING FROM MOLIER
345 
346  // below 1 keV it can give bananas so prevet (1 keV)
347  if(efEnergy<0.001) efEnergy = 0.001;
348 
349  // total mometum square in Geant4 internal energy2 units which is MeV2
350  G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2);
351  G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2);
352  G4int matindx = mat->GetIndex();
353  G4double bc = fgGSTable->GetMoliereBc(matindx);
354  scrA = fgGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc);
355  // total elastic mean free path in Geant4 internal lenght units
356  lambda0 = beta2/bc;//*(1.+scrA);
357  g1 = 2.0*scrA*((1.0+scrA)*G4Log(1.0/scrA+1.0)-1.0);
358  lambda1 = lambda0/g1;
359  }
360 
361  return lambda1;
362 }
std::vector< G4Element * > G4ElementVector
const G4Material * GetMaterial() const
static G4PWATotalXsecTable * fgPWAXsecTable
size_t GetIndex() const
Definition: G4Material.hh:262
int G4int
Definition: G4Types.hh:78
Float_t mat
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
float electron_mass_c2
Definition: hepunit.py:274
static G4GoudsmitSaundersonTable * fgGSTable
const G4PWATotalXsecZ * GetPWATotalXsecForZet(G4int Z) const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4int GetPWATotalXsecEnergyBinIndex(G4double energy) const
G4double GetMoliereXc2(G4int matindx)
int G4lrint(double ad)
Definition: templates.hh:163
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4double GetMoliereBc(G4int matindx)
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
double G4double
Definition: G4Types.hh:76
const G4MaterialCutsCouple * currentCouple
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Initialise()

void G4GoudsmitSaundersonMscModel::Initialise ( const G4ParticleDefinition p,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 205 of file G4GoudsmitSaundersonMscModel.cc.

206  {
207  SetParticle(p);
209  //
210  // -create static GoudsmitSaundersonTable and PWATotalXsecTable is necessary
211  //
212  if(IsMaster()) {
213 
214  // Could delete as well if they are exist but then the same data are reloaded
215  // in case of e+ for example.
216 /* if(fgGSTable) {
217  delete fgGSTable;
218  fgGSTable = 0;
219  }
220  if(fgPWAXsecTable) {
221  delete fgPWAXsecTable;
222  fgPWAXsecTable = 0;
223  }
224 */
225  if(!fgGSTable)
227 
228  // G4GSTable will be initialised (data are loaded) only if they are not there yet.
230 
232  if(!fgPWAXsecTable)
234 
235  // Makes sure that all (and only those) atomic PWA data are in memory that
236  // are in the current MaterilaTable.
238  }
239  }
240 }
static G4PWATotalXsecTable * fgPWAXsecTable
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:90
static G4GoudsmitSaundersonTable * fgGSTable
void SetParticle(const G4ParticleDefinition *p)
Here is the call graph for this function:

◆ operator=()

G4GoudsmitSaundersonMscModel& G4GoudsmitSaundersonMscModel::operator= ( const G4GoudsmitSaundersonMscModel right)
private
Here is the caller graph for this function:

◆ Randomizetlimit()

G4double G4GoudsmitSaundersonMscModel::Randomizetlimit ( )
inlineprivate

Definition at line 221 of file G4GoudsmitSaundersonMscModel.hh.

222 {
223  G4double temptlimit = tlimit;
224  do {
225  temptlimit = G4RandGauss::shoot(rndmEngineMod,tlimit,0.3*tlimit);
226  } while ( (temptlimit<0.) || (temptlimit > 2.*tlimit));
227 
228  return temptlimit;
229 }
ThreeVector shoot(const G4int Ap, const G4int Af)
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleMSC()

void G4GoudsmitSaundersonMscModel::SampleMSC ( )

Definition at line 874 of file G4GoudsmitSaundersonMscModel.cc.

874  {
875  fIsNoScatteringInMSC = false;
876  // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
877  G4double kineticEnergy = currentKinEnergy;
878 
880  // Energy loss correction: 2 version
881  G4double eloss = 0.0;
882 // if (fTheTrueStepLenght > currentRange*dtrl) {
883  eloss = kineticEnergy -
885 // } else {
886 // eloss = fTheTrueStepLenght*GetDEDX(particle,kineticEnergy,currentCouple);
887 // }
888 
889 
890  G4double tau = 0.;// = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy
891  G4double tau2 = 0.;// = tau*tau;
892  G4double eps0 = 0.;// = eloss/kineticEnergy0; // energy loss fraction to the begin step energy
893  G4double epsm = 0.;// = eloss/kineticEnergy; // energy loss fraction to the mean step energy
894 
895 
896  // - init.
897  G4double efEnergy = kineticEnergy;
898  G4double efStep = fTheTrueStepLenght;
899 
900  G4double kineticEnergy0 = kineticEnergy;
901  if(fgIsUseAccurate) { // - use accurate energy loss correction
902  kineticEnergy -= 0.5*eloss; // mean energy along the full step
903 
904  // other parameters for energy loss corrections
905  tau = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy
906  tau2 = tau*tau;
907  eps0 = eloss/kineticEnergy0; // energy loss fraction to the begin step energy
908  epsm = eloss/kineticEnergy; // energy loss fraction to the mean step energy
909 
910  efEnergy = kineticEnergy * (1. - epsm*epsm*(6.+10.*tau+5.*tau2)/(24.*tau2+48.*tau+72.));
911  G4double dum = 0.166666*(4.+tau*(6.+tau*(7.+tau*(4.+tau))))*
912  (epsm/((tau+1.)*(tau+2.)))*(epsm/((tau+1.)*(tau+2.)));
913  efStep =fTheTrueStepLenght*(1.-dum);
914  } else { // - take only mean energy
915  kineticEnergy -= 0.5*eloss; // mean energy along the full step
916  efEnergy = kineticEnergy;
917  G4double factor = 1./(1. + 0.9784671*kineticEnergy); //0.9784671 = 1/(2*rm)
918  eps0= eloss/kineticEnergy0;
919  epsm= eps0/(1.-0.5*eps0);
920  G4double temp = 0.3*(1. - factor*(1. - 0.333333*factor))*eps0*eps0;
921  efStep = fTheTrueStepLenght*(1. + temp);
922  }
923 
925  // Compute elastic mfp, first transport mfp, screening parameter, and G1
927 
928  G4double lambdan=0.;
929 
930  if(fLambda0>0.0) { lambdan=efStep/fLambda0; }
931  if(lambdan<=1.0e-12) {
934  fIsNoScatteringInMSC = true;
935  return;
936  }
937 
938  // correction from neglected term 1+A (only for Moliere parameter)
940  lambdan = lambdan/(1.+fScrA);
941 
942  G4double Qn1 = lambdan *fG1;//2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
943 
945  // Sample scattering angles
946  // new direction, relative to the orriginal one is in {us,vs,ws}
947  G4double cosTheta1 =1.0,sinTheta1 =0.0, cosTheta2 =1.0, sinTheta2 =0.0;
948  G4double cosPhi1=1.0,sinPhi1=0.0, cosPhi2 =1.0, sinPhi2 =0.0;
949  G4double us=0.0,vs=0.0,ws=1.0,x_coord=0.0,y_coord=0.0,z_coord=1.0;
950  G4double u2=0.0,v2=0.0;
951 
952  // if we are above the upper grid limit with lambdaxG1=true-length/first-trans-mfp
953  // => izotropic distribution: lambG1_max =7.992 but set it to 7
954  if(0.5*Qn1 > 7.0){
955  G4double vrand[2];
956  rndmEngineMod->flatArray(2,vrand); //get 2 random number
957  cosTheta1 = 1.-2.*vrand[0];
958  sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+cosTheta1));
959  cosTheta2 = 1.-2.*vrand[1];
960  sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+cosTheta2));
961  } else {
962  // sample 2 scattering cost1, sint1, cost2 and sint2 for half path
963  fgGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta1, sinTheta1);
964  fgGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta2, sinTheta2);
965  if(cosTheta1 + cosTheta2 == 2.) { // no scattering happened
968  fIsNoScatteringInMSC = true;
969  return;
970  }
971  }
972  // sample 2 azimuthal angles
973  G4double vrand[2];
974  rndmEngineMod->flatArray(2,vrand); //get 2 random number
975  G4double phi1 = CLHEP::twopi*vrand[0];
976  sinPhi1 = std::sin(phi1);
977  cosPhi1 = std::cos(phi1);
978  G4double phi2 = CLHEP::twopi*vrand[1];
979  sinPhi2 = std::sin(phi2);
980  cosPhi2 = std::cos(phi2);
981 
982  // compute final direction realtive to z-dir
983  u2 = sinTheta2*cosPhi2;
984  v2 = sinTheta2*sinPhi2;
985  G4double u2p = cosTheta1*u2 + sinTheta1*cosTheta2;
986  us = u2p*cosPhi1 - v2*sinPhi1;
987  vs = u2p*sinPhi1 + v2*cosPhi1;
988  ws = cosTheta1*cosTheta2 - sinTheta1*u2;
989 
991  fTheNewDirection.set(us,vs,ws);
992 
993  // set the fTheZPathLenght if we don't sample displacement and
994  // we should do everything at the step-limit-phase before we return
997 
998  // in optimized-mode if the current-safety > current-range we do not use dispalcement
999  if(fIsNoDisplace)
1000  return;
1001 
1003  // Compute final position
1004  if(fgIsUseAccurate) {
1005  // correction parameter
1006  G4double par =1.;
1007  if(Qn1<0.7) par = 1.;
1008  else if (Qn1<7.0) par = -0.031376*Qn1+1.01356;
1009  else par = 0.79;
1010 
1011  // Moments with energy loss correction
1012  // --first the uncorrected (for energy loss) values of gamma, eta, a1=a2=0.5*(1-eta), delta
1013  // gamma = G_2/G_1 based on G2 computed from A by using the Wentzel DCS form of G2
1014  G4double loga = G4Log(1.0+1.0/fScrA);
1015  G4double gamma = 6.0*fScrA*(1.0 + fScrA)*(loga*(1.0 + 2.0*fScrA) - 2.0)/fG1;
1016  // sample eta from p(eta)=2*eta i.e. P(eta) = eta_square ;-> P(eta) = rand --> eta = sqrt(rand)
1017  G4double eta = std::sqrt(rndmEngineMod->flat());
1018  G4double eta1 = 0.5*(1 - eta); // used more than once
1019  // 0.5 +sqrt(6)/6 = 0.9082483;
1020  // 1/(4*sqrt(6)) = 0.1020621;
1021  // (4-sqrt(6)/(24*sqrt(6))) = 0.026374715
1022  // delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1 without energy loss cor.
1023  G4double delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1;
1024 
1025  // compute alpha1 and alpha2 for energy loss correction
1026  G4double temp1 = 2.0 + tau;
1027  G4double temp = (2.0+tau*temp1)/((tau+1.0)*temp1);
1028  //Take logarithmic dependence
1029  temp = temp - (tau+1.0)/((tau+2.0)*(loga*(1.0+fScrA)-1.0));
1030  temp = temp * epsm;
1031  temp1 = 1.0 - temp;
1032  delta = delta + 0.40824829*(eps0*(tau+1.0)/((tau+2.0)*
1033  (loga*(1.0+fScrA)-1.0)*(loga*(1.0+2.0*fScrA)-2.0)) - 0.25*temp*temp);
1034  G4double b = eta*delta;
1035  G4double c = eta*(1.0-delta);
1036 
1037  //Calculate transport direction cosines:
1038  // ut,vt,wt is the final position divided by the true step length
1039  G4double w1v2 = cosTheta1*v2;
1040  G4double ut = b*sinTheta1*cosPhi1 + c*(cosPhi1*u2 - sinPhi1*w1v2) + eta1*us*temp1;
1041  G4double vt = b*sinTheta1*sinPhi1 + c*(sinPhi1*u2 + cosPhi1*w1v2) + eta1*vs*temp1;
1042  G4double wt = eta1*(1+temp) + b*cosTheta1 + c*cosTheta2 + eta1*ws*temp1;
1043 
1044  // long step correction
1045  ut *=par;
1046  vt *=par;
1047  wt *=par;
1048 
1049  // final position relative to the pre-step point in the scattering frame
1050  // ut = x_f/s so needs to multiply by s
1051  x_coord = ut*fTheTrueStepLenght;
1052  y_coord = vt*fTheTrueStepLenght;
1053  z_coord = wt*fTheTrueStepLenght;
1054 
1056  // We sample in the step limit so set fTheZPathLenght = transportDistance
1057  // and lateral displacement (x_coord,y_coord,z_coord-transportDistance)
1058  //Calculate transport distance
1059  G4double transportDistance = std::sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
1060  // protection
1061  if(transportDistance>fTheTrueStepLenght)
1062  transportDistance = fTheTrueStepLenght;
1063  fTheZPathLenght = transportDistance;
1064  }
1065  // else:: we sample in the DoIt so
1066  // the fTheZPathLenght was already set and was taken as transport along zet
1067  fTheDisplacementVector.set(x_coord,y_coord,z_coord-fTheZPathLenght);
1068  } else {
1069  // compute zz = <z>/tPathLength
1070  // s -> true-path-length
1071  // z -> geom-path-length:: when PRESTA is used z =(def.) <z>
1072  // r -> lateral displacement = s/2 sin(theta) => x_f = r cos(phi); y_f = r sin(phi)
1073  G4double zz = 0.0;
1075  // We sample in the step limit so set fTheZPathLenght = transportDistance
1076  // and lateral displacement (x_coord,y_coord,z_coord-transportDistance)
1077  if(Qn1<0.1) { // use 3-order Taylor approximation of (1-exp(-x))/x around x=0
1078  zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666667 - 0.041666667*Qn1)); // 1/6 =0.166..7 ; 1/24=0.041..
1079  } else {
1080  zz = (1.-G4Exp(-Qn1))/Qn1;
1081  }
1082  } else {
1083  // we sample in the DoIt so
1084  // the fTheZPathLenght was already set and was taken as transport along zet
1086  }
1087 
1088  G4double rr = (1.-zz*zz)/(1.-ws*ws); // s^2 >= <z>^2+r^2 :: where r^2 = s^2/4 sin^2(theta)
1089  if(rr >= 0.25) rr = 0.25; // (1-<z>^2/s^2)/sin^2(theta) >= r^2/(s^2 sin^2(theta)) = 1/4 must hold
1090  G4double rperp = fTheTrueStepLenght*std::sqrt(rr); // this is r/sint
1091  x_coord = rperp*us;
1092  y_coord = rperp*vs;
1093  z_coord = zz*fTheTrueStepLenght;
1094 
1096  G4double transportDistance = std::sqrt(x_coord*x_coord + y_coord*y_coord + z_coord*z_coord);
1097  fTheZPathLenght = transportDistance;
1098  }
1099 
1100  fTheDisplacementVector.set(x_coord,y_coord,z_coord- fTheZPathLenght);
1101  }
1102 }
void set(double x, double y, double z)
const G4ParticleDefinition * particle
virtual double flat()=0
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:315
float electron_mass_c2
Definition: hepunit.py:274
void Sampling(G4double, G4double, G4double, G4double &, G4double &)
static G4GoudsmitSaundersonTable * fgGSTable
G4double G4Log(G4double x)
Definition: G4Log.hh:230
Double_t zz
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const G4double factor
G4double GetTransportMeanFreePath(const G4ParticleDefinition *, G4double)
double G4double
Definition: G4Types.hh:76
static const double twopi
Definition: SystemOfUnits.h:54
virtual void flatArray(const int size, double *vect)=0
const G4MaterialCutsCouple * currentCouple
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleScattering()

G4ThreeVector & G4GoudsmitSaundersonMscModel::SampleScattering ( const G4ThreeVector oldDirection,
G4double  safety 
)
virtual

Reimplemented from G4VMscModel.

Definition at line 811 of file G4GoudsmitSaundersonMscModel.cc.

811  {
812  if(steppingAlgorithm == fUseDistanceToBoundary && fIsEverythingWasDone && fIsSingleScattering){ // ONLY single scattering is done in advance
813  // single scattering was and scattering happend
814  fTheNewDirection.rotateUz(oldDirection);
815 // fTheDisplacementVector.set(0.,0.,0.);
816  fParticleChange->ProposeMomentumDirection(fTheNewDirection);
817  return fTheDisplacementVector;
818  } else if(steppingAlgorithm == fUseSafety) {
819  if(fIsEndedUpOnBoundary) {// do nothing
820 // fTheDisplacementVector.set(0.,0.,0.);
821  return fTheDisplacementVector;
822  } else if(fIsEverythingWasDone) {// evrything is done if not optimizations case !!!
823  // check single scattering and see if it happened
824  if(fIsSingleScattering) {
825  fTheNewDirection.rotateUz(oldDirection);
826  // fTheDisplacementVector.set(0.,0.,0.);
827  fParticleChange->ProposeMomentumDirection(fTheNewDirection);
828  return fTheDisplacementVector;
829  }
830 
831  // check if multiple scattering happened and do things only if scattering was really happening
833  fTheNewDirection.rotateUz(oldDirection);
834  fTheDisplacementVector.rotateUz(oldDirection);
835  fParticleChange->ProposeMomentumDirection(fTheNewDirection);
836  }
837  // The only thing that could happen if we are here (fUseSafety and fIsEverythingWasDone)
838  // is that single scattering was tried but did not win so scattering did not happen.
839  // So no displacement and no scattering
840  return fTheDisplacementVector;
841  }
842  //
843  // The only thing that could still happen with fUseSafety is that we are in the
844  // optimization branch: so sample MSC angle here (no displacement)
845  }
846 
847  //else
848  // MSC needs to be done here
849  SampleMSC();
850  if(!fIsNoScatteringInMSC) {
851  fTheNewDirection.rotateUz(oldDirection);
852  fParticleChange->ProposeMomentumDirection(fTheNewDirection);
853  if(!fIsNoDisplace)
854  fTheDisplacementVector.rotateUz(oldDirection);
855  }
856  return fTheDisplacementVector;
857 }
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:187
Here is the call graph for this function:

◆ SetOptionPWAScreening()

void G4GoudsmitSaundersonMscModel::SetOptionPWAScreening ( G4bool  opt)
inline

Definition at line 137 of file G4GoudsmitSaundersonMscModel.hh.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetParticle()

void G4GoudsmitSaundersonMscModel::SetParticle ( const G4ParticleDefinition p)
inlineprivate

Definition at line 209 of file G4GoudsmitSaundersonMscModel.hh.

210 {
211  if (p != particle) {
212  particle = p;
214  mass = p->GetPDGMass();
215  }
216 }
const G4ParticleDefinition * particle
int G4int
Definition: G4Types.hh:78
static const double eplus
G4double GetPDGCharge() const
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SingleScattering()

void G4GoudsmitSaundersonMscModel::SingleScattering ( G4double cost,
G4double sint 
)

Definition at line 860 of file G4GoudsmitSaundersonMscModel.cc.

860  {
861  G4double rand1 = rndmEngineMod->flat();
862  // sampling 1-cos(theta)
863  G4double dum0 = 2.0*fScrA*rand1/(1.0 - rand1 + fScrA);
864  // add protections
865  if(dum0 < 0.0) dum0 = 0.0;
866  else if(dum0 > 2.0 ) dum0 = 2.0;
867  // compute cos(theta) and sin(theta)
868  cost = 1.0 - dum0;
869  sint = std::sqrt(dum0*(2.0 - dum0));
870 }
virtual double flat()=0
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ StartTracking()

void G4GoudsmitSaundersonMscModel::StartTracking ( G4Track *  track)
virtual

Reimplemented from G4VEmModel.

Definition at line 365 of file G4GoudsmitSaundersonMscModel.cc.

Here is the call graph for this function:

Member Data Documentation

◆ charge

G4int G4GoudsmitSaundersonMscModel::charge
private

Definition at line 162 of file G4GoudsmitSaundersonMscModel.hh.

◆ currentCouple

const G4MaterialCutsCouple* G4GoudsmitSaundersonMscModel::currentCouple
private

Definition at line 171 of file G4GoudsmitSaundersonMscModel.hh.

◆ currentKinEnergy

G4double G4GoudsmitSaundersonMscModel::currentKinEnergy
private

Definition at line 157 of file G4GoudsmitSaundersonMscModel.hh.

◆ currentMaterialIndex

G4int G4GoudsmitSaundersonMscModel::currentMaterialIndex
private

Definition at line 162 of file G4GoudsmitSaundersonMscModel.hh.

◆ currentRange

G4double G4GoudsmitSaundersonMscModel::currentRange
private

Definition at line 158 of file G4GoudsmitSaundersonMscModel.hh.

◆ fG1

G4double G4GoudsmitSaundersonMscModel::fG1
private

Definition at line 185 of file G4GoudsmitSaundersonMscModel.hh.

◆ fgGSTable

G4GoudsmitSaundersonTable * G4GoudsmitSaundersonMscModel::fgGSTable = 0
staticprivate

Definition at line 173 of file G4GoudsmitSaundersonMscModel.hh.

◆ fgIsOptimizationOn

G4bool G4GoudsmitSaundersonMscModel::fgIsOptimizationOn = true
staticprivate

Definition at line 204 of file G4GoudsmitSaundersonMscModel.hh.

◆ fgIsUseAccurate

G4bool G4GoudsmitSaundersonMscModel::fgIsUseAccurate = true
staticprivate

Definition at line 203 of file G4GoudsmitSaundersonMscModel.hh.

◆ fgPWAXsecTable

G4PWATotalXsecTable * G4GoudsmitSaundersonMscModel::fgPWAXsecTable = 0
staticprivate

Definition at line 174 of file G4GoudsmitSaundersonMscModel.hh.

◆ firstStep

G4bool G4GoudsmitSaundersonMscModel::firstStep
private

Definition at line 164 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsEndedUpOnBoundary

G4bool G4GoudsmitSaundersonMscModel::fIsEndedUpOnBoundary
private

Definition at line 193 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsEverythingWasDone

G4bool G4GoudsmitSaundersonMscModel::fIsEverythingWasDone
private

Definition at line 196 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsFirstRealStep

G4bool G4GoudsmitSaundersonMscModel::fIsFirstRealStep
private

Definition at line 201 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsInsideSkin

G4bool G4GoudsmitSaundersonMscModel::fIsInsideSkin
private

Definition at line 199 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsMultipleSacettring

G4bool G4GoudsmitSaundersonMscModel::fIsMultipleSacettring
private

Definition at line 194 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsNoDisplace

G4bool G4GoudsmitSaundersonMscModel::fIsNoDisplace
private

Definition at line 198 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsNoScatteringInMSC

G4bool G4GoudsmitSaundersonMscModel::fIsNoScatteringInMSC
private

Definition at line 197 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsSingleScattering

G4bool G4GoudsmitSaundersonMscModel::fIsSingleScattering
private

Definition at line 195 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsUsePWATotalXsecData

G4bool G4GoudsmitSaundersonMscModel::fIsUsePWATotalXsecData
private

Definition at line 176 of file G4GoudsmitSaundersonMscModel.hh.

◆ fIsWasOnBoundary

G4bool G4GoudsmitSaundersonMscModel::fIsWasOnBoundary
private

Definition at line 200 of file G4GoudsmitSaundersonMscModel.hh.

◆ fLambda0

G4double G4GoudsmitSaundersonMscModel::fLambda0
private

Definition at line 182 of file G4GoudsmitSaundersonMscModel.hh.

◆ fLambda1

G4double G4GoudsmitSaundersonMscModel::fLambda1
private

Definition at line 183 of file G4GoudsmitSaundersonMscModel.hh.

◆ fParticleChange

G4ParticleChangeForMSC* G4GoudsmitSaundersonMscModel::fParticleChange
private

Definition at line 170 of file G4GoudsmitSaundersonMscModel.hh.

◆ fr

G4double G4GoudsmitSaundersonMscModel::fr
private

Definition at line 160 of file G4GoudsmitSaundersonMscModel.hh.

◆ fScrA

G4double G4GoudsmitSaundersonMscModel::fScrA
private

Definition at line 184 of file G4GoudsmitSaundersonMscModel.hh.

◆ fTheDisplacementVector

G4ThreeVector G4GoudsmitSaundersonMscModel::fTheDisplacementVector
private

Definition at line 190 of file G4GoudsmitSaundersonMscModel.hh.

◆ fTheNewDirection

G4ThreeVector G4GoudsmitSaundersonMscModel::fTheNewDirection
private

Definition at line 191 of file G4GoudsmitSaundersonMscModel.hh.

◆ fTheTransportDistance

G4double G4GoudsmitSaundersonMscModel::fTheTransportDistance
private

Definition at line 188 of file G4GoudsmitSaundersonMscModel.hh.

◆ fTheTrueStepLenght

G4double G4GoudsmitSaundersonMscModel::fTheTrueStepLenght
private

Definition at line 187 of file G4GoudsmitSaundersonMscModel.hh.

◆ fTheZPathLenght

G4double G4GoudsmitSaundersonMscModel::fTheZPathLenght
private

Definition at line 189 of file G4GoudsmitSaundersonMscModel.hh.

◆ fZeff

G4double G4GoudsmitSaundersonMscModel::fZeff
private

Definition at line 179 of file G4GoudsmitSaundersonMscModel.hh.

◆ geombig

G4double G4GoudsmitSaundersonMscModel::geombig
private

Definition at line 160 of file G4GoudsmitSaundersonMscModel.hh.

◆ geomlimit

G4double G4GoudsmitSaundersonMscModel::geomlimit
private

Definition at line 160 of file G4GoudsmitSaundersonMscModel.hh.

◆ highKEnergy

G4double G4GoudsmitSaundersonMscModel::highKEnergy
private

Definition at line 156 of file G4GoudsmitSaundersonMscModel.hh.

◆ lambdalimit

G4double G4GoudsmitSaundersonMscModel::lambdalimit
private

Definition at line 161 of file G4GoudsmitSaundersonMscModel.hh.

◆ lowKEnergy

G4double G4GoudsmitSaundersonMscModel::lowKEnergy
private

Definition at line 155 of file G4GoudsmitSaundersonMscModel.hh.

◆ mass

G4double G4GoudsmitSaundersonMscModel::mass
private

Definition at line 166 of file G4GoudsmitSaundersonMscModel.hh.

◆ par1

G4double G4GoudsmitSaundersonMscModel::par1
private

Definition at line 166 of file G4GoudsmitSaundersonMscModel.hh.

◆ par2

G4double G4GoudsmitSaundersonMscModel::par2
private

Definition at line 166 of file G4GoudsmitSaundersonMscModel.hh.

◆ par3

G4double G4GoudsmitSaundersonMscModel::par3
private

Definition at line 166 of file G4GoudsmitSaundersonMscModel.hh.

◆ particle

const G4ParticleDefinition* G4GoudsmitSaundersonMscModel::particle
private

Definition at line 169 of file G4GoudsmitSaundersonMscModel.hh.

◆ presafety

G4double G4GoudsmitSaundersonMscModel::presafety
private

Definition at line 178 of file G4GoudsmitSaundersonMscModel.hh.

◆ rangeinit

G4double G4GoudsmitSaundersonMscModel::rangeinit
private

Definition at line 160 of file G4GoudsmitSaundersonMscModel.hh.

◆ rndmEngineMod

CLHEP::HepRandomEngine* G4GoudsmitSaundersonMscModel::rndmEngineMod
private

Definition at line 152 of file G4GoudsmitSaundersonMscModel.hh.

◆ taulim

G4double G4GoudsmitSaundersonMscModel::taulim
private

Definition at line 166 of file G4GoudsmitSaundersonMscModel.hh.

◆ tausmall

G4double G4GoudsmitSaundersonMscModel::tausmall
private

Definition at line 166 of file G4GoudsmitSaundersonMscModel.hh.

◆ tgeom

G4double G4GoudsmitSaundersonMscModel::tgeom
private

Definition at line 161 of file G4GoudsmitSaundersonMscModel.hh.

◆ theManager

G4LossTableManager* G4GoudsmitSaundersonMscModel::theManager
private

Definition at line 168 of file G4GoudsmitSaundersonMscModel.hh.

◆ tlimit

G4double G4GoudsmitSaundersonMscModel::tlimit
private

Definition at line 161 of file G4GoudsmitSaundersonMscModel.hh.

◆ tlimitminfix2

G4double G4GoudsmitSaundersonMscModel::tlimitminfix2
private

Definition at line 166 of file G4GoudsmitSaundersonMscModel.hh.


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