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
G4WentzelVIRelModel Class Reference

#include <G4WentzelVIRelModel.hh>

+ Inheritance diagram for G4WentzelVIRelModel:

Public Member Functions

 G4WentzelVIRelModel (const G4String &nam="WentzelVIUni")
 
virtual ~G4WentzelVIRelModel ()
 
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 69 of file G4WentzelVIRelModel.hh.

Constructor & Destructor Documentation

◆ G4WentzelVIRelModel()

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

Definition at line 71 of file G4WentzelVIRelModel.cc.

71 :
72 G4VMscModel(nam),
73 numlimit(0.1),
74 currentCouple(0),
75 cosThetaMin(1.0),
76 inside(false),
77 singleScatteringMode(false)
78{
79 invsqrt12 = 1./sqrt(12.);
80 tlimitminfix = 1.e-6*mm;
81 lowEnergyLimit = 1.0*eV;
82 particle = 0;
83 nelments = 5;
84 xsecn.resize(nelments);
85 prob.resize(nelments);
86 theManager = G4LossTableManager::Instance();
87 fNistManager = G4NistManager::Instance();
88 fG4pow = G4Pow::GetInstance();
89 wokvi = new G4WentzelVIRelXSection();
90
91 preKinEnergy = tPathLength = zPathLength = lambdaeff = currentRange = xtsec = 0;
92 currentMaterialIndex = 0;
93 cosThetaMax = cosTetMaxNuc = 1.0;
94
95 fParticleChange = 0;
96 currentCuts = 0;
97 currentMaterial = 0;
98}
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

◆ ~G4WentzelVIRelModel()

G4WentzelVIRelModel::~G4WentzelVIRelModel ( )
virtual

Definition at line 102 of file G4WentzelVIRelModel.cc.

103{
104 delete wokvi;
105}

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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

136{
137 G4double cross = 0.0;
138 if(p != particle) { SetupParticle(p); }
139 if(kinEnergy < lowEnergyLimit) { return cross; }
140 if(!CurrentCouple()) {
141 G4Exception("G4WentzelVIRelModel::ComputeCrossSectionPerAtom", "em0011",
142 FatalException, " G4MaterialCutsCouple is not defined");
143 return 0.0;
144 }
145 DefineMaterial(CurrentCouple());
146 cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
147 if(cosTetMaxNuc < 1.0) {
148 cosTetMaxNuc = wokvi->SetupTarget(G4lrint(Z), cutEnergy);
149 cross = wokvi->ComputeTransportCrossSectionPerAtom(cosTetMaxNuc);
150 /*
151 if(p->GetParticleName() == "e-")
152 G4cout << "G4WentzelVIRelModel::CS: Z= " << G4int(Z) << " e(MeV)= " << kinEnergy
153 << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
154 << " " << particle->GetParticleName() << G4endl;
155 */
156 }
157 return cross;
158}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4double
double G4double
Definition: G4Types.hh:64
G4VEmModel::CurrentCouple
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:377
G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
Definition: G4WentzelVIRelXSection.cc:192
G4WentzelVIRelXSection::SetupTarget
G4double SetupTarget(G4int Z, G4double cut=DBL_MAX)
Definition: G4WentzelVIRelXSection.cc:147
G4WentzelVIRelXSection::SetupKinematic
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
Definition: G4WentzelVIRelXSection.hh:179
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 G4WentzelVIRelModel::ComputeGeomPathLength ( G4double  truePathLength)
virtual

Reimplemented from G4VMscModel.

Definition at line 266 of file G4WentzelVIRelModel.cc.

267{
268 tPathLength = truelength;
269 zPathLength = tPathLength;
270
271 if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) {
272 G4double tau = tPathLength/lambdaeff;
273 //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
274 // << " Leff= " << lambdaeff << " tau= " << tau << G4endl;
275 // small step
276 if(tau < numlimit) {
277 zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
278
279 // medium step
280 } else {
281 G4double e1 = 0.0;
282 if(currentRange > tPathLength) {
283 e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
284 }
285 e1 = 0.5*(e1 + preKinEnergy);
286 cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
287 lambdaeff = GetTransportMeanFreePath(particle,e1);
288 zPathLength = lambdaeff*(1.0 - exp(-tPathLength/lambdaeff));
289 }
290 } else { lambdaeff = DBL_MAX; }
291 //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
292 return zPathLength;
293}
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 G4WentzelVIRelModel::ComputeTruePathLengthLimit ( const G4Track &  track,
G4double &  currentMinimalStep 
)
virtual

Reimplemented from G4VMscModel.

Definition at line 170 of file G4WentzelVIRelModel.cc.

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

Reimplemented from G4VMscModel.

Definition at line 297 of file G4WentzelVIRelModel.cc.

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

◆ Initialise()

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

Implements G4VEmModel.

Definition at line 109 of file G4WentzelVIRelModel.cc.

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

◆ SampleScattering()

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

Reimplemented from G4VMscModel.

Definition at line 388 of file G4WentzelVIRelModel.cc.

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

◆ StartTracking()

void G4WentzelVIRelModel::StartTracking ( G4Track *  track)
virtual

Reimplemented from G4VEmModel.

Definition at line 162 of file G4WentzelVIRelModel.cc.

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