Geant4_10
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 74543 2013-10-14 12:38:02Z 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 "G4LossTableManager.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  numlimit(0.1),
77  currentCouple(0),
78  cosThetaMin(1.0),
79  inside(false),
80  singleScatteringMode(false)
81 {
82  invsqrt12 = 1./sqrt(12.);
83  tlimitminfix = 1.e-6*mm;
84  lowEnergyLimit = 1.0*eV;
85  particle = 0;
86  nelments = 5;
87  xsecn.resize(nelments);
88  prob.resize(nelments);
89  theManager = G4LossTableManager::Instance();
90  wokvi = new G4WentzelOKandVIxSection();
91  fixedCut = -1.0;
92 
93  preKinEnergy = tPathLength = zPathLength = lambdaeff = currentRange
94  = xtsec = 0;
95  currentMaterialIndex = 0;
96  cosThetaMax = cosTetMaxNuc = 1.0;
97 
98  fParticleChange = 0;
99  currentCuts = 0;
100  currentMaterial = 0;
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104 
106 {
107  delete wokvi;
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111 
113  const G4DataVector& cuts)
114 {
115  // reset parameters
116  SetupParticle(p);
117  currentRange = 0.0;
118 
119  cosThetaMax = cos(PolarAngleLimit());
120  //G4cout << "G4WentzelVIModel::Initialise " << p->GetParticleName() << G4endl;
121  wokvi->Initialise(p, cosThetaMax);
122  /*
123  G4cout << "G4WentzelVIModel: " << particle->GetParticleName()
124  << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax
125  << G4endl;
126  */
127  currentCuts = &cuts;
128 
129  // set values of some data members
130  fParticleChange = GetParticleChangeForMSC(p);
131 }
132 
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 
136  const G4ParticleDefinition* p,
137  G4double kinEnergy,
139  G4double cutEnergy, G4double)
140 {
141  G4double cross = 0.0;
142  if(p != particle) { SetupParticle(p); }
143  if(kinEnergy < lowEnergyLimit) { return cross; }
144  if(!CurrentCouple()) {
145  G4Exception("G4WentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
146  FatalException, " G4MaterialCutsCouple is not defined");
147  return 0.0;
148  }
149  DefineMaterial(CurrentCouple());
150  cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
151  if(cosTetMaxNuc < 1.0) {
152  G4double cut = cutEnergy;
153  if(fixedCut > 0.0) { cut = fixedCut; }
154  cosTetMaxNuc = wokvi->SetupTarget(G4lrint(Z), cut);
155  cross = wokvi->ComputeTransportCrossSectionPerAtom(cosTetMaxNuc);
156  /*
157  if(p->GetParticleName() == "e-")
158  G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= " << kinEnergy
159  << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
160  << " " << particle->GetParticleName() << G4endl;
161  */
162  }
163  return cross;
164 }
165 
166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
167 
169 {
170  SetupParticle(track->GetDynamicParticle()->GetDefinition());
171  inside = false;
172 }
173 
174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
175 
177  const G4Track& track,
178  G4double& currentMinimalStep)
179 {
180  G4double tlimit = currentMinimalStep;
181  const G4DynamicParticle* dp = track.GetDynamicParticle();
182  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
183  G4StepStatus stepStatus = sp->GetStepStatus();
184  singleScatteringMode = false;
185  //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
186  // << stepStatus << " " << track.GetDefinition()->GetParticleName()
187  // << G4endl;
188 
189  // initialisation for each step, lambda may be computed from scratch
190  preKinEnergy = dp->GetKineticEnergy();
191  DefineMaterial(track.GetMaterialCutsCouple());
192  lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy);
193  currentRange = GetRange(particle,preKinEnergy,currentCouple);
194  cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial);
195 
196  //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
197  // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
198 
199  // extra check for abnormal situation
200  // this check needed to run MSC with eIoni and eBrem inactivated
201  if(tlimit > currentRange) { tlimit = currentRange; }
202 
203  // stop here if small range particle
204  if(inside || tlimit < tlimitminfix) {
205  return ConvertTrueToGeom(tlimit, currentMinimalStep);
206  }
207 
208  // pre step
209  G4double presafety = sp->GetSafety();
210  // far from geometry boundary
211  if(currentRange < presafety) {
212  inside = true;
213  return ConvertTrueToGeom(tlimit, currentMinimalStep);
214  }
215 
216  // compute presafety again if presafety <= 0 and no boundary
217  // i.e. when it is needed for optimization purposes
218  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
219  presafety = ComputeSafety(sp->GetPosition(), tlimit);
220  if(currentRange < presafety) {
221  inside = true;
222  return ConvertTrueToGeom(tlimit, currentMinimalStep);
223  }
224  }
225  /*
226  G4cout << "e(MeV)= " << preKinEnergy/MeV
227  << " " << particle->GetParticleName()
228  << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
229  << " R(mm)= " <<currentRange/mm
230  << " L0(mm^-1)= " << lambdaeff*mm
231  <<G4endl;
232  */
233 
234  // natural limit for high energy
235  G4double rlimit = std::max(facrange*currentRange,
236  0.7*(1.0 - cosTetMaxNuc)*lambdaeff);
237 
238  // low-energy e-
239  if(cosThetaMax > cosTetMaxNuc) {
240  rlimit = std::min(rlimit, facsafety*presafety);
241  }
242 
243  // cut correction
244  G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
245  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= "
246  // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax
247  //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
248  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
249 
250  if(rlimit < tlimit) { tlimit = rlimit; }
251 
252  tlimit = std::max(tlimit, tlimitminfix);
253 
254  // step limit in infinite media
255  tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
256 
257  //compute geomlimit and force few steps within a volume
259  && stepStatus == fGeomBoundary) {
260 
261  G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
262  tlimit = std::min(tlimit, geomlimit/facgeom);
263  }
264 
265  /*
266  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
267  << " L0= " << lambdaeff << " R= " << currentRange
268  << "tlimit= " << tlimit
269  << " currentMinimalStep= " << currentMinimalStep << G4endl;
270  */
271  return ConvertTrueToGeom(tlimit, currentMinimalStep);
272 }
273 
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
275 
277 {
278  tPathLength = truelength;
279  zPathLength = tPathLength;
280 
281  if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) {
282  G4double tau = tPathLength/lambdaeff;
283  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
284  // << " Leff= " << lambdaeff << " tau= " << tau << G4endl;
285  // small step
286  if(tau < numlimit) {
287  zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
288 
289  // medium step
290  } else {
291  G4double e1 = 0.0;
292  if(currentRange > tPathLength) {
293  e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
294  }
295  e1 = 0.5*(e1 + preKinEnergy);
296  cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
297  lambdaeff = GetTransportMeanFreePath(particle,e1);
298  zPathLength = lambdaeff*(1.0 - G4Exp(-tPathLength/lambdaeff));
299  }
300  } else { lambdaeff = DBL_MAX; }
301  //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
302  return zPathLength;
303 }
304 
305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
306 
308 {
309  // initialisation of single scattering x-section
310  xtsec = 0.0;
311  cosThetaMin = cosTetMaxNuc;
312  /*
313  G4cout << "ComputeTrueStepLength: Step= " << geomStepLength
314  << " Lambda= " << lambdaeff
315  << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
316  */
317  // pathalogical case
318  if(lambdaeff == DBL_MAX) {
319  singleScatteringMode = true;
320  zPathLength = geomStepLength;
321  tPathLength = geomStepLength;
322  cosThetaMin = 1.0;
323 
324  // normal case
325  } else {
326 
327  // small step use only single scattering
328  static const G4double singleScatLimit = 1.0e-7;
329  if(geomStepLength < lambdaeff*singleScatLimit*(1.0 - cosTetMaxNuc)) {
330  singleScatteringMode = true;
331  cosThetaMin = 1.0;
332  zPathLength = geomStepLength;
333  tPathLength = geomStepLength;
334 
335  // step defined by transportation
336  } else if(geomStepLength != zPathLength) {
337 
338  // step defined by transportation
339  zPathLength = geomStepLength;
340  G4double tau = geomStepLength/lambdaeff;
341  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
342 
343  // energy correction for a big step
344  if(tau > numlimit) {
345  G4double e1 = 0.0;
346  if(currentRange > tPathLength) {
347  e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
348  }
349  e1 = 0.5*(e1 + preKinEnergy);
350  cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
351  lambdaeff = GetTransportMeanFreePath(particle,e1);
352  tau = zPathLength/lambdaeff;
353 
354  if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
355  else { tPathLength = currentRange; }
356  }
357  }
358  }
359 
360  // check of step length
361  // define threshold angle between single and multiple scattering
362  if(!singleScatteringMode) { cosThetaMin = 1.0 - 1.5*tPathLength/lambdaeff; }
363 
364  // recompute transport cross section - do not change energy
365  // anymore - cannot be applied for big steps
366  if(cosThetaMin > cosTetMaxNuc) {
367 
368  // new computation
369  G4double cross = ComputeXSectionPerVolume();
370  //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec << G4endl;
371  if(cross <= 0.0) {
372  singleScatteringMode = true;
373  tPathLength = zPathLength;
374  lambdaeff = DBL_MAX;
375  cosThetaMin = 1.0;
376  } else if(xtsec > 0.0) {
377 
378  lambdaeff = 1./cross;
379  G4double tau = zPathLength*cross;
380  if(tau < numlimit) {
381  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
382  }
383  else if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
384  else { tPathLength = currentRange; }
385 
386  if(tPathLength > currentRange) { tPathLength = currentRange; }
387  }
388  }
389 
390  /*
391  G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
392  <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
393  G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
394  << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
395  << " e(MeV)= " << preKinEnergy/MeV << " "
396  << singleScatteringMode << G4endl;
397  */
398  return tPathLength;
399 }
400 
401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
402 
405  G4double safety)
406 {
407  fDisplacement.set(0.0,0.0,0.0);
408  //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
409  // << particle->GetParticleName() << G4endl;
410 
411  // ignore scattering for zero step length and energy below the limit
412  if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0)
413  { return fDisplacement; }
414 
415  G4double invlambda = 0.0;
416  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
417 
418  // use average kinetic energy over the step
419  G4double cut = (*currentCuts)[currentMaterialIndex];
420  if(fixedCut > 0.0) { cut = fixedCut; }
421  /*
422  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
423  << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
424  << " xmsc= " << tPathLength*invlambda
425  << " safety= " << safety << G4endl;
426  */
427 
428  // step limit due msc
429  G4int nMscSteps = 1;
430  G4double x0 = tPathLength;
431  G4double z0 = x0*invlambda;
432  G4double zzz = 0.0;
433 
434  // large scattering angle case - two step approach
435 
436  if(!singleScatteringMode) {
437  static const G4double zzmin = 0.05;
438  if(z0 > zzmin && safety > tlimitminfix) {
439  x0 *= 0.5;
440  z0 *= 0.5;
441  nMscSteps = 2;
442  }
443  if(z0 > zzmin) { zzz = G4Exp(-1.0/z0); }
444  }
445 
446  // step limit due to single scattering
447  G4double x1 = 2*tPathLength;
448  if(0.0 < xtsec) { x1 = -G4Log(G4UniformRand())/xtsec; }
449 
450  // no scattering case
451  if(singleScatteringMode && x1 > tPathLength)
452  { return fDisplacement; }
453 
454  const G4ElementVector* theElementVector =
455  currentMaterial->GetElementVector();
456  G4int nelm = currentMaterial->GetNumberOfElements();
457 
458  // geometry
459  G4double sint, cost, phi;
460  G4ThreeVector temp(0.0,0.0,1.0);
461 
462  // current position and direction relative to the end point
463  // because of magnetic field geometry is computed relatively to the
464  // end point of the step
465  G4ThreeVector dir(0.0,0.0,1.0);
466  fDisplacement.set(0.0,0.0,-zPathLength);
467 
468  G4double mscfac = zPathLength/tPathLength;
469 
470  // start a loop
471  G4double x2 = x0;
472  G4double step, z;
473  G4bool singleScat;
474  /*
475  G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
476  << " 1-cost1= " << 1 - cosThetaMin << " " << singleScatteringMode
477  << " xtsec= " << xtsec << " " << nMscSteps << G4endl;
478  */
479  do {
480 
481  //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl;
482  // single scattering case
483  if(singleScatteringMode && x1 > x2) { break; }
484 
485  // what is next single of multiple?
486  if(x1 <= x2) {
487  step = x1;
488  singleScat = true;
489  } else {
490  step = x2;
491  singleScat = false;
492  }
493 
494  //G4cout << "# step(mm)= "<< step<< " singlScat= "<< singleScat << G4endl;
495 
496  // new position
497  fDisplacement += step*mscfac*dir;
498 
499  if(singleScat) {
500 
501  // select element
502  G4int i = 0;
503  if(nelm > 1) {
504  G4double qsec = G4UniformRand()*xtsec;
505  for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
506  }
507  G4double cosTetM =
508  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
509  //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " "
510  // << prob[i] << G4endl;
511  temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
512 
513  // direction is changed
514  temp.rotateUz(dir);
515  dir = temp;
516 
517  // new proposed step length
518  x2 -= step;
519  x1 = -G4Log(G4UniformRand())/xtsec;
520 
521  // multiple scattering
522  } else {
523  --nMscSteps;
524  x1 -= step;
525  x2 = x0;
526 
527  // sample z in interval 0 - 1
528  do {
529  z = -z0*G4Log(1.0 - (1.0 - zzz)*G4UniformRand());
530  } while(z > 1.0);
531  cost = 1.0 - 2.0*z/*factCM*/;
532  if(cost > 1.0) { cost = 1.0; }
533  else if(cost < -1.0) { cost =-1.0; }
534  sint = sqrt((1.0 - cost)*(1.0 + cost));
535  phi = twopi*G4UniformRand();
536  G4double vx1 = sint*cos(phi);
537  G4double vy1 = sint*sin(phi);
538 
539  // change direction
540  temp.set(vx1,vy1,cost);
541  temp.rotateUz(dir);
542  dir = temp;
543 
544  // lateral displacement
545  if (latDisplasment && safety > tlimitminfix) {
546  G4double rms = invsqrt12*sqrt(2*z0);
547  G4double r = x0*mscfac;
548  G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0));
549  G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0));
550  G4double dz;
551  G4double d = r*r - dx*dx - dy*dy;
552 
553  // change position
554  if(d >= 0.0) {
555  dz = sqrt(d) - r;
556  temp.set(dx,dy,dz);
557  temp.rotateUz(dir);
558  fDisplacement += temp;
559  }
560  }
561  }
562  } while (0 < nMscSteps);
563 
564  dir.rotateUz(oldDirection);
565 
566  //G4cout<<"G4WentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl;
567  // end of sampling -------------------------------
568 
569  fParticleChange->ProposeMomentumDirection(dir);
570 
571  // lateral displacement
572  fDisplacement.rotateUz(oldDirection);
573 
574  /*
575  G4cout << " r(mm)= " << fDisplacement.mag()
576  << " safety= " << safety
577  << " trueStep(mm)= " << tPathLength
578  << " geomStep(mm)= " << zPathLength
579  << " x= " << fDisplacement.x()
580  << " y= " << fDisplacement.y()
581  << " z= " << fDisplacement.z()
582  << G4endl;
583  */
584 
585  //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
586  return fDisplacement;
587 }
588 
589 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
590 
591 G4double G4WentzelVIModel::ComputeXSectionPerVolume()
592 {
593  // prepare recomputation of x-sections
594  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
595  const G4double* theAtomNumDensityVector =
596  currentMaterial->GetVecNbOfAtomsPerVolume();
597  G4int nelm = currentMaterial->GetNumberOfElements();
598  if(nelm > nelments) {
599  nelments = nelm;
600  xsecn.resize(nelm);
601  prob.resize(nelm);
602  }
603  G4double cut = (*currentCuts)[currentMaterialIndex];
604  if(fixedCut > 0.0) { cut = fixedCut; }
605  // cosTetMaxNuc = wokvi->GetCosThetaNuc();
606 
607  // check consistency
608  xtsec = 0.0;
609  if(cosTetMaxNuc > cosThetaMin) { return 0.0; }
610 
611  // loop over elements
612  G4double xs = 0.0;
613  for (G4int i=0; i<nelm; ++i) {
614  G4double costm =
615  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
616  G4double density = theAtomNumDensityVector[i];
617 
618  G4double esec = 0.0;
619  if(costm < cosThetaMin) {
620 
621  // recompute the transport x-section
622  if(1.0 > cosThetaMin) {
623  xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosThetaMin);
624  }
625  // recompute the total x-section
626  G4double nucsec = wokvi->ComputeNuclearCrossSection(cosThetaMin, costm);
627  esec = wokvi->ComputeElectronCrossSection(cosThetaMin, costm);
628  nucsec += esec;
629  if(nucsec > 0.0) { esec /= nucsec; }
630  xtsec += nucsec*density;
631  }
632  xsecn[i] = xtsec;
633  prob[i] = esec;
634  //G4cout << i << " xs= " << xs << " xtsec= " << xtsec
635  // << " 1-cosThetaMin= " << 1-cosThetaMin
636  // << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
637  }
638 
639  //G4cout << "ComputeXS result: xsec(1/mm)= " << xs
640  // << " txsec(1/mm)= " << xtsec <<G4endl;
641  return xs;
642 }
643 
644 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void set(double x, double y, double z)
virtual ~G4WentzelVIModel()
G4double facgeom
Definition: G4VMscModel.hh:177
ThreeVector shoot(const G4int Ap, const G4int Af)
static G4LossTableManager * Instance()
G4WentzelVIModel(const G4String &nam="WentzelVIUni")
std::vector< G4Element * > G4ElementVector
Double_t x2[nxs]
Definition: Style.C:19
G4double GetKineticEnergy() const
Float_t d
Definition: plot.C:237
const G4DynamicParticle * GetDynamicParticle() const
G4double facrange
Definition: G4VMscModel.hh:176
G4ThreeVector SampleSingleScattering(G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio=0.0)
G4double GetProductionCut(G4int index) const
const char * p
Definition: xmltok.h:285
G4StepStatus GetStepStatus() const
virtual G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety)
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:246
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4bool latDisplasment
Definition: G4VMscModel.hh:189
G4ParticleDefinition * GetDefinition() const
const G4Step * GetStep() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
int G4int
Definition: G4Types.hh:78
G4double SetupTarget(G4int Z, G4double cut=DBL_MAX)
virtual G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep)
G4StepStatus
Definition: G4StepStatus.hh:49
G4StepPoint * GetPreStepPoint() const
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:308
G4double density
Definition: TRTMaterials.hh:39
G4double ComputeTransportCrossSectionPerAtom(G4double CosThetaMax)
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:204
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:426
#define G4UniformRand()
Definition: Randomize.hh:87
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:288
Float_t Z
Definition: plot.C:39
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:89
bool G4bool
Definition: G4Types.hh:79
G4double SetupKinematic(G4double kinEnergy, const G4Material *mat)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
Double_t x1[nxs]
Definition: Style.C:18
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:256
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:185
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:345
G4double GetRadlen() const
Definition: G4Material.hh:218
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double ComputeSafety(const G4ThreeVector &position, G4double limit)
Definition: G4VMscModel.hh:238
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
jump r
Definition: plot.C:36
G4double ComputeElectronCrossSection(G4double CosThetaMin, G4double CosThetaMax)
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:178
void StartTracking(G4Track *)
G4double ComputeNuclearCrossSection(G4double CosThetaMin, G4double CosThetaMax)
G4double PolarAngleLimit() const
Definition: G4VEmModel.hh:620
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetSafety() const
tuple z
Definition: test.py:28
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:186
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
double G4double
Definition: G4Types.hh:76
virtual G4double ComputeTrueStepLength(G4double geomStepLength)
G4ProductionCuts * GetProductionCuts() const
#define DBL_MAX
Definition: templates.hh:83
TDirectory * dir
Definition: macro.C:5
void Initialise(const G4ParticleDefinition *, G4double CosThetaLim)
virtual G4double ComputeGeomPathLength(G4double truePathLength)