ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4WentzelVIModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4WentzelVIModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4WentzelVIModel
33 //
34 // Author: V.Ivanchenko
35 //
36 // Creation date: 09.04.2008 from G4MuMscModel
37 //
38 // Modifications:
39 // 27-05-2010 V.Ivanchenko added G4WentzelOKandVIxSection class to
40 // compute cross sections and sample scattering angle
41 //
42 //
43 // Class Description:
44 //
45 // Implementation of the model of multiple scattering based on
46 // G.Wentzel, Z. Phys. 40 (1927) 590.
47 // H.W.Lewis, Phys Rev 78 (1950) 526.
48 // J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
49 // L.Urban, CERN-OPEN-2006-077.
50 
51 // -------------------------------------------------------------------
52 //
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 #include "G4WentzelVIModel.hh"
58 #include "G4PhysicalConstants.hh"
59 #include "G4SystemOfUnits.hh"
60 #include "Randomize.hh"
62 #include "G4PhysicsTableHelper.hh"
63 #include "G4ElementVector.hh"
64 #include "G4ProductionCutsTable.hh"
65 #include "G4EmParameters.hh"
66 #include "G4Log.hh"
67 #include "G4Exp.hh"
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 
71 using namespace std;
72 
74  : G4VMscModel(nam),
75  ssFactor(1.05),
76  invssFactor(1.0),
77  currentCouple(nullptr),
78  cosThetaMin(1.0),
79  cosThetaMax(-1.0),
80  fSecondMoments(nullptr),
81  idx2(0),
82  numlimit(0.1),
83  singleScatteringMode(false),
84  isCombined(comb),
85  useSecondMoment(false)
86 {
88  invsqrt12 = 1./sqrt(12.);
89  tlimitminfix = 1.e-6*mm;
90  lowEnergyLimit = 1.0*eV;
91  particle = nullptr;
92  nelments = 5;
93  xsecn.resize(nelments);
94  prob.resize(nelments);
96  fixedCut = -1.0;
97 
98  minNCollisions = 10;
99 
101  = currentRange = xtsec = cosTetMaxNuc = 0.0;
103 
104  fParticleChange = nullptr;
105  currentCuts = nullptr;
106  currentMaterial = nullptr;
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
112 {
113  delete wokvi;
114  if(fSecondMoments && IsMaster()) {
115  delete fSecondMoments;
116  fSecondMoments = nullptr;
117  }
118 }
119 
120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121 
123  const G4DataVector& cuts)
124 {
125  // reset parameters
126  SetupParticle(p);
128  currentRange = 0.0;
129 
130  if(isCombined) {
131  G4double tet = PolarAngleLimit();
132  if(tet <= 0.0) { cosThetaMax = 1.0; }
133  else if(tet < CLHEP::pi) { cosThetaMax = cos(tet); }
134  }
135  //G4cout << "G4WentzelVIModel::Initialise " << p->GetParticleName()
136  // << " " << this << " " << wokvi << G4endl;
137 
139  /*
140  G4cout << "G4WentzelVIModel: " << particle->GetParticleName()
141  << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax
142  << " SingScatFactor= " << ssFactor
143  << G4endl;
144  */
145  currentCuts = &cuts;
146 
147  // set values of some data members
149 
150  // build second moment table only if transport table is build
152  if(useSecondMoment && IsMaster() && table) {
153 
154  //G4cout << "### G4WentzelVIModel::Initialise: build 2nd moment table "
155  // << table << G4endl;
156  fSecondMoments =
158  // Access to materials
159  const G4ProductionCutsTable* theCoupleTable =
161  size_t numOfCouples = theCoupleTable->GetTableSize();
162 
163  G4bool splineFlag = true;
164  G4PhysicsVector* aVector = nullptr;
165  G4PhysicsVector* bVector = nullptr;
168  if(emin < emax) {
170  *G4lrint(std::log10(emax/emin));
171  if(n < 3) { n = 3; }
172 
173  for(size_t i=0; i<numOfCouples; ++i) {
174 
175  //G4cout<< "i= " << i << " Flag= " << fSecondMoments->GetFlag(i)
176  // << G4endl;
177  if(fSecondMoments->GetFlag(i)) {
178  DefineMaterial(theCoupleTable->GetMaterialCutsCouple(i));
179 
180  delete (*fSecondMoments)[i];
181  if(!aVector) {
182  aVector = new G4PhysicsLogVector(emin, emax, n);
183  bVector = aVector;
184  } else {
185  bVector = new G4PhysicsVector(*aVector);
186  }
187  for(size_t j=0; j<n; ++j) {
188  G4double e = bVector->Energy(j);
189  bVector->PutValue(j, ComputeSecondMoment(p, e)*e*e);
190  }
191  if(splineFlag) { bVector->FillSecondDerivatives(); }
192  (*fSecondMoments)[i] = bVector;
193  }
194  }
195  }
196  //G4cout << *fSecondMoments << G4endl;
197  }
198 }
199 
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
201 
203  G4VEmModel* masterModel)
204 {
205  fSecondMoments = static_cast<G4WentzelVIModel*>(masterModel)
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210 
212 {
213  if(cup != currentCouple) {
214  currentCouple = cup;
215  SetCurrentCouple(cup);
216  currentMaterial = cup->GetMaterial();
218  }
219 }
220 
221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
222 
224  const G4ParticleDefinition* p,
225  G4double kinEnergy,
227  G4double cutEnergy, G4double)
228 {
229  G4double cross = 0.0;
230  SetupParticle(p);
231  if(kinEnergy < lowEnergyLimit) { return cross; }
232  if(!CurrentCouple()) {
233  G4Exception("G4WentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
234  FatalException, " G4MaterialCutsCouple is not defined");
235  return 0.0;
236  }
239  if(cosTetMaxNuc < 1.0) {
240  G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy;
241  G4double cost = wokvi->SetupTarget(G4lrint(Z), cut);
243  /*
244  if(p->GetParticleName() == "e-")
245  G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= "<<kinEnergy
246  << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
247  << " " << particle->GetParticleName() << G4endl;
248  */
249  }
250  return cross;
251 }
252 
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254 
256 {
257  /*
258  G4cout << "G4WentzelVIModel::StartTracking " << track << " " << this << " "
259  << track->GetParticleDefinition()->GetParticleName()
260  << " workvi: " << wokvi << G4endl;
261  */
263 }
264 
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266 
268  const G4Track& track,
269  G4double& currentMinimalStep)
270 {
271  G4double tlimit = currentMinimalStep;
272  const G4DynamicParticle* dp = track.GetDynamicParticle();
273  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
274  G4StepStatus stepStatus = sp->GetStepStatus();
275  singleScatteringMode = false;
276 
277  //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
278  // << stepStatus << " " << track.GetDefinition()->GetParticleName()
279  // << G4endl;
280 
281  // initialisation for each step, lambda may be computed from scratch
285  const G4double logPreKinEnergy = dp->GetLogKineticEnergy();
289 
290  //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
291  // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
292 
293  // extra check for abnormal situation
294  // this check needed to run MSC with eIoni and eBrem inactivated
295  if(tlimit > currentRange) { tlimit = currentRange; }
296 
297  // stop here if small range particle
298  if(tlimit < tlimitminfix) {
299  return ConvertTrueToGeom(tlimit, currentMinimalStep);
300  }
301 
302  // pre step
303  G4double presafety = sp->GetSafety();
304  // far from geometry boundary
305  if(currentRange < presafety) {
306  return ConvertTrueToGeom(tlimit, currentMinimalStep);
307  }
308 
309  // compute presafety again if presafety <= 0 and no boundary
310  // i.e. when it is needed for optimization purposes
311  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
312  presafety = ComputeSafety(sp->GetPosition(), tlimit);
313  if(currentRange < presafety) {
314  return ConvertTrueToGeom(tlimit, currentMinimalStep);
315  }
316  }
317  /*
318  G4cout << "e(MeV)= " << preKinEnergy/MeV
319  << " " << particle->GetParticleName()
320  << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
321  << " R(mm)= " <<currentRange/mm
322  << " L0(mm^-1)= " << lambdaeff*mm
323  << G4endl;
324  */
325  // natural limit for high energy
328 
329  // low-energy e-
330  if(cosThetaMax > cosTetMaxNuc) {
331  rlimit = std::min(rlimit, facsafety*presafety);
332  }
333 
334  // cut correction
336  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= "
337  // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax
338  //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
339  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
340 
341  tlimit = std::min(tlimit, rlimit);
342  tlimit = std::max(tlimit, tlimitminfix);
343 
344  // step limit in infinite media
345  tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
346 
347  //compute geomlimit and force few steps within a volume
349  && stepStatus == fGeomBoundary) {
350 
351  G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
352  tlimit = std::min(tlimit, geomlimit/facgeom);
353  }
354  /*
355  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
356  << " L0= " << lambdaeff << " R= " << currentRange
357  << " tlimit= " << tlimit
358  << " currentMinimalStep= " << currentMinimalStep << G4endl;
359  */
360  return ConvertTrueToGeom(tlimit, currentMinimalStep);
361 }
362 
363 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
364 
366 {
367  zPathLength = tPathLength = truelength;
368 
369  // small step use only single scattering
370  cosThetaMin = 1.0;
372  //G4cout << "xtsec= " << xtsec << " Nav= "
373  // << zPathLength*xtsec << G4endl;
374  if(0.0 >= lambdaeff || G4int(zPathLength*xtsec) < minNCollisions) {
375  singleScatteringMode = true;
376  lambdaeff = DBL_MAX;
377 
378  } else {
379  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
380  // << " Leff= " << lambdaeff << G4endl;
381  // small step
384  zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
385 
386  // medium step
387  } else {
388  G4double e1 = 0.0;
389  if(currentRange > tPathLength) {
391  }
392  effKinEnergy = 0.5*(e1 + preKinEnergy);
395  //G4cout << " tLength= "<< tPathLength<< " Leff= " << lambdaeff << G4endl;
397  if(tPathLength*numlimit < lambdaeff) {
398  zPathLength *= (1.0 - G4Exp(-tPathLength/lambdaeff));
399  }
400  }
401  }
402  //G4cout << "Comp.geom: zLength= "<<zPathLength<<" tLength= "
403  // << tPathLength<< " Leff= " << lambdaeff << G4endl;
404  return zPathLength;
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408 
410 {
411  // initialisation of single scattering x-section
412  /*
413  G4cout << "ComputeTrueStepLength: Step= " << geomStepLength
414  << " geomL= " << zPathLength
415  << " Lambda= " << lambdaeff
416  << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
417  */
419  zPathLength = tPathLength = geomStepLength;
420 
421  } else {
422 
423  // step defined by transportation
424  // change both geom and true step lengths
425  if(geomStepLength < zPathLength) {
426 
427  // single scattering
428  if(G4int(geomStepLength*xtsec) < minNCollisions) {
429  zPathLength = tPathLength = geomStepLength;
430  lambdaeff = DBL_MAX;
431  singleScatteringMode = true;
432 
433  // multiple scattering
434  } else {
435  // small step
436  if(geomStepLength < numlimit*lambdaeff) {
437  G4double tau = geomStepLength/lambdaeff;
438  tPathLength = geomStepLength*(1.0 + 0.5*tau + tau*tau/3.0);
439 
440  // energy correction for a big step
441  } else {
442  tPathLength *= geomStepLength/zPathLength;
443  G4double e1 = 0.0;
444  if(currentRange > tPathLength) {
446  }
447  effKinEnergy = 0.5*(e1 + preKinEnergy);
450  G4double tau = geomStepLength/lambdaeff;
451 
452  if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
453  else { tPathLength = currentRange; }
454  }
455  zPathLength = geomStepLength;
456  }
457  }
458  }
459  // check of step length
460  // define threshold angle between single and multiple scattering
461  if(!singleScatteringMode) {
463  xtsec = 0.0;
464 
465  // recompute transport cross section - do not change energy
466  // anymore - cannot be applied for big steps
467  if(cosThetaMin > cosTetMaxNuc) {
468  // new computation
470  //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec
471  // << " 1-cosTMin= " << 1.0 - cosThetaMin << G4endl;
472  if(cross <= 0.0) {
473  singleScatteringMode = true;
475  lambdaeff = DBL_MAX;
476  cosThetaMin = 1.0;
477  } else if(xtsec > 0.0) {
478 
479  lambdaeff = 1./cross;
480  G4double tau = zPathLength*cross;
481  if(tau < numlimit) {
482  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
483  } else if(tau < 0.999999) {
484  tPathLength = -lambdaeff*G4Log(1.0 - tau);
485  } else {
487  }
488  }
489  }
490  }
492  /*
493  G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
494  <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
495  G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
496  << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
497  << " e(MeV)= " << preKinEnergy/MeV << " "
498  << " SSmode= " << singleScatteringMode << G4endl;
499  */
500  return tPathLength;
501 }
502 
503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
504 
507  G4double /*safety*/)
508 {
509  fDisplacement.set(0.0,0.0,0.0);
510  //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
511  // << particle->GetParticleName() << G4endl;
512 
513  // ignore scattering for zero step length and energy below the limit
514  if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0)
515  { return fDisplacement; }
516 
517  G4double invlambda = 0.0;
518  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
519 
520  // use average kinetic energy over the step
521  G4double cut = (*currentCuts)[currentMaterialIndex];
522  if(fixedCut > 0.0) { cut = fixedCut; }
523  /*
524  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
525  << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
526  << " xmsc= " << tPathLength*invlambda
527  << " safety= " << safety << G4endl;
528  */
529  // step limit due msc
530  G4int nMscSteps = 1;
531  G4double x0 = tPathLength;
532  G4double z0 = x0*invlambda;
533  //G4double zzz = 0.0;
534  G4double prob2 = 0.0;
535 
536  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
537 
538  // large scattering angle case - two step approach
539  if(!singleScatteringMode) {
540  static const G4double zzmin = 0.05;
541  if(useSecondMoment) {
542  G4double z1 = invlambda*invlambda;
544  prob2 = (z2 - z1)/(1.5*z1 - z2);
545  }
546  // if(z0 > zzmin && safety > tlimitminfix) {
547  if(z0 > zzmin) {
548  x0 *= 0.5;
549  z0 *= 0.5;
550  nMscSteps = 2;
551  }
552  //if(z0 > zzmin) { zzz = G4Exp(-1.0/z0); }
553  G4double zzz = 0.0;
554  if(z0 > zzmin) {
555  zzz = G4Exp(-1.0/z0);
556  z0 += zzz;
557  prob2 *= (1 + zzz);
558  }
559  prob2 /= (1 + prob2);
560  }
561 
562  // step limit due to single scattering
563  G4double x1 = 2*tPathLength;
564  if(0.0 < xtsec) { x1 = -G4Log(rndmEngine->flat())/xtsec; }
565 
566  // no scattering case
567  if(singleScatteringMode && x1 > tPathLength)
568  { return fDisplacement; }
569 
570  const G4ElementVector* theElementVector =
573 
574  // geometry
575  G4double sint, cost, phi;
576  G4ThreeVector temp(0.0,0.0,1.0);
577 
578  // current position and direction relative to the end point
579  // because of magnetic field geometry is computed relatively to the
580  // end point of the step
581  G4ThreeVector dir(0.0,0.0,1.0);
582  fDisplacement.set(0.0,0.0,-zPathLength);
583 
585 
586  // start a loop
587  G4double x2 = x0;
588  G4double step, z;
589  G4bool singleScat;
590  /*
591  G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
592  << " 1-cost1= " << 1 - cosThetaMin << " SSmode= " << singleScatteringMode
593  << " xtsec= " << xtsec << " Nst= " << nMscSteps << G4endl;
594  */
595  do {
596 
597  //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl;
598  // single scattering case
599  if(singleScatteringMode && x1 > x2) {
600  fDisplacement += x2*mscfac*dir;
601  break;
602  }
603 
604  // what is next single of multiple?
605  if(x1 <= x2) {
606  step = x1;
607  singleScat = true;
608  } else {
609  step = x2;
610  singleScat = false;
611  }
612 
613  //G4cout << "# step(mm)= "<< step<< " singlScat= "<< singleScat << G4endl;
614 
615  // new position
616  fDisplacement += step*mscfac*dir;
617 
618  if(singleScat) {
619 
620  // select element
621  G4int i = 0;
622  if(nelm > 1) {
623  G4double qsec = rndmEngine->flat()*xtsec;
624  for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
625  }
626  G4double cosTetM =
627  wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
628  //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " "
629  // << prob[i] << G4endl;
630  temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
631 
632  // direction is changed
633  temp.rotateUz(dir);
634  dir = temp;
635  //G4cout << dir << G4endl;
636 
637  // new proposed step length
638  x2 -= step;
639  x1 = -G4Log(rndmEngine->flat())/xtsec;
640 
641  // multiple scattering
642  } else {
643  --nMscSteps;
644  x1 -= step;
645  x2 = x0;
646 
647  // sample z in interval 0 - 1
648  G4bool isFirst = true;
649  if(prob2 > 0.0 && rndmEngine->flat() < prob2) { isFirst = false; }
650  do {
651  //z = -z0*G4Log(1.0 - (1.0 - zzz)*rndmEngine->flat());
652  if(isFirst) { z = -G4Log(rndmEngine->flat()); }
653  else { z = G4RandGamma::shoot(rndmEngine, 2.0, 2.0); }
654  z *= z0;
655  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
656  } while(z > 1.0);
657 
658  cost = 1.0 - 2.0*z/*factCM*/;
659  if(cost > 1.0) { cost = 1.0; }
660  else if(cost < -1.0) { cost =-1.0; }
661  sint = sqrt((1.0 - cost)*(1.0 + cost));
662  phi = twopi*rndmEngine->flat();
663  G4double vx1 = sint*cos(phi);
664  G4double vy1 = sint*sin(phi);
665 
666  // lateral displacement
667  if (latDisplasment) {
668  G4double rms = invsqrt12*sqrt(2*z0);
669  G4double r = x0*mscfac;
670  G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0));
671  G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0));
672  G4double d = r*r - dx*dx - dy*dy;
673 
674  // change position
675  if(d >= 0.0) {
676  temp.set(dx,dy,sqrt(d) - r);
677  temp.rotateUz(dir);
678  fDisplacement += temp;
679  }
680  }
681  // change direction
682  temp.set(vx1,vy1,cost);
683  temp.rotateUz(dir);
684  dir = temp;
685  }
686  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
687  } while (0 < nMscSteps);
688 
689  dir.rotateUz(oldDirection);
690 
691  //G4cout<<"G4WentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl;
692  // end of sampling -------------------------------
693 
695 
696  // lateral displacement
697  fDisplacement.rotateUz(oldDirection);
698 
699  /*
700  G4cout << " r(mm)= " << fDisplacement.mag()
701  << " safety= " << safety
702  << " trueStep(mm)= " << tPathLength
703  << " geomStep(mm)= " << zPathLength
704  << " x= " << fDisplacement.x()
705  << " y= " << fDisplacement.y()
706  << " z= " << fDisplacement.z()
707  << G4endl;
708  */
709 
710  //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
711  return fDisplacement;
712 }
713 
714 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
715 
717 {
718  // prepare recomputation of x-sections
719  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
720  const G4double* theAtomNumDensityVector =
723  if(nelm > nelments) {
724  nelments = nelm;
725  xsecn.resize(nelm);
726  prob.resize(nelm);
727  }
728 
729  // check consistency
730  xtsec = 0.0;
731  if(cosTetMaxNuc >= cosTheta) { return 0.0; }
732 
733  G4double cut = (*currentCuts)[currentMaterialIndex];
734  if(fixedCut > 0.0) { cut = fixedCut; }
735 
736  // loop over elements
737  G4double xs = 0.0;
738  for (G4int i=0; i<nelm; ++i) {
739  G4double costm =
740  wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
741  G4double density = theAtomNumDensityVector[i];
742 
743  G4double esec = 0.0;
744  if(costm < cosTheta) {
745 
746  // recompute the transport x-section
747  if(1.0 > cosTheta) {
748  xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosTheta);
749  }
750  // recompute the total x-section
751  G4double nucsec = wokvi->ComputeNuclearCrossSection(cosTheta, costm);
752  esec = wokvi->ComputeElectronCrossSection(cosTheta, costm);
753  nucsec += esec;
754  if(nucsec > 0.0) { esec /= nucsec; }
755  xtsec += nucsec*density;
756  }
757  xsecn[i] = xtsec;
758  prob[i] = esec;
759  //G4cout << i << " xs= " << xs << " xtsec= " << xtsec
760  // << " 1-cosTheta= " << 1-cosTheta
761  // << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
762  }
763 
764  //G4cout << "ComputeXS result: xsec(1/mm)= " << xs
765  // << " txsec(1/mm)= " << xtsec <<G4endl;
766  return xs;
767 }
768 
769 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
770 
772  G4double kinEnergy)
773 {
774  G4double xs = 0.0;
775 
776  SetupParticle(p);
778 
779  if(cosTetMaxNuc >= 1.0) { return xs; }
780 
781  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
782  const G4double* theAtomNumDensityVector =
785 
786  G4double cut = (*currentCuts)[currentMaterialIndex];
787  if(fixedCut > 0.0) { cut = fixedCut; }
788 
789  // loop over elements
790  for (G4int i=0; i<nelm; ++i) {
791  G4double costm =
792  wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
793  xs += theAtomNumDensityVector[i]
795  }
796  return xs;
797 }
798 
799 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
800 
802 {
803  if(val > 0.05) {
804  ssFactor = val;
805  invssFactor = 1.0/(val - 0.05);
806  }
807 }
808 
809 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......