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