Geant4.10
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4WentzelVIModel Class Reference

#include <G4WentzelVIModel.hh>

Inheritance diagram for G4WentzelVIModel:
G4VMscModel G4VEmModel

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 G4ThreeVector &, G4double safety)
 
virtual G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep)
 
virtual G4double ComputeGeomPathLength (G4double truePathLength)
 
virtual G4double ComputeTrueStepLength (G4double geomStepLength)
 
void SetFixedCut (G4double)
 
G4double GetFixedCut () const
 
G4WentzelOKandVIxSection * GetWVICrossSection ()
 
- 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)
 
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 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 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 *)
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector
< G4EmElementSelector * > * 		GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
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)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
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
 
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 SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
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
G4ElementData * fElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 

Detailed Description

Definition at line 69 of file G4WentzelVIModel.hh.

Constructor & Destructor Documentation

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

Definition at line 74 of file G4WentzelVIModel.cc.

References python.hepunit::eV, G4LossTableManager::Instance(), and python.hepunit::mm.

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  wokvi = new G4WentzelOKandVIxSection();
91  fixedCut = -1.0;
92 
93  preKinEnergy = tPathLength = zPathLength = lambdaeff = currentRange
94  = 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:112
G4VMscModel::G4VMscModel
G4VMscModel(const G4String &nam)
Definition: G4VMscModel.cc:58
G4WentzelVIModel::~G4WentzelVIModel ( )
virtual

Definition at line 105 of file G4WentzelVIModel.cc.

106 {
107  delete wokvi;
108 }

Member Function Documentation

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 135 of file G4WentzelVIModel.cc.

References G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom(), G4VEmModel::CurrentCouple(), FatalException, G4Exception(), G4lrint(), G4WentzelOKandVIxSection::SetupKinematic(), and G4WentzelOKandVIxSection::SetupTarget().

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

Reimplemented from G4VMscModel.

Definition at line 276 of file G4WentzelVIModel.cc.

References DBL_MAX, G4Exp(), G4VMscModel::GetEnergy(), G4VMscModel::GetTransportMeanFreePath(), and G4WentzelOKandVIxSection::SetupKinematic().

277 {
278  tPathLength = truelength;
279  zPathLength = tPathLength;
280 
281  if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) {
282  G4double tau = tPathLength/lambdaeff;
283  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
284  // << " Leff= " << lambdaeff << " tau= " << tau << G4endl;
285  // small step
286  if(tau < numlimit) {
287  zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
288 
289  // medium step
290  } else {
291  G4double e1 = 0.0;
292  if(currentRange > tPathLength) {
293  e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
294  }
295  e1 = 0.5*(e1 + preKinEnergy);
296  cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
297  lambdaeff = GetTransportMeanFreePath(particle,e1);
298  zPathLength = lambdaeff*(1.0 - G4Exp(-tPathLength/lambdaeff));
299  }
300  } else { lambdaeff = DBL_MAX; }
301  //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
302  return zPathLength;
303 }
G4VMscModel::GetEnergy
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:308
G4WentzelOKandVIxSection::SetupKinematic
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
Definition: G4WentzelOKandVIxSection.hh:177
G4VMscModel::GetTransportMeanFreePath
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:345
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4double
double G4double
Definition: G4Types.hh:76
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4double G4WentzelVIModel::ComputeTruePathLengthLimit ( const G4Track &  track,
G4double &  currentMinimalStep 
)
virtual

Reimplemented from G4VMscModel.

Definition at line 176 of file G4WentzelVIModel.cc.

References G4VMscModel::ComputeGeomLimit(), G4VMscModel::ComputeSafety(), G4VMscModel::ConvertTrueToGeom(), G4VMscModel::facgeom, G4VMscModel::facrange, G4VMscModel::facsafety, fGeomBoundary, fUseDistanceToBoundary, G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterialCutsCouple(), G4StepPoint::GetPosition(), G4Step::GetPreStepPoint(), G4ProductionCuts::GetProductionCut(), G4MaterialCutsCouple::GetProductionCuts(), G4Material::GetRadlen(), G4VMscModel::GetRange(), G4StepPoint::GetSafety(), G4Track::GetStep(), G4StepPoint::GetStepStatus(), G4VMscModel::GetTransportMeanFreePath(), G4INCL::Math::max(), G4INCL::Math::min(), G4WentzelOKandVIxSection::SetupKinematic(), G4InuclParticleNames::sp, and G4VMscModel::steppingAlgorithm.

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

Reimplemented from G4VMscModel.

Definition at line 307 of file G4WentzelVIModel.cc.

References DBL_MAX, G4Log(), G4VMscModel::GetEnergy(), G4VMscModel::GetTransportMeanFreePath(), and G4WentzelOKandVIxSection::SetupKinematic().

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

Definition at line 198 of file G4WentzelVIModel.hh.

199 {
200  return fixedCut;
201 }
G4WentzelOKandVIxSection * G4WentzelVIModel::GetWVICrossSection ( )
inline

Definition at line 205 of file G4WentzelVIModel.hh.

206 {
207  return wokvi;
208 }
void G4WentzelVIModel::Initialise ( const G4ParticleDefinition *  p,
const G4DataVector &  cuts 
)
virtual

Implements G4VEmModel.

Definition at line 112 of file G4WentzelVIModel.cc.

References G4VMscModel::GetParticleChangeForMSC(), G4WentzelOKandVIxSection::Initialise(), and G4VEmModel::PolarAngleLimit().

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

Reimplemented from G4VMscModel.

Definition at line 404 of file G4WentzelVIModel.cc.

References DBL_MAX, G4VMscModel::fDisplacement, G4Exp(), G4Log(), G4lrint(), G4UniformRand, G4Material::GetElementVector(), G4Material::GetNumberOfElements(), G4VMscModel::latDisplasment, G4ParticleChangeForMSC::ProposeMomentumDirection(), CLHEP::Hep3Vector::rotateUz(), G4WentzelOKandVIxSection::SampleSingleScattering(), CLHEP::Hep3Vector::set(), G4WentzelOKandVIxSection::SetupTarget(), G4INCL::DeJongSpin::shoot(), python.hepunit::twopi, z, and G4InuclParticleNames::z0.

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

Definition at line 191 of file G4WentzelVIModel.hh.

192 {
193  fixedCut = val;
194 }
void G4WentzelVIModel::StartTracking ( G4Track *  track)
virtual

Reimplemented from G4VEmModel.

Definition at line 168 of file G4WentzelVIModel.cc.

References G4DynamicParticle::GetDefinition(), and G4Track::GetDynamicParticle().

169 {
170  SetupParticle(track->GetDynamicParticle()->GetDefinition());
171  inside = false;
172 }
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:57:21 for Geant4.10 by   doxygen 1.8.7