Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4WentzelVIModel Class Reference

#include <G4WentzelVIModel.hh>

Inheritance diagram for G4WentzelVIModel:
Collaboration diagram for G4WentzelVIModel:

Public Member Functions

 G4WentzelVIModel (G4bool comb=true, const G4String &nam="WentzelVIUni")
 
virtual ~G4WentzelVIModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel) override
 
virtual void StartTracking (G4Track *) override
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=DBL_MAX, G4double emax=DBL_MAX) override
 
virtual G4ThreeVectorSampleScattering (const G4ThreeVector &, G4double safety) override
 
virtual G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep) override
 
virtual G4double ComputeGeomPathLength (G4double truePathLength) override
 
virtual G4double ComputeTrueStepLength (G4double geomStepLength) override
 
void SetFixedCut (G4double)
 
G4double GetFixedCut () const
 
G4WentzelOKandVIxSectionGetWVICrossSection ()
 
void SetUseSecondMoment (G4bool)
 
G4bool UseSecondMoment () const
 
G4PhysicsTableGetSecondMomentTable ()
 
G4double SecondMoment (const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double kineticEnergy)
 
void SetSingleScatteringFactor (G4double)
 
- 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 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 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 SetFluctuationFlag (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)
 

Protected Member Functions

void DefineMaterial (const G4MaterialCutsCouple *)
 
- 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

G4WentzelOKandVIxSectionwokvi
 
G4double tlimitminfix
 
G4double ssFactor
 
G4double invssFactor
 
G4double preKinEnergy
 
G4double tPathLength
 
G4double zPathLength
 
G4double lambdaeff
 
G4double currentRange
 
G4double cosTetMaxNuc
 
G4int currentMaterialIndex
 
const G4MaterialCutsCouplecurrentCouple
 
const G4MaterialcurrentMaterial
 
const G4ParticleDefinitionparticle
 
G4bool singleScatteringMode
 
- 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
 
G4bool lossFlucFlag
 

Additional Inherited Members

- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 66 of file G4WentzelVIModel.hh.

Constructor & Destructor Documentation

G4WentzelVIModel::G4WentzelVIModel ( G4bool  comb = true,
const G4String nam = "WentzelVIUni" 
)
explicit

Definition at line 74 of file G4WentzelVIModel.cc.

74  :
75  G4VMscModel(nam),
76  ssFactor(1.05),
77  invssFactor(1.0),
78  currentCouple(nullptr),
79  singleScatteringMode(false),
80  cosThetaMin(1.0),
81  cosThetaMax(-1.0),
82  fSecondMoments(nullptr),
83  idx2(0),
84  numlimit(0.1),
85  isCombined(combined),
86  useSecondMoment(false)
87 {
89  invsqrt12 = 1./sqrt(12.);
90  tlimitminfix = 1.e-6*mm;
91  lowEnergyLimit = 1.0*eV;
92  particle = 0;
93  nelments = 5;
94  xsecn.resize(nelments);
95  prob.resize(nelments);
96  wokvi = new G4WentzelOKandVIxSection(combined);
97  fixedCut = -1.0;
98 
99  minNCollisions = 10;
100 
101  preKinEnergy = effKinEnergy = tPathLength = zPathLength = lambdaeff
102  = currentRange = xtsec = cosTetMaxNuc = 0.0;
104 
105  fParticleChange = nullptr;
106  currentCuts = nullptr;
107  currentMaterial = nullptr;
108 }
void SetSingleScatteringFactor(G4double)
G4WentzelOKandVIxSection * wokvi
static constexpr double mm
Definition: G4SIunits.hh:115
const G4MaterialCutsCouple * currentCouple
const G4Material * currentMaterial
static constexpr double eV
Definition: G4SIunits.hh:215
const G4ParticleDefinition * particle
G4VMscModel(const G4String &nam)
Definition: G4VMscModel.cc:60

Here is the call graph for this function:

G4WentzelVIModel::~G4WentzelVIModel ( )
virtual

Definition at line 112 of file G4WentzelVIModel.cc.

113 {
114  delete wokvi;
115  if(fSecondMoments && IsMaster()) {
116  delete fSecondMoments;
117  fSecondMoments = 0;
118  }
119 }
G4WentzelOKandVIxSection * wokvi
G4bool IsMaster() const
Definition: G4VEmModel.hh:725

Here is the call graph for this function:

Member Function Documentation

G4double G4WentzelVIModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition p,
G4double  KineticEnergy,
G4double  AtomicNumber,
G4double  AtomicWeight = 0.,
G4double  cut = DBL_MAX,
G4double  emax = DBL_MAX 
)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 210 of file G4WentzelVIModel.cc.

215 {
216  G4double cross = 0.0;
217  if(p != particle) { SetupParticle(p); }
218  if(kinEnergy < lowEnergyLimit) { return cross; }
219  if(!CurrentCouple()) {
220  G4Exception("G4WentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
221  FatalException, " G4MaterialCutsCouple is not defined");
222  return 0.0;
223  }
226  if(cosTetMaxNuc < 1.0) {
227  G4double cut = cutEnergy;
228  if(fixedCut > 0.0) { cut = fixedCut; }
229  G4double cost = wokvi->SetupTarget(G4lrint(Z), cut);
231  /*
232  if(p->GetParticleName() == "e-")
233  G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= "<<kinEnergy
234  << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
235  << " " << particle->GetParticleName() << G4endl;
236  */
237  }
238  return cross;
239 }
G4WentzelOKandVIxSection * wokvi
void DefineMaterial(const G4MaterialCutsCouple *)
const G4Material * currentMaterial
G4double SetupTarget(G4int Z, G4double cut=DBL_MAX)
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:454
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
int G4lrint(double ad)
Definition: templates.hh:163
const G4ParticleDefinition * particle
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4WentzelVIModel::ComputeGeomPathLength ( G4double  truePathLength)
overridevirtual

Reimplemented from G4VMscModel.

Definition at line 347 of file G4WentzelVIModel.cc.

348 {
349  zPathLength = tPathLength = truelength;
350 
351  // small step use only single scattering
352  cosThetaMin = 1.0;
353  ComputeTransportXSectionPerVolume(cosThetaMin);
354  //G4cout << "xtsec= " << xtsec << " Nav= "
355  // << zPathLength*xtsec << G4endl;
356  if(0.0 >= lambdaeff || G4int(zPathLength*xtsec) < minNCollisions) {
357  singleScatteringMode = true;
358  lambdaeff = DBL_MAX;
359 
360  } else {
361  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
362  // << " Leff= " << lambdaeff << G4endl;
363  // small step
364  if(tPathLength < numlimit*lambdaeff) {
366  zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
367 
368  // medium step
369  } else {
370  G4double e1 = 0.0;
371  if(currentRange > tPathLength) {
373  }
374  effKinEnergy = 0.5*(e1 + preKinEnergy);
376  lambdaeff = GetTransportMeanFreePath(particle, effKinEnergy);
377  //G4cout << " tLength= "<< tPathLength<< " Leff= " << lambdaeff << G4endl;
379  if(tPathLength*numlimit < lambdaeff) {
380  zPathLength *= (1.0 - G4Exp(-tPathLength/lambdaeff));
381  }
382  }
383  }
384  //G4cout << "Comp.geom: zLength= "<<zPathLength<<" tLength= "
385  // << tPathLength<< " Leff= " << lambdaeff << G4endl;
386  return zPathLength;
387 }
G4WentzelOKandVIxSection * wokvi
const G4MaterialCutsCouple * currentCouple
const G4Material * currentMaterial
int G4int
Definition: G4Types.hh:78
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:303
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:340
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
const G4ParticleDefinition * particle
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

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

Reimplemented from G4VMscModel.

Reimplemented in G4LowEWentzelVIModel.

Definition at line 250 of file G4WentzelVIModel.cc.

253 {
254  G4double tlimit = currentMinimalStep;
255  const G4DynamicParticle* dp = track.GetDynamicParticle();
256  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
257  G4StepStatus stepStatus = sp->GetStepStatus();
258  singleScatteringMode = false;
259 
260  //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
261  // << stepStatus << " " << track.GetDefinition()->GetParticleName()
262  // << G4endl;
263 
264  // initialisation for each step, lambda may be computed from scratch
266  effKinEnergy = preKinEnergy;
271 
272  //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
273  // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
274 
275  // extra check for abnormal situation
276  // this check needed to run MSC with eIoni and eBrem inactivated
277  if(tlimit > currentRange) { tlimit = currentRange; }
278 
279  // stop here if small range particle
280  if(tlimit < tlimitminfix) {
281  return ConvertTrueToGeom(tlimit, currentMinimalStep);
282  }
283 
284  // pre step
285  G4double presafety = sp->GetSafety();
286  // far from geometry boundary
287  if(currentRange < presafety) {
288  return ConvertTrueToGeom(tlimit, currentMinimalStep);
289  }
290 
291  // compute presafety again if presafety <= 0 and no boundary
292  // i.e. when it is needed for optimization purposes
293  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
294  presafety = ComputeSafety(sp->GetPosition(), tlimit);
295  if(currentRange < presafety) {
296  return ConvertTrueToGeom(tlimit, currentMinimalStep);
297  }
298  }
299  /*
300  G4cout << "e(MeV)= " << preKinEnergy/MeV
301  << " " << particle->GetParticleName()
302  << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
303  << " R(mm)= " <<currentRange/mm
304  << " L0(mm^-1)= " << lambdaeff*mm
305  << G4endl;
306  */
307  // natural limit for high energy
309  (1.0 - cosTetMaxNuc)*lambdaeff*invssFactor);
310 
311  // low-energy e-
312  if(cosThetaMax > cosTetMaxNuc) {
313  rlimit = std::min(rlimit, facsafety*presafety);
314  }
315 
316  // cut correction
318  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= "
319  // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax
320  //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
321  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
322 
323  tlimit = std::min(tlimit, rlimit);
324  tlimit = std::max(tlimit, tlimitminfix);
325 
326  // step limit in infinite media
327  tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
328 
329  //compute geomlimit and force few steps within a volume
331  && stepStatus == fGeomBoundary) {
332 
333  G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
334  tlimit = std::min(tlimit, geomlimit/facgeom);
335  }
336  /*
337  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
338  << " L0= " << lambdaeff << " R= " << currentRange
339  << " tlimit= " << tlimit
340  << " currentMinimalStep= " << currentMinimalStep << G4endl;
341  */
342  return ConvertTrueToGeom(tlimit, currentMinimalStep);
343 }
G4double facgeom
Definition: G4VMscModel.hh:176
G4WentzelOKandVIxSection * wokvi
G4double GetKineticEnergy() const
const G4DynamicParticle * GetDynamicParticle() const
G4double facrange
Definition: G4VMscModel.hh:175
G4double GetProductionCut(G4int index) const
G4StepStatus GetStepStatus() const
void DefineMaterial(const G4MaterialCutsCouple *)
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:245
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
const G4MaterialCutsCouple * currentCouple
const G4Material * currentMaterial
const G4Step * GetStep() const
G4double ComputeSafety(const G4ThreeVector &position, G4double limit=DBL_MAX)
Definition: G4VMscModel.hh:237
G4StepStatus
Definition: G4StepStatus.hh:49
G4StepPoint * GetPreStepPoint() const
const G4ThreeVector & GetPosition() const
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:283
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:255
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:340
G4double GetRadlen() const
Definition: G4Material.hh:220
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double facsafety
Definition: G4VMscModel.hh:177
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
const G4ParticleDefinition * particle
G4double GetSafety() const
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:185
double G4double
Definition: G4Types.hh:76
G4ProductionCuts * GetProductionCuts() const

Here is the call graph for this function:

G4double G4WentzelVIModel::ComputeTrueStepLength ( G4double  geomStepLength)
overridevirtual

Reimplemented from G4VMscModel.

Definition at line 391 of file G4WentzelVIModel.cc.

392 {
393  // initialisation of single scattering x-section
394  /*
395  G4cout << "ComputeTrueStepLength: Step= " << geomStepLength
396  << " geomL= " << zPathLength
397  << " Lambda= " << lambdaeff
398  << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
399  */
401  zPathLength = tPathLength = geomStepLength;
402 
403  } else {
404 
405  // step defined by transportation
406  // change both geom and true step lengths
407  if(geomStepLength < zPathLength) {
408 
409  // single scattering
410  if(G4int(geomStepLength*xtsec) < minNCollisions) {
411  zPathLength = tPathLength = geomStepLength;
412  lambdaeff = DBL_MAX;
413  singleScatteringMode = true;
414 
415  // multiple scattering
416  } else {
417  // small step
418  if(geomStepLength < numlimit*lambdaeff) {
419  G4double tau = geomStepLength/lambdaeff;
420  tPathLength = geomStepLength*(1.0 + 0.5*tau + tau*tau/3.0);
421 
422  // energy correction for a big step
423  } else {
424  tPathLength *= geomStepLength/zPathLength;
425  G4double e1 = 0.0;
426  if(currentRange > tPathLength) {
428  }
429  effKinEnergy = 0.5*(e1 + preKinEnergy);
432  G4double tau = geomStepLength/lambdaeff;
433 
434  if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
435  else { tPathLength = currentRange; }
436  }
437  zPathLength = geomStepLength;
438  }
439  }
440  }
441  // check of step length
442  // define threshold angle between single and multiple scattering
443  if(!singleScatteringMode) {
444  cosThetaMin -= ssFactor*tPathLength/lambdaeff;
445  xtsec = 0.0;
446 
447  // recompute transport cross section - do not change energy
448  // anymore - cannot be applied for big steps
449  if(cosThetaMin > cosTetMaxNuc) {
450  // new computation
451  G4double cross = ComputeTransportXSectionPerVolume(cosThetaMin);
452  //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec
453  // << " 1-cosTMin= " << 1.0 - cosThetaMin << G4endl;
454  if(cross <= 0.0) {
455  singleScatteringMode = true;
457  lambdaeff = DBL_MAX;
458  cosThetaMin = 1.0;
459  } else if(xtsec > 0.0) {
460 
461  lambdaeff = 1./cross;
462  G4double tau = zPathLength*cross;
463  if(tau < numlimit) {
464  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
465  } else if(tau < 0.999999) {
466  tPathLength = -lambdaeff*G4Log(1.0 - tau);
467  } else {
469  }
470  }
471  }
472  }
474  /*
475  G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
476  <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
477  G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
478  << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
479  << " e(MeV)= " << preKinEnergy/MeV << " "
480  << " SSmode= " << singleScatteringMode << G4endl;
481  */
482  return tPathLength;
483 }
G4WentzelOKandVIxSection * wokvi
const G4MaterialCutsCouple * currentCouple
const G4Material * currentMaterial
int G4int
Definition: G4Types.hh:78
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:303
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:340
G4double G4Log(G4double x)
Definition: G4Log.hh:230
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
const G4ParticleDefinition * particle
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

void G4WentzelVIModel::DefineMaterial ( const G4MaterialCutsCouple cup)
inlineprotected

Definition at line 205 of file G4WentzelVIModel.hh.

206 {
207  if(cup != currentCouple) {
208  currentCouple = cup;
209  SetCurrentCouple(cup);
210  currentMaterial = cup->GetMaterial();
212  }
213 }
const G4MaterialCutsCouple * currentCouple
const G4Material * currentMaterial
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:447
const G4Material * GetMaterial() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4WentzelVIModel::GetFixedCut ( ) const
inline

Definition at line 235 of file G4WentzelVIModel.hh.

236 {
237  return fixedCut;
238 }
G4PhysicsTable * G4WentzelVIModel::GetSecondMomentTable ( )
inline

Definition at line 263 of file G4WentzelVIModel.hh.

264 {
265  return fSecondMoments;
266 }

Here is the caller graph for this function:

G4WentzelOKandVIxSection * G4WentzelVIModel::GetWVICrossSection ( )
inline

Definition at line 242 of file G4WentzelVIModel.hh.

243 {
244  return wokvi;
245 }
G4WentzelOKandVIxSection * wokvi
void G4WentzelVIModel::Initialise ( const G4ParticleDefinition p,
const G4DataVector cuts 
)
overridevirtual

Implements G4VEmModel.

Definition at line 123 of file G4WentzelVIModel.cc.

125 {
126  // reset parameters
127  SetupParticle(p);
128  currentRange = 0.0;
129 
130  if(isCombined) {
132  if(tet <= 0.0) { cosThetaMax = 1.0; }
133  else if(tet < CLHEP::pi) { cosThetaMax = cos(tet); }
134  }
135 
136  //G4cout << "G4WentzelVIModel::Initialise " << p->GetParticleName() << G4endl;
137  wokvi->Initialise(p, cosThetaMax);
138  /*
139  G4cout << "G4WentzelVIModel: " << particle->GetParticleName()
140  << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax
141  << " SingScatFactor= " << ssFactor
142  << G4endl;
143  */
144  currentCuts = &cuts;
145 
146  // set values of some data members
147  fParticleChange = GetParticleChangeForMSC(p);
148 
149  // build second moment table only if transport table is build
151  if(useSecondMoment && IsMaster() && table) {
152 
153  //G4cout << "### G4WentzelVIModel::Initialise: build 2nd moment table "
154  // << table << G4endl;
155  fSecondMoments =
157  // Access to materials
158  const G4ProductionCutsTable* theCoupleTable =
160  size_t numOfCouples = theCoupleTable->GetTableSize();
161 
162  G4bool splineFlag = true;
163  G4PhysicsVector* aVector = nullptr;
164  G4PhysicsVector* bVector = nullptr;
167  if(emin < emax) {
169  *G4lrint(std::log10(emax/emin));
170  if(n < 3) { n = 3; }
171 
172  for(size_t i=0; i<numOfCouples; ++i) {
173 
174  //G4cout<< "i= " << i << " Flag= " << fSecondMoments->GetFlag(i)
175  // << G4endl;
176  if(fSecondMoments->GetFlag(i)) {
177  DefineMaterial(theCoupleTable->GetMaterialCutsCouple(i));
178 
179  delete (*fSecondMoments)[i];
180  if(!aVector) {
181  aVector = new G4PhysicsLogVector(emin, emax, n);
182  bVector = aVector;
183  } else {
184  bVector = new G4PhysicsVector(*aVector);
185  }
186  for(size_t j=0; j<n; ++j) {
187  G4double e = bVector->Energy(j);
188  bVector->PutValue(j, ComputeSecondMoment(p, e)*e*e);
189  }
190  if(splineFlag) { bVector->FillSecondDerivatives(); }
191  (*fSecondMoments)[i] = bVector;
192  }
193  }
194  }
195  //G4cout << *fSecondMoments << G4endl;
196  }
197 }
G4double LowEnergyActivationLimit() const
Definition: G4VEmModel.hh:657
G4double HighEnergyActivationLimit() const
Definition: G4VEmModel.hh:650
G4int NumberOfBinsPerDecade() const
G4WentzelOKandVIxSection * wokvi
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
void DefineMaterial(const G4MaterialCutsCouple *)
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
G4PhysicsTable * GetCrossSectionTable()
Definition: G4VEmModel.hh:833
void FillSecondDerivatives()
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
static G4double tet[DIM]
bool G4bool
Definition: G4Types.hh:79
void PutValue(size_t index, G4double theValue)
const G4int n
G4double Energy(size_t index) const
static const G4double emax
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=nullptr)
Definition: G4VMscModel.cc:91
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4EmParameters * Instance()
G4double PolarAngleLimit() const
Definition: G4VEmModel.hh:664
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4bool GetFlag(size_t i) const
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54
void Initialise(const G4ParticleDefinition *, G4double CosThetaLim)

Here is the call graph for this function:

void G4WentzelVIModel::InitialiseLocal ( const G4ParticleDefinition ,
G4VEmModel masterModel 
)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 201 of file G4WentzelVIModel.cc.

203 {
204  fSecondMoments = static_cast<G4WentzelVIModel*>(masterModel)
206 }
G4PhysicsTable * GetSecondMomentTable()

Here is the call graph for this function:

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

Reimplemented from G4VMscModel.

Definition at line 488 of file G4WentzelVIModel.cc.

490 {
491  fDisplacement.set(0.0,0.0,0.0);
492  //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
493  // << particle->GetParticleName() << G4endl;
494 
495  // ignore scattering for zero step length and energy below the limit
496  if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0)
497  { return fDisplacement; }
498 
499  G4double invlambda = 0.0;
500  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
501 
502  // use average kinetic energy over the step
503  G4double cut = (*currentCuts)[currentMaterialIndex];
504  if(fixedCut > 0.0) { cut = fixedCut; }
505  /*
506  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
507  << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
508  << " xmsc= " << tPathLength*invlambda
509  << " safety= " << safety << G4endl;
510  */
511  // step limit due msc
512  G4int nMscSteps = 1;
513  G4double x0 = tPathLength;
514  G4double z0 = x0*invlambda;
515  //G4double zzz = 0.0;
516  G4double prob2 = 0.0;
517 
518  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
519 
520  // large scattering angle case - two step approach
521  if(!singleScatteringMode) {
522  static const G4double zzmin = 0.05;
523  if(useSecondMoment) {
524  G4double z1 = invlambda*invlambda;
525  G4double z2 = SecondMoment(particle, currentCouple, effKinEnergy);
526  prob2 = (z2 - z1)/(1.5*z1 - z2);
527  }
528  // if(z0 > zzmin && safety > tlimitminfix) {
529  if(z0 > zzmin) {
530  x0 *= 0.5;
531  z0 *= 0.5;
532  nMscSteps = 2;
533  }
534  //if(z0 > zzmin) { zzz = G4Exp(-1.0/z0); }
535  G4double zzz = 0.0;
536  if(z0 > zzmin) {
537  zzz = G4Exp(-1.0/z0);
538  z0 += zzz;
539  prob2 *= (1 + zzz);
540  }
541  prob2 /= (1 + prob2);
542  }
543 
544  // step limit due to single scattering
545  G4double x1 = 2*tPathLength;
546  if(0.0 < xtsec) { x1 = -G4Log(rndmEngine->flat())/xtsec; }
547 
548  // no scattering case
549  if(singleScatteringMode && x1 > tPathLength)
550  { return fDisplacement; }
551 
552  const G4ElementVector* theElementVector =
555 
556  // geometry
557  G4double sint, cost, phi;
558  G4ThreeVector temp(0.0,0.0,1.0);
559 
560  // current position and direction relative to the end point
561  // because of magnetic field geometry is computed relatively to the
562  // end point of the step
563  G4ThreeVector dir(0.0,0.0,1.0);
564  fDisplacement.set(0.0,0.0,-zPathLength);
565 
567 
568  // start a loop
569  G4double x2 = x0;
570  G4double step, z;
571  G4bool singleScat;
572  /*
573  G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
574  << " 1-cost1= " << 1 - cosThetaMin << " SSmode= " << singleScatteringMode
575  << " xtsec= " << xtsec << " Nst= " << nMscSteps << G4endl;
576  */
577  do {
578 
579  //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl;
580  // single scattering case
581  if(singleScatteringMode && x1 > x2) {
582  fDisplacement += x2*mscfac*dir;
583  break;
584  }
585 
586  // what is next single of multiple?
587  if(x1 <= x2) {
588  step = x1;
589  singleScat = true;
590  } else {
591  step = x2;
592  singleScat = false;
593  }
594 
595  //G4cout << "# step(mm)= "<< step<< " singlScat= "<< singleScat << G4endl;
596 
597  // new position
598  fDisplacement += step*mscfac*dir;
599 
600  if(singleScat) {
601 
602  // select element
603  G4int i = 0;
604  if(nelm > 1) {
605  G4double qsec = rndmEngine->flat()*xtsec;
606  for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
607  }
608  G4double cosTetM =
609  wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
610  //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " "
611  // << prob[i] << G4endl;
612  temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
613 
614  // direction is changed
615  temp.rotateUz(dir);
616  dir = temp;
617  //G4cout << dir << G4endl;
618 
619  // new proposed step length
620  x2 -= step;
621  x1 = -G4Log(rndmEngine->flat())/xtsec;
622 
623  // multiple scattering
624  } else {
625  --nMscSteps;
626  x1 -= step;
627  x2 = x0;
628 
629  // sample z in interval 0 - 1
630  G4bool isFirst = true;
631  if(prob2 > 0.0 && rndmEngine->flat() < prob2) { isFirst = false; }
632  do {
633  //z = -z0*G4Log(1.0 - (1.0 - zzz)*rndmEngine->flat());
634  if(isFirst) { z = -G4Log(rndmEngine->flat()); }
635  else { z = G4RandGamma::shoot(rndmEngine, 2.0, 2.0); }
636  z *= z0;
637  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
638  } while(z > 1.0);
639 
640  cost = 1.0 - 2.0*z/*factCM*/;
641  if(cost > 1.0) { cost = 1.0; }
642  else if(cost < -1.0) { cost =-1.0; }
643  sint = sqrt((1.0 - cost)*(1.0 + cost));
644  phi = twopi*rndmEngine->flat();
645  G4double vx1 = sint*cos(phi);
646  G4double vy1 = sint*sin(phi);
647 
648  // lateral displacement
649  if (latDisplasment) {
650  G4double rms = invsqrt12*sqrt(2*z0);
651  G4double r = x0*mscfac;
652  G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0));
653  G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0));
654  G4double d = r*r - dx*dx - dy*dy;
655 
656  // change position
657  if(d >= 0.0) {
658  temp.set(dx,dy,sqrt(d) - r);
659  temp.rotateUz(dir);
660  fDisplacement += temp;
661  }
662  }
663  // change direction
664  temp.set(vx1,vy1,cost);
665  temp.rotateUz(dir);
666  dir = temp;
667  }
668  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
669  } while (0 < nMscSteps);
670 
671  dir.rotateUz(oldDirection);
672 
673  //G4cout<<"G4WentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl;
674  // end of sampling -------------------------------
675 
676  fParticleChange->ProposeMomentumDirection(dir);
677 
678  // lateral displacement
679  fDisplacement.rotateUz(oldDirection);
680 
681  /*
682  G4cout << " r(mm)= " << fDisplacement.mag()
683  << " safety= " << safety
684  << " trueStep(mm)= " << tPathLength
685  << " geomStep(mm)= " << zPathLength
686  << " x= " << fDisplacement.x()
687  << " y= " << fDisplacement.y()
688  << " z= " << fDisplacement.z()
689  << G4endl;
690  */
691 
692  //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
693  return fDisplacement;
694 }
void set(double x, double y, double z)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4WentzelOKandVIxSection * wokvi
std::vector< G4Element * > G4ElementVector
G4bool latDisplasment
Definition: G4VMscModel.hh:188
const G4MaterialCutsCouple * currentCouple
virtual double flat()=0
const G4Material * currentMaterial
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
G4double SetupTarget(G4int Z, G4double cut=DBL_MAX)
static constexpr double twopi
Definition: G4SIunits.hh:76
bool G4bool
Definition: G4Types.hh:79
G4double SecondMoment(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double kineticEnergy)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:184
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
const G4ParticleDefinition * particle
tuple z
Definition: test.py:28
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
G4ThreeVector & SampleSingleScattering(G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio=0.0)

Here is the call graph for this function:

G4double G4WentzelVIModel::SecondMoment ( const G4ParticleDefinition part,
const G4MaterialCutsCouple couple,
G4double  kineticEnergy 
)
inline

Definition at line 271 of file G4WentzelVIModel.hh.

274 {
275  G4double x = 0.0;
276  if(useSecondMoment) {
277  DefineMaterial(couple);
278  if(fSecondMoments) {
279  x = (*fSecondMoments)[(*theDensityIdx)[currentMaterialIndex]]
280  ->Value(ekin, idx2)
281  *(*theDensityFactor)[currentMaterialIndex]/(ekin*ekin);
282  } else {
283  x = ComputeSecondMoment(part, ekin);
284  }
285  }
286  return x;
287 }
void DefineMaterial(const G4MaterialCutsCouple *)
tuple x
Definition: test.py:50
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4WentzelVIModel::SetFixedCut ( G4double  val)
inline

Definition at line 228 of file G4WentzelVIModel.hh.

229 {
230  fixedCut = val;
231 }
void G4WentzelVIModel::SetSingleScatteringFactor ( G4double  val)

Definition at line 783 of file G4WentzelVIModel.cc.

784 {
785  if(val > 0.05) {
786  ssFactor = val;
787  invssFactor = 1.0/(val - 0.05);
788  }
789 }

Here is the caller graph for this function:

void G4WentzelVIModel::SetUseSecondMoment ( G4bool  val)
inline

Definition at line 249 of file G4WentzelVIModel.hh.

250 {
251  useSecondMoment = val;
252 }
void G4WentzelVIModel::StartTracking ( G4Track track)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 243 of file G4WentzelVIModel.cc.

244 {
245  SetupParticle(track->GetDynamicParticle()->GetDefinition());
246 }
const G4DynamicParticle * GetDynamicParticle() const
G4ParticleDefinition * GetDefinition() const

Here is the call graph for this function:

G4bool G4WentzelVIModel::UseSecondMoment ( ) const
inline

Definition at line 256 of file G4WentzelVIModel.hh.

257 {
258  return useSecondMoment;
259 }

Member Data Documentation

G4double G4WentzelVIModel::cosTetMaxNuc
protected

Definition at line 154 of file G4WentzelVIModel.hh.

const G4MaterialCutsCouple* G4WentzelVIModel::currentCouple
protected

Definition at line 158 of file G4WentzelVIModel.hh.

const G4Material* G4WentzelVIModel::currentMaterial
protected

Definition at line 159 of file G4WentzelVIModel.hh.

G4int G4WentzelVIModel::currentMaterialIndex
protected

Definition at line 157 of file G4WentzelVIModel.hh.

G4double G4WentzelVIModel::currentRange
protected

Definition at line 153 of file G4WentzelVIModel.hh.

G4double G4WentzelVIModel::invssFactor
protected

Definition at line 146 of file G4WentzelVIModel.hh.

G4double G4WentzelVIModel::lambdaeff
protected

Definition at line 152 of file G4WentzelVIModel.hh.

const G4ParticleDefinition* G4WentzelVIModel::particle
protected

Definition at line 161 of file G4WentzelVIModel.hh.

G4double G4WentzelVIModel::preKinEnergy
protected

Definition at line 149 of file G4WentzelVIModel.hh.

G4bool G4WentzelVIModel::singleScatteringMode
protected

Definition at line 164 of file G4WentzelVIModel.hh.

G4double G4WentzelVIModel::ssFactor
protected

Definition at line 145 of file G4WentzelVIModel.hh.

G4double G4WentzelVIModel::tlimitminfix
protected

Definition at line 144 of file G4WentzelVIModel.hh.

G4double G4WentzelVIModel::tPathLength
protected

Definition at line 150 of file G4WentzelVIModel.hh.

G4WentzelOKandVIxSection* G4WentzelVIModel::wokvi
protected

Definition at line 142 of file G4WentzelVIModel.hh.

G4double G4WentzelVIModel::zPathLength
protected

Definition at line 151 of file G4WentzelVIModel.hh.


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