Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4WentzelVIModel.cc 105734 2017-08-16 12:58:28Z 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 
75  : G4VMscModel(nam),
76  ssFactor(1.05),
77  invssFactor(1.0),
78  currentCouple(nullptr),
79  cosThetaMin(1.0),
80  cosThetaMax(-1.0),
81  fSecondMoments(nullptr),
82  idx2(0),
83  numlimit(0.1),
84  singleScatteringMode(false),
85  isCombined(comb),
86  useSecondMoment(false)
87 {
89  invsqrt12 = 1./sqrt(12.);
90  tlimitminfix = 1.e-6*mm;
91  lowEnergyLimit = 1.0*eV;
92  particle = nullptr;
93  nelments = 5;
94  xsecn.resize(nelments);
95  prob.resize(nelments);
97  fixedCut = -1.0;
98 
99  minNCollisions = 10;
100 
102  = currentRange = xtsec = cosTetMaxNuc = 0.0;
104 
105  fParticleChange = nullptr;
106  currentCuts = nullptr;
107  currentMaterial = nullptr;
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111 
113 {
114  delete wokvi;
115  if(fSecondMoments && IsMaster()) {
116  delete fSecondMoments;
117  fSecondMoments = nullptr;
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) {
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
288 
289  //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
290  // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
291 
292  // extra check for abnormal situation
293  // this check needed to run MSC with eIoni and eBrem inactivated
294  if(tlimit > currentRange) { tlimit = currentRange; }
295 
296  // stop here if small range particle
297  if(tlimit < tlimitminfix) {
298  return ConvertTrueToGeom(tlimit, currentMinimalStep);
299  }
300 
301  // pre step
302  G4double presafety = sp->GetSafety();
303  // far from geometry boundary
304  if(currentRange < presafety) {
305  return ConvertTrueToGeom(tlimit, currentMinimalStep);
306  }
307 
308  // compute presafety again if presafety <= 0 and no boundary
309  // i.e. when it is needed for optimization purposes
310  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
311  presafety = ComputeSafety(sp->GetPosition(), tlimit);
312  if(currentRange < presafety) {
313  return ConvertTrueToGeom(tlimit, currentMinimalStep);
314  }
315  }
316  /*
317  G4cout << "e(MeV)= " << preKinEnergy/MeV
318  << " " << particle->GetParticleName()
319  << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
320  << " R(mm)= " <<currentRange/mm
321  << " L0(mm^-1)= " << lambdaeff*mm
322  << G4endl;
323  */
324  // natural limit for high energy
326  (1.0 - cosTetMaxNuc)*lambdaeff*invssFactor);
327 
328  // low-energy e-
329  if(cosThetaMax > cosTetMaxNuc) {
330  rlimit = std::min(rlimit, facsafety*presafety);
331  }
332 
333  // cut correction
335  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= "
336  // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax
337  //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
338  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
339 
340  tlimit = std::min(tlimit, rlimit);
341  tlimit = std::max(tlimit, tlimitminfix);
342 
343  // step limit in infinite media
344  tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
345 
346  //compute geomlimit and force few steps within a volume
348  && stepStatus == fGeomBoundary) {
349 
350  G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
351  tlimit = std::min(tlimit, geomlimit/facgeom);
352  }
353  /*
354  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
355  << " L0= " << lambdaeff << " R= " << currentRange
356  << " tlimit= " << tlimit
357  << " currentMinimalStep= " << currentMinimalStep << G4endl;
358  */
359  return ConvertTrueToGeom(tlimit, currentMinimalStep);
360 }
361 
362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
363 
365 {
366  zPathLength = tPathLength = truelength;
367 
368  // small step use only single scattering
369  cosThetaMin = 1.0;
371  //G4cout << "xtsec= " << xtsec << " Nav= "
372  // << zPathLength*xtsec << G4endl;
373  if(0.0 >= lambdaeff || G4int(zPathLength*xtsec) < minNCollisions) {
374  singleScatteringMode = true;
375  lambdaeff = DBL_MAX;
376 
377  } else {
378  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
379  // << " Leff= " << lambdaeff << G4endl;
380  // small step
383  zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
384 
385  // medium step
386  } else {
387  G4double e1 = 0.0;
388  if(currentRange > tPathLength) {
390  }
391  effKinEnergy = 0.5*(e1 + preKinEnergy);
394  //G4cout << " tLength= "<< tPathLength<< " Leff= " << lambdaeff << G4endl;
396  if(tPathLength*numlimit < lambdaeff) {
397  zPathLength *= (1.0 - G4Exp(-tPathLength/lambdaeff));
398  }
399  }
400  }
401  //G4cout << "Comp.geom: zLength= "<<zPathLength<<" tLength= "
402  // << tPathLength<< " Leff= " << lambdaeff << G4endl;
403  return zPathLength;
404 }
405 
406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
407 
409 {
410  // initialisation of single scattering x-section
411  /*
412  G4cout << "ComputeTrueStepLength: Step= " << geomStepLength
413  << " geomL= " << zPathLength
414  << " Lambda= " << lambdaeff
415  << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
416  */
418  zPathLength = tPathLength = geomStepLength;
419 
420  } else {
421 
422  // step defined by transportation
423  // change both geom and true step lengths
424  if(geomStepLength < zPathLength) {
425 
426  // single scattering
427  if(G4int(geomStepLength*xtsec) < minNCollisions) {
428  zPathLength = tPathLength = geomStepLength;
429  lambdaeff = DBL_MAX;
430  singleScatteringMode = true;
431 
432  // multiple scattering
433  } else {
434  // small step
435  if(geomStepLength < numlimit*lambdaeff) {
436  G4double tau = geomStepLength/lambdaeff;
437  tPathLength = geomStepLength*(1.0 + 0.5*tau + tau*tau/3.0);
438 
439  // energy correction for a big step
440  } else {
441  tPathLength *= geomStepLength/zPathLength;
442  G4double e1 = 0.0;
443  if(currentRange > tPathLength) {
445  }
446  effKinEnergy = 0.5*(e1 + preKinEnergy);
449  G4double tau = geomStepLength/lambdaeff;
450 
451  if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
452  else { tPathLength = currentRange; }
453  }
454  zPathLength = geomStepLength;
455  }
456  }
457  }
458  // check of step length
459  // define threshold angle between single and multiple scattering
460  if(!singleScatteringMode) {
462  xtsec = 0.0;
463 
464  // recompute transport cross section - do not change energy
465  // anymore - cannot be applied for big steps
466  if(cosThetaMin > cosTetMaxNuc) {
467  // new computation
469  //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec
470  // << " 1-cosTMin= " << 1.0 - cosThetaMin << G4endl;
471  if(cross <= 0.0) {
472  singleScatteringMode = true;
474  lambdaeff = DBL_MAX;
475  cosThetaMin = 1.0;
476  } else if(xtsec > 0.0) {
477 
478  lambdaeff = 1./cross;
479  G4double tau = zPathLength*cross;
480  if(tau < numlimit) {
481  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
482  } else if(tau < 0.999999) {
483  tPathLength = -lambdaeff*G4Log(1.0 - tau);
484  } else {
486  }
487  }
488  }
489  }
491  /*
492  G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
493  <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
494  G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
495  << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
496  << " e(MeV)= " << preKinEnergy/MeV << " "
497  << " SSmode= " << singleScatteringMode << G4endl;
498  */
499  return tPathLength;
500 }
501 
502 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
503 
506  G4double /*safety*/)
507 {
508  fDisplacement.set(0.0,0.0,0.0);
509  //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
510  // << particle->GetParticleName() << G4endl;
511 
512  // ignore scattering for zero step length and energy below the limit
513  if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0)
514  { return fDisplacement; }
515 
516  G4double invlambda = 0.0;
517  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
518 
519  // use average kinetic energy over the step
520  G4double cut = (*currentCuts)[currentMaterialIndex];
521  if(fixedCut > 0.0) { cut = fixedCut; }
522  /*
523  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
524  << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
525  << " xmsc= " << tPathLength*invlambda
526  << " safety= " << safety << G4endl;
527  */
528  // step limit due msc
529  G4int nMscSteps = 1;
530  G4double x0 = tPathLength;
531  G4double z0 = x0*invlambda;
532  //G4double zzz = 0.0;
533  G4double prob2 = 0.0;
534 
535  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
536 
537  // large scattering angle case - two step approach
538  if(!singleScatteringMode) {
539  static const G4double zzmin = 0.05;
540  if(useSecondMoment) {
541  G4double z1 = invlambda*invlambda;
543  prob2 = (z2 - z1)/(1.5*z1 - z2);
544  }
545  // if(z0 > zzmin && safety > tlimitminfix) {
546  if(z0 > zzmin) {
547  x0 *= 0.5;
548  z0 *= 0.5;
549  nMscSteps = 2;
550  }
551  //if(z0 > zzmin) { zzz = G4Exp(-1.0/z0); }
552  G4double zzz = 0.0;
553  if(z0 > zzmin) {
554  zzz = G4Exp(-1.0/z0);
555  z0 += zzz;
556  prob2 *= (1 + zzz);
557  }
558  prob2 /= (1 + prob2);
559  }
560 
561  // step limit due to single scattering
562  G4double x1 = 2*tPathLength;
563  if(0.0 < xtsec) { x1 = -G4Log(rndmEngine->flat())/xtsec; }
564 
565  // no scattering case
566  if(singleScatteringMode && x1 > tPathLength)
567  { return fDisplacement; }
568 
569  const G4ElementVector* theElementVector =
572 
573  // geometry
574  G4double sint, cost, phi;
575  G4ThreeVector temp(0.0,0.0,1.0);
576 
577  // current position and direction relative to the end point
578  // because of magnetic field geometry is computed relatively to the
579  // end point of the step
580  G4ThreeVector dir(0.0,0.0,1.0);
581  fDisplacement.set(0.0,0.0,-zPathLength);
582 
584 
585  // start a loop
586  G4double x2 = x0;
587  G4double step, z;
588  G4bool singleScat;
589  /*
590  G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
591  << " 1-cost1= " << 1 - cosThetaMin << " SSmode= " << singleScatteringMode
592  << " xtsec= " << xtsec << " Nst= " << nMscSteps << G4endl;
593  */
594  do {
595 
596  //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl;
597  // single scattering case
598  if(singleScatteringMode && x1 > x2) {
599  fDisplacement += x2*mscfac*dir;
600  break;
601  }
602 
603  // what is next single of multiple?
604  if(x1 <= x2) {
605  step = x1;
606  singleScat = true;
607  } else {
608  step = x2;
609  singleScat = false;
610  }
611 
612  //G4cout << "# step(mm)= "<< step<< " singlScat= "<< singleScat << G4endl;
613 
614  // new position
615  fDisplacement += step*mscfac*dir;
616 
617  if(singleScat) {
618 
619  // select element
620  G4int i = 0;
621  if(nelm > 1) {
622  G4double qsec = rndmEngine->flat()*xtsec;
623  for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
624  }
625  G4double cosTetM =
626  wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
627  //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " "
628  // << prob[i] << G4endl;
629  temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
630 
631  // direction is changed
632  temp.rotateUz(dir);
633  dir = temp;
634  //G4cout << dir << G4endl;
635 
636  // new proposed step length
637  x2 -= step;
638  x1 = -G4Log(rndmEngine->flat())/xtsec;
639 
640  // multiple scattering
641  } else {
642  --nMscSteps;
643  x1 -= step;
644  x2 = x0;
645 
646  // sample z in interval 0 - 1
647  G4bool isFirst = true;
648  if(prob2 > 0.0 && rndmEngine->flat() < prob2) { isFirst = false; }
649  do {
650  //z = -z0*G4Log(1.0 - (1.0 - zzz)*rndmEngine->flat());
651  if(isFirst) { z = -G4Log(rndmEngine->flat()); }
652  else { z = G4RandGamma::shoot(rndmEngine, 2.0, 2.0); }
653  z *= z0;
654  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
655  } while(z > 1.0);
656 
657  cost = 1.0 - 2.0*z/*factCM*/;
658  if(cost > 1.0) { cost = 1.0; }
659  else if(cost < -1.0) { cost =-1.0; }
660  sint = sqrt((1.0 - cost)*(1.0 + cost));
661  phi = twopi*rndmEngine->flat();
662  G4double vx1 = sint*cos(phi);
663  G4double vy1 = sint*sin(phi);
664 
665  // lateral displacement
666  if (latDisplasment) {
667  G4double rms = invsqrt12*sqrt(2*z0);
668  G4double r = x0*mscfac;
669  G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0));
670  G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0));
671  G4double d = r*r - dx*dx - dy*dy;
672 
673  // change position
674  if(d >= 0.0) {
675  temp.set(dx,dy,sqrt(d) - r);
676  temp.rotateUz(dir);
677  fDisplacement += temp;
678  }
679  }
680  // change direction
681  temp.set(vx1,vy1,cost);
682  temp.rotateUz(dir);
683  dir = temp;
684  }
685  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
686  } while (0 < nMscSteps);
687 
688  dir.rotateUz(oldDirection);
689 
690  //G4cout<<"G4WentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl;
691  // end of sampling -------------------------------
692 
694 
695  // lateral displacement
696  fDisplacement.rotateUz(oldDirection);
697 
698  /*
699  G4cout << " r(mm)= " << fDisplacement.mag()
700  << " safety= " << safety
701  << " trueStep(mm)= " << tPathLength
702  << " geomStep(mm)= " << zPathLength
703  << " x= " << fDisplacement.x()
704  << " y= " << fDisplacement.y()
705  << " z= " << fDisplacement.z()
706  << G4endl;
707  */
708 
709  //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
710  return fDisplacement;
711 }
712 
713 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
714 
716 {
717  // prepare recomputation of x-sections
718  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
719  const G4double* theAtomNumDensityVector =
722  if(nelm > nelments) {
723  nelments = nelm;
724  xsecn.resize(nelm);
725  prob.resize(nelm);
726  }
727 
728  // check consistency
729  xtsec = 0.0;
730  if(cosTetMaxNuc >= cosTheta) { return 0.0; }
731 
732  G4double cut = (*currentCuts)[currentMaterialIndex];
733  if(fixedCut > 0.0) { cut = fixedCut; }
734 
735  // loop over elements
736  G4double xs = 0.0;
737  for (G4int i=0; i<nelm; ++i) {
738  G4double costm =
739  wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
740  G4double density = theAtomNumDensityVector[i];
741 
742  G4double esec = 0.0;
743  if(costm < cosTheta) {
744 
745  // recompute the transport x-section
746  if(1.0 > cosTheta) {
747  xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosTheta);
748  }
749  // recompute the total x-section
750  G4double nucsec = wokvi->ComputeNuclearCrossSection(cosTheta, costm);
751  esec = wokvi->ComputeElectronCrossSection(cosTheta, costm);
752  nucsec += esec;
753  if(nucsec > 0.0) { esec /= nucsec; }
754  xtsec += nucsec*density;
755  }
756  xsecn[i] = xtsec;
757  prob[i] = esec;
758  //G4cout << i << " xs= " << xs << " xtsec= " << xtsec
759  // << " 1-cosTheta= " << 1-cosTheta
760  // << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
761  }
762 
763  //G4cout << "ComputeXS result: xsec(1/mm)= " << xs
764  // << " txsec(1/mm)= " << xtsec <<G4endl;
765  return xs;
766 }
767 
768 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
769 
771  G4double kinEnergy)
772 {
773  G4double xs = 0.0;
774 
775  SetupParticle(p);
777 
778  if(cosTetMaxNuc >= 1.0) { return xs; }
779 
780  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
781  const G4double* theAtomNumDensityVector =
784 
785  G4double cut = (*currentCuts)[currentMaterialIndex];
786  if(fixedCut > 0.0) { cut = fixedCut; }
787 
788  // loop over elements
789  for (G4int i=0; i<nelm; ++i) {
790  G4double costm =
791  wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
792  xs += theAtomNumDensityVector[i]
794  }
795  return xs;
796 }
797 
798 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
799 
801 {
802  if(val > 0.05) {
803  ssFactor = val;
804  invssFactor = 1.0/(val - 0.05);
805  }
806 }
807 
808 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double Energy(size_t index) const
void set(double x, double y, double z)
G4double ComputeSecondTransportMoment(G4double CosThetaMax)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=DBL_MAX, G4double emax=DBL_MAX) override
G4ThreeVector & SampleSingleScattering(G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio)
G4double facsafety
Definition: G4VMscModel.hh:177
G4int NumberOfBinsPerDecade() const
void DefineMaterial(const G4MaterialCutsCouple *)
virtual ~G4WentzelVIModel()
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:184
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
G4WentzelVIModel(G4bool comb=true, const G4String &nam="WentzelVIUni")
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:245
G4StepPoint * GetPreStepPoint() const
static constexpr double mm
Definition: G4SIunits.hh:115
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:255
G4PhysicsTable * GetSecondMomentTable()
G4double ComputeSafety(const G4ThreeVector &position, G4double limit=DBL_MAX)
Definition: G4VMscModel.hh:237
Float_t x1[n_points_granero]
Definition: compare.C:5
virtual G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep) override
virtual G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety) override
G4double GetRadlen() const
Definition: G4Material.hh:221
const char * p
Definition: xmltok.h:285
G4ParticleChangeForMSC * fParticleChange
Double_t z
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
G4double facgeom
Definition: G4VMscModel.hh:176
std::vector< G4double > prob
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
void SetSingleScatteringFactor(G4double)
static const G4double emax
const G4ParticleDefinition * GetParticleDefinition() const
G4double SecondMoment(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double kineticEnergy)
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:455
G4double LowEnergyActivationLimit() const
Definition: G4VEmModel.hh:630
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double GetProductionCut(G4int index) const
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:616
Float_t Z
G4double facrange
Definition: G4VMscModel.hh:175
G4double ComputeTransportXSectionPerVolume(G4double cosTheta)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4StepStatus GetStepStatus() const
G4PhysicsTable * GetCrossSectionTable()
Definition: G4VEmModel.hh:808
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
virtual double flat()=0
G4double PolarAngleLimit() const
Definition: G4VEmModel.hh:637
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) override
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:283
G4bool latDisplasment
Definition: G4VMscModel.hh:188
const G4ThreeVector & GetPosition() const
virtual G4double ComputeTrueStepLength(G4double geomStepLength) override
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double ComputeNuclearCrossSection(G4double CosThetaMin, G4double CosThetaMax)
virtual G4double ComputeGeomPathLength(G4double truePathLength) override
void FillSecondDerivatives()
Float_t d
G4WentzelOKandVIxSection * wokvi
G4StepStatus
Definition: G4StepStatus.hh:49
static constexpr double eV
Definition: G4SIunits.hh:215
const G4Step * GetStep() const
const G4MaterialCutsCouple * currentCouple
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double ComputeSecondMoment(const G4ParticleDefinition *, G4double kineticEnergy)
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:462
G4double SetupTarget(G4int Z, G4double cut)
static G4ProductionCutsTable * GetProductionCutsTable()
const G4ParticleDefinition * particle
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:340
void Initialise(const G4ParticleDefinition *, G4double CosThetaLim)
int G4lrint(double ad)
Definition: templates.hh:151
std::vector< G4Element * > G4ElementVector
int G4int
Definition: G4Types.hh:78
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=nullptr)
Definition: G4VMscModel.cc:91
virtual void StartTracking(G4Track *) override
TDirectory * dir
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:185
G4double GetSafety() const
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4double GetKineticEnergy() const
void SetupParticle(const G4ParticleDefinition *)
G4bool IsMaster() const
Definition: G4VEmModel.hh:700
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4ProductionCuts * GetProductionCuts() const
G4double HighEnergyActivationLimit() const
Definition: G4VEmModel.hh:623
G4PhysicsTable * fSecondMoments
Char_t n[5]
const G4Material * GetMaterial() const
const G4DataVector * currentCuts
G4bool GetFlag(size_t i) const
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:303
const G4Material * currentMaterial
Float_t x2[n_points_geant4]
Definition: compare.C:26
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
#define DBL_MAX
Definition: templates.hh:83
const G4DynamicParticle * GetDynamicParticle() const
static G4EmParameters * Instance()
void PutValue(size_t index, G4double theValue)
static constexpr double pi
Definition: SystemOfUnits.h:54
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double ComputeElectronCrossSection(G4double CosThetaMin, G4double CosThetaMax)
static G4double tet[DIM]
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
std::vector< G4double > xsecn