Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4GoudsmitSaundersonMscModel Class Reference

#include <G4GoudsmitSaundersonMscModel.hh>

+ Inheritance diagram for G4GoudsmitSaundersonMscModel:

Public Member Functions

 G4GoudsmitSaundersonMscModel (const G4String &nam="GoudsmitSaunderson")
 
 ~G4GoudsmitSaundersonMscModel () override
 
void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
void InitialiseLocal (const G4ParticleDefinition *p, G4VEmModel *masterModel) override
 
G4ThreeVectorSampleScattering (const G4ThreeVector &, G4double safety) override
 
G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep) override
 
G4double ComputeGeomPathLength (G4double truePathLength) override
 
G4double ComputeTrueStepLength (G4double geomStepLength) override
 
G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX) override
 
void StartTracking (G4Track *) override
 
void SampleMSC ()
 
G4double GetTransportMeanFreePath (const G4ParticleDefinition *, G4double)
 
void SetOptionPWACorrection (G4bool opt)
 
G4bool GetOptionPWACorrection () const
 
void SetOptionMottCorrection (G4bool opt)
 
G4bool GetOptionMottCorrection () const
 
G4GoudsmitSaundersonTableGetGSTable ()
 
G4GSPWACorrectionsGetPWACorrection ()
 
G4GoudsmitSaundersonMscModeloperator= (const G4GoudsmitSaundersonMscModel &right)=delete
 
 G4GoudsmitSaundersonMscModel (const G4GoudsmitSaundersonMscModel &)=delete
 
- Public Member Functions inherited from G4VMscModel
 G4VMscModel (const G4String &nam)
 
 ~G4VMscModel () override
 
virtual G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &stepLimit)=0
 
virtual G4double ComputeGeomPathLength (G4double truePathLength)=0
 
virtual G4double ComputeTrueStepLength (G4double geomPathLength)=0
 
virtual G4ThreeVectorSampleScattering (const G4ThreeVector &, G4double safety)=0
 
void InitialiseParameters (const G4ParticleDefinition *)
 
void DumpParameters (std::ostream &out) const
 
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 SetLambdaLimit (G4double)
 
void SetSafetyFactor (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 GetDEDX (const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple, G4double logKineticEnergy)
 
G4double GetRange (const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
 
G4double GetRange (const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple, G4double logKineticEnergy)
 
G4double GetEnergy (const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
 
G4double GetTransportMeanFreePath (const G4ParticleDefinition *part, G4double kinEnergy)
 
G4double GetTransportMeanFreePath (const G4ParticleDefinition *part, G4double kinEnergy, G4double logKinEnergy)
 
G4VMscModeloperator= (const G4VMscModel &right)=delete
 
 G4VMscModel (const G4VMscModel &)=delete
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)=0
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
 
virtual 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 StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, const G4double &length, G4double &eloss)
 
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 FillNumberOfSecondaries (G4int &numberOfTriplets, G4int &numberOfRecoil)
 
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)
 
const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectTargetAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementGetCurrentElement (const G4Material *mat=nullptr) const
 
G4int SelectRandomAtomNumber (const G4Material *) const
 
const G4IsotopeGetCurrentIsotope (const G4Element *elm=nullptr) const
 
G4int SelectIsotopeNumber (const G4Element *) const
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
G4VEmModelGetTripletModel ()
 
void SetTripletModel (G4VEmModel *)
 
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) const
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
void SetUseBaseMaterials (G4bool val)
 
G4bool UseBaseMaterials () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 
G4VEmModeloperator= (const G4VEmModel &right)=delete
 
 G4VEmModel (const G4VEmModel &)=delete
 

Additional Inherited Members

- Protected Member Functions inherited from G4VMscModel
G4ParticleChangeForMSCGetParticleChangeForMSC (const G4ParticleDefinition *p=nullptr)
 
G4double ConvertTrueToGeom (G4double &tLength, G4double &gLength)
 
void SetUseSplineForMSC (G4bool val)
 
- 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 = 0.04
 
G4double facgeom = 2.5
 
G4double facsafety = 0.6
 
G4double skin = 1.0
 
G4double dtrl = 0.05
 
G4double lambdalimit
 
G4double geomMin
 
G4double geomMax
 
G4ThreeVector fDisplacement
 
G4MscStepLimitType steppingAlgorithm
 
G4bool samplez = false
 
G4bool latDisplasment = true
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData = nullptr
 
G4VParticleChangepParticleChange = nullptr
 
G4PhysicsTablexSectionTable = nullptr
 
const G4MaterialpBaseMaterial = nullptr
 
const std::vector< G4double > * theDensityFactor = nullptr
 
const std::vector< G4int > * theDensityIdx = nullptr
 
G4double inveplus
 
G4double pFactor = 1.0
 
size_t currentCoupleIndex = 0
 
size_t basedCoupleIndex = 0
 
G4bool lossFlucFlag = true
 

Detailed Description

Definition at line 138 of file G4GoudsmitSaundersonMscModel.hh.

Constructor & Destructor Documentation

◆ G4GoudsmitSaundersonMscModel() [1/2]

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

Definition at line 163 of file G4GoudsmitSaundersonMscModel.cc.

164 : G4VMscModel(nam) {
165 charge = 0;
166 currentMaterialIndex = -1;
167 //
168 fr = 0.1;
169 rangeinit = 1.e+21;
170 geombig = 1.e+50*mm;
171 geomlimit = geombig;
172 tgeom = geombig;
173 tlimit = 1.e+10*mm;
174 presafety = 0.*mm;
175 //
176 particle = nullptr;
177 theManager = G4LossTableManager::Instance();
178 firstStep = true;
179 currentKinEnergy = 0.0;
180 currentRange = 0.0;
181 //
182 tlimitminfix2 = 1.*nm;
183 tausmall = 1.e-16;
184 mass = electron_mass_c2;
185 taulim = 1.e-6;
186 //
187 currentCouple = nullptr;
188 fParticleChange = nullptr;
189 //
190 fZeff = 1.;
191 //
192 par1 = 0.;
193 par2 = 0.;
194 par3 = 0.;
195 //
196 // Moliere screeing parameter will be used and (by default) corrections are
197 // appalied to the integrated quantities (screeing parameter, elastic mfp, first
198 // and second moments) derived from the corresponding PWA quantities
199 // this PWA correction is ignored if Mott-correction is set to true because
200 // Mott-correction contains all these corrections as well
201 fIsUsePWACorrection = true;
202 //
203 fIsUseMottCorrection = false;
204 //
205 fLambda0 = 0.0; // elastic mean free path
206 fLambda1 = 0.0; // first transport mean free path
207 fScrA = 0.0; // screening parameter
208 fG1 = 0.0; // first transport coef.
209 //
210 fMCtoScrA = 1.0;
211 fMCtoQ1 = 1.0;
212 fMCtoG2PerG1 = 1.0;
213 //
214 fTheTrueStepLenght = 0.;
215 fTheTransportDistance = 0.;
216 fTheZPathLenght = 0.;
217 //
218 fTheDisplacementVector.set(0.,0.,0.);
219 fTheNewDirection.set(0.,0.,1.);
220 //
221 fIsEverythingWasDone = false;
222 fIsMultipleSacettring = false;
223 fIsSingleScattering = false;
224 fIsEndedUpOnBoundary = false;
225 fIsNoScatteringInMSC = false;
226 fIsNoDisplace = false;
227 fIsInsideSkin = false;
228 fIsWasOnBoundary = false;
229 fIsFirstRealStep = false;
230 rndmEngineMod = G4Random::getTheEngine();
231 //
232 fGSTable = nullptr;
233 fPWACorrection = nullptr;
234}
void set(double x, double y, double z)
static G4LossTableManager * Instance()

◆ ~G4GoudsmitSaundersonMscModel()

G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel ( )
override

Definition at line 237 of file G4GoudsmitSaundersonMscModel.cc.

237 {
238 if (IsMaster()) {
239 if (fGSTable) {
240 delete fGSTable;
241 fGSTable = nullptr;
242 }
243 if (fPWACorrection) {
244 delete fPWACorrection;
245 fPWACorrection = nullptr;
246 }
247 }
248}
G4bool IsMaster() const
Definition: G4VEmModel.hh:725

◆ G4GoudsmitSaundersonMscModel() [2/2]

G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel ( const G4GoudsmitSaundersonMscModel )
delete

Member Function Documentation

◆ ComputeGeomPathLength()

G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength ( G4double  truePathLength)
overridevirtual

Implements G4VMscModel.

Definition at line 781 of file G4GoudsmitSaundersonMscModel.cc.

781 {
782 // convert true ->geom
783 // It is called from the step limitation ComputeTruePathLengthLimit if
784 // !fIsEverythingWasDone but protect:
785 par1 = -1.;
786 par2 = par3 = 0.;
787 // if fIsEverythingWasDone = TRUE => fTheZPathLenght is already set
788 // so return with the already known value
789 // Otherwise:
790 if (!fIsEverythingWasDone) {
791 // this correction needed to run MSC with eIoni and eBrem inactivated
792 // and makes no harm for a normal run
793 fTheTrueStepLenght = std::min(fTheTrueStepLenght,currentRange);
794 // do the true -> geom transformation
795 fTheZPathLenght = fTheTrueStepLenght;
796 // z = t for very small true-path-length
797 if (fTheTrueStepLenght<tlimitminfix2) {
798 return fTheZPathLenght;
799 }
800 G4double tau = fTheTrueStepLenght/fLambda1;
801 if (tau<=tausmall) {
802 fTheZPathLenght = std::min(fTheTrueStepLenght, fLambda1);
803 } else if (fTheTrueStepLenght<currentRange*dtrl) {
804 if (tau<taulim) fTheZPathLenght = fTheTrueStepLenght*(1.-0.5*tau) ;
805 else fTheZPathLenght = fLambda1*(1.-G4Exp(-tau));
806 } else if (currentKinEnergy<mass || fTheTrueStepLenght==currentRange) {
807 par1 = 1./currentRange ; // alpha =1/range_init for Ekin<mass
808 par2 = 1./(par1*fLambda1) ; // 1/(alphaxlambda01)
809 par3 = 1.+par2 ; // 1+1/
810 if (fTheTrueStepLenght<currentRange) {
811 fTheZPathLenght = 1./(par1*par3) * (1.-std::pow(1.-par1*fTheTrueStepLenght,par3));
812 } else {
813 fTheZPathLenght = 1./(par1*par3);
814 }
815 } else {
816 G4double rfin = std::max(currentRange-fTheTrueStepLenght, 0.01*currentRange);
817 G4double T1 = GetEnergy(particle,rfin,currentCouple);
818 G4double lambda1 = GetTransportMeanFreePathOnly(particle,T1);
819 //
820 par1 = (fLambda1-lambda1)/(fLambda1*fTheTrueStepLenght); // alpha
821 par2 = 1./(par1*fLambda1);
822 par3 = 1.+par2 ;
823 G4Pow *g4calc = G4Pow::GetInstance();
824 fTheZPathLenght = 1./(par1*par3) * (1.-g4calc->powA(1.-par1*fTheTrueStepLenght,par3));
825 }
826 }
827 fTheZPathLenght = std::min(fTheZPathLenght, fLambda1);
828 //
829 return fTheZPathLenght;
830}
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
double G4double
Definition: G4Types.hh:83
Definition: G4Pow.hh:49
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
G4double dtrl
Definition: G4VMscModel.hh:203
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.cc:223

◆ ComputeTruePathLengthLimit()

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

Implements G4VMscModel.

Definition at line 469 of file G4GoudsmitSaundersonMscModel.cc.

470 {
471 G4double skindepth = 0.;
472 //
473 const G4DynamicParticle* dp = track.GetDynamicParticle();
475 G4StepStatus stepStatus = sp->GetStepStatus();
476 currentCouple = track.GetMaterialCutsCouple();
477 SetCurrentCouple(currentCouple);
478 currentMaterialIndex = (G4int)currentCouple->GetMaterial()->GetIndex();
479 currentKinEnergy = dp->GetKineticEnergy();
480 currentRange = GetRange(particle,currentKinEnergy,currentCouple,
481 dp->GetLogKineticEnergy());
482 // elastic and first transport mfp, screening parameter and G1 are also set
483 // (Mott-correction will be used if it was requested by the user)
484 fLambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
485 // Set initial values:
486 // : lengths are initialised to currentMinimalStep which is the true, minimum
487 // step length from all other physics
488 fTheTrueStepLenght = currentMinimalStep;
489 fTheTransportDistance = currentMinimalStep;
490 fTheZPathLenght = currentMinimalStep; // will need to be converted
491 fTheDisplacementVector.set(0.,0.,0.);
492 fTheNewDirection.set(0.,0.,1.);
493
494 // Can everything be done in the step limit phase ?
495 fIsEverythingWasDone = false;
496 // Multiple scattering needs to be sample ?
497 fIsMultipleSacettring = false;
498 // Single scattering needs to be sample ?
499 fIsSingleScattering = false;
500 // Was zero deflection in multiple scattering sampling ?
501 fIsNoScatteringInMSC = false;
502 // Do not care about displacement in MSC sampling
503 // ( used only in the case of gIsOptimizationOn = true)
504 fIsNoDisplace = false;
505 // get pre-step point safety
506 presafety = sp->GetSafety();
507 //
508 fZeff = currentCouple->GetMaterial()->GetIonisation()->GetZeffective();
509 // distance will take into account max-fluct.
510 G4double distance = currentRange;
511 distance *= (1.20-fZeff*(1.62e-2-9.22e-5*fZeff));
512 //
513 // Possible optimization : if the distance is samller than the safety -> the
514 // particle will never leave this volume -> dispalcement
515 // as the effect of multiple elastic scattering can be skipped
516 // Important : this optimization can cause problems if one does scoring
517 // in a bigger volume since MSC won't be done deep inside the volume when
518 // distance < safety so don't use optimized-mode in such case.
519 if (gIsOptimizationOn && (distance<presafety)) {
520 // Indicate that we need to do MSC after transportation and no dispalcement.
521 fIsMultipleSacettring = true;
522 fIsNoDisplace = true;
524 //Compute geomlimit (and presafety) :
525 // - geomlimit will be:
526 // == the straight line distance to the boundary if currentRange is
527 // longer than that
528 // == a big value [geombig = 1.e50*mm] if currentRange is shorter than
529 // the straight line distance to the boundary
530 // - presafety will be updated as well
531 // So the particle can travell 'gemlimit' distance (along a straight
532 // line!) in its current direction:
533 // (1) before reaching a boundary (geomlimit < geombig) OR
534 // (2) before reaching its current range (geomlimit == geombig)
535 geomlimit = ComputeGeomLimit(track, presafety, currentRange);
536 // Record that the particle is on a boundary
537 if ( (stepStatus==fGeomBoundary) || (stepStatus==fUndefined && presafety==0.0)) {
538 fIsWasOnBoundary = true;
539 }
540 // Set skin depth = skin x elastic_mean_free_path
541 skindepth = skin*fLambda0;
542 // Init the flag that indicates that the particle are within a skindepth
543 // distance from a boundary
544 fIsInsideSkin = false;
545 // Check if we can try Single Scattering because we are within skindepth
546 // distance from/to a boundary OR the current minimum true-step-length is
547 // shorter than skindepth. NOTICE: the latest has only efficieny reasons
548 // because the MSC angular sampling is fine for any short steps but much
549 // faster to try single scattering in case of short steps.
550 if ((stepStatus==fGeomBoundary) || (presafety<skindepth) || (fTheTrueStepLenght<skindepth)) {
551 // check if we are within skindepth distance from a boundary
552 if ((stepStatus == fGeomBoundary) || (presafety < skindepth)) {
553 fIsInsideSkin = true;
554 fIsWasOnBoundary = true;
555 }
556 //Try single scattering:
557 // - sample distance to next single scattering interaction (sslimit)
558 // - compare to current minimum length
559 // == if sslimit is the shorter:
560 // - set the step length to sslimit
561 // - indicate that single scattering needs to be done
562 // == else : nothing to do
563 //- in both cases, the step length was very short so geometrical and
564 // true path length are the same
565 G4double sslimit = -1.*fLambda0*G4Log(G4UniformRand());
566 // compare to current minimum step length
567 if (sslimit<fTheTrueStepLenght) {
568 fTheTrueStepLenght = sslimit;
569 fIsSingleScattering = true;
570 }
571 // short step -> true step length equal to geometrical path length
572 fTheZPathLenght = fTheTrueStepLenght;
573 // Set taht everything is done in step-limit phase so no MSC call
574 // We will check if we need to perform the single-scattering angular
575 // sampling i.e. if single elastic scattering was the winer!
576 fIsEverythingWasDone = true;
577 } else {
578 // After checking we know that we cannot try single scattering so we will
579 // need to make an MSC step
580 // Indicate that we need to make and MSC step. We do not check if we can
581 // do it now i.e. if presafety>final_true_step_length so we let the
582 // fIsEverythingWasDone = false which indicates that we will perform
583 // MSC after transportation.
584 fIsMultipleSacettring = true;
585 // Init the first-real-step falg: it will indicate if we do the first
586 // non-single scattering step in this volume with this particle
587 fIsFirstRealStep = false;
588 // If previously the partcile was on boundary it was within skin as
589 // well. When it is not within skin anymore it has just left the skin
590 // so we make the first real MSC step with the particle.
591 if (fIsWasOnBoundary && !fIsInsideSkin) {
592 // reset the 'was on boundary' indicator flag
593 fIsWasOnBoundary = false;
594 fIsFirstRealStep = true;
595 }
596 // If this is the first-real msc step (the partcile has just left the
597 // skin) or this is the first step with the particle (was born or
598 // primary):
599 // - set the initial range that will be used later to limit its step
600 // (only in this volume, because after boundary crossing at the
601 // first-real MSC step we will reset)
602 // - don't let the partcile to cross the volume just in one step
603 if (firstStep || fIsFirstRealStep || rangeinit>1.e+20) {
604 rangeinit = currentRange;
605 // If geomlimit < geombig than the particle might reach the boundary
606 // along its initial direction before losing its energy (in this step)
607 // Otherwise we can be sure that the particle will lose it energy
608 // before reaching the boundary along a starigth line so there is no
609 // geometrical limit appalied. [However, tgeom is set only in the
610 // first or the first-real MSC step. After the first or first real
611 // MSC step the direction will change tgeom won't guaranty anything!
612 // But we will try to end up within skindepth from the boundary using
613 // the actual value of geomlimit(See later at step reduction close to
614 // boundary).]
615 if (geomlimit<geombig) {
616 // transfrom straight line distance to the boundary to real step
617 // length based on the mean values (using the prestep point
618 // first-transport mean free path i.e. no energy loss correction)
619 if ((1.-geomlimit/fLambda1)> 0.) {
620 geomlimit = -fLambda1*G4Log(1.-geomlimit/fLambda1);
621 }
622 // the 2-different case that could lead us here
623 if (firstStep) {
624 tgeom = 2.*geomlimit/facgeom;
625 } else {
626 tgeom = geomlimit/facgeom;
627 }
628 } else {
629 tgeom = geombig;
630 }
631 }
632 // True step length limit from range factor. Noteice, that the initial
633 // range is used that was set at the first step or first-real MSC step
634 // in this volume with this particle.
635 tlimit = facrange*rangeinit;
636 // Take the minimum of the true step length limits coming from
637 // geometrical constraint or range-factor limitation
638 tlimit = std::min(tlimit,tgeom);
639 // Step reduction close to boundary: we try to end up within skindepth
640 // from the boundary ( Notice: in case of mag. field it might not work
641 // because geomlimit is the straigth line distance to the boundary in
642 // the currect direction (if geomlimit<geombig) and mag. field can
643 // change the initial direction. So te particle might hit some boundary
644 // before in a different direction. However, here we restrict the true
645 // path length to this (straight line) lenght so the corresponding
646 // transport distance (straight line) will be even shorter than
647 // geomlimit-0.999*skindepth after the change of true->geom.
648 if (geomlimit<geombig) {
649 tlimit = std::min(tlimit, geomlimit-0.999*skindepth);
650 }
651 // randomize 1st step or 1st 'normal' step in volume
652 if (firstStep || fIsFirstRealStep) {
653 fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit());
654 } else {
655 fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit);
656 }
657 }
658 } else if (steppingAlgorithm==fUseSafetyPlus) { // THE ERROR_FREE stepping alg.
659 presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
660 geomlimit = presafety;
661 // Set skin depth = skin x elastic_mean_free_path
662 skindepth = skin*fLambda0;
663 // Check if we can try Single Scattering because we are within skindepth
664 // distance from/to a boundary OR the current minimum true-step-length is
665 // shorter than skindepth. NOTICE: the latest has only efficieny reasons
666 // because the MSC angular sampling is fine for any short steps but much
667 // faster to try single scattering in case of short steps.
668 if ((stepStatus==fGeomBoundary) || (presafety<skindepth) || (fTheTrueStepLenght<skindepth)) {
669 //Try single scattering:
670 // - sample distance to next single scattering interaction (sslimit)
671 // - compare to current minimum length
672 // == if sslimit is the shorter:
673 // - set the step length to sslimit
674 // - indicate that single scattering needs to be done
675 // == else : nothing to do
676 //- in both cases, the step length was very short so geometrical and
677 // true path length are the same
678 G4double sslimit = -1.*fLambda0*G4Log(G4UniformRand());
679 // compare to current minimum step length
680 if (sslimit<fTheTrueStepLenght) {
681 fTheTrueStepLenght = sslimit;
682 fIsSingleScattering = true;
683 }
684 // short step -> true step length equal to geometrical path length
685 fTheZPathLenght = fTheTrueStepLenght;
686 // Set taht everything is done in step-limit phase so no MSC call
687 // We will check if we need to perform the single-scattering angular
688 // sampling i.e. if single elastic scattering was the winer!
689 fIsEverythingWasDone = true;
690 } else {
691 // After checking we know that we cannot try single scattering so we will
692 // need to make an MSC step
693 // Indicate that we need to make and MSC step.
694 fIsMultipleSacettring = true;
695 fIsEverythingWasDone = true;
696 // limit from range factor
697 fTheTrueStepLenght = std::min(fTheTrueStepLenght, facrange*currentRange);
698 // never let the particle go further than the safety if we are out of the skin
699 // if we are here we are out of the skin, presafety > 0.
700 if (fTheTrueStepLenght>presafety) {
701 fTheTrueStepLenght = std::min(fTheTrueStepLenght, presafety);
702 }
703 // make sure that we are still within the aplicability of condensed histry model
704 // i.e. true step length is not longer than first transport mean free path.
705 // We schould take into account energy loss along 0.5x lambda_transport1
706 // step length as well. So let it 0.5 x lambda_transport1
707 fTheTrueStepLenght = std::min(fTheTrueStepLenght, fLambda1*0.5);
708 }
709 } else {
710 // This is the default stepping algorithm: the fastest but the least
711 // accurate that corresponds to fUseSafety in Urban model. Note, that GS
712 // model can handle any short steps so we do not need the minimum limits
713 //
714 // NO single scattering in case of skin or short steps (by defult the MSC
715 // model will be single or even no scattering in case of short steps
716 // compared to the elastic mean free path.)
717 //
718 // indicate that MSC needs to be done (always and always after transportation)
719 fIsMultipleSacettring = true;
720 if (stepStatus!=fGeomBoundary) {
721 presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
722 }
723 // Far from boundary-> in optimized mode do not sample dispalcement.
724 if ((distance<presafety) && (gIsOptimizationOn)) {
725 fIsNoDisplace = true;
726 } else {
727 // Urban like
728 if (firstStep || (stepStatus==fGeomBoundary) || rangeinit>1.e+20) {
729 rangeinit = currentRange;
730 fr = facrange;
731// We don't use this: we won't converge to the single scattering results with
732// decreasing range-factor.
733// rangeinit = std::max(rangeinit, fLambda1);
734// if(fLambda1 > lambdalimit) {
735// fr *= (0.75+0.25*fLambda1/lambdalimit);
736// }
737
738 }
739 //step limit
740 tlimit = std::max(fr*rangeinit, facsafety*presafety);
741 // first step randomization
742 if (firstStep || stepStatus==fGeomBoundary) {
743 fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit());
744 } else {
745 fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit);
746 }
747 }
748 }
749 //
750 // unset first-step
751 firstStep =false;
752 // performe single scattering, multiple scattering if this later can be done safely here
753 if (fIsEverythingWasDone) {
754 if (fIsSingleScattering) {
755 // sample single scattering
756 //G4double ekin = 0.5*(currentKinEnergy + GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple));
757 G4double lekin = G4Log(currentKinEnergy);
758 G4double pt2 = currentKinEnergy*(currentKinEnergy+2.0*CLHEP::electron_mass_c2);
759 G4double beta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
760 G4double cost = fGSTable->SingleScattering(1., fScrA, lekin, beta2, currentMaterialIndex);
761 // protection
762 if (cost<-1.) cost = -1.;
763 if (cost> 1.) cost = 1.;
764 // compute sint
765 G4double dum = 1.-cost;
766 G4double sint = std::sqrt(dum*(2.-dum));
767 G4double phi = CLHEP::twopi*G4UniformRand();
768 G4double sinPhi = std::sin(phi);
769 G4double cosPhi = std::cos(phi);
770 fTheNewDirection.set(sint*cosPhi,sint*sinPhi,cost);
771 } else if (fIsMultipleSacettring) {
772 // sample multiple scattering
773 SampleMSC(); // fTheZPathLenght, fTheDisplacementVector and fTheNewDirection will be set
774 } // and if single scattering but it was longer => nothing to do
775 } //else { do nothing here but after transportation
776 //
777 return ConvertTrueToGeom(fTheTrueStepLenght,currentMinimalStep);
778}
G4double G4Log(G4double x)
Definition: G4Log.hh:227
@ fUseSafetyPlus
@ fUseDistanceToBoundary
G4StepStatus
Definition: G4StepStatus.hh:40
@ fGeomBoundary
Definition: G4StepStatus.hh:43
@ fUndefined
Definition: G4StepStatus.hh:55
int G4int
Definition: G4Types.hh:85
#define G4UniformRand()
Definition: Randomize.hh:52
G4double GetLogKineticEnergy() const
G4double GetKineticEnergy() const
G4double GetTransportMeanFreePath(const G4ParticleDefinition *, G4double)
G4double SingleScattering(G4double lambdaval, G4double scra, G4double lekin, G4double beta2, G4int matindx)
G4double GetZeffective() const
const G4Material * GetMaterial() const
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:221
size_t GetIndex() const
Definition: G4Material.hh:255
G4StepPoint * GetPreStepPoint() const
const G4DynamicParticle * GetDynamicParticle() const
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
const G4Step * GetStep() const
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:468
G4double facrange
Definition: G4VMscModel.hh:199
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:296
G4double skin
Definition: G4VMscModel.hh:202
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.cc:188
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:209
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:286
G4double ComputeSafety(const G4ThreeVector &position, G4double limit=DBL_MAX)
Definition: G4VMscModel.hh:278
G4double facsafety
Definition: G4VMscModel.hh:201
G4double facgeom
Definition: G4VMscModel.hh:200

◆ ComputeTrueStepLength()

G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength ( G4double  geomStepLength)
overridevirtual

Implements G4VMscModel.

Definition at line 833 of file G4GoudsmitSaundersonMscModel.cc.

833 {
834 // init
835 fIsEndedUpOnBoundary = false;
836 // step defined other than transportation
837 if (geomStepLength==fTheZPathLenght) {
838 return fTheTrueStepLenght;
839 }
840 // else ::
841 // - set the flag that transportation was the winer so DoNothin in DOIT !!
842 // - convert geom -> true by using the mean value
843 fIsEndedUpOnBoundary = true; // OR LAST STEP
844 fTheZPathLenght = geomStepLength;
845 // was a short single scattering step
846 if (fIsEverythingWasDone && !fIsMultipleSacettring) {
847 fTheTrueStepLenght = geomStepLength;
848 return fTheTrueStepLenght;
849 }
850 // t = z for very small step
851 if (geomStepLength<tlimitminfix2) {
852 fTheTrueStepLenght = geomStepLength;
853 // recalculation
854 } else {
855 G4double tlength = geomStepLength;
856 if (geomStepLength>fLambda1*tausmall) {
857 if (par1< 0.) {
858 tlength = -fLambda1*G4Log(1.-geomStepLength/fLambda1) ;
859 } else {
860 if (par1*par3*geomStepLength<1.) {
861 G4Pow *g4calc = G4Pow::GetInstance();
862 tlength = (1.-g4calc->powA( 1.-par1*par3*geomStepLength,1./par3))/par1;
863 } else {
864 tlength = currentRange;
865 }
866 }
867 if (tlength<geomStepLength || tlength>fTheTrueStepLenght) {
868 tlength = geomStepLength;
869 }
870 }
871 fTheTrueStepLenght = tlength;
872 }
873 //
874 return fTheTrueStepLenght;
875}

◆ CrossSectionPerVolume()

G4double G4GoudsmitSaundersonMscModel::CrossSectionPerVolume ( const G4Material mat,
const G4ParticleDefinition ,
G4double  kineticEnergy,
G4double  cutEnergy = 0.0,
G4double  maxEnergy = DBL_MAX 
)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 309 of file G4GoudsmitSaundersonMscModel.cc.

313 {
314 G4double xsecTr1 = 0.; // cross section per volume i.e. macroscopic 1st transport cross section
315 //
316 fLambda0 = 0.0; // elastic mean free path
317 fLambda1 = 0.0; // first transport mean free path
318 fScrA = 0.0; // screening parameter
319 fG1 = 0.0; // first transport coef.
320 // use Moliere's screening (with Mott-corretion if it was requested)
321 G4double efEnergy = std::max(kineticEnergy, 10.*CLHEP::eV);
322 // total mometum square
323 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2);
324 // beta square
325 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2);
326 // current material index
327 G4int matindx = (G4int)mat->GetIndex();
328 // Moliere's b_c
329 G4double bc = fGSTable->GetMoliereBc(matindx);
330 // get the Mott-correcton factors if Mott-correcton was requested by the user
331 fMCtoScrA = 1.0;
332 fMCtoQ1 = 1.0;
333 fMCtoG2PerG1 = 1.0;
334 G4double scpCor = 1.0;
335 if (fIsUseMottCorrection) {
336 fGSTable->GetMottCorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1);
337 // ! no scattering power correction since the current couple is not set before this interface method is called
338 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
339 } else if (fIsUsePWACorrection) {
340 fPWACorrection->GetPWACorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1);
341 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
342 }
343 // screening parameter:
344 // - if Mott-corretioncorrection: the Screened-Rutherford times Mott-corretion DCS with this
345 // screening parameter gives back the (elsepa) PWA first transport cross section
346 // - if PWA correction: he Screened-Rutherford DCS with this screening parameter
347 // gives back the (elsepa) PWA first transport cross section
348 fScrA = fGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc)*fMCtoScrA;
349 // elastic mean free path in Geant4 internal lenght units: the neglected (1+screening parameter) term is corrected
350 // (if Mott-corretion: the corrected screening parameter is used for this (1+A) correction + Moliere b_c is also
351 // corrected with the screening parameter correction)
352 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scpCor;
353 // first transport coefficient (if Mott-corretion: the corrected screening parameter is used (it will be fully
354 // consistent with the one used during the pre-computation of the Mott-correted GS angular distributions))
355 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/fScrA+1.0)-1.0);
356 // first transport mean free path
357 fLambda1 = fLambda0/fG1;
358 xsecTr1 = 1./fLambda1;
359 return xsecTr1;
360}
void GetPWACorrectionFactors(G4double logekin, G4double beta2, G4int matindx, G4double &corToScr, G4double &corToQ1, G4double &corToG2PerG1)
void GetMottCorrectionFactors(G4double logekin, G4double beta2, G4int matindx, G4double &mcToScr, G4double &mcToQ1, G4double &mcToG2PerG1)
G4double GetMoliereBc(G4int matindx)
G4double GetMoliereXc2(G4int matindx)

◆ GetGSTable()

G4GoudsmitSaundersonTable * G4GoudsmitSaundersonMscModel::GetGSTable ( )
inline

Definition at line 176 of file G4GoudsmitSaundersonMscModel.hh.

176{ return fGSTable; }

Referenced by InitialiseLocal().

◆ GetOptionMottCorrection()

G4bool G4GoudsmitSaundersonMscModel::GetOptionMottCorrection ( ) const
inline

Definition at line 174 of file G4GoudsmitSaundersonMscModel.hh.

174{ return fIsUseMottCorrection; }

Referenced by InitialiseLocal().

◆ GetOptionPWACorrection()

G4bool G4GoudsmitSaundersonMscModel::GetOptionPWACorrection ( ) const
inline

Definition at line 170 of file G4GoudsmitSaundersonMscModel.hh.

170{ return fIsUsePWACorrection; }

Referenced by InitialiseLocal().

◆ GetPWACorrection()

G4GSPWACorrections * G4GoudsmitSaundersonMscModel::GetPWACorrection ( )
inline

Definition at line 178 of file G4GoudsmitSaundersonMscModel.hh.

178{ return fPWACorrection; }

Referenced by InitialiseLocal().

◆ GetTransportMeanFreePath()

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

Definition at line 365 of file G4GoudsmitSaundersonMscModel.cc.

366 {
367 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
368 G4double efEnergy = kineticEnergy;
369 //
370 const G4Material* mat = currentCouple->GetMaterial();
371 //
372 fLambda0 = 0.0; // elastic mean free path
373 fLambda1 = 0.0; // first transport mean free path
374 fScrA = 0.0; // screening parameter
375 fG1 = 0.0; // first transport coef.
376
377 // use Moliere's screening (with Mott-corretion if it was requested)
378 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.*CLHEP::eV;
379 // total mometum square
380 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2);
381 // beta square
382 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2);
383 // current material index
384 G4int matindx = (G4int)mat->GetIndex();
385 // Moliere's b_c
386 G4double bc = fGSTable->GetMoliereBc(matindx);
387 // get the Mott-correcton factors if Mott-correcton was requested by the user
388 fMCtoScrA = 1.0;
389 fMCtoQ1 = 1.0;
390 fMCtoG2PerG1 = 1.0;
391 G4double scpCor = 1.0;
392 if (fIsUseMottCorrection) {
393 fGSTable->GetMottCorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1);
394 scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
395 } else if (fIsUsePWACorrection) {
396 fPWACorrection->GetPWACorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1);
397 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
398 }
399 // screening parameter:
400 // - if Mott-corretioncorrection: the Screened-Rutherford times Mott-corretion DCS with this
401 // screening parameter gives back the (elsepa) PWA first transport cross section
402 // - if PWA correction: he Screened-Rutherford DCS with this screening parameter
403 // gives back the (elsepa) PWA first transport cross section
404 fScrA = fGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc)*fMCtoScrA;
405 // elastic mean free path in Geant4 internal lenght units: the neglected (1+screening parameter) term is corrected
406 // (if Mott-corretion: the corrected screening parameter is used for this (1+A) correction + Moliere b_c is also
407 // corrected with the screening parameter correction)
408 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scpCor;
409 // first transport coefficient (if Mott-corretion: the corrected screening parameter is used (it will be fully
410 // consistent with the one used during the pre-computation of the Mott-correted GS angular distributions))
411 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/fScrA+1.0)-1.0);
412 // first transport mean free path
413 fLambda1 = fLambda0/fG1;
414
415 return fLambda1;
416}
G4double ComputeScatteringPowerCorrection(const G4MaterialCutsCouple *matcut, G4double ekin)

Referenced by ComputeTruePathLengthLimit(), and SampleMSC().

◆ Initialise()

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

Implements G4VEmModel.

Definition at line 251 of file G4GoudsmitSaundersonMscModel.cc.

251 {
252 SetParticle(p);
254 // -create GoudsmitSaundersonTable and init its Mott-correction member if
255 // Mott-correction was required
256 if (IsMaster()) {
257 // get the Mott-correction flag from EmParameters
258 if (G4EmParameters::Instance()->UseMottCorrection()) {
259 fIsUseMottCorrection = true;
260 }
261 // Mott-correction includes other way of PWA x-section corrections so deactivate it even if it was true
262 // when Mott-correction is activated by the user
263 if (fIsUseMottCorrection) {
264 fIsUsePWACorrection = false;
265 }
266 // clear GS-table
267 if (fGSTable) {
268 delete fGSTable;
269 fGSTable = nullptr;
270 }
271 // clear PWA corrections table if any
272 if (fPWACorrection) {
273 delete fPWACorrection;
274 fPWACorrection = nullptr;
275 }
276 // create GS-table
277 G4bool isElectron = true;
278 if (p->GetPDGCharge()>0.) {
279 isElectron = false;
280 }
281 fGSTable = new G4GoudsmitSaundersonTable(isElectron);
282 // G4GSTable will be initialised:
283 // - Screened-Rutherford DCS based GS angular distributions will be loaded only if they are not there yet
284 // - Mott-correction will be initialised if Mott-correction was requested to be used
285 fGSTable->SetOptionMottCorrection(fIsUseMottCorrection);
286 // - set PWA correction (correction to integrated quantites from Dirac-PWA)
287 fGSTable->SetOptionPWACorrection(fIsUsePWACorrection);
288 // init
290 // create PWA corrections table if it was requested (and not disactivated because active Mott-correction)
291 if (fIsUsePWACorrection) {
292 fPWACorrection = new G4GSPWACorrections(isElectron);
293 fPWACorrection->Initialise();
294 }
295 }
296 fParticleChange = GetParticleChangeForMSC(p);
297}
bool G4bool
Definition: G4Types.hh:86
static G4EmParameters * Instance()
void Initialise(G4double lownergylimit, G4double highenergylimit)
G4double GetPDGCharge() const
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=nullptr)
Definition: G4VMscModel.cc:77
void InitialiseParameters(const G4ParticleDefinition *)
Definition: G4VMscModel.cc:115
G4bool isElectron(G4int ityp)

◆ InitialiseLocal()

void G4GoudsmitSaundersonMscModel::InitialiseLocal ( const G4ParticleDefinition p,
G4VEmModel masterModel 
)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 300 of file G4GoudsmitSaundersonMscModel.cc.

300 {
301 fGSTable = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetGSTable();
302 fIsUseMottCorrection = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetOptionMottCorrection();
303 fIsUsePWACorrection = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetOptionPWACorrection();
304 fPWACorrection = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetPWACorrection();
305}
G4GoudsmitSaundersonTable * GetGSTable()

◆ operator=()

G4GoudsmitSaundersonMscModel & G4GoudsmitSaundersonMscModel::operator= ( const G4GoudsmitSaundersonMscModel right)
delete

◆ SampleMSC()

void G4GoudsmitSaundersonMscModel::SampleMSC ( )

Definition at line 923 of file G4GoudsmitSaundersonMscModel.cc.

923 {
924 fIsNoScatteringInMSC = false;
925 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
926 G4double kineticEnergy = currentKinEnergy;
927 //
928 // Energy loss correction: 2 version
929 G4double eloss = 0.0;
930// if (fTheTrueStepLenght > currentRange*dtrl) {
931 eloss = kineticEnergy - GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple);
932// } else {
933// eloss = fTheTrueStepLenght*GetDEDX(particle,kineticEnergy,currentCouple);
934// }
935
936 G4double tau = 0.;// = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy
937 G4double tau2 = 0.;// = tau*tau;
938 G4double eps0 = 0.;// = eloss/kineticEnergy0; // energy loss fraction to the begin step energy
939 G4double epsm = 0.;// = eloss/kineticEnergy; // energy loss fraction to the mean step energy
940
941 // - init.
942 G4double efEnergy = kineticEnergy;
943 G4double efStep = fTheTrueStepLenght;
944
945 G4double kineticEnergy0 = kineticEnergy;
946 if (gIsUseAccurate) { // - use accurate energy loss correction
947 kineticEnergy -= 0.5*eloss; // mean energy along the full step
948 // other parameters for energy loss corrections
949 tau = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy
950 tau2 = tau*tau;
951 eps0 = eloss/kineticEnergy0; // energy loss fraction to the begin step energy
952 epsm = eloss/kineticEnergy; // energy loss fraction to the mean step energy
953
954 efEnergy = kineticEnergy * (1.-epsm*epsm*(6.+10.*tau+5.*tau2)/(24.*tau2+48.*tau+72.));
955 G4double dum = 0.166666*(4.+tau*(6.+tau*(7.+tau*(4.+tau))))*(epsm/((tau+1.)*(tau+2.)))*(epsm/((tau+1.)*(tau+2.)));
956 efStep = fTheTrueStepLenght*(1.-dum);
957 } else { // - take only mean energy
958 kineticEnergy -= 0.5*eloss; // mean energy along the full step
959 efEnergy = kineticEnergy;
960 G4double factor = 1./(1.+0.9784671*kineticEnergy); //0.9784671 = 1/(2*m_e)
961 eps0 = eloss/kineticEnergy0;
962 epsm = eps0/(1.-0.5*eps0);
963 G4double temp = 0.3*(1 -factor*(1.-0.333333*factor))*eps0*eps0;
964 efStep = fTheTrueStepLenght*(1.+temp);
965 }
966 //
967 // compute elastic mfp, first transport mfp, screening parameter, and G1 (with Mott-correction
968 // if it was requested by the user)
969 fLambda1 = GetTransportMeanFreePath(particle, efEnergy);
970 // s/lambda_el
971 G4double lambdan=0.;
972 if (fLambda0>0.0) {
973 lambdan=efStep/fLambda0;
974 }
975 if (lambdan<=1.0e-12) {
976 if (fIsEverythingWasDone) {
977 fTheZPathLenght = fTheTrueStepLenght;
978 }
979 fIsNoScatteringInMSC = true;
980 return;
981 }
982 // first moment: 2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
983 G4double Qn1 = lambdan *fG1;
984 // sample scattering angles
985 // new direction, relative to the orriginal one is in {uss,vss,wss}
986 G4double cosTheta1 = 1.0, sinTheta1 = 0.0, cosTheta2 = 1.0, sinTheta2 = 0.0;
987 G4double cosPhi1 = 1.0, sinPhi1 = 0.0, cosPhi2 = 1.0, sinPhi2 = 0.0;
988 G4double uss = 0.0, vss = 0.0, wss = 1.0;
989 G4double x_coord = 0.0, y_coord = 0.0, z_coord = 1.0;
990 G4double u2 = 0.0, v2 = 0.0;
991 // if we are above the upper grid limit with lambdaxG1=true-length/first-trans-mfp
992 // => izotropic distribution: lambG1_max =7.992 but set it to 7
993 if (0.5*Qn1 > 7.0){
994 cosTheta1 = 1.-2.*G4UniformRand();
995 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+cosTheta1));
996 cosTheta2 = 1.-2.*G4UniformRand();
997 sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+cosTheta2));
998 } else {
999 // sample 2 scattering cost1, sint1, cost2 and sint2 for half path
1000 G4double lekin = G4Log(efEnergy);
1001 G4double pt2 = efEnergy*(efEnergy+2.0*CLHEP::electron_mass_c2);
1002 G4double beta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
1003 // backup GS angular dtr pointer (kinetic energy and delta index in case of Mott-correction)
1004 // if the first was an msc sampling (the same will be used if the second is also an msc step)
1006 G4int mcEkinIdx = -1;
1007 G4int mcDeltIdx = -1;
1008 G4double transfPar = 0.;
1009 G4bool isMsc = fGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta1, sinTheta1, lekin, beta2,
1010 currentMaterialIndex, &gsDtr, mcEkinIdx, mcDeltIdx, transfPar,
1011 true);
1012 fGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta2, sinTheta2, lekin, beta2,
1013 currentMaterialIndex, &gsDtr, mcEkinIdx, mcDeltIdx, transfPar, !isMsc);
1014 if (cosTheta1+cosTheta2==2.) { // no scattering happened
1015 if (fIsEverythingWasDone)
1016 fTheZPathLenght = fTheTrueStepLenght;
1017 fIsNoScatteringInMSC = true;
1018 return;
1019 }
1020 }
1021 // sample 2 azimuthal angles
1022 G4double phi1 = CLHEP::twopi*G4UniformRand();
1023 sinPhi1 = std::sin(phi1);
1024 cosPhi1 = std::cos(phi1);
1025 G4double phi2 = CLHEP::twopi*G4UniformRand();
1026 sinPhi2 = std::sin(phi2);
1027 cosPhi2 = std::cos(phi2);
1028
1029 // compute final direction realtive to z-dir
1030 u2 = sinTheta2*cosPhi2;
1031 v2 = sinTheta2*sinPhi2;
1032 G4double u2p = cosTheta1*u2 + sinTheta1*cosTheta2;
1033 uss = u2p*cosPhi1 - v2*sinPhi1;
1034 vss = u2p*sinPhi1 + v2*cosPhi1;
1035 wss = cosTheta1*cosTheta2 - sinTheta1*u2;
1036
1037 // set new direction (is scattering frame)
1038 fTheNewDirection.set(uss,vss,wss);
1039
1040 // set the fTheZPathLenght if we don't sample displacement and
1041 // we should do everything at the step-limit-phase before we return
1042 if(fIsNoDisplace && fIsEverythingWasDone)
1043 fTheZPathLenght = fTheTrueStepLenght;
1044
1045 // in optimized-mode if the current-safety > current-range we do not use dispalcement
1046 if(fIsNoDisplace)
1047 return;
1048
1049 //////////////////////////////////////////////////////////////////////
1050 // Compute final position
1051 Qn1 *= fMCtoQ1;
1052 if (gIsUseAccurate) {
1053 // correction parameter
1054 G4double par =1.;
1055 if(Qn1<0.7) par = 1.;
1056 else if (Qn1<7.0) par = -0.031376*Qn1+1.01356;
1057 else par = 0.79;
1058
1059 // Moments with energy loss correction
1060 // --first the uncorrected (for energy loss) values of gamma, eta, a1=a2=0.5*(1-eta), delta
1061 // gamma = G_2/G_1 based on G2 computed from A by using the Wentzel DCS form of G2
1062 G4double loga = G4Log(1.0+1.0/fScrA);
1063 G4double gamma = 6.0*fScrA*(1.0 + fScrA)*(loga*(1.0 + 2.0*fScrA) - 2.0)/fG1;
1064 gamma *= fMCtoG2PerG1;
1065 // sample eta from p(eta)=2*eta i.e. P(eta) = eta_square ;-> P(eta) = rand --> eta = sqrt(rand)
1066 G4double eta = std::sqrt(G4UniformRand());
1067 G4double eta1 = 0.5*(1 - eta); // used more than once
1068 // 0.5 +sqrt(6)/6 = 0.9082483;
1069 // 1/(4*sqrt(6)) = 0.1020621;
1070 // (4-sqrt(6)/(24*sqrt(6))) = 0.026374715
1071 // delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1 without energy loss cor.
1072 G4double delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1;
1073
1074 // compute alpha1 and alpha2 for energy loss correction
1075 G4double temp1 = 2.0 + tau;
1076 G4double temp = (2.0+tau*temp1)/((tau+1.0)*temp1);
1077 //Take logarithmic dependence
1078 temp = temp - (tau+1.0)/((tau+2.0)*(loga*(1.0+fScrA)-1.0));
1079 temp = temp * epsm;
1080 temp1 = 1.0 - temp;
1081 delta = delta + 0.40824829*(eps0*(tau+1.0)/((tau+2.0)*
1082 (loga*(1.0+fScrA)-1.0)*(loga*(1.0+2.0*fScrA)-2.0)) - 0.25*temp*temp);
1083 G4double b = eta*delta;
1084 G4double c = eta*(1.0-delta);
1085
1086 //Calculate transport direction cosines:
1087 // ut,vt,wt is the final position divided by the true step length
1088 G4double w1v2 = cosTheta1*v2;
1089 G4double ut = b*sinTheta1*cosPhi1 + c*(cosPhi1*u2 - sinPhi1*w1v2) + eta1*uss*temp1;
1090 G4double vt = b*sinTheta1*sinPhi1 + c*(sinPhi1*u2 + cosPhi1*w1v2) + eta1*vss*temp1;
1091 G4double wt = eta1*(1+temp) + b*cosTheta1 + c*cosTheta2 + eta1*wss*temp1;
1092
1093 // long step correction
1094 ut *=par;
1095 vt *=par;
1096 wt *=par;
1097
1098 // final position relative to the pre-step point in the scattering frame
1099 // ut = x_f/s so needs to multiply by s
1100 x_coord = ut*fTheTrueStepLenght;
1101 y_coord = vt*fTheTrueStepLenght;
1102 z_coord = wt*fTheTrueStepLenght;
1103
1104 if(fIsEverythingWasDone){
1105 // We sample in the step limit so set fTheZPathLenght = transportDistance
1106 // and lateral displacement (x_coord,y_coord,z_coord-transportDistance)
1107 //Calculate transport distance
1108 G4double transportDistance = std::sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
1109 // protection
1110 if(transportDistance>fTheTrueStepLenght)
1111 transportDistance = fTheTrueStepLenght;
1112 fTheZPathLenght = transportDistance;
1113 }
1114 // else:: we sample in the DoIt so
1115 // the fTheZPathLenght was already set and was taken as transport along zet
1116 fTheDisplacementVector.set(x_coord,y_coord,z_coord-fTheZPathLenght);
1117 } else {
1118 // compute zz = <z>/tPathLength
1119 // s -> true-path-length
1120 // z -> geom-path-length:: when PRESTA is used z =(def.) <z>
1121 // r -> lateral displacement = s/2 sin(theta) => x_f = r cos(phi); y_f = r sin(phi)
1122 G4double zz = 0.0;
1123 if(fIsEverythingWasDone){
1124 // We sample in the step limit so set fTheZPathLenght = transportDistance
1125 // and lateral displacement (x_coord,y_coord,z_coord-transportDistance)
1126 if(Qn1<0.1) { // use 3-order Taylor approximation of (1-exp(-x))/x around x=0
1127 zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666667 - 0.041666667*Qn1)); // 1/6 =0.166..7 ; 1/24=0.041..
1128 } else {
1129 zz = (1.-G4Exp(-Qn1))/Qn1;
1130 }
1131 } else {
1132 // we sample in the DoIt so
1133 // the fTheZPathLenght was already set and was taken as transport along zet
1134 zz = fTheZPathLenght/fTheTrueStepLenght;
1135 }
1136
1137 G4double rr = (1.-zz*zz)/(1.-wss*wss); // s^2 >= <z>^2+r^2 :: where r^2 = s^2/4 sin^2(theta)
1138 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
1139 G4double rperp = fTheTrueStepLenght*std::sqrt(rr); // this is r/sint
1140 x_coord = rperp*uss;
1141 y_coord = rperp*vss;
1142 z_coord = zz*fTheTrueStepLenght;
1143
1144 if(fIsEverythingWasDone){
1145 G4double transportDistance = std::sqrt(x_coord*x_coord + y_coord*y_coord + z_coord*z_coord);
1146 fTheZPathLenght = transportDistance;
1147 }
1148
1149 fTheDisplacementVector.set(x_coord,y_coord,z_coord- fTheZPathLenght);
1150 }
1151}
G4bool Sampling(G4double lambdaval, G4double qval, G4double scra, G4double &cost, G4double &sint, G4double lekin, G4double beta2, G4int matindx, GSMSCAngularDtr **gsDtr, G4int &mcekini, G4int &mcdelti, G4double &transfPar, G4bool isfirst)

Referenced by ComputeTruePathLengthLimit(), and SampleScattering().

◆ SampleScattering()

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

Implements G4VMscModel.

Definition at line 878 of file G4GoudsmitSaundersonMscModel.cc.

878 {
879 if (steppingAlgorithm==fUseDistanceToBoundary && fIsEverythingWasDone && fIsSingleScattering) {
880 // single scattering was and scattering happend
881 fTheNewDirection.rotateUz(oldDirection);
882 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
883 return fTheDisplacementVector;
884 } else if (steppingAlgorithm==fUseSafetyPlus) { // error-free stepping
885 if (fIsEndedUpOnBoundary) { // do nothing on the boundary
886 return fTheDisplacementVector;
887 } else if (fIsEverythingWasDone) { // evrything is done if not optimizations case !!!
888 // check single scattering and see if it happened
889 if (fIsSingleScattering) {
890 fTheNewDirection.rotateUz(oldDirection);
891 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
892 return fTheDisplacementVector;
893 }
894 // check if multiple scattering happened and do things only if scattering was really happening
895 if (fIsMultipleSacettring && !fIsNoScatteringInMSC) {
896 fTheNewDirection.rotateUz(oldDirection);
897 fTheDisplacementVector.rotateUz(oldDirection);
898 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
899 }
900 // The only thing that could happen if we are here (fUseSafety and fIsEverythingWasDone)
901 // is that single scattering was tried but did not win so scattering did not happen.
902 // So no displacement and no scattering
903 return fTheDisplacementVector;
904 }
905 //
906 // The only thing that could still happen with fUseSafetyPlus is that we are in the
907 // optimization branch: so sample MSC angle here (no displacement)
908 }
909 //else MSC needs to be done here
910 SampleMSC();
911 if (!fIsNoScatteringInMSC) {
912 fTheNewDirection.rotateUz(oldDirection);
913 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
914 if (!fIsNoDisplace) {
915 fTheDisplacementVector.rotateUz(oldDirection);
916 }
917 }
918 //
919 return fTheDisplacementVector;
920}
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:33
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)

◆ SetOptionMottCorrection()

void G4GoudsmitSaundersonMscModel::SetOptionMottCorrection ( G4bool  opt)
inline

Definition at line 172 of file G4GoudsmitSaundersonMscModel.hh.

172{ fIsUseMottCorrection = opt; }

◆ SetOptionPWACorrection()

void G4GoudsmitSaundersonMscModel::SetOptionPWACorrection ( G4bool  opt)
inline

Definition at line 168 of file G4GoudsmitSaundersonMscModel.hh.

168{ fIsUsePWACorrection = opt; }

◆ StartTracking()

void G4GoudsmitSaundersonMscModel::StartTracking ( G4Track track)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 461 of file G4GoudsmitSaundersonMscModel.cc.

461 {
462 SetParticle(track->GetDynamicParticle()->GetDefinition());
463 firstStep = true;
464 tlimit = tgeom = rangeinit = geombig;
465 rangeinit = 1.e+21;
466}
G4ParticleDefinition * GetDefinition() const

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