Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4WentzelVIModel Class Reference

#include <G4WentzelVIModel.hh>

+ Inheritance diagram for G4WentzelVIModel:

Public Member Functions

 G4WentzelVIModel (const G4String &nam="WentzelVIUni")
 
virtual ~G4WentzelVIModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
void StartTracking (G4Track *)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=DBL_MAX, G4double emax=DBL_MAX)
 
virtual G4ThreeVector & SampleScattering (const G4DynamicParticle *, G4double safety)
 
virtual G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep)
 
virtual G4double ComputeGeomPathLength (G4double truePathLength)
 
virtual G4double ComputeTrueStepLength (G4double geomStepLength)
 
- Public Member Functions inherited from G4VMscModel
 G4VMscModel (const G4String &nam)
 
virtual ~G4VMscModel ()
 
virtual G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &stepLimit)
 
virtual G4double ComputeGeomPathLength (G4double truePathLength)
 
virtual G4double ComputeTrueStepLength (G4double geomPathLength)
 
virtual G4ThreeVector & SampleScattering (const G4DynamicParticle *, G4double safety)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double tmax)
 
void SetStepLimitType (G4MscStepLimitType)
 
void SetLateralDisplasmentFlag (G4bool val)
 
void SetRangeFactor (G4double)
 
void SetGeomFactor (G4double)
 
void SetSkin (G4double)
 
void SetSampleZ (G4bool)
 
G4VEnergyLossProcess * GetIonisation () const
 
void SetIonisation (G4VEnergyLossProcess *, const G4ParticleDefinition *part)
 
G4double ComputeSafety (const G4ThreeVector &position, G4double limit)
 
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 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 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 ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., 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 *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
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 G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *)
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
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
 
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 ForceBuildTable (G4bool val)
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4Element * GetCurrentElement () const
 

Additional Inherited Members

- Protected Member Functions inherited from G4VMscModel
G4ParticleChangeForMSC * GetParticleChangeForMSC (const G4ParticleDefinition *p=0)
 
G4double ConvertTrueToGeom (G4double &tLength, G4double &gLength)
 
- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCouple * CurrentCouple () 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
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 

Detailed Description

Definition at line 71 of file G4WentzelVIModel.hh.

Constructor & Destructor Documentation

◆ G4WentzelVIModel()

G4WentzelVIModel::G4WentzelVIModel ( const G4String &  nam = "WentzelVIUni")

Definition at line 74 of file G4WentzelVIModel.cc.

74 :
75 G4VMscModel(nam),
76 numlimit(0.1),
77 currentCouple(0),
78 cosThetaMin(1.0),
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 G4WentzelOKandVIxSection();
93
94 preKinEnergy = tPathLength = zPathLength = lambdaeff = currentRange = xtsec = 0;
95 currentMaterialIndex = 0;
96 cosThetaMax = cosTetMaxNuc = 1.0;
97
98 fParticleChange = 0;
99 currentCuts = 0;
100 currentMaterial = 0;
101}
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:107
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:50

◆ ~G4WentzelVIModel()

G4WentzelVIModel::~G4WentzelVIModel ( )
virtual

Definition at line 105 of file G4WentzelVIModel.cc.

106{
107 delete wokvi;
108}

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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

Reimplemented from G4VEmModel.

Definition at line 134 of file G4WentzelVIModel.cc.

139{
140 G4double cross = 0.0;
141 if(p != particle) { SetupParticle(p); }
142 if(kinEnergy < lowEnergyLimit) { return cross; }
143 if(!CurrentCouple()) {
144 G4Exception("G4WentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
145 FatalException, " G4MaterialCutsCouple is not defined");
146 return 0.0;
147 }
148 DefineMaterial(CurrentCouple());
149 cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
150 if(cosTetMaxNuc < 1.0) {
151 cosTetMaxNuc = wokvi->SetupTarget(G4lrint(Z), cutEnergy);
152 cross = wokvi->ComputeTransportCrossSectionPerAtom(cosTetMaxNuc);
153 /*
154 if(p->GetParticleName() == "e-")
155 G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= " << kinEnergy
156 << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
157 << " " << particle->GetParticleName() << G4endl;
158 */
159 }
160 return cross;
161}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4double
double G4double
Definition: G4Types.hh:64
G4VEmModel::CurrentCouple
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:377
G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
Definition: G4WentzelOKandVIxSection.cc:192
G4WentzelOKandVIxSection::SetupKinematic
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
Definition: G4WentzelOKandVIxSection.hh:179
G4WentzelOKandVIxSection::SetupTarget
G4double SetupTarget(G4int Z, G4double cut=DBL_MAX)
Definition: G4WentzelOKandVIxSection.cc:151
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163

◆ ComputeGeomPathLength()

G4double G4WentzelVIModel::ComputeGeomPathLength ( G4double  truePathLength)
virtual

Reimplemented from G4VMscModel.

Definition at line 272 of file G4WentzelVIModel.cc.

273{
274 tPathLength = truelength;
275 zPathLength = tPathLength;
276
277 if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) {
278 G4double tau = tPathLength/lambdaeff;
279 //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
280 // << " Leff= " << lambdaeff << " tau= " << tau << G4endl;
281 // small step
282 if(tau < numlimit) {
283 zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
284
285 // medium step
286 } else {
287 G4double e1 = 0.0;
288 if(currentRange > tPathLength) {
289 e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
290 }
291 e1 = 0.5*(e1 + preKinEnergy);
292 cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
293 lambdaeff = GetTransportMeanFreePath(particle,e1);
294 zPathLength = lambdaeff*(1.0 - exp(-tPathLength/lambdaeff));
295 }
296 } else { lambdaeff = DBL_MAX; }
297 //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
298 return zPathLength;
299}
G4VMscModel::GetTransportMeanFreePath
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:332
G4VMscModel::GetEnergy
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:304
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83

◆ ComputeTruePathLengthLimit()

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

Reimplemented from G4VMscModel.

Definition at line 173 of file G4WentzelVIModel.cc.

176{
177 G4double tlimit = currentMinimalStep;
178 const G4DynamicParticle* dp = track.GetDynamicParticle();
179 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
180 G4StepStatus stepStatus = sp->GetStepStatus();
181 singleScatteringMode = false;
182 //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
183 // << stepStatus << G4endl;
184
185 // initialisation for each step, lambda may be computed from scratch
186 preKinEnergy = dp->GetKineticEnergy();
187 DefineMaterial(track.GetMaterialCutsCouple());
188 lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy);
189 currentRange = GetRange(particle,preKinEnergy,currentCouple);
190 cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial);
191
192 //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
193 // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
194
195 // extra check for abnormal situation
196 // this check needed to run MSC with eIoni and eBrem inactivated
197 if(tlimit > currentRange) { tlimit = currentRange; }
198
199 // stop here if small range particle
200 if(inside || tlimit < tlimitminfix) {
201 return ConvertTrueToGeom(tlimit, currentMinimalStep);
202 }
203
204 // pre step
205 G4double presafety = sp->GetSafety();
206 // far from geometry boundary
207 if(currentRange < presafety) {
208 inside = true;
209 return ConvertTrueToGeom(tlimit, currentMinimalStep);
210 }
211
212 // compute presafety again if presafety <= 0 and no boundary
213 // i.e. when it is needed for optimization purposes
214 if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
215 presafety = ComputeSafety(sp->GetPosition(), tlimit);
216 if(currentRange < presafety) {
217 inside = true;
218 return ConvertTrueToGeom(tlimit, currentMinimalStep);
219 }
220 }
221 /*
222 G4cout << "e(MeV)= " << preKinEnergy/MeV
223 << " " << particle->GetParticleName()
224 << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
225 << " R(mm)= " <<currentRange/mm
226 << " L0(mm^-1)= " << lambdaeff*mm
227 <<G4endl;
228 */
229
230 // natural limit for high energy
231 G4double rlimit = std::max(facrange*currentRange,
232 0.7*(1.0 - cosTetMaxNuc)*lambdaeff);
233
234 // low-energy e-
235 if(cosThetaMax > cosTetMaxNuc) {
236 rlimit = std::min(rlimit, facsafety*presafety);
237 }
238
239 // cut correction
240 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
241 //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= "
242 // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax
243 //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
244 if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
245
246 if(rlimit < tlimit) { tlimit = rlimit; }
247
248 tlimit = std::max(tlimit, tlimitminfix);
249
250 // step limit in infinite media
251 tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
252
253 //compute geomlimit and force few steps within a volume
254 if (steppingAlgorithm == fUseDistanceToBoundary
255 && stepStatus == fGeomBoundary) {
256
257 G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
258 tlimit = std::min(tlimit, geomlimit/facgeom);
259 }
260
261 /*
262 G4cout << particle->GetParticleName() << " e= " << preKinEnergy
263 << " L0= " << lambdaeff << " R= " << currentRange
264 << "tlimit= " << tlimit
265 << " currentMinimalStep= " << currentMinimalStep << G4endl;
266 */
267 return ConvertTrueToGeom(tlimit, currentMinimalStep);
268}
fUseDistanceToBoundary
@ fUseDistanceToBoundary
Definition: G4MscStepLimitType.hh:51
G4StepStatus
G4StepStatus
Definition: G4StepStatus.hh:51
fGeomBoundary
@ fGeomBoundary
Definition: G4StepStatus.hh:54
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4MaterialCutsCouple::GetProductionCuts
G4ProductionCuts * GetProductionCuts() const
Definition: G4MaterialCutsCouple.hh:127
G4Material::GetRadlen
G4double GetRadlen() const
Definition: G4Material.hh:219
G4ProductionCuts::GetProductionCut
G4double GetProductionCut(G4int index) const
Definition: G4ProductionCuts.hh:157
G4Step::GetPreStepPoint
G4StepPoint * GetPreStepPoint() const
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4Track::GetStep
const G4Step * GetStep() const
G4VMscModel::facrange
G4double facrange
Definition: G4VMscModel.hh:176
G4VMscModel::ComputeGeomLimit
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:256
G4VMscModel::GetRange
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:288
G4VMscModel::ComputeSafety
G4double ComputeSafety(const G4ThreeVector &position, G4double limit)
Definition: G4VMscModel.hh:238
G4VMscModel::steppingAlgorithm
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:186
G4VMscModel::ConvertTrueToGeom
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:246
G4VMscModel::facsafety
G4double facsafety
Definition: G4VMscModel.hh:178
G4VMscModel::facgeom
G4double facgeom
Definition: G4VMscModel.hh:177

◆ ComputeTrueStepLength()

G4double G4WentzelVIModel::ComputeTrueStepLength ( G4double  geomStepLength)
virtual

Reimplemented from G4VMscModel.

Definition at line 303 of file G4WentzelVIModel.cc.

304{
305 // initialisation of single scattering x-section
306 xtsec = 0.0;
307 cosThetaMin = cosTetMaxNuc;
308 /*
309 G4cout << "ComputeTrueStepLength: Step= " << geomStepLength
310 << " Lambda= " << lambdaeff
311 << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
312 */
313 // pathalogical case
314 if(lambdaeff == DBL_MAX) {
315 singleScatteringMode = true;
316 zPathLength = geomStepLength;
317 tPathLength = geomStepLength;
318 cosThetaMin = 1.0;
319
320 // normal case
321 } else {
322
323 // small step use only single scattering
324 const G4double singleScatLimit = 1.0e-7;
325 if(geomStepLength < lambdaeff*singleScatLimit*(1.0 - cosTetMaxNuc)) {
326 singleScatteringMode = true;
327 cosThetaMin = 1.0;
328 zPathLength = geomStepLength;
329 tPathLength = geomStepLength;
330
331 // step defined by transportation
332 } else if(geomStepLength != zPathLength) {
333
334 // step defined by transportation
335 zPathLength = geomStepLength;
336 G4double tau = geomStepLength/lambdaeff;
337 tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
338
339 // energy correction for a big step
340 if(tau > numlimit) {
341 G4double e1 = 0.0;
342 if(currentRange > tPathLength) {
343 e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
344 }
345 e1 = 0.5*(e1 + preKinEnergy);
346 cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
347 lambdaeff = GetTransportMeanFreePath(particle,e1);
348 tau = zPathLength/lambdaeff;
349
350 if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); }
351 else { tPathLength = currentRange; }
352 }
353 }
354 }
355
356 // check of step length
357 // define threshold angle between single and multiple scattering
358 if(!singleScatteringMode) { cosThetaMin = 1.0 - 1.5*tPathLength/lambdaeff; }
359
360 // recompute transport cross section - do not change energy
361 // anymore - cannot be applied for big steps
362 if(cosThetaMin > cosTetMaxNuc) {
363
364 // new computation
365 G4double cross = ComputeXSectionPerVolume();
366 //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec << G4endl;
367 if(cross <= 0.0) {
368 singleScatteringMode = true;
369 tPathLength = zPathLength;
370 lambdaeff = DBL_MAX;
371 cosThetaMin = 1.0;
372 } else if(xtsec > 0.0) {
373
374 lambdaeff = 1./cross;
375 G4double tau = zPathLength*cross;
376 if(tau < numlimit) { tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); }
377 else if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); }
378 else { tPathLength = currentRange; }
379
380 if(tPathLength > currentRange) { tPathLength = currentRange; }
381 }
382 }
383
384 /*
385 G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
386 <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
387 G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
388 << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
389 << " e(MeV)= " << preKinEnergy/MeV << " "
390 << singleScatteringMode << G4endl;
391 */
392 return tPathLength;
393}

◆ Initialise()

void G4WentzelVIModel::Initialise ( const G4ParticleDefinition *  p,
const G4DataVector &  cuts 
)
virtual

Implements G4VEmModel.

Definition at line 112 of file G4WentzelVIModel.cc.

114{
115 // reset parameters
116 SetupParticle(p);
117 currentRange = 0.0;
118
119 cosThetaMax = cos(PolarAngleLimit());
120 wokvi->Initialise(p, cosThetaMax);
121 /*
122 G4cout << "G4WentzelVIModel: " << particle->GetParticleName()
123 << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax
124 << G4endl;
125 */
126 currentCuts = &cuts;
127
128 // set values of some data members
129 fParticleChange = GetParticleChangeForMSC(p);
130}
G4VEmModel::PolarAngleLimit
G4double PolarAngleLimit() const
Definition: G4VEmModel.hh:550
G4VMscModel::GetParticleChangeForMSC
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:89
G4WentzelOKandVIxSection::Initialise
void Initialise(const G4ParticleDefinition *, G4double CosThetaLim)
Definition: G4WentzelOKandVIxSection.cc:118

◆ SampleScattering()

G4ThreeVector & G4WentzelVIModel::SampleScattering ( const G4DynamicParticle *  dynParticle,
G4double  safety 
)
virtual

Reimplemented from G4VMscModel.

Definition at line 398 of file G4WentzelVIModel.cc.

400{
401 fDisplacement.set(0.0,0.0,0.0);
402 //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
403 // << particle->GetParticleName() << G4endl;
404
405 // ignore scattering for zero step length and energy below the limit
406 G4double tkin = dynParticle->GetKineticEnergy();
407 if(tkin < lowEnergyLimit || tPathLength <= 0.0)
408 { return fDisplacement; }
409
410 G4double invlambda = 0.0;
411 if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
412
413 // use average kinetic energy over the step
414 G4double cut = (*currentCuts)[currentMaterialIndex];
415 /*
416 G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
417 << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
418 << " xmsc= " << tPathLength*invlambda
419 << " safety= " << safety << G4endl;
420 */
421
422 // step limit due msc
423 G4double x0 = tPathLength;
424 // const G4double thinlimit = 0.5;
425 const G4double thinlimit = 0.1;
426 G4int nMscSteps = 1;
427 // large scattering angle case - two step approach
428 if(!singleScatteringMode && tPathLength*invlambda > thinlimit
429 && safety > tlimitminfix) {
430 x0 *= 0.5;
431 nMscSteps = 2;
432 }
433
434 // step limit due to single scattering
435 G4double x1 = 2*tPathLength;
436 if(0.0 < xtsec) { x1 = -log(G4UniformRand())/xtsec; }
437
438 // no scattering case
439 if(singleScatteringMode && x1 > tPathLength)
440 { return fDisplacement; }
441
442 const G4ElementVector* theElementVector =
443 currentMaterial->GetElementVector();
444 G4int nelm = currentMaterial->GetNumberOfElements();
445
446 // geometry
447 G4double sint, cost, phi;
448 G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
449 G4ThreeVector temp(0.0,0.0,1.0);
450
451 // current position and direction relative to the end point
452 // because of magnetic field geometry is computed relatively to the
453 // end point of the step
454 G4ThreeVector dir(0.0,0.0,1.0);
455 fDisplacement.set(0.0,0.0,-zPathLength);
456
457 G4double mscfac = zPathLength/tPathLength;
458
459 // start a loop
460 G4double x2 = x0;
461 G4double step, z;
462 G4bool singleScat;
463 /*
464 G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
465 << " 1-cost1= " << 1 - cosThetaMin << " " << singleScatteringMode
466 << " xtsec= " << xtsec << G4endl;
467 */
468 do {
469
470 // single scattering case
471 if(singleScatteringMode || x1 <= x2) {
472 step = x1;
473 singleScat = true;
474 } else {
475 step = x2;
476 singleScat = false;
477 }
478
479 // new position
480 fDisplacement += step*mscfac*dir;
481
482 if(singleScat) {
483
484 // select element
485 G4int i = 0;
486 if(nelm > 1) {
487 G4double qsec = G4UniformRand()*xtsec;
488 for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
489 }
490 G4double cosTetM =
491 wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
492 //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " " << prob[i] << G4endl;
493 temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
494
495 // direction is changed
496 temp.rotateUz(dir);
497 dir = temp;
498
499 // new proposed step length
500 x2 -= step;
501 x1 = -log(G4UniformRand())/xtsec;
502 if(singleScatteringMode && x1 > x2) { break; }
503
504 // multiple scattering
505 } else {
506 --nMscSteps;
507 x1 -= step;
508 x2 = x0;
509
510 // for multiple scattering x0 should be used as a step size
511 G4double z0 = x0*invlambda;
512
513 // sample z in interval 0 - 1
514 do {
515 G4double zzz = 0.0;
516 if(z0 > 0.1) { zzz = exp(-1.0/z0); }
517 z = -z0*log(1.0 - (1.0 - zzz)*G4UniformRand());
518 } while(z > 1.0);
519 cost = 1.0 - 2.0*z/*factCM*/;
520 if(cost > 1.0) { cost = 1.0; }
521 else if(cost < -1.0) { cost =-1.0; }
522 sint = sqrt((1.0 - cost)*(1.0 + cost));
523 phi = twopi*G4UniformRand();
524 G4double vx1 = sint*cos(phi);
525 G4double vy1 = sint*sin(phi);
526
527 // change direction
528 temp.set(vx1,vy1,cost);
529 temp.rotateUz(dir);
530 dir = temp;
531
532 // lateral displacement
533 if (latDisplasment && safety > tlimitminfix) {
534 G4double rms = invsqrt12*sqrt(2*z0);
535 G4double r = x0*mscfac;
536 G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0));
537 G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0));
538 G4double dz;
539 G4double d = r*r - dx*dx - dy*dy;
540 if(d >= 0.0) { dz = sqrt(d) - r; }
541 else { dx = dy = dz = 0.0; }
542
543 // change position
544 temp.set(dx,dy,dz);
545 temp.rotateUz(dir);
546 fDisplacement += temp;
547 }
548 }
549 } while (0 < nMscSteps);
550
551 dir.rotateUz(oldDirection);
552
553 //G4cout << "G4WentzelVIModel sampling of scattering is done" << G4endl;
554 // end of sampling -------------------------------
555
556 fParticleChange->ProposeMomentumDirection(dir);
557
558 // lateral displacement
559 fDisplacement.rotateUz(oldDirection);
560
561 /*
562 G4cout << " r(mm)= " << fDisplacement.mag()
563 << " safety= " << safety
564 << " trueStep(mm)= " << tPathLength
565 << " geomStep(mm)= " << zPathLength
566 << " x= " << fDisplacement.x()
567 << " y= " << fDisplacement.y()
568 << " z= " << fDisplacement.z()
569 << G4endl;
570 */
571
572 //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
573 return fDisplacement;
574}
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:189
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
G4ParticleChangeForMSC::ProposeMomentumDirection
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4VMscModel::latDisplasment
G4bool latDisplasment
Definition: G4VMscModel.hh:189
G4VMscModel::fDisplacement
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:185
G4WentzelOKandVIxSection::SampleSingleScattering
G4ThreeVector SampleSingleScattering(G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio=0.0)
Definition: G4WentzelOKandVIxSection.cc:289

◆ StartTracking()

void G4WentzelVIModel::StartTracking ( G4Track *  track)
virtual

Reimplemented from G4VEmModel.

Definition at line 165 of file G4WentzelVIModel.cc.

166{
167 SetupParticle(track->GetDynamicParticle()->GetDefinition());
168 inside = false;
169}
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
The documentation for this class was generated from the following files: