QwGeant4
QweakSimWentzelVIModel Class Reference

#include <QweakSimWentzelVIModel.hh>

Inherits G4VMscModel.

Public Member Functions

 QweakSimWentzelVIModel (const G4String &nam="WentzelVIUni")
 
virtual ~QweakSimWentzelVIModel ()
 
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 ()
 

Private Member Functions

G4double ComputeXSectionPerVolume ()
 
void SetupParticle (const G4ParticleDefinition *)
 
void DefineMaterial (const G4MaterialCutsCouple *)
 
QweakSimWentzelVIModeloperator= (const QweakSimWentzelVIModel &right)
 
 QweakSimWentzelVIModel (const QweakSimWentzelVIModel &)
 

Private Attributes

G4LossTableManager * theManager
 
G4ParticleChangeForMSC * fParticleChange
 
G4WentzelOKandVIxSection * wokvi
 
const G4DataVector * currentCuts
 
G4double tlimitminfix
 
G4double invsqrt12
 
G4double fixedCut
 
G4double preKinEnergy
 
G4double tPathLength
 
G4double zPathLength
 
G4double lambdaeff
 
G4double currentRange
 
G4double xtsec
 
std::vector< G4double > xsecn
 
std::vector< G4double > prob
 
G4int nelments
 
G4double numlimit
 
G4int currentMaterialIndex
 
const G4MaterialCutsCouple * currentCouple
 
const G4Material * currentMaterial
 
G4double cosThetaMin
 
G4double cosThetaMax
 
G4double cosTetMaxNuc
 
const G4ParticleDefinition * particle
 
G4double lowEnergyLimit
 
G4bool inside
 
G4bool singleScatteringMode
 
G4bool ePolarized
 
G4ThreeVector polarization
 
G4double eEnergy
 
G4bool debugPrint
 

Detailed Description

Definition at line 69 of file QweakSimWentzelVIModel.hh.

Constructor & Destructor Documentation

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

Definition at line 75 of file QweakSimWentzelVIModel.cc.

References cosTetMaxNuc, cosThetaMax, currentCuts, currentMaterial, currentMaterialIndex, currentRange, fixedCut, fParticleChange, invsqrt12, lambdaeff, lowEnergyLimit, nelments, particle, preKinEnergy, prob, theManager, tlimitminfix, tPathLength, wokvi, xsecn, xtsec, and zPathLength.

75  :
76  G4VMscModel(nam),
77  numlimit(0.1),
78  currentCouple(0),
79  cosThetaMin(1.0),
80  inside(false),
82 {
83  invsqrt12 = 1./sqrt(12.);
84  tlimitminfix = 1.e-6*mm;
85  lowEnergyLimit = 1.0*eV;
86  particle = 0;
87  nelments = 5;
88  xsecn.resize(nelments);
89  prob.resize(nelments);
90  theManager = G4LossTableManager::Instance();
91  wokvi = new G4WentzelOKandVIxSection();
92  fixedCut = -1.0;
93 
95  = xtsec = 0;
97  cosThetaMax = cosTetMaxNuc = 1.0;
98 
99  fParticleChange = 0;
100  currentCuts = 0;
101  currentMaterial = 0;
102 }
std::vector< G4double > prob
std::vector< G4double > xsecn
const G4Material * currentMaterial
const G4DataVector * currentCuts
G4LossTableManager * theManager
const G4MaterialCutsCouple * currentCouple
const G4ParticleDefinition * particle
G4ParticleChangeForMSC * fParticleChange
G4WentzelOKandVIxSection * wokvi
QweakSimWentzelVIModel::~QweakSimWentzelVIModel ( )
virtual

Definition at line 106 of file QweakSimWentzelVIModel.cc.

References wokvi.

107 {
108  delete wokvi;
109 }
G4WentzelOKandVIxSection * wokvi
QweakSimWentzelVIModel::QweakSimWentzelVIModel ( const QweakSimWentzelVIModel )
private

Member Function Documentation

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

Definition at line 136 of file QweakSimWentzelVIModel.cc.

References cosTetMaxNuc, currentMaterial, DefineMaterial(), fixedCut, lowEnergyLimit, particle, SetupParticle(), and wokvi.

141 {
142  G4double cross = 0.0;
143  if(p != particle) { SetupParticle(p); }
144  if(kinEnergy < lowEnergyLimit) { return cross; }
145  if(!CurrentCouple()) {
146  G4Exception("QweakSimWentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
147  FatalException, " G4MaterialCutsCouple is not defined");
148  return 0.0;
149  }
150  DefineMaterial(CurrentCouple());
151  cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
152  if(cosTetMaxNuc < 1.0) {
153  G4double cut = cutEnergy;
154  if(fixedCut > 0.0) { cut = fixedCut; }
155  cosTetMaxNuc = wokvi->SetupTarget(G4lrint(Z), cut);
156  cross = wokvi->ComputeTransportCrossSectionPerAtom(cosTetMaxNuc);
157  /*
158  if(p->GetParticleName() == "e-")
159  G4cout << "QweakSimWentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= " << kinEnergy
160  << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
161  << " " << particle->GetParticleName() << G4endl;
162  */
163  }
164  return cross;
165 }
void SetupParticle(const G4ParticleDefinition *)
const G4Material * currentMaterial
void DefineMaterial(const G4MaterialCutsCouple *)
const G4ParticleDefinition * particle
G4WentzelOKandVIxSection * wokvi

+ Here is the call graph for this function:

G4double QweakSimWentzelVIModel::ComputeGeomPathLength ( G4double  truePathLength)
virtual

Definition at line 293 of file QweakSimWentzelVIModel.cc.

References cosTetMaxNuc, currentCouple, currentMaterial, currentRange, lambdaeff, numlimit, particle, preKinEnergy, tPathLength, wokvi, and zPathLength.

294 {
295  tPathLength = truelength;
297 
298  if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) {
299  G4double tau = tPathLength/lambdaeff;
300  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
301  // << " Leff= " << lambdaeff << " tau= " << tau << G4endl;
302  // small step
303  if(tau < numlimit) {
304  zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
305 
306  // medium step
307  } else {
308  G4double e1 = 0.0;
309  if(currentRange > tPathLength) {
311  }
312  e1 = 0.5*(e1 + preKinEnergy);
313  cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
314  lambdaeff = GetTransportMeanFreePath(particle,e1);
315  zPathLength = lambdaeff*(1.0 - G4Exp(-tPathLength/lambdaeff));
316  }
317  } else { lambdaeff = DBL_MAX; }
318  //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
319  return zPathLength;
320 }
const G4Material * currentMaterial
const G4MaterialCutsCouple * currentCouple
const G4ParticleDefinition * particle
G4WentzelOKandVIxSection * wokvi
G4double QweakSimWentzelVIModel::ComputeTruePathLengthLimit ( const G4Track &  track,
G4double &  currentMinimalStep 
)
virtual

Definition at line 177 of file QweakSimWentzelVIModel.cc.

References cosTetMaxNuc, cosThetaMax, currentCouple, currentMaterial, currentRange, debugPrint, DefineMaterial(), eEnergy, ePolarized, inside, lambdaeff, particle, polarization, preKinEnergy, singleScatteringMode, tlimitminfix, and wokvi.

180 {
181  G4double tlimit = currentMinimalStep;
182  const G4DynamicParticle* dp = track.GetDynamicParticle();
183  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
184  G4StepStatus stepStatus = sp->GetStepStatus();
185  singleScatteringMode = false;
186 
187  // FIXME
188  ePolarized=false;
189  debugPrint=false;
190  if(strcmp(track.GetParticleDefinition()->GetParticleName().data() , "e-") == 0)
191  if(strcmp(track.GetMaterial()->GetName(),"PBA") == 0){
192  if(track.GetPolarization().getR() >= 0.1) debugPrint=true;
193  if(sqrt(pow(track.GetPolarization().getX(),2)+pow(track.GetPolarization().getY(),2))>0.01){
194  ePolarized=true;
195  polarization=track.GetPolarization();
196  eEnergy=track.GetTotalEnergy();
197  }
198  }
199  debugPrint=false;
200  // FIXME
201 
202  //G4cout << "QweakSimWentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
203  // << stepStatus << " " << track.GetDefinition()->GetParticleName()
204  // << G4endl;
205 
206  // initialisation for each step, lambda may be computed from scratch
207  preKinEnergy = dp->GetKineticEnergy();
208  DefineMaterial(track.GetMaterialCutsCouple());
209  lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy);
211  cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial);
212 
213  //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
214  // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
215 
216  // extra check for abnormal situation
217  // this check needed to run MSC with eIoni and eBrem inactivated
218  if(tlimit > currentRange) { tlimit = currentRange; }
219 
220  // stop here if small range particle
221  if(inside || tlimit < tlimitminfix) {
222  return ConvertTrueToGeom(tlimit, currentMinimalStep);
223  }
224 
225  // pre step
226  G4double presafety = sp->GetSafety();
227  // far from geometry boundary
228  if(currentRange < presafety) {
229  inside = true;
230  return ConvertTrueToGeom(tlimit, currentMinimalStep);
231  }
232 
233  // compute presafety again if presafety <= 0 and no boundary
234  // i.e. when it is needed for optimization purposes
235  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
236  presafety = ComputeSafety(sp->GetPosition(), tlimit);
237  if(currentRange < presafety) {
238  inside = true;
239  return ConvertTrueToGeom(tlimit, currentMinimalStep);
240  }
241  }
242  /*
243  G4cout << "e(MeV)= " << preKinEnergy/MeV
244  << " " << particle->GetParticleName()
245  << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
246  << " R(mm)= " <<currentRange/mm
247  << " L0(mm^-1)= " << lambdaeff*mm
248  <<G4endl;
249  */
250 
251  // natural limit for high energy
252  G4double rlimit = std::max(facrange*currentRange,
253  0.7*(1.0 - cosTetMaxNuc)*lambdaeff);
254 
255  // low-energy e-
256  if(cosThetaMax > cosTetMaxNuc) {
257  rlimit = std::min(rlimit, facsafety*presafety);
258  }
259 
260  // cut correction
261  G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
262  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= "
263  // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax
264  //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
265  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
266 
267  if(rlimit < tlimit) { tlimit = rlimit; }
268 
269  tlimit = std::max(tlimit, tlimitminfix);
270 
271  // step limit in infinite media
272  tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
273 
274  //compute geomlimit and force few steps within a volume
275  if (steppingAlgorithm == fUseDistanceToBoundary
276  && stepStatus == fGeomBoundary) {
277 
278  G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
279  tlimit = std::min(tlimit, geomlimit/facgeom);
280  }
281 
282  /*
283  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
284  << " L0= " << lambdaeff << " R= " << currentRange
285  << "tlimit= " << tlimit
286  << " currentMinimalStep= " << currentMinimalStep << G4endl;
287  */
288  return ConvertTrueToGeom(tlimit, currentMinimalStep);
289 }
const G4Material * currentMaterial
const G4MaterialCutsCouple * currentCouple
void DefineMaterial(const G4MaterialCutsCouple *)
const G4ParticleDefinition * particle
G4WentzelOKandVIxSection * wokvi

+ Here is the call graph for this function:

G4double QweakSimWentzelVIModel::ComputeTrueStepLength ( G4double  geomStepLength)
virtual

Definition at line 324 of file QweakSimWentzelVIModel.cc.

References ComputeXSectionPerVolume(), cosTetMaxNuc, cosThetaMin, currentCouple, currentMaterial, currentRange, lambdaeff, numlimit, particle, preKinEnergy, singleScatteringMode, tPathLength, wokvi, xtsec, and zPathLength.

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

+ Here is the call graph for this function:

G4double QweakSimWentzelVIModel::ComputeXSectionPerVolume ( )
private

Definition at line 650 of file QweakSimWentzelVIModel.cc.

References cosTetMaxNuc, cosThetaMin, currentMaterial, currentMaterialIndex, fixedCut, nelments, prob, wokvi, xsecn, and xtsec.

Referenced by ComputeTrueStepLength().

651 {
652  // prepare recomputation of x-sections
653  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
654  const G4double* theAtomNumDensityVector =
655  currentMaterial->GetVecNbOfAtomsPerVolume();
656  G4int nelm = currentMaterial->GetNumberOfElements();
657  if(nelm > nelments) {
658  nelments = nelm;
659  xsecn.resize(nelm);
660  prob.resize(nelm);
661  }
662  G4double cut = (*currentCuts)[currentMaterialIndex];
663  if(fixedCut > 0.0) { cut = fixedCut; }
664  // cosTetMaxNuc = wokvi->GetCosThetaNuc();
665 
666  // check consistency
667  xtsec = 0.0;
668  if(cosTetMaxNuc > cosThetaMin) { return 0.0; }
669 
670  // loop over elements
671  G4double xs = 0.0;
672  for (G4int i=0; i<nelm; ++i) {
673  G4double costm =
674  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
675  G4double density = theAtomNumDensityVector[i];
676 
677  G4double esec = 0.0;
678  if(costm < cosThetaMin) {
679 
680  // recompute the transport x-section
681  if(1.0 > cosThetaMin) {
682  xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosThetaMin);
683  }
684  // recompute the total x-section
685  G4double nucsec = wokvi->ComputeNuclearCrossSection(cosThetaMin, costm);
686  esec = wokvi->ComputeElectronCrossSection(cosThetaMin, costm);
687  nucsec += esec;
688  if(nucsec > 0.0) { esec /= nucsec; }
689  xtsec += nucsec*density;
690  }
691  xsecn[i] = xtsec;
692  prob[i] = esec;
693  //G4cout << i << " xs= " << xs << " xtsec= " << xtsec
694  // << " 1-cosThetaMin= " << 1-cosThetaMin
695  // << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
696  }
697 
698  //G4cout << "ComputeXS result: xsec(1/mm)= " << xs
699  // << " txsec(1/mm)= " << xtsec <<G4endl;
700  return xs;
701 }
std::vector< G4double > prob
std::vector< G4double > xsecn
const G4Material * currentMaterial
G4WentzelOKandVIxSection * wokvi

+ Here is the caller graph for this function:

void QweakSimWentzelVIModel::DefineMaterial ( const G4MaterialCutsCouple *  cup)
inlineprivate

Definition at line 174 of file QweakSimWentzelVIModel.hh.

References currentCouple, currentMaterial, and currentMaterialIndex.

Referenced by ComputeCrossSectionPerAtom(), and ComputeTruePathLengthLimit().

175 {
176  if(cup != currentCouple) {
177  currentCouple = cup;
178  SetCurrentCouple(cup);
179  currentMaterial = cup->GetMaterial();
180  currentMaterialIndex = currentCouple->GetIndex();
181  }
182 }
const G4Material * currentMaterial
const G4MaterialCutsCouple * currentCouple

+ Here is the caller graph for this function:

G4double QweakSimWentzelVIModel::GetFixedCut ( ) const
inline

Definition at line 204 of file QweakSimWentzelVIModel.hh.

References fixedCut.

205 {
206  return fixedCut;
207 }
G4WentzelOKandVIxSection * QweakSimWentzelVIModel::GetWVICrossSection ( )
inline

Definition at line 211 of file QweakSimWentzelVIModel.hh.

References wokvi.

212 {
213  return wokvi;
214 }
G4WentzelOKandVIxSection * wokvi
void QweakSimWentzelVIModel::Initialise ( const G4ParticleDefinition *  p,
const G4DataVector &  cuts 
)
virtual

Definition at line 113 of file QweakSimWentzelVIModel.cc.

References cosThetaMax, currentCuts, currentRange, fParticleChange, SetupParticle(), and wokvi.

115 {
116  // reset parameters
117  SetupParticle(p);
118  currentRange = 0.0;
119 
120  cosThetaMax = cos(PolarAngleLimit());
121  //G4cout << "QweakSimWentzelVIModel::Initialise " << p->GetParticleName() << G4endl;
122  wokvi->Initialise(p, cosThetaMax);
123  /*
124  G4cout << "QweakSimWentzelVIModel: " << particle->GetParticleName()
125  << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax
126  << G4endl;
127  */
128  currentCuts = &cuts;
129 
130  // set values of some data members
131  fParticleChange = GetParticleChangeForMSC(p);
132 }
void SetupParticle(const G4ParticleDefinition *)
const G4DataVector * currentCuts
G4ParticleChangeForMSC * fParticleChange
G4WentzelOKandVIxSection * wokvi

+ Here is the call graph for this function:

QweakSimWentzelVIModel& QweakSimWentzelVIModel::operator= ( const QweakSimWentzelVIModel right)
private
G4ThreeVector & QweakSimWentzelVIModel::SampleScattering ( const G4ThreeVector &  oldDirection,
G4double  safety 
)
virtual

Definition at line 421 of file QweakSimWentzelVIModel.cc.

References cosThetaMin, currentMaterial, currentMaterialIndex, debugPrint, eEnergy, ePolarized, fixedCut, fParticleChange, invsqrt12, lambdaeff, lowEnergyLimit, polarization, preKinEnergy, prob, singleScatteringMode, tlimitminfix, tPathLength, wokvi, xsecn, xtsec, and zPathLength.

423 {
424  fDisplacement.set(0.0,0.0,0.0);
425  //G4cout << "!##! QweakSimWentzelVIModel::SampleScattering for "
426  // << particle->GetParticleName() << G4endl;
427 
428  // ignore scattering for zero step length and energy below the limit
429  if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0)
430  { return fDisplacement; }
431 
432  G4double invlambda = 0.0;
433  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
434 
435  // use average kinetic energy over the step
436  G4double cut = (*currentCuts)[currentMaterialIndex];
437  if(fixedCut > 0.0) { cut = fixedCut; }
438  /*
439  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
440  << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
441  << " xmsc= " << tPathLength*invlambda
442  << " safety= " << safety << G4endl;
443  */
444 
445  // step limit due msc
446  G4int nMscSteps = 1;
447  G4double x0 = tPathLength;
448  G4double z0 = x0*invlambda;
449  G4double zzz = 0.0;
450 
451  // large scattering angle case - two step approach
452 
453  if(!singleScatteringMode) {
454  static const G4double zzmin = 0.05;
455  if(z0 > zzmin && safety > tlimitminfix) {
456  x0 *= 0.5;
457  z0 *= 0.5;
458  nMscSteps = 2;
459  }
460  if(z0 > zzmin) { zzz = G4Exp(-1.0/z0); }
461  }
462 
463  // step limit due to single scattering
464  G4double x1 = 2*tPathLength;
465  if(0.0 < xtsec) { x1 = -G4Log(G4UniformRand())/xtsec; }
466 
467  // no scattering case
468  if(singleScatteringMode && x1 > tPathLength)
469  { return fDisplacement; }
470 
471  const G4ElementVector* theElementVector =
472  currentMaterial->GetElementVector();
473  G4int nelm = currentMaterial->GetNumberOfElements();
474 
475  // geometry
476  G4double sint, cost, phi;
477  G4ThreeVector temp(0.0,0.0,1.0);
478 
479  // current position and direction relative to the end point
480  // because of magnetic field geometry is computed relatively to the
481  // end point of the step
482  G4ThreeVector dir(0.0,0.0,1.0);
483  fDisplacement.set(0.0,0.0,-zPathLength);
484 
485  G4double mscfac = zPathLength/tPathLength;
486 
487  // start a loop
488  G4double x2 = x0;
489  G4double step, z;
490  G4bool singleScat;
491  /*
492  G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
493  << " 1-cost1= " << 1 - cosThetaMin << " " << singleScatteringMode
494  << " xtsec= " << xtsec << " " << nMscSteps << G4endl;
495  */
496  do {
497 
498  //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl;
499  // single scattering case
500  if(singleScatteringMode && x1 > x2) { break; }
501 
502  // what is next single of multiple?
503  if(x1 <= x2) {
504  step = x1;
505  singleScat = true;
506  } else {
507  step = x2;
508  singleScat = false;
509  }
510 
511  //G4cout << "# step(mm)= "<< step<< " singlScat= "<< singleScat << G4endl;
512 
513  // new position
514  fDisplacement += step*mscfac*dir;
515 
516  if(singleScat) {
517 
518  // select element
519  G4int i = 0;
520  if(nelm > 1) {
521  G4double qsec = G4UniformRand()*xtsec;
522  for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
523  }
524  G4double cosTetM =
525  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
526  //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " "
527  // << prob[i] << G4endl;
528  temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
529 
530  //FIXME
531  cost=cos(temp.getTheta());
532  sint = sqrt((1.0 - cost)*(1.0 + cost));
533  phi=temp.getPhi();
534  if(ePolarized){
535  if(debugPrint) G4cout<<" ~~sc~~ "<<nMscSteps<<G4endl;
536  G4double _prob=G4UniformRand();
537  G4double _amplitude=1.0/eEnergy * sint *
538  sqrt(pow(polarization.getX(),2)+pow(polarization.getY(),2));//scale by transvers polarization
539  if(_amplitude > 1 ) _amplitude=1;
540  if( _prob < _amplitude * sin(phi-pi) )
541  phi-=pi;
542  phi+= polarization.getPhi() - oldDirection.getPhi();
543  if(phi<0) phi+=twopi;
544  else if(phi>twopi) phi=fmod(phi,twopi);
545  G4double vx1 = sint*cos(phi);
546  G4double vy1 = sint*sin(phi);
547  temp.set(vx1,vy1,cost);
548  }
549  //FIXME
550 
551  // direction is changed
552  temp.rotateUz(dir);
553  dir = temp;
554 
555  // new proposed step length
556  x2 -= step;
557  x1 = -G4Log(G4UniformRand())/xtsec;
558 
559  // multiple scattering
560  } else {
561  --nMscSteps;
562  x1 -= step;
563  x2 = x0;
564 
565  // sample z in interval 0 - 1
566  do {
567  z = -z0*G4Log(1.0 - (1.0 - zzz)*G4UniformRand());
568  } while(z > 1.0);
569  cost = 1.0 - 2.0*z/*factCM*/;
570  if(cost > 1.0) { cost = 1.0; }
571  else if(cost < -1.0) { cost =-1.0; }
572  sint = sqrt((1.0 - cost)*(1.0 + cost));
573  phi = twopi*G4UniformRand();
574  //FIXME
575  if(ePolarized){
576  if(debugPrint) G4cout<<" ~~msc~~ "<<nMscSteps<<G4endl;
577  G4double _prob=G4UniformRand();
578  G4double _amplitude=1.0/eEnergy * sint*
579  sqrt(pow(polarization.getX(),2)+pow(polarization.getY(),2));//scale by transvers polarization
580  if(_amplitude > 1 ) _amplitude=1;
581  if( _prob < _amplitude * sin(phi-pi) )
582  phi-=pi;
583  phi+= polarization.getPhi() - oldDirection.getPhi();
584  if(phi<0) phi+=twopi;
585  else if(phi>twopi) phi=fmod(phi,twopi);
586  }
587  //FIXME
588  G4double vx1 = sint*cos(phi);
589  G4double vy1 = sint*sin(phi);
590 
591  // lateral displacement
592  if (latDisplasment && safety > tlimitminfix) {
593  G4double rms = invsqrt12*sqrt(2*z0);
594  G4double r = x0*mscfac;
595  G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0));
596  G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0));
597  G4double dz;
598  G4double d = r*r - dx*dx - dy*dy;
599 
600  // change position
601  if(d >= 0.0) {
602  dz = sqrt(d) - r;
603  temp.set(dx,dy,dz);
604  temp.rotateUz(dir);
605  fDisplacement += temp;
606  }
607  }
608  // change direction
609  temp.set(vx1,vy1,cost);
610  temp.rotateUz(dir);
611  dir = temp;
612  }
613  } while (0 < nMscSteps);
614 
615  dir.rotateUz(oldDirection);
616 
617  //FIXME
618  if(debugPrint){
619  G4cout<<" WentzelMS cth, th, phi old.angle(new)" << cost << " " << acos(cost) << " " << phi << " " <<oldDirection.angle(dir) << G4endl;
620  G4cout<<" WentzelMS old dir: R th phi "<<oldDirection.getR()<<" "<<oldDirection.getTheta()<<" "<<oldDirection.getPhi()<<G4endl;
621  G4cout<<" WentzelMS new dir: R th phi "<<dir.getR()<<" "<<dir.getTheta()<<" "<<dir.getPhi()<<G4endl;
622  }
623  //FIXME
624 
625  //G4cout<<"QweakSimWentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl;
626  // end of sampling -------------------------------
627 
628  fParticleChange->ProposeMomentumDirection(dir);
629 
630  // lateral displacement
631  fDisplacement.rotateUz(oldDirection);
632 
633  /*
634  G4cout << " r(mm)= " << fDisplacement.mag()
635  << " safety= " << safety
636  << " trueStep(mm)= " << tPathLength
637  << " geomStep(mm)= " << zPathLength
638  << " x= " << fDisplacement.x()
639  << " y= " << fDisplacement.y()
640  << " z= " << fDisplacement.z()
641  << G4endl;
642  */
643 
644  //G4cout<< "QweakSimWentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
645  return fDisplacement;
646 }
std::vector< G4double > prob
std::vector< G4double > xsecn
const G4Material * currentMaterial
G4ParticleChangeForMSC * fParticleChange
G4WentzelOKandVIxSection * wokvi
void QweakSimWentzelVIModel::SetFixedCut ( G4double  val)
inline

Definition at line 197 of file QweakSimWentzelVIModel.hh.

References fixedCut.

198 {
199  fixedCut = val;
200 }
void QweakSimWentzelVIModel::SetupParticle ( const G4ParticleDefinition *  p)
inlineprivate

Definition at line 186 of file QweakSimWentzelVIModel.hh.

References particle, and wokvi.

Referenced by ComputeCrossSectionPerAtom(), Initialise(), and StartTracking().

187 {
188  // Initialise mass and charge
189  if(p != particle) {
190  particle = p;
191  wokvi->SetupParticle(p);
192  }
193 }
const G4ParticleDefinition * particle
G4WentzelOKandVIxSection * wokvi

+ Here is the caller graph for this function:

void QweakSimWentzelVIModel::StartTracking ( G4Track *  track)

Definition at line 169 of file QweakSimWentzelVIModel.cc.

References inside, and SetupParticle().

170 {
171  SetupParticle(track->GetDynamicParticle()->GetDefinition());
172  inside = false;
173 }
void SetupParticle(const G4ParticleDefinition *)

+ Here is the call graph for this function:

Field Documentation

G4double QweakSimWentzelVIModel::cosThetaMax
private
G4double QweakSimWentzelVIModel::cosThetaMin
private
const G4MaterialCutsCouple* QweakSimWentzelVIModel::currentCouple
private
const G4DataVector* QweakSimWentzelVIModel::currentCuts
private

Definition at line 124 of file QweakSimWentzelVIModel.hh.

Referenced by Initialise(), and QweakSimWentzelVIModel().

G4int QweakSimWentzelVIModel::currentMaterialIndex
private
G4double QweakSimWentzelVIModel::currentRange
private
G4bool QweakSimWentzelVIModel::debugPrint
private

Definition at line 166 of file QweakSimWentzelVIModel.hh.

Referenced by ComputeTruePathLengthLimit(), and SampleScattering().

G4double QweakSimWentzelVIModel::eEnergy
private

Definition at line 165 of file QweakSimWentzelVIModel.hh.

Referenced by ComputeTruePathLengthLimit(), and SampleScattering().

G4bool QweakSimWentzelVIModel::ePolarized
private

Definition at line 163 of file QweakSimWentzelVIModel.hh.

Referenced by ComputeTruePathLengthLimit(), and SampleScattering().

G4double QweakSimWentzelVIModel::fixedCut
private
G4ParticleChangeForMSC* QweakSimWentzelVIModel::fParticleChange
private

Definition at line 121 of file QweakSimWentzelVIModel.hh.

Referenced by Initialise(), QweakSimWentzelVIModel(), and SampleScattering().

G4bool QweakSimWentzelVIModel::inside
private

Definition at line 160 of file QweakSimWentzelVIModel.hh.

Referenced by ComputeTruePathLengthLimit(), and StartTracking().

G4double QweakSimWentzelVIModel::invsqrt12
private

Definition at line 127 of file QweakSimWentzelVIModel.hh.

Referenced by QweakSimWentzelVIModel(), and SampleScattering().

G4double QweakSimWentzelVIModel::lambdaeff
private
G4double QweakSimWentzelVIModel::lowEnergyLimit
private
G4int QweakSimWentzelVIModel::nelments
private

Definition at line 141 of file QweakSimWentzelVIModel.hh.

Referenced by ComputeXSectionPerVolume(), and QweakSimWentzelVIModel().

G4double QweakSimWentzelVIModel::numlimit
private

Definition at line 143 of file QweakSimWentzelVIModel.hh.

Referenced by ComputeGeomPathLength(), and ComputeTrueStepLength().

const G4ParticleDefinition* QweakSimWentzelVIModel::particle
private
G4ThreeVector QweakSimWentzelVIModel::polarization
private

Definition at line 164 of file QweakSimWentzelVIModel.hh.

Referenced by ComputeTruePathLengthLimit(), and SampleScattering().

G4double QweakSimWentzelVIModel::preKinEnergy
private
std::vector<G4double> QweakSimWentzelVIModel::prob
private
G4bool QweakSimWentzelVIModel::singleScatteringMode
private
G4LossTableManager* QweakSimWentzelVIModel::theManager
private

Definition at line 120 of file QweakSimWentzelVIModel.hh.

Referenced by QweakSimWentzelVIModel().

G4double QweakSimWentzelVIModel::tlimitminfix
private
G4double QweakSimWentzelVIModel::tPathLength
private
std::vector<G4double> QweakSimWentzelVIModel::xsecn
private
G4double QweakSimWentzelVIModel::xtsec
private
G4double QweakSimWentzelVIModel::zPathLength
private

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