Geant4  10.02.p01
G4WentzelVIModel.cc
Go to the documentation of this file.
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 // $Id: G4WentzelVIModel.cc 91726 2015-08-03 15:41:36Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4WentzelVIModel
34 //
35 // Author: V.Ivanchenko
36 //
37 // Creation date: 09.04.2008 from G4MuMscModel
38 //
39 // Modifications:
40 // 27-05-2010 V.Ivanchenko added G4WentzelOKandVIxSection class to
41 // compute cross sections and sample scattering angle
42 //
43 //
44 // Class Description:
45 //
46 // Implementation of the model of multiple scattering based on
47 // G.Wentzel, Z. Phys. 40 (1927) 590.
48 // H.W.Lewis, Phys Rev 78 (1950) 526.
49 // J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
50 // L.Urban, CERN-OPEN-2006-077.
51 
52 // -------------------------------------------------------------------
53 //
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 
58 #include "G4WentzelVIModel.hh"
59 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "Randomize.hh"
63 #include "G4PhysicsTableHelper.hh"
64 #include "G4ElementVector.hh"
65 #include "G4ProductionCutsTable.hh"
66 #include "G4EmParameters.hh"
67 #include "G4Log.hh"
68 #include "G4Exp.hh"
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 using namespace std;
73 
74 const G4int minNCollisions = 10;
75 
77  G4VMscModel(nam),
78  ssFactor(1.05),
79  invssFactor(1.0),
80  currentCouple(0),
81  singleScatteringMode(false),
82  cosThetaMin(1.0),
83  cosThetaMax(-1.0),
84  fSecondMoments(0),
85  idx2(0),
86  numlimit(0.1),
87  isCombined(combined),
88  useSecondMoment(false)
89 {
91  invsqrt12 = 1./sqrt(12.);
92  tlimitminfix = 1.e-6*mm;
93  lowEnergyLimit = 1.0*eV;
94  particle = 0;
95  nelments = 5;
96  xsecn.resize(nelments);
97  prob.resize(nelments);
98  wokvi = new G4WentzelOKandVIxSection(combined);
99  fixedCut = -1.0;
100 
102  = currentRange = xtsec = cosTetMaxNuc = 0.0;
104 
105  fParticleChange = 0;
106  currentCuts = 0;
107  currentMaterial = 0;
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111 
113 {
114  delete wokvi;
115  if(fSecondMoments && IsMaster()) {
116  delete fSecondMoments;
117  fSecondMoments = 0;
118  }
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
124  const G4DataVector& cuts)
125 {
126  // reset parameters
127  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 
136  //G4cout << "G4WentzelVIModel::Initialise " << p->GetParticleName() << G4endl;
138  /*
139  G4cout << "G4WentzelVIModel: " << particle->GetParticleName()
140  << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax
141  << " SingScatFactor= " << ssFactor
142  << G4endl;
143  */
144  currentCuts = &cuts;
145 
146  // set values of some data members
148 
149  // build second moment table only if transport table is build
151  if(useSecondMoment && IsMaster() && table) {
152 
153  //G4cout << "### G4WentzelVIModel::Initialise: build 2nd moment table "
154  // << table << G4endl;
155  fSecondMoments =
157  // Access to materials
158  const G4ProductionCutsTable* theCoupleTable =
160  size_t numOfCouples = theCoupleTable->GetTableSize();
161 
162  G4bool splineFlag = true;
163  G4PhysicsVector* aVector = 0;
164  G4PhysicsVector* bVector = 0;
167  if(emin < emax) {
169  *G4lrint(std::log10(emax/emin));
170  if(n < 3) { n = 3; }
171 
172  for(size_t i=0; i<numOfCouples; ++i) {
173 
174  //G4cout<< "i= " << i << " Flag= " << fSecondMoments->GetFlag(i)
175  // << G4endl;
176  if(fSecondMoments->GetFlag(i)) {
177  DefineMaterial(theCoupleTable->GetMaterialCutsCouple(i));
178 
179  delete (*fSecondMoments)[i];
180  if(!aVector) {
181  aVector = new G4PhysicsLogVector(emin, emax, n);
182  bVector = aVector;
183  } else {
184  bVector = new G4PhysicsVector(*aVector);
185  }
186  for(size_t j=0; j<n; ++j) {
187  G4double e = bVector->Energy(j);
188  bVector->PutValue(j, ComputeSecondMoment(p, e)*e*e);
189  }
190  if(splineFlag) { bVector->FillSecondDerivatives(); }
191  (*fSecondMoments)[i] = bVector;
192  }
193  }
194  }
195  //G4cout << *fSecondMoments << G4endl;
196  }
197 }
198 
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
200 
202  G4VEmModel* masterModel)
203 {
204  fSecondMoments = static_cast<G4WentzelVIModel*>(masterModel)
206 }
207 
208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209 
211  const G4ParticleDefinition* p,
212  G4double kinEnergy,
213  G4double Z, G4double,
214  G4double cutEnergy, G4double)
215 {
216  G4double cross = 0.0;
217  if(p != particle) { SetupParticle(p); }
218  if(kinEnergy < lowEnergyLimit) { return cross; }
219  if(!CurrentCouple()) {
220  G4Exception("G4WentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
221  FatalException, " G4MaterialCutsCouple is not defined");
222  return 0.0;
223  }
226  if(cosTetMaxNuc < 1.0) {
227  G4double cut = cutEnergy;
228  if(fixedCut > 0.0) { cut = fixedCut; }
229  G4double cost = wokvi->SetupTarget(G4lrint(Z), cut);
231  /*
232  if(p->GetParticleName() == "e-")
233  G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= "<<kinEnergy
234  << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
235  << " " << particle->GetParticleName() << G4endl;
236  */
237  }
238  return cross;
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 
244 {
246 }
247 
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249 
251  const G4Track& track,
252  G4double& currentMinimalStep)
253 {
254  G4double tlimit = currentMinimalStep;
255  const G4DynamicParticle* dp = track.GetDynamicParticle();
256  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
257  G4StepStatus stepStatus = sp->GetStepStatus();
258  singleScatteringMode = false;
259 
260  //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
261  // << stepStatus << " " << track.GetDefinition()->GetParticleName()
262  // << G4endl;
263 
264  // initialisation for each step, lambda may be computed from scratch
271 
272  //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
273  // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
274 
275  // extra check for abnormal situation
276  // this check needed to run MSC with eIoni and eBrem inactivated
277  if(tlimit > currentRange) { tlimit = currentRange; }
278 
279  // stop here if small range particle
280  if(tlimit < tlimitminfix) {
281  return ConvertTrueToGeom(tlimit, currentMinimalStep);
282  }
283 
284  // pre step
285  G4double presafety = sp->GetSafety();
286  // far from geometry boundary
287  if(currentRange < presafety) {
288  return ConvertTrueToGeom(tlimit, currentMinimalStep);
289  }
290 
291  // compute presafety again if presafety <= 0 and no boundary
292  // i.e. when it is needed for optimization purposes
293  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
294  presafety = ComputeSafety(sp->GetPosition(), tlimit);
295  if(currentRange < presafety) {
296  return ConvertTrueToGeom(tlimit, currentMinimalStep);
297  }
298  }
299  /*
300  G4cout << "e(MeV)= " << preKinEnergy/MeV
301  << " " << particle->GetParticleName()
302  << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
303  << " R(mm)= " <<currentRange/mm
304  << " L0(mm^-1)= " << lambdaeff*mm
305  << G4endl;
306  */
307  // natural limit for high energy
309  (1.0 - cosTetMaxNuc)*lambdaeff*invssFactor);
310 
311  // low-energy e-
312  if(cosThetaMax > cosTetMaxNuc) {
313  rlimit = std::min(rlimit, facsafety*presafety);
314  }
315 
316  // cut correction
318  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= "
319  // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax
320  //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
321  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
322 
323  tlimit = std::min(tlimit, rlimit);
324  tlimit = std::max(tlimit, tlimitminfix);
325 
326  // step limit in infinite media
327  tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
328 
329  //compute geomlimit and force few steps within a volume
331  && stepStatus == fGeomBoundary) {
332 
333  G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
334  tlimit = std::min(tlimit, geomlimit/facgeom);
335  }
336  /*
337  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
338  << " L0= " << lambdaeff << " R= " << currentRange
339  << " tlimit= " << tlimit
340  << " currentMinimalStep= " << currentMinimalStep << G4endl;
341  */
342  return ConvertTrueToGeom(tlimit, currentMinimalStep);
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
346 
348 {
349  zPathLength = tPathLength = truelength;
350 
351  // small step use only single scattering
352  cosThetaMin = 1.0;
354  //G4cout << "xtsec= " << xtsec << " Nav= "
355  // << zPathLength*xtsec << G4endl;
356  if(0.0 >= lambdaeff || G4int(zPathLength*xtsec) < minNCollisions) {
357  singleScatteringMode = true;
358  lambdaeff = DBL_MAX;
359 
360  } else {
361  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
362  // << " Leff= " << lambdaeff << G4endl;
363  // small step
366  zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
367 
368  // medium step
369  } else {
370  G4double e1 = 0.0;
371  if(currentRange > tPathLength) {
373  }
374  effKinEnergy = 0.5*(e1 + preKinEnergy);
377  //G4cout << " tLength= "<< tPathLength<< " Leff= " << lambdaeff << G4endl;
379  if(tPathLength*numlimit < lambdaeff) {
380  zPathLength *= (1.0 - G4Exp(-tPathLength/lambdaeff));
381  }
382  }
383  }
384  //G4cout << "Comp.geom: zLength= "<<zPathLength<<" tLength= "
385  // << tPathLength<< " Leff= " << lambdaeff << G4endl;
386  return zPathLength;
387 }
388 
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390 
392 {
393  // initialisation of single scattering x-section
394  /*
395  G4cout << "ComputeTrueStepLength: Step= " << geomStepLength
396  << " geomL= " << zPathLength
397  << " Lambda= " << lambdaeff
398  << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
399  */
401  zPathLength = tPathLength = geomStepLength;
402 
403  } else {
404 
405  // step defined by transportation
406  // change both geom and true step lengths
407  if(geomStepLength < zPathLength) {
408 
409  // single scattering
410  if(G4int(geomStepLength*xtsec) < minNCollisions) {
411  zPathLength = tPathLength = geomStepLength;
412  lambdaeff = DBL_MAX;
413  singleScatteringMode = true;
414 
415  // multiple scattering
416  } else {
417  // small step
418  if(geomStepLength < numlimit*lambdaeff) {
419  G4double tau = geomStepLength/lambdaeff;
420  tPathLength = geomStepLength*(1.0 + 0.5*tau + tau*tau/3.0);
421 
422  // energy correction for a big step
423  } else {
424  tPathLength *= geomStepLength/zPathLength;
425  G4double e1 = 0.0;
426  if(currentRange > tPathLength) {
428  }
429  effKinEnergy = 0.5*(e1 + preKinEnergy);
432  G4double tau = geomStepLength/lambdaeff;
433 
434  if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
435  else { tPathLength = currentRange; }
436  }
437  zPathLength = geomStepLength;
438  }
439  }
440  }
441  // check of step length
442  // define threshold angle between single and multiple scattering
443  if(!singleScatteringMode) {
445  xtsec = 0.0;
446 
447  // recompute transport cross section - do not change energy
448  // anymore - cannot be applied for big steps
449  if(cosThetaMin > cosTetMaxNuc) {
450  // new computation
452  //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec
453  // << " 1-cosTMin= " << 1.0 - cosThetaMin << G4endl;
454  if(cross <= 0.0) {
455  singleScatteringMode = true;
457  lambdaeff = DBL_MAX;
458  cosThetaMin = 1.0;
459  } else if(xtsec > 0.0) {
460 
461  lambdaeff = 1./cross;
462  G4double tau = zPathLength*cross;
463  if(tau < numlimit) {
464  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
465  } else if(tau < 0.999999) {
466  tPathLength = -lambdaeff*G4Log(1.0 - tau);
467  } else {
469  }
470  }
471  }
472  }
474  /*
475  G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
476  <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
477  G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
478  << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
479  << " e(MeV)= " << preKinEnergy/MeV << " "
480  << " SSmode= " << singleScatteringMode << G4endl;
481  */
482  return tPathLength;
483 }
484 
485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
486 
489  G4double /*safety*/)
490 {
491  fDisplacement.set(0.0,0.0,0.0);
492  //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
493  // << particle->GetParticleName() << G4endl;
494 
495  // ignore scattering for zero step length and energy below the limit
496  if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0)
497  { return fDisplacement; }
498 
499  G4double invlambda = 0.0;
500  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
501 
502  // use average kinetic energy over the step
503  G4double cut = (*currentCuts)[currentMaterialIndex];
504  if(fixedCut > 0.0) { cut = fixedCut; }
505  /*
506  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
507  << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
508  << " xmsc= " << tPathLength*invlambda
509  << " safety= " << safety << G4endl;
510  */
511  // step limit due msc
512  G4int nMscSteps = 1;
513  G4double x0 = tPathLength;
514  G4double z0 = x0*invlambda;
515  //G4double zzz = 0.0;
516  G4double prob2 = 0.0;
517 
518  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
519 
520  // large scattering angle case - two step approach
521  if(!singleScatteringMode) {
522  static const G4double zzmin = 0.05;
523  if(useSecondMoment) {
524  G4double z1 = invlambda*invlambda;
526  prob2 = (z2 - z1)/(1.5*z1 - z2);
527  }
528  // if(z0 > zzmin && safety > tlimitminfix) {
529  if(z0 > zzmin) {
530  x0 *= 0.5;
531  z0 *= 0.5;
532  nMscSteps = 2;
533  }
534  //if(z0 > zzmin) { zzz = G4Exp(-1.0/z0); }
535  G4double zzz = 0.0;
536  if(z0 > zzmin) {
537  zzz = G4Exp(-1.0/z0);
538  z0 += zzz;
539  prob2 *= (1 + zzz);
540  }
541  prob2 /= (1 + prob2);
542  }
543 
544  // step limit due to single scattering
545  G4double x1 = 2*tPathLength;
546  if(0.0 < xtsec) { x1 = -G4Log(rndmEngine->flat())/xtsec; }
547 
548  // no scattering case
549  if(singleScatteringMode && x1 > tPathLength)
550  { return fDisplacement; }
551 
552  const G4ElementVector* theElementVector =
555 
556  // geometry
557  G4double sint, cost, phi;
558  G4ThreeVector temp(0.0,0.0,1.0);
559 
560  // current position and direction relative to the end point
561  // because of magnetic field geometry is computed relatively to the
562  // end point of the step
563  G4ThreeVector dir(0.0,0.0,1.0);
564  fDisplacement.set(0.0,0.0,-zPathLength);
565 
567 
568  // start a loop
569  G4double x2 = x0;
570  G4double step, z;
571  G4bool singleScat;
572  /*
573  G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
574  << " 1-cost1= " << 1 - cosThetaMin << " SSmode= " << singleScatteringMode
575  << " xtsec= " << xtsec << " Nst= " << nMscSteps << G4endl;
576  */
577  do {
578 
579  //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl;
580  // single scattering case
581  if(singleScatteringMode && x1 > x2) {
582  fDisplacement += x2*mscfac*dir;
583  break;
584  }
585 
586  // what is next single of multiple?
587  if(x1 <= x2) {
588  step = x1;
589  singleScat = true;
590  } else {
591  step = x2;
592  singleScat = false;
593  }
594 
595  //G4cout << "# step(mm)= "<< step<< " singlScat= "<< singleScat << G4endl;
596 
597  // new position
598  fDisplacement += step*mscfac*dir;
599 
600  if(singleScat) {
601 
602  // select element
603  G4int i = 0;
604  if(nelm > 1) {
605  G4double qsec = rndmEngine->flat()*xtsec;
606  for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
607  }
608  G4double cosTetM =
609  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
610  //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " "
611  // << prob[i] << G4endl;
612  temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
613 
614  // direction is changed
615  temp.rotateUz(dir);
616  dir = temp;
617  //G4cout << dir << G4endl;
618 
619  // new proposed step length
620  x2 -= step;
621  x1 = -G4Log(rndmEngine->flat())/xtsec;
622 
623  // multiple scattering
624  } else {
625  --nMscSteps;
626  x1 -= step;
627  x2 = x0;
628 
629  // sample z in interval 0 - 1
630  G4bool isFirst = true;
631  if(prob2 > 0.0 && rndmEngine->flat() < prob2) { isFirst = false; }
632  do {
633  //z = -z0*G4Log(1.0 - (1.0 - zzz)*rndmEngine->flat());
634  if(isFirst) { z = -G4Log(rndmEngine->flat()); }
635  else { z = G4RandGamma::shoot(rndmEngine, 2.0, 2.0); }
636  z *= z0;
637  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
638  } while(z > 1.0);
639 
640  cost = 1.0 - 2.0*z/*factCM*/;
641  if(cost > 1.0) { cost = 1.0; }
642  else if(cost < -1.0) { cost =-1.0; }
643  sint = sqrt((1.0 - cost)*(1.0 + cost));
644  phi = twopi*rndmEngine->flat();
645  G4double vx1 = sint*cos(phi);
646  G4double vy1 = sint*sin(phi);
647 
648  // lateral displacement
649  if (latDisplasment) {
650  G4double rms = invsqrt12*sqrt(2*z0);
651  G4double r = x0*mscfac;
652  G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(rndmEngine, 0.0, 1.0));
653  G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(rndmEngine, 0.0, 1.0));
654  G4double dz;
655  G4double d = r*r - dx*dx - dy*dy;
656 
657  // change position
658  if(d >= 0.0) {
659  dz = sqrt(d) - r;
660  temp.set(dx,dy,dz);
661  temp.rotateUz(dir);
662  fDisplacement += temp;
663  }
664  }
665  // change direction
666  temp.set(vx1,vy1,cost);
667  temp.rotateUz(dir);
668  dir = temp;
669  }
670  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
671  } while (0 < nMscSteps);
672 
673  dir.rotateUz(oldDirection);
674 
675  //G4cout<<"G4WentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl;
676  // end of sampling -------------------------------
677 
679 
680  // lateral displacement
681  fDisplacement.rotateUz(oldDirection);
682 
683  /*
684  G4cout << " r(mm)= " << fDisplacement.mag()
685  << " safety= " << safety
686  << " trueStep(mm)= " << tPathLength
687  << " geomStep(mm)= " << zPathLength
688  << " x= " << fDisplacement.x()
689  << " y= " << fDisplacement.y()
690  << " z= " << fDisplacement.z()
691  << G4endl;
692  */
693 
694  //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
695  return fDisplacement;
696 }
697 
698 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
699 
701 {
702  // prepare recomputation of x-sections
703  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
704  const G4double* theAtomNumDensityVector =
707  if(nelm > nelments) {
708  nelments = nelm;
709  xsecn.resize(nelm);
710  prob.resize(nelm);
711  }
712 
713  // check consistency
714  xtsec = 0.0;
715  if(cosTetMaxNuc >= cosTheta) { return 0.0; }
716 
717  G4double cut = (*currentCuts)[currentMaterialIndex];
718  if(fixedCut > 0.0) { cut = fixedCut; }
719 
720  // loop over elements
721  G4double xs = 0.0;
722  for (G4int i=0; i<nelm; ++i) {
723  G4double costm =
724  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
725  G4double density = theAtomNumDensityVector[i];
726 
727  G4double esec = 0.0;
728  if(costm < cosTheta) {
729 
730  // recompute the transport x-section
731  if(1.0 > cosTheta) {
732  xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosTheta);
733  }
734  // recompute the total x-section
735  G4double nucsec = wokvi->ComputeNuclearCrossSection(cosTheta, costm);
736  esec = wokvi->ComputeElectronCrossSection(cosTheta, costm);
737  nucsec += esec;
738  if(nucsec > 0.0) { esec /= nucsec; }
739  xtsec += nucsec*density;
740  }
741  xsecn[i] = xtsec;
742  prob[i] = esec;
743  //G4cout << i << " xs= " << xs << " xtsec= " << xtsec
744  // << " 1-cosTheta= " << 1-cosTheta
745  // << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
746  }
747 
748  //G4cout << "ComputeXS result: xsec(1/mm)= " << xs
749  // << " txsec(1/mm)= " << xtsec <<G4endl;
750  return xs;
751 }
752 
753 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
754 
756  G4double kinEnergy)
757 {
758  G4double xs = 0.0;
759 
760  SetupParticle(p);
762 
763  if(cosTetMaxNuc >= 1.0) { return xs; }
764 
765  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
766  const G4double* theAtomNumDensityVector =
769 
770  G4double cut = (*currentCuts)[currentMaterialIndex];
771  if(fixedCut > 0.0) { cut = fixedCut; }
772 
773  // loop over elements
774  for (G4int i=0; i<nelm; ++i) {
775  G4double costm =
776  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
777  xs += theAtomNumDensityVector[i]
779  }
780  return xs;
781 }
782 
783 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
784 
786 {
787  if(val > 0.05) {
788  ssFactor = val;
789  invssFactor = 1.0/(val - 0.05);
790  }
791 }
792 
793 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double LowEnergyActivationLimit() const
Definition: G4VEmModel.hh:655
G4double HighEnergyActivationLimit() const
Definition: G4VEmModel.hh:648
virtual ~G4WentzelVIModel()
G4int NumberOfBinsPerDecade() const
void SetSingleScatteringFactor(G4double)
G4double facgeom
Definition: G4VMscModel.hh:178
ThreeVector shoot(const G4int Ap, const G4int Af)
G4WentzelOKandVIxSection * wokvi
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
std::vector< G4double > prob
std::vector< G4Element * > G4ElementVector
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
const G4DynamicParticle * GetDynamicParticle() const
G4double facrange
Definition: G4VMscModel.hh:177
G4double z
Definition: TRTMaterials.hh:39
G4double GetProductionCut(G4int index) const
G4StepStatus GetStepStatus() const
virtual G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety)
G4double ComputeSecondTransportMoment(G4double CosThetaMax)
void DefineMaterial(const G4MaterialCutsCouple *)
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:247
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4bool latDisplasment
Definition: G4VMscModel.hh:190
const G4MaterialCutsCouple * currentCouple
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
G4PhysicsTable * GetCrossSectionTable()
Definition: G4VEmModel.hh:826
G4ParticleChangeForMSC * fParticleChange
G4ParticleDefinition * GetDefinition() const
const G4Material * currentMaterial
G4WentzelVIModel(G4bool comb=true, const G4String &nam="WentzelVIUni")
const G4Step * GetStep() const
G4double ComputeSafety(const G4ThreeVector &position, G4double limit=DBL_MAX)
Definition: G4VMscModel.hh:239
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
void FillSecondDerivatives()
G4double SetupTarget(G4int Z, G4double cut=DBL_MAX)
virtual G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep)
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
G4StepStatus
Definition: G4StepStatus.hh:49
G4StepPoint * GetPreStepPoint() const
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:315
G4double density
Definition: TRTMaterials.hh:39
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
std::vector< G4double > xsecn
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:452
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=DBL_MAX, G4double emax=DBL_MAX)
const G4ThreeVector & GetPosition() const
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:295
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:90
bool G4bool
Definition: G4Types.hh:79
G4double SecondMoment(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double kineticEnergy)
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
static const double twopi
Definition: G4SIunits.hh:75
void PutValue(size_t index, G4double theValue)
G4double ComputeTransportXSectionPerVolume(G4double cosTheta)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4int n
G4double Energy(size_t index) const
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:257
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:186
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:352
G4double GetRadlen() const
Definition: G4Material.hh:220
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4double emax
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:212
G4double ComputeElectronCrossSection(G4double CosThetaMin, G4double CosThetaMax)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static const double pi
Definition: G4SIunits.hh:74
const G4int minNCollisions
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double facsafety
Definition: G4VMscModel.hh:179
void SetupParticle(const G4ParticleDefinition *)
const G4DataVector * currentCuts
void StartTracking(G4Track *)
G4PhysicsTable * GetSecondMomentTable()
G4double ComputeNuclearCrossSection(G4double CosThetaMin, G4double CosThetaMax)
static G4EmParameters * Instance()
G4double PolarAngleLimit() const
Definition: G4VEmModel.hh:662
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
const G4ParticleDefinition * particle
G4double GetSafety() const
G4double ComputeSecondMoment(const G4ParticleDefinition *, G4double kineticEnergy)
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:187
G4PhysicsTable * fSecondMoments
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4bool GetFlag(size_t i) const
double G4double
Definition: G4Types.hh:76
virtual G4double ComputeTrueStepLength(G4double geomStepLength)
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
G4ProductionCuts * GetProductionCuts() const
#define DBL_MAX
Definition: templates.hh:83
static const double mm
Definition: G4SIunits.hh:114
void Initialise(const G4ParticleDefinition *, G4double CosThetaLim)
virtual G4double ComputeGeomPathLength(G4double truePathLength)
G4ThreeVector & SampleSingleScattering(G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio=0.0)