Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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) override
 
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 level, const G4ParticleDefinition *, G4double kineticEnergy)
 
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 * > *)
 
virtual 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=nullptr)
 
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)
 

Additional Inherited Members

- Protected Member Functions inherited from G4VMscModel
G4ParticleChangeForMSCGetParticleChangeForMSC (const G4ParticleDefinition *p=nullptr)
 
G4double ConvertTrueToGeom (G4double &tLength, G4double &gLength)
 
- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
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
 
G4VParticleChangepParticleChange
 
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::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;
128  currentMaterialIndex = -1;
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;
140  theManager = G4LossTableManager::Instance();
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;
148  mass = electron_mass_c2;
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.;
165  fTheTransportDistance = 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()
static constexpr double mm
Definition: G4SIunits.hh:115
static constexpr double electron_mass_c2
G4double facsafety
Definition: G4VMscModel.hh:177
static constexpr double nm
Definition: G4SIunits.hh:112
static constexpr double GeV
Definition: G4SIunits.hh:217
G4VMscModel(const G4String &nam)
Definition: G4VMscModel.cc:60
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

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 }
G4bool IsMaster() const
Definition: G4VEmModel.hh:717

Here is the call graph for this function:

Member Function Documentation

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
718  fTheTrueStepLenght = std::min(fTheTrueStepLenght,currentRange);
719 
720  // do the true -> geom transformation
721  fTheZPathLenght = fTheTrueStepLenght;
722 
723  // z = t for very small true-path-length
724  if(fTheTrueStepLenght < tlimitminfix2) return fTheZPathLenght;
725 
726  G4double tau = fTheTrueStepLenght/fLambda1;
727 
728  if (tau <= tausmall) {
729  fTheZPathLenght = std::min(fTheTrueStepLenght, fLambda1);
730  } else if (fTheTrueStepLenght < currentRange*dtrl) {
731  if(tau < taulim) fTheZPathLenght = fTheTrueStepLenght*(1.-0.5*tau) ;
732  else fTheZPathLenght = fLambda1*(1.-G4Exp(-tau));
733  } else if(currentKinEnergy < mass || fTheTrueStepLenght == currentRange) {
734  par1 = 1./currentRange ; // alpha =1/range_init for Ekin<mass
735  par2 = 1./(par1*fLambda1) ; // 1/(alphaxlambda01)
736  par3 = 1.+par2 ; // 1+1/
737  if(fTheTrueStepLenght < currentRange) {
738  fTheZPathLenght = 1./(par1*par3) * (1.-std::pow(1.-par1*fTheTrueStepLenght,par3));
739  } else {
740  fTheZPathLenght = 1./(par1*par3);
741  }
742 
743  } else {
744  G4double rfin = std::max(currentRange-fTheTrueStepLenght, 0.01*currentRange);
745  G4double T1 = GetEnergy(particle,rfin,currentCouple);
746  G4double lambda1 = GetTransportMeanFreePathOnly(particle,T1);
747 
748  par1 = (fLambda1-lambda1)/(fLambda1*fTheTrueStepLenght); // alpha
749  par2 = 1./(par1*fLambda1);
750  par3 = 1.+par2 ;
751  G4Pow *g4calc = G4Pow::GetInstance();
752  fTheZPathLenght = 1./(par1*par3) * (1.-g4calc->powA(1.-par1*fTheTrueStepLenght,par3));
753  }
754  }
755  fTheZPathLenght = std::min(fTheZPathLenght, fLambda1);
756 
757  return fTheZPathLenght;
758 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4double dtrl
Definition: G4VMscModel.hh:179
Definition: G4Pow.hh:56
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:303
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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();
385  SetCurrentCouple(currentCouple);
386  currentMaterialIndex = currentCouple->GetIndex();
387  currentKinEnergy = dp->GetKineticEnergy();
388  currentRange = GetRange(particle,currentKinEnergy,currentCouple);
389 
390 
391 
392  // *** Elastic and first transport mfp, fScrA and fG1 are also set in this !!!
393  fLambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
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
421  fZeff = currentCouple->GetMaterial()->GetIonisation()->GetZeffective();
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)
452  geomlimit = ComputeGeomLimit(track, presafety, currentRange);
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
489  fTheZPathLenght = fTheTrueStepLenght;
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.
509  if(fIsWasOnBoundary && !fIsInsideSkin){
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
522  if(firstStep || fIsFirstRealStep){
523  rangeinit = currentRange;
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.)
539  geomlimit = -fLambda1*G4Log(1.-geomlimit/fLambda1);
540  // the 2-different case that could lead us here
541  if(firstStep)
542  tgeom = 2.*geomlimit/facgeom;
543  else
544  tgeom = geomlimit/facgeom;
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.
553  tlimit = facrange*rangeinit;
554  // Take the minimum of the true step length limits coming from
555  // geometrical constraint or range-factor limitation
556  tlimit = std::min(tlimit,tgeom);
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) {
572  fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit());
573  } else {
574  fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit);
575  }
576  }
577  } else if (steppingAlgorithm == fUseSafety) { // THE ERROR_FREE stepping alg.
578  presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
579  geomlimit = presafety;
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
606  fTheZPathLenght = fTheTrueStepLenght;
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
619  fTheTrueStepLenght = std::min(fTheTrueStepLenght, facrange*currentRange);
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.
623  if(fTheTrueStepLenght > presafety)
624  fTheTrueStepLenght = std::min(fTheTrueStepLenght, presafety);
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
630  fTheTrueStepLenght = std::min(fTheTrueStepLenght, fLambda1*0.4);
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) {
671  fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit());
672  } else {
673  fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit);
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
682  if(fIsEverythingWasDone){
683  if(fIsSingleScattering){
684  // sample single scattering
685  G4double cost,sint,cosPhi,sinPhi;
686  SingleScattering(cost, sint);
687  G4double phi = CLHEP::twopi*rndmEngineMod->flat();
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)
G4double facgeom
Definition: G4VMscModel.hh:176
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
G4double GetKineticEnergy() const
const G4DynamicParticle * GetDynamicParticle() const
G4double facrange
Definition: G4VMscModel.hh:175
G4StepStatus GetStepStatus() const
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:245
G4double skin
Definition: G4VMscModel.hh:178
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
virtual double flat()=0
const G4Step * GetStep() const
G4double ComputeSafety(const G4ThreeVector &position, G4double limit=DBL_MAX)
Definition: G4VMscModel.hh:237
G4StepStatus
Definition: G4StepStatus.hh:49
G4double GetZeffective() const
G4StepPoint * GetPreStepPoint() const
void SingleScattering(G4double &cost, G4double &sint)
const G4ThreeVector & GetPosition() const
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:283
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:255
G4double G4Log(G4double x)
Definition: G4Log.hh:230
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double facsafety
Definition: G4VMscModel.hh:177
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:444
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetSafety() const
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:185
G4double GetTransportMeanFreePath(const G4ParticleDefinition *, G4double)
double G4double
Definition: G4Types.hh:76
static constexpr double twopi
Definition: SystemOfUnits.h:55
const G4Material * GetMaterial() const

Here is the call graph for this function:

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
779  if(fIsEverythingWasDone && !fIsMultipleSacettring) {
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 *g4calc = G4Pow::GetInstance();
796  tlength = (1.-g4calc->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
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
Definition: G4Pow.hh:56
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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 
254  const G4Material* mat = currentCouple->GetMaterial();
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 
261  if(fIsUsePWATotalXsecData) {
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 
287  fScrA = fgGSTable->GetScreeningParam(fG1);
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);
303  fLambda1 = fLambda0/fG1;
304  }
305 
306  return fLambda1;
307 }
std::vector< G4Element * > G4ElementVector
size_t GetIndex() const
Definition: G4Material.hh:262
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
static constexpr double electron_mass_c2
G4int GetPWATotalXsecEnergyBinIndex(G4double energy) const
const G4PWATotalXsecZ * GetPWATotalXsecForZet(G4int Z) const
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double GetMoliereXc2(G4int matindx)
int G4lrint(double ad)
Definition: templates.hh:163
G4double GetMoliereBc(G4int matindx)
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76
const G4Material * GetMaterial() const

Here is the call graph for this function:

Here is the caller graph for this function:

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

Implements G4VEmModel.

Definition at line 205 of file G4GoudsmitSaundersonMscModel.cc.

206  {
207  SetParticle(p);
208  fParticleChange = GetParticleChangeForMSC(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)
226  fgGSTable = new G4GoudsmitSaundersonTable();
227 
228  // G4GSTable will be initialised (data are loaded) only if they are not there yet.
229  fgGSTable->Initialise();
230 
231  if(fIsUsePWATotalXsecData) {
232  if(!fgPWAXsecTable)
233  fgPWAXsecTable = new G4PWATotalXsecTable();
234 
235  // Makes sure that all (and only those) atomic PWA data are in memory that
236  // are in the current MaterilaTable.
237  fgPWAXsecTable->Initialise();
238  }
239  }
240 }
G4bool IsMaster() const
Definition: G4VEmModel.hh:717
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=nullptr)
Definition: G4VMscModel.cc:91

Here is the call graph for this function:

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 -
884  GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple);
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
926  fLambda1 = GetTransportMeanFreePath(particle,efEnergy);
927 
928  G4double lambdan=0.;
929 
930  if(fLambda0>0.0) { lambdan=efStep/fLambda0; }
931  if(lambdan<=1.0e-12) {
932  if(fIsEverythingWasDone)
933  fTheZPathLenght = fTheTrueStepLenght;
934  fIsNoScatteringInMSC = true;
935  return;
936  }
937 
938  // correction from neglected term 1+A (only for Moliere parameter)
939  if(!fIsUsePWATotalXsecData)
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 {uss,vss,wss}
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 uss=0.0,vss=0.0,wss=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
966  if(fIsEverythingWasDone)
967  fTheZPathLenght = fTheTrueStepLenght;
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  uss = u2p*cosPhi1 - v2*sinPhi1;
987  vss = u2p*sinPhi1 + v2*cosPhi1;
988  wss = cosTheta1*cosTheta2 - sinTheta1*u2;
989 
991  fTheNewDirection.set(uss,vss,wss);
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
995  if(fIsNoDisplace && fIsEverythingWasDone)
996  fTheZPathLenght = fTheTrueStepLenght;
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*uss*temp1;
1041  G4double vt = b*sinTheta1*sinPhi1 + c*(sinPhi1*u2 + cosPhi1*w1v2) + eta1*vss*temp1;
1042  G4double wt = eta1*(1+temp) + b*cosTheta1 + c*cosTheta2 + eta1*wss*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 
1055  if(fIsEverythingWasDone){
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;
1074  if(fIsEverythingWasDone){
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
1085  zz = fTheZPathLenght/fTheTrueStepLenght;
1086  }
1087 
1088  G4double rr = (1.-zz*zz)/(1.-wss*wss); // 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*uss;
1092  y_coord = rperp*vss;
1093  z_coord = zz*fTheTrueStepLenght;
1094 
1095  if(fIsEverythingWasDone){
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)
virtual double flat()=0
static constexpr double electron_mass_c2
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:303
void Sampling(G4double, G4double, G4double, G4double &, G4double &)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetTransportMeanFreePath(const G4ParticleDefinition *, G4double)
double G4double
Definition: G4Types.hh:76
virtual void flatArray(const int size, double *vect)=0
static constexpr double twopi
Definition: SystemOfUnits.h:55

Here is the call graph for this function:

Here is the caller graph for this function:

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
832  if(fIsMultipleSacettring && !fIsNoScatteringInMSC){
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
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:185

Here is the call graph for this function:

void G4GoudsmitSaundersonMscModel::SetOptionPWAScreening ( G4bool  opt)
inline

Definition at line 137 of file G4GoudsmitSaundersonMscModel.hh.

137 {fIsUsePWATotalXsecData=opt;}

Here is the caller graph for this function:

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:

void G4GoudsmitSaundersonMscModel::StartTracking ( G4Track track)
virtual

Reimplemented from G4VEmModel.

Definition at line 365 of file G4GoudsmitSaundersonMscModel.cc.

366 {
367  SetParticle(track->GetDynamicParticle()->GetDefinition());
368  firstStep = true;
369  tlimit = tgeom = rangeinit = geombig;
370 }
const G4DynamicParticle * GetDynamicParticle() const
G4ParticleDefinition * GetDefinition() const

Here is the call graph for this function:


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