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

#include <G4WentzelVIRelModel.hh>

Inheritance diagram for G4WentzelVIRelModel:
Collaboration diagram for G4WentzelVIRelModel:

Public Member Functions

 G4WentzelVIRelModel (G4bool combined=true)
 
virtual ~G4WentzelVIRelModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &) 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
 
- Public Member Functions inherited from G4VMscModel
 G4VMscModel (const G4String &nam)
 
virtual ~G4VMscModel ()
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double tmax) override
 
void SetStepLimitType (G4MscStepLimitType)
 
void SetLateralDisplasmentFlag (G4bool val)
 
void SetRangeFactor (G4double)
 
void SetGeomFactor (G4double)
 
void SetSkin (G4double)
 
void SetSampleZ (G4bool)
 
G4VEnergyLossProcessGetIonisation () const
 
void SetIonisation (G4VEnergyLossProcess *, const G4ParticleDefinition *part)
 
G4double ComputeSafety (const G4ThreeVector &position, G4double limit=DBL_MAX)
 
G4double ComputeGeomLimit (const G4Track &, G4double &presafety, G4double limit)
 
G4double GetDEDX (const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
 
G4double GetRange (const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
 
G4double GetEnergy (const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
 
G4double GetTransportMeanFreePath (const G4ParticleDefinition *part, G4double kinEnergy)
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double 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)
 

Additional Inherited Members

- Protected Member Functions inherited from G4VMscModel
G4ParticleChangeForMSCGetParticleChangeForMSC (const G4ParticleDefinition *p=nullptr)
 
G4double ConvertTrueToGeom (G4double &tLength, G4double &gLength)
 
- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Protected Attributes inherited from G4VMscModel
G4double facrange
 
G4double facgeom
 
G4double facsafety
 
G4double skin
 
G4double dtrl
 
G4double lambdalimit
 
G4double geomMin
 
G4double geomMax
 
G4ThreeVector fDisplacement
 
G4MscStepLimitType steppingAlgorithm
 
G4bool samplez
 
G4bool latDisplasment
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChangepParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 69 of file G4WentzelVIRelModel.hh.

Constructor & Destructor Documentation

G4WentzelVIRelModel::G4WentzelVIRelModel ( G4bool  combined = true)
explicit

Definition at line 73 of file G4WentzelVIRelModel.cc.

73  :
74  G4VMscModel("WentzelVIRel"),
75  numlimit(0.1),
76  currentCouple(nullptr),
77  cosThetaMin(1.0),
78  isCombined(combined),
79  inside(false),
80  singleScatteringMode(false)
81 {
82  invsqrt12 = 1./sqrt(12.);
83  tlimitminfix = 1.e-6*mm;
84  lowEnergyLimit = 1.0*eV;
85  particle = 0;
86  nelments = 5;
87  xsecn.resize(nelments);
88  prob.resize(nelments);
89  theManager = G4LossTableManager::Instance();
90  fNistManager = G4NistManager::Instance();
91  fG4pow = G4Pow::GetInstance();
92  wokvi = new G4WentzelVIRelXSection(combined);
93 
94  preKinEnergy = tPathLength = zPathLength = lambdaeff = currentRange =
95  xtsec = 0;
96  currentMaterialIndex = 0;
97  cosThetaMax = cosTetMaxNuc = 1.0;
98 
99  fParticleChange = nullptr;
100  currentCuts = nullptr;
101  currentMaterial = nullptr;
102 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static G4LossTableManager * Instance()
static constexpr double mm
Definition: G4SIunits.hh:115
static G4NistManager * Instance()
static constexpr double eV
Definition: G4SIunits.hh:215
G4VMscModel(const G4String &nam)
Definition: G4VMscModel.cc:60

Here is the call graph for this function:

G4WentzelVIRelModel::~G4WentzelVIRelModel ( )
virtual

Definition at line 106 of file G4WentzelVIRelModel.cc.

107 {
108  delete wokvi;
109 }

Member Function Documentation

G4double G4WentzelVIRelModel::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 140 of file G4WentzelVIRelModel.cc.

145 {
146  G4double cross = 0.0;
147  if(p != particle) { SetupParticle(p); }
148  if(kinEnergy < lowEnergyLimit) { return cross; }
149  if(!CurrentCouple()) {
150  G4Exception("G4WentzelVIRelModel::ComputeCrossSectionPerAtom", "em0011",
151  FatalException, " G4MaterialCutsCouple is not defined");
152  return cross;
153  }
154  DefineMaterial(CurrentCouple());
155  G4int iz = G4lrint(Z);
156  G4double tmass = proton_mass_c2;
157  if(1 < iz) {
158  tmass = fNistManager->GetAtomicMassAmu(iz)*amu_c2;
159  }
160  cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial,
161  cutEnergy, tmass);
162  if(cosTetMaxNuc < 1.0) {
163  G4double cost = wokvi->SetupTarget(iz, cutEnergy);
164  cross = wokvi->ComputeTransportCrossSectionPerAtom(cost);
165  /*
166  if(p->GetParticleName() == "e-")
167  G4cout << "G4WentzelVIRelModel::CS: Z= " << G4int(Z)
168  << " e(MeV)= " << kinEnergy
169  << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
170  << " " << particle->GetParticleName() << G4endl;
171  */
172  }
173  return cross;
174 }
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
int G4int
Definition: G4Types.hh:78
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:454
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat, G4double cut, G4double tmass)
G4double SetupTarget(G4int Z, G4double cut)
float proton_mass_c2
Definition: hepunit.py:275
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
G4double GetAtomicMassAmu(const G4String &symb) const
double G4double
Definition: G4Types.hh:76
float amu_c2
Definition: hepunit.py:277

Here is the call graph for this function:

G4double G4WentzelVIRelModel::ComputeGeomPathLength ( G4double  truePathLength)
overridevirtual

Reimplemented from G4VMscModel.

Definition at line 287 of file G4WentzelVIRelModel.cc.

288 {
289  tPathLength = truelength;
290  zPathLength = tPathLength;
291 
292  if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) {
293  G4double tau = tPathLength/lambdaeff;
294  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
295  // << " Leff= " << lambdaeff << " tau= " << tau << G4endl;
296  // small step
297  if(tau < numlimit) {
298  zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
299 
300  // medium step
301  } else {
302  G4double e1 = 0.0;
303  if(currentRange > tPathLength) {
304  e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
305  }
306  e1 = 0.5*(e1 + preKinEnergy);
307  cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial, 0.0,
309  lambdaeff = GetTransportMeanFreePath(particle,e1);
310  zPathLength = lambdaeff*(1.0 - G4Exp(-tPathLength/lambdaeff));
311  }
312  } else { lambdaeff = DBL_MAX; }
313  //G4cout<<"Comp.geom: zLength= "<<zPathLength
314  // <<" tLength= "<<tPathLength<<G4endl;
315  return zPathLength;
316 }
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:303
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat, G4double cut, G4double tmass)
float proton_mass_c2
Definition: hepunit.py:275
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:340
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

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

Reimplemented from G4VMscModel.

Definition at line 187 of file G4WentzelVIRelModel.cc.

190 {
191  G4double tlimit = currentMinimalStep;
192  const G4DynamicParticle* dp = track.GetDynamicParticle();
193  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
194  G4StepStatus stepStatus = sp->GetStepStatus();
195  singleScatteringMode = false;
196  //G4cout << "G4WentzelVIRelModel::ComputeTruePathLengthLimit stepStatus= "
197  // << stepStatus << G4endl;
198 
199 
200  // initialisation for each step, lambda may be computed from scratch
201  preKinEnergy = dp->GetKineticEnergy();
202  DefineMaterial(track.GetMaterialCutsCouple());
203  lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy);
204  currentRange = GetRange(particle,preKinEnergy,currentCouple);
205 
206  G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
207  cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial,
208  rcut, proton_mass_c2);
209 
210  // extra check for abnormal situation
211  // this check needed to run MSC with eIoni and eBrem inactivated
212  if(tlimit > currentRange) { tlimit = currentRange; }
213 
214  // stop here if small range particle
215  if(inside || tlimit < tlimitminfix) {
216  return ConvertTrueToGeom(tlimit, currentMinimalStep);
217  }
218 
219  // pre step
220  G4double presafety = sp->GetSafety();
221  // far from geometry boundary
222  if(currentRange < presafety) {
223  inside = true;
224  return ConvertTrueToGeom(tlimit, currentMinimalStep);
225  }
226 
227  // compute presafety again if presafety <= 0 and no boundary
228  // i.e. when it is needed for optimization purposes
229  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
230  presafety = ComputeSafety(sp->GetPosition(), tlimit);
231  if(currentRange < presafety) {
232  inside = true;
233  return ConvertTrueToGeom(tlimit, currentMinimalStep);
234  }
235  }
236  /*
237  G4cout << "e(MeV)= " << preKinEnergy/MeV
238  << " " << particle->GetParticleName()
239  << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
240  << " R(mm)= " <<currentRange/mm
241  << " L0(mm^-1)= " << lambdaeff*mm
242  <<G4endl;
243  */
244 
245  // natural limit for high energy
246  G4double rlimit = std::max(facrange*currentRange,
247  0.7*(1.0 - cosTetMaxNuc)*lambdaeff);
248 
249  // low-energy e-
250  if(cosThetaMax > cosTetMaxNuc) {
251  rlimit = std::min(rlimit, facsafety*presafety);
252  }
253 
254  // cut correction
255  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit
256  // << " presafety= " << presafety
257  // << " 1-cosThetaMax= " <<1-cosThetaMax
258  // << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
259  // << G4endl;
260  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
261 
262  if(rlimit < tlimit) { tlimit = rlimit; }
263 
264  tlimit = std::max(tlimit, tlimitminfix);
265 
266  // step limit in infinite media
267  tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
268 
269  //compute geomlimit and force few steps within a volume
271  {
272  G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
273  tlimit = std::min(tlimit, geomlimit/facgeom);
274  }
275 
276  /*
277  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
278  << " L0= " << lambdaeff << " R= " << currentRange
279  << "tlimit= " << tlimit
280  << " currentMinimalStep= " << currentMinimalStep << G4endl;
281  */
282  return ConvertTrueToGeom(tlimit, currentMinimalStep);
283 }
G4double facgeom
Definition: G4VMscModel.hh:176
G4double GetKineticEnergy() const
const G4DynamicParticle * GetDynamicParticle() const
G4double facrange
Definition: G4VMscModel.hh:175
G4double GetProductionCut(G4int index) const
G4StepStatus GetStepStatus() const
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:245
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
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 cut, G4double tmass)
float proton_mass_c2
Definition: hepunit.py:275
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
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 G4WentzelVIRelModel::ComputeTrueStepLength ( G4double  geomStepLength)
overridevirtual

Reimplemented from G4VMscModel.

Definition at line 320 of file G4WentzelVIRelModel.cc.

321 {
322  // initialisation of single scattering x-section
323  xtsec = 0.0;
324  cosThetaMin = cosTetMaxNuc;
325 
326  //G4cout << "Step= " << geomStepLength << " Lambda= " << lambdaeff
327  // << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
328  // pathalogical case
329  if(lambdaeff == DBL_MAX) {
330  singleScatteringMode = true;
331  zPathLength = geomStepLength;
332  tPathLength = geomStepLength;
333  cosThetaMin = 1.0;
334 
335  // normal case
336  } else {
337 
338  // small step use only single scattering
339  static const G4double singleScatLimit = 1.0e-7;
340  if(geomStepLength < lambdaeff*singleScatLimit*(1.0 - cosTetMaxNuc)) {
341  singleScatteringMode = true;
342  cosThetaMin = 1.0;
343  zPathLength = geomStepLength;
344  tPathLength = geomStepLength;
345 
346  // step defined by transportation
347  } else if(geomStepLength != zPathLength) {
348 
349  // step defined by transportation
350  zPathLength = geomStepLength;
351  G4double tau = geomStepLength/lambdaeff;
352  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
353 
354  // energy correction for a big step
355  if(tau > numlimit) {
356  G4double e1 = 0.0;
357  if(currentRange > tPathLength) {
358  e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
359  }
360  e1 = 0.5*(e1 + preKinEnergy);
361  cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial, 0.0,
363  lambdaeff = GetTransportMeanFreePath(particle,e1);
364  tau = zPathLength/lambdaeff;
365 
366  if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
367  else { tPathLength = currentRange; }
368  }
369  }
370  }
371 
372  // check of step length
373  // define threshold angle between single and multiple scattering
374  if(!singleScatteringMode) { cosThetaMin = 1.0 - 1.5*tPathLength/lambdaeff; }
375 
376  // recompute transport cross section - do not change energy
377  // anymore - cannot be applied for big steps
378  if(cosThetaMin > cosTetMaxNuc) {
379 
380  // new computation
381  G4double cross = ComputeXSectionPerVolume();
382  //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec << G4endl;
383  if(cross <= 0.0) {
384  singleScatteringMode = true;
385  tPathLength = zPathLength;
386  lambdaeff = DBL_MAX;
387  } else if(xtsec > 0.0) {
388 
389  lambdaeff = 1./cross;
390  G4double tau = zPathLength*cross;
391  if(tau < numlimit) {
392  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
393  }
394  else if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
395  else { tPathLength = currentRange; }
396 
397  if(tPathLength > currentRange) { tPathLength = currentRange; }
398  }
399  }
400 
401  /*
402  G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
403  <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
404  G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
405  << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
406  << " e(MeV)= " << preKinEnergy/MeV << " " << singleScatteringMode
407  << G4endl;
408  */
409  return tPathLength;
410 }
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:303
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat, G4double cut, G4double tmass)
float proton_mass_c2
Definition: hepunit.py:275
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:340
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

void G4WentzelVIRelModel::Initialise ( const G4ParticleDefinition p,
const G4DataVector cuts 
)
overridevirtual

Implements G4VEmModel.

Definition at line 113 of file G4WentzelVIRelModel.cc.

115 {
116  // reset parameters
117  SetupParticle(p);
118  currentRange = 0.0;
119 
120  if(isCombined) {
122  if(tet >= pi) { cosThetaMax = -1.0; }
123  else if(tet > 0.0) { cosThetaMax = cos(tet); }
124  }
125 
126  wokvi->Initialise(p, cosThetaMax);
127  /*
128  G4cout << "G4WentzelVIRelModel: " << particle->GetParticleName()
129  << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax
130  << G4endl;
131  */
132  currentCuts = &cuts;
133 
134  // set values of some data members
135  fParticleChange = GetParticleChangeForMSC(p);
136 }
static G4double tet[DIM]
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=nullptr)
Definition: G4VMscModel.cc:91
G4double PolarAngleLimit() const
Definition: G4VEmModel.hh:664
static constexpr double pi
Definition: G4SIunits.hh:75
void Initialise(const G4ParticleDefinition *, G4double CosThetaLim)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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

Reimplemented from G4VMscModel.

Definition at line 415 of file G4WentzelVIRelModel.cc.

417 {
418  fDisplacement.set(0.0,0.0,0.0);
419  //G4cout << "!##! G4WentzelVIRelModel::SampleScattering for "
420  // << particle->GetParticleName() << G4endl;
421 
422  // ignore scattering for zero step length and energy below the limit
423  if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0)
424  { return fDisplacement; }
425 
426  G4double invlambda = 0.0;
427  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
428 
429  // use average kinetic energy over the step
430  G4double cut = (*currentCuts)[currentMaterialIndex];
431  /*
432  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
433  << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
434  << " xmsc= " << tPathLength*invlambda
435  << " safety= " << safety << G4endl;
436  */
437 
438  // step limit due msc
439  G4double x0 = tPathLength;
440  // const G4double thinlimit = 0.5;
441  static const G4double thinlimit = 0.1;
442  G4int nMscSteps = 1;
443  // large scattering angle case - two step approach
444  if(tPathLength*invlambda > thinlimit && safety > tlimitminfix) {
445  x0 *= 0.5;
446  nMscSteps = 2;
447  }
448 
449  // step limit due to single scattering
450  G4double x1 = 2*tPathLength;
451  if(0.0 < xtsec) { x1 = -G4Log(G4UniformRand())/xtsec; }
452 
453  const G4ElementVector* theElementVector =
454  currentMaterial->GetElementVector();
455  G4int nelm = currentMaterial->GetNumberOfElements();
456 
457  // geometry
458  G4double sint, cost, phi;
459  G4ThreeVector temp(0.0,0.0,1.0);
460 
461  // current position and direction relative to the end point
462  // because of magnetic field geometry is computed relatively to the
463  // end point of the step
464  G4ThreeVector dir(0.0,0.0,1.0);
465  fDisplacement.set(0.0,0.0,-zPathLength);
466  G4double mscfac = zPathLength/tPathLength;
467 
468  // start a loop
469  G4double x2 = x0;
470  G4double step, z;
471  G4bool singleScat;
472  /*
473  G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
474  << " 1-cost1= " << 1 - cosThetaMin << " " << singleScatteringMode
475  << " xtsec= " << xtsec << G4endl;
476  */
477  do {
478 
479  // single scattering case
480  if(x1 < x2) {
481  step = x1;
482  singleScat = true;
483  } else {
484  step = x2;
485  singleScat = false;
486  }
487 
488  // new position
489  fDisplacement += step*mscfac*dir;
490 
491  if(singleScat) {
492 
493  // select element
494  G4int i = 0;
495  if(nelm > 1) {
496  G4double qsec = G4UniformRand()*xtsec;
497  for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
498  }
499  G4double cosTetM =
500  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
501  //G4cout << "!!! " << cosThetaMin << " " << cosTetM
502  // << " " << prob[i] << G4endl;
503  temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
504 
505  // direction is changed
506  temp.rotateUz(dir);
507  dir = temp;
508 
509  // new proposed step length
510  x1 = -G4Log(G4UniformRand())/xtsec;
511  x2 -= step;
512  if(x2 <= 0.0) { --nMscSteps; }
513 
514  // multiple scattering
515  } else {
516  --nMscSteps;
517  x1 -= step;
518  x2 = x0;
519 
520  // for multiple scattering x0 should be used as a step size
521  if(!singleScatteringMode) {
522  G4double z0 = x0*invlambda;
523 
524  // correction to keep first moment
525 
526  // sample z in interval 0 - 1
527  if(z0 > 5.0) { z = G4UniformRand(); }
528  else {
529  G4double zzz = 0.0;
530  if(z0 > 0.01) { zzz = G4Exp(-1.0/z0); }
531  z = -z0*G4Log(1.0 - (1.0 - zzz)*G4UniformRand());
532  // /(1.0 - (1.0/z0 + 1.0)*zzz);
533  }
534 
535  cost = 1.0 - 2.0*z/*factCM*/;
536  if(cost > 1.0) { cost = 1.0; }
537  else if(cost < -1.0) { cost =-1.0; }
538  sint = sqrt((1.0 - cost)*(1.0 + cost));
539  phi = twopi*G4UniformRand();
540  G4double vx1 = sint*cos(phi);
541  G4double vy1 = sint*sin(phi);
542 
543  // lateral displacement
544  if (latDisplasment && x0 > tlimitminfix) {
545  G4double rms = invsqrt12*sqrt(2*z0);
546  G4double r = x0*mscfac;
547  G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0));
548  G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0));
549  G4double dz;
550  G4double d = r*r - dx*dx - dy*dy;
551  if(d >= 0.0) { dz = sqrt(d) - r; }
552  else { dx = dy = dz = 0.0; }
553 
554  // change position
555  temp.set(dx,dy,dz);
556  temp.rotateUz(dir);
557  fDisplacement += temp;
558  }
559  // change direction
560  temp.set(vx1,vy1,cost);
561  temp.rotateUz(dir);
562  dir = temp;
563  }
564  }
565  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
566  } while (0 < nMscSteps);
567 
568  dir.rotateUz(oldDirection);
569 
570  //G4cout << "G4WentzelVIRelModel sampling of scattering is done" << G4endl;
571  // end of sampling -------------------------------
572 
573  fParticleChange->ProposeMomentumDirection(dir);
574 
575  // lateral displacement
576  fDisplacement.rotateUz(oldDirection);
577 
578  /*
579  G4cout << " r(mm)= " << fDisplacement.mag()
580  << " safety= " << safety
581  << " trueStep(mm)= " << tPathLength
582  << " geomStep(mm)= " << zPathLength
583  << " x= " << fDisplacement.x()
584  << " y= " << fDisplacement.y()
585  << " z= " << fDisplacement.z()
586  << G4endl;
587  */
588 
589  //G4cout<< "G4WentzelVIRelModel::SampleScattering end NewDir= "
590  // << dir<< G4endl;
591  return fDisplacement;
592 }
void set(double x, double y, double z)
ThreeVector shoot(const G4int Ap, const G4int Af)
std::vector< G4Element * > G4ElementVector
G4bool latDisplasment
Definition: G4VMscModel.hh:188
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
static constexpr double twopi
Definition: G4SIunits.hh:76
#define G4UniformRand()
Definition: Randomize.hh:97
G4double SetupTarget(G4int Z, G4double cut)
bool G4bool
Definition: G4Types.hh:79
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)
int G4lrint(double ad)
Definition: templates.hh:163
tuple z
Definition: test.py:28
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76
G4ThreeVector & SampleSingleScattering(G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio)
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

void G4WentzelVIRelModel::StartTracking ( G4Track track)
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 178 of file G4WentzelVIRelModel.cc.

179 {
180  SetupParticle(track->GetDynamicParticle()->GetDefinition());
181  inside = false;
183 }
virtual void StartTracking(G4Track *)
Definition: G4VEmModel.cc:292
const G4DynamicParticle * GetDynamicParticle() const
G4ParticleDefinition * GetDefinition() const

Here is the call graph for this function:


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