Geant4  10.01.p03
G4GoudsmitSaundersonMscModel.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: G4GoudsmitSaundersonMscModel.cc 88979 2015-03-17 10:10:21Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 // File name: G4GoudsmitSaundersonMscModel
33 //
34 // Author: Omrane Kadri
35 //
36 // Creation date: 20.02.2009
37 //
38 // Modifications:
39 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
40 //
41 // 15.04.2009 O.Kadri: cleanup: discard no scattering and single scattering theta
42 // sampling from SampleCosineTheta() which means the splitting
43 // step into two sub-steps occur only for msc regime
44 //
45 // 12.06.2009 O.Kadri: linear log-log extrapolation of lambda0 & lambda1 between 1 GeV - 100 TeV
46 // adding a theta min limit due to screening effect of the atomic nucleus
47 // 26.08.2009 O.Kadri: Cubic Spline interpolation was replaced with polynomial method
48 // within CalculateIntegrals method
49 // 05.10.2009 O.Kadri: tuning small angle theta distributions
50 // assuming the case of lambdan<1 as single scattering regime
51 // tuning theta sampling for theta below the screening angle
52 // 08.02.2010 O.Kadri: bugfix in compound xsection calculation and small angle computation
53 // adding a rejection condition to hard collision angular sampling
54 // ComputeTruePathLengthLimit was taken from G4WentzelVIModel
55 // 26.03.2010 O.Kadri: direct xsection calculation not inverse of the inverse
56 // angular sampling without large angle rejection method
57 // longitudinal displacement is computed exactly from <z>
58 // 12.05.2010 O.Kadri: exchange between target and projectile has as a condition the particle type (e-/e-)
59 // some cleanup to minimize time consuming (adding lamdan12 & Qn12, changing the error to 1.0e-12 for scrA)
60 //
61 //
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 //REFERENCES:
64 //Ref.1:E. Benedito et al.,"Mixed simulation ... cross-sections", NIMB 174 (2001) pp 91-110;
65 //Ref.2:I. Kawrakow et al.,"On the condensed ... transport",NIMB 142 (1998) pp 253-280;
66 //Ref.3:I. Kawrakow et al.,"On the representation ... calculations",NIMB 134 (1998) pp 325-336;
67 //Ref.4:Bielajew et al.,".....", NIMB 173 (2001) 332-343;
68 //Ref.5:F. Salvat et al.,"ELSEPA--Dirac partial ...molecules", Comp.Phys.Comm.165 (2005) pp 157-190;
69 //Ref.6:G4UrbanMscModel G4 9.2;
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 
74 #include "G4PhysicalConstants.hh"
75 #include "G4SystemOfUnits.hh"
76 
78 #include "G4MaterialCutsCouple.hh"
79 #include "G4DynamicParticle.hh"
80 #include "G4Electron.hh"
81 #include "G4Positron.hh"
82 
83 #include "G4LossTableManager.hh"
84 #include "G4Track.hh"
85 #include "G4PhysicsTable.hh"
86 #include "Randomize.hh"
87 #include "G4Log.hh"
88 #include "G4Exp.hh"
89 #include <fstream>
90 
91 using namespace std;
92 
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101  : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV)
102 {
106  =lambda11=mass=0.0;
108 
109  fr=0.02,rangeinit=0.,masslimite=0.6*MeV,
110  particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
111  tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
114  inside=false;insideskin=false;
115  samplez=false;
116  firstStep = true;
117 
118  currentCouple = 0;
119  fParticleChange = 0;
120 
122 
123  if(ener[0] < 0.0){
124  G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
126  }
127 }
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 {
131  delete GSTable;
132 }
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135  const G4DataVector&)
136 {
138  SetParticle(p);
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143 
144 G4double
146  G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
147 {
148  G4double kinEnergy = kineticEnergy;
149  if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
150  if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
151 
152  G4double cs(0.0), cs0(0.0);
153  CalculateIntegrals(p,Z,kinEnergy,cs0,cs);
154 
155  return cs;
156 }
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
158 
161  G4double)
162 {
163  fDisplacement.set(0.0,0.0,0.0);
164  G4double kineticEnergy = currentKinEnergy;
165  //dynParticle->GetKineticEnergy();
166  if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)||
167  (tPathLength/tausmall < lambda1)) { return fDisplacement; }
168 
170  // Effective energy
171  G4double eloss = 0.0;
172  if (tPathLength > currentRange*dtrl) {
173  eloss = kineticEnergy -
175  } else {
176  eloss = tPathLength*GetDEDX(particle,kineticEnergy,currentCouple);
177  }
178  /*
179  G4double ttau = kineticEnergy/electron_mass_c2;
180  G4double ttau2 = ttau*ttau;
181  G4double epsilonpp = eloss/kineticEnergy;
182  G4double cst1 = epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
183  kineticEnergy *= (1 - cst1);
184  */
185  kineticEnergy -= 0.5*eloss;
186 
188  // additivity rule for mixture and compound xsection's
189  const G4Material* mat = currentCouple->GetMaterial();
190  const G4ElementVector* theElementVector = mat->GetElementVector();
191  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
192  G4int nelm = mat->GetNumberOfElements();
193  G4double s0(0.0), s1(0.0);
194  lambda0 = 0.0;
195  for(G4int i=0;i<nelm;i++)
196  {
197  CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,s0,s1);
198  lambda0 += (theAtomNumDensityVector[i]*s0);
199  }
200  if(lambda0>0.0) { lambda0 =1./lambda0; }
201 
202  // Newton-Raphson root's finding method of scrA from:
203  // Sig1(PWA)/Sig0(PWA)=g1=2*scrA*((1+scrA)*log(1+1/scrA)-1)
204  G4double g1=0.0;
205  if(lambda1>0.0) { g1 = lambda0/lambda1; }
206 
207  G4double logx0,x1,delta;
208  G4double x0=g1*0.5;
209  // V.Ivanchenko added limit of the loop
210  for(G4int i=0;i<1000;++i)
211  {
212  logx0 = G4Log(1.+1./x0);
213  x1 = x0-(x0*((1.+x0)*logx0-1.0)-g1*0.5)/( (1.+2.*x0)*logx0-2.0);
214 
215  // V.Ivanchenko cut step size of iterative procedure
216  if(x1 < 0.0) { x1 = 0.5*x0; }
217  else if(x1 > 2*x0) { x1 = 2*x0; }
218  else if(x1 < 0.5*x0) { x1 = 0.5*x0; }
219  delta = std::fabs( x1 - x0 );
220  x0 = x1;
221  if(delta < 1.0e-3*x1) { break;}
222  }
223  G4double scrA = x1;
224 
225  G4double lambdan=0.;
226 
227  if(lambda0>0.0) { lambdan=tPathLength/lambda0; }
228  if(lambdan<=1.0e-12) { return fDisplacement; }
229 
230  //G4cout << "E(eV)= " << kineticEnergy/eV << " L0= " << lambda0
231  // << " L1= " << lambda1 << G4endl;
232 
233  G4double Qn1 = lambdan *g1;//2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
234  G4double Qn12 = 0.5*Qn1;
235 
236  G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
237  G4double cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
238  G4double us=0.0,vs=0.0,ws=1.0,wss=0.,x_coord=0.0,y_coord=0.0,z_coord=1.0;
239 
240  G4double epsilon1=G4UniformRand();
241  G4double expn = G4Exp(-lambdan);
242 
243  if(epsilon1<expn)// no scattering
244  { return fDisplacement; }
245  else if((epsilon1<((1.+lambdan)*expn))||(lambdan<1.))//single or plural scattering (Rutherford DCS's)
246  {
247  G4double xi=G4UniformRand();
248  xi= 2.*scrA*xi/(1.-xi + scrA);
249  if(xi<0.)xi=0.;
250  else if(xi>2.)xi=2.;
251  ws=(1. - xi);
252  wss=std::sqrt(xi*(2.-xi));
253  G4double phi0=CLHEP::twopi*G4UniformRand();
254  us=wss*cos(phi0);
255  vs=wss*sin(phi0);
256  }
257  else // multiple scattering
258  {
259  // Ref.2 subsection 4.4 "The best solution found"
260  // Sample first substep scattering angle
261  SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
262  G4double phi1 = CLHEP::twopi*G4UniformRand();
263  cosPhi1 = cos(phi1);
264  sinPhi1 = sin(phi1);
265 
266  // Sample second substep scattering angle
267  SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
268  G4double phi2 = CLHEP::twopi*G4UniformRand();
269  cosPhi2 = cos(phi2);
270  sinPhi2 = sin(phi2);
271 
272  // Overall scattering direction
273  us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2)
274  + cosTheta2*sinTheta1*cosPhi1;
275  vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2)
276  + cosTheta2*sinTheta1*sinPhi1;
277  ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2;
278 
279  //small angle approximation for theta less than screening angle
280  if(1. - ws < 0.5*scrA)
281  {
282  G4int i=0;
283  do{i++;
284  ws=1.+Qn12*G4Log(G4UniformRand());
285  }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
286  if(i>=19)ws=cos(sqrt(scrA));
287  wss=std::sqrt((1. - ws)*(1. + ws));
288  us=wss*std::cos(phi1);
289  vs=wss*std::sin(phi1);
290  }
291  }
292 
293  //G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
294  G4ThreeVector newDirection(us,vs,ws);
295  newDirection.rotateUz(oldDirection);
297 
298  // corresponding to error less than 1% in the exact formula of <z>
299  if(Qn1<0.02) { z_coord = 1.0 - Qn1*(0.5 - Qn1/6.); }
300  else { z_coord = (1.-G4Exp(-Qn1))/Qn1; }
301  G4double rr = zPathLength*std::sqrt((1.- z_coord*z_coord)/(1.-ws*ws));
302  x_coord = rr*us;
303  y_coord = rr*vs;
304 
305  // displacement is computed relatively to the end point
306  z_coord -= 1.0;
307  z_coord *= zPathLength;
308  /*
309  G4cout << "G4GS::SampleSecondaries: e(MeV)= " << kineticEnergy
310  << " sinTheta= " << sqrt(1.0 - ws*ws)
311  << " trueStep(mm)= " << tPathLength
312  << " geomStep(mm)= " << zPathLength
313  << G4endl;
314  */
315 
316  fDisplacement.set(x_coord,y_coord,z_coord);
317  fDisplacement.rotateUz(oldDirection);
318 
319  return fDisplacement;
320  }
321 
322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
323 
324 void
326  G4double &cost, G4double &sint)
327 {
328  G4double r1,tet,xi=0.;
329  G4double Qn1 = 2.* lambdan;
330  if(scrA < 10.) { Qn1 *= scrA*((1.+scrA)*G4Log(1.+1./scrA)-1.); }
331  else { Qn1*= (1.0 - 0.5/scrA - 0.5/(scrA*scrA)) ; }
332  if (Qn1<0.001)
333  {
334  do{
335  r1=G4UniformRand();
336  xi=-0.5*Qn1*G4Log(G4UniformRand());
337  tet=acos(1.-xi);
338  }while(tet*r1*r1>sin(tet));
339  }
340  else if(Qn1>0.5) { xi=2.*G4UniformRand(); }//isotropic distribution
341  else{ xi=2.*(GSTable->SampleTheta(lambdan,scrA,G4UniformRand()));}
342 
343 
344  if(xi<0.)xi=0.;
345  else if(xi>2.)xi=2.;
346 
347  cost=(1. - xi);
348  sint=sqrt(xi*(2.-xi));
349 
350 }
351 
352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
353 // Polynomial log-log interpolation of Lambda0 and Lambda1 between 100 eV - 1 GeV
354 // linear log-log extrapolation between 1 GeV - 100 TeV
355 
356 void
358  G4double kinEnergy,G4double &Sig0,
359  G4double &Sig1)
360 {
361  G4double x1,x2,y1,y2,acoeff,bcoeff;
362  G4double kineticE = kinEnergy;
363  if(kineticE<lowKEnergy)kineticE=lowKEnergy;
364  if(kineticE>highKEnergy)kineticE=highKEnergy;
365  kineticE /= eV;
366  G4double logE=G4Log(kineticE);
367 
368  G4int iZ = G4int(Z);
369  if(iZ > 103) iZ = 103;
370 
371  G4int enerInd=0;
372  for(G4int i=0;i<105;i++)
373  {
374  if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;}
375  }
376 
377  if(p==G4Electron::Electron())
378  {
379  if(kineticE<=1.0e+9)//Interpolation of the form y=ax²+b
380  {
381  x1=ener[enerInd];
382  x2=ener[enerInd+1];
383  y1=TCSE[iZ-1][enerInd];
384  y2=TCSE[iZ-1][enerInd+1];
385  acoeff=(y2-y1)/(x2*x2-x1*x1);
386  bcoeff=y2-acoeff*x2*x2;
387  Sig0=acoeff*logE*logE+bcoeff;
388  Sig0 =G4Exp(Sig0);
389  y1=FTCSE[iZ-1][enerInd];
390  y2=FTCSE[iZ-1][enerInd+1];
391  acoeff=(y2-y1)/(x2*x2-x1*x1);
392  bcoeff=y2-acoeff*x2*x2;
393  Sig1=acoeff*logE*logE+bcoeff;
394  Sig1=G4Exp(Sig1);
395  }
396  else //Interpolation of the form y=ax+b
397  {
398  x1=ener[104];
399  x2=ener[105];
400  y1=TCSE[iZ-1][104];
401  y2=TCSE[iZ-1][105];
402  Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
403  Sig0=G4Exp(Sig0);
404  y1=FTCSE[iZ-1][104];
405  y2=FTCSE[iZ-1][105];
406  Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
407  Sig1=G4Exp(Sig1);
408  }
409  }
410  if(p==G4Positron::Positron())
411  {
412  if(kinEnergy<=1.0e+9)
413  {
414  x1=ener[enerInd];
415  x2=ener[enerInd+1];
416  y1=TCSP[iZ-1][enerInd];
417  y2=TCSP[iZ-1][enerInd+1];
418  acoeff=(y2-y1)/(x2*x2-x1*x1);
419  bcoeff=y2-acoeff*x2*x2;
420  Sig0=acoeff*logE*logE+bcoeff;
421  Sig0 =G4Exp(Sig0);
422  y1=FTCSP[iZ-1][enerInd];
423  y2=FTCSP[iZ-1][enerInd+1];
424  acoeff=(y2-y1)/(x2*x2-x1*x1);
425  bcoeff=y2-acoeff*x2*x2;
426  Sig1=acoeff*logE*logE+bcoeff;
427  Sig1=G4Exp(Sig1);
428  }
429  else
430  {
431  x1=ener[104];
432  x2=ener[105];
433  y1=TCSP[iZ-1][104];
434  y2=TCSP[iZ-1][105];
435  Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
436  Sig0 =G4Exp(Sig0);
437  y1=FTCSP[iZ-1][104];
438  y2=FTCSP[iZ-1][105];
439  Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
440  Sig1=G4Exp(Sig1);
441  }
442  }
443 
444  Sig0 *= barn;
445  Sig1 *= barn;
446 
447 }
448 
449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
450 
452 {
454  firstStep = true;
455  inside = false;
456  insideskin = false;
457  tlimit = geombig;
459 }
460 
461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
462 //t->g->t step transformations taken from Ref.6
463 
464 G4double
466  G4double& currentMinimalStep)
467 {
468  tPathLength = currentMinimalStep;
469  const G4DynamicParticle* dp = track.GetDynamicParticle();
470  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
471  G4StepStatus stepStatus = sp->GetStepStatus();
477 
479 
480  // stop here if small range particle
481  if(inside || tPathLength < tlimitminfix) {
482  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
483  }
485 
486  G4double presafety = sp->GetSafety();
487 
488  //G4cout << "G4GS::StepLimit tPathLength= "
489  // <<tPathLength<<" safety= " << presafety
490  // << " range= " <<currentRange<< " lambda= "<<lambda1
491  // << " Alg: " << steppingAlgorithm <<G4endl;
492 
493  // far from geometry boundary
494  if(currentRange < presafety)
495  {
496  inside = true;
497  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
498  }
499 
500  // standard version
501  //
503  {
504  //compute geomlimit and presafety
505  G4double geomlimit = ComputeGeomLimit(track, presafety, tPathLength);
506 
507  // is far from boundary
508  if(currentRange <= presafety)
509  {
510  inside = true;
511  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
512  }
513 
514  smallstep += 1.;
515  insideskin = false;
516 
517  if(firstStep || stepStatus == fGeomBoundary)
518  {
520  if(firstStep) smallstep = 1.e10;
521  else smallstep = 1.;
522 
523  //define stepmin here (it depends on lambda!)
524  //rough estimation of lambda_elastic/lambda_transport
526  rat = 1.e-3/(rat*(10.+rat)) ;
527  //stepmin ~ lambda_elastic
528  stepmin = rat*lambda1;
530  //define tlimitmin
531  tlimitmin = 10.*stepmin;
533 
534  //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
535  // << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
536  // constraint from the geometry
537  if((geomlimit < geombig) && (geomlimit > geommin))
538  {
539  if(stepStatus == fGeomBoundary)
540  tgeom = geomlimit/facgeom;
541  else
542  tgeom = 2.*geomlimit/facgeom;
543  }
544  else
545  tgeom = geombig;
546 
547  }
548 
549  //step limit
551  if(tlimit < facsafety*presafety)
552  tlimit = facsafety*presafety;
553 
554  //lower limit for tlimit
556 
557  if(tlimit > tgeom) tlimit = tgeom;
558 
559  //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
560  // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
561 
562  // shortcut
563  if((tPathLength < tlimit) && (tPathLength < presafety) &&
564  (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
565  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
566 
567  // step reduction near to boundary
568  if(smallstep < skin)
569  {
570  tlimit = stepmin;
571  insideskin = true;
572  }
573  else if(geomlimit < geombig)
574  {
575  if(geomlimit > skindepth)
576  {
577  if(tlimit > geomlimit-0.999*skindepth)
578  tlimit = geomlimit-0.999*skindepth;
579  }
580  else
581  {
582  insideskin = true;
583  if(tlimit > stepmin) tlimit = stepmin;
584  }
585  }
586 
587  if(tlimit < stepmin) tlimit = stepmin;
588 
590 
591  }
592  // for 'normal' simulation with or without magnetic field
593  // there no small step/single scattering at boundaries
594  else if(steppingAlgorithm == fUseSafety)
595  {
596  // compute presafety again if presafety <= 0 and no boundary
597  // i.e. when it is needed for optimization purposes
598  if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
599  presafety = ComputeSafety(sp->GetPosition(),tPathLength);
600 
601  // is far from boundary
602  if(currentRange < presafety)
603  {
604  inside = true;
605  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
606  }
607 
608  if(firstStep || stepStatus == fGeomBoundary)
609  {
611  fr = facrange;
612  // 9.1 like stepping for e+/e- only (not for muons,hadrons)
613  if(mass < masslimite)
614  {
615  if(lambda1 > currentRange)
616  rangeinit = lambda1;
617  if(lambda1 > lambdalimit)
618  fr *= 0.75+0.25*lambda1/lambdalimit;
619  }
620 
621  //lower limit for tlimit
623  rat = 1.e-3/(rat*(10.+rat)) ;
624  tlimitmin = 10.*lambda1*rat;
626  }
627  //step limit
628  tlimit = fr*rangeinit;
629 
630  if(tlimit < facsafety*presafety)
631  tlimit = facsafety*presafety;
632 
633  //lower limit for tlimit
635 
637  }
638 
639  // version similar to 7.1 (needed for some experiments)
640  else
641  {
642  if (stepStatus == fGeomBoundary)
643  {
645  else tlimit = facrange*lambda1;
646 
649  }
650  }
651  //G4cout << "tPathLength= " << tPathLength
652  // << " currentMinimalStep= " << currentMinimalStep << G4endl;
653  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
654 }
655 
656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
657 // taken from Ref.6
659 {
660  firstStep = false;
661  par1 = -1. ;
662  par2 = par3 = 0. ;
663 
664  // do the true -> geom transformation
666 
667  // z = t for very small tPathLength
668  if(tPathLength < tlimitminfix) { return zPathLength; }
669 
670  // this correction needed to run MSC with eIoni and eBrem inactivated
671  // and makes no harm for a normal run
674 
676 
677  if ((tau <= tausmall) || insideskin) {
680  return zPathLength;
681  }
682 
683  G4double zmean = tPathLength;
684  if (tPathLength < currentRange*dtrl) {
685  if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
686  else zmean = lambda1*(1.-G4Exp(-tau));
687  } else if(currentKinEnergy < mass || tPathLength == currentRange) {
688  par1 = 1./currentRange ;
689  par2 = 1./(par1*lambda1) ;
690  par3 = 1.+par2 ;
692  zmean = (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3) ;
693  else
694  zmean = 1./(par1*par3) ;
695  } else {
697 
699 
701  par2 = 1./(par1*lambda1) ;
702  par3 = 1.+par2 ;
703  zmean = (1.-G4Exp(par3*G4Log(lambda11/lambda1)))/(par1*par3) ;
704  }
705 
706  zPathLength = zmean ;
707  // sample z
708  if(samplez) {
709 
710  const G4double ztmax = 0.99;
711  G4double zt = zmean/tPathLength ;
712 
713  if (tPathLength > stepmin && zt < ztmax) {
714 
715  G4double u,cz1;
716  if(zt >= 0.333333333) {
717 
718  G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
719  cz1 = 1.+cz ;
720  G4double u0 = cz/cz1 ;
721  G4double grej ;
722  do {
723  u = G4Exp(G4Log(G4UniformRand())/cz1) ;
724  grej = exp(cz*G4Log(u/u0))*(1.-u)/(1.-u0) ;
725  } while (grej < G4UniformRand()) ;
726 
727  } else {
728  cz1 = 1./zt-1.;
729  u = 1.-G4Exp(G4Log(G4UniformRand())/cz1) ;
730  }
732  }
733  }
735  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda1 << G4endl;
736 
737  return zPathLength;
738 }
739 
740 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
741 // taken from Ref.6
742 G4double
744 {
745  // step defined other than transportation
746  if(geomStepLength == zPathLength && tPathLength <= currentRange)
747  return tPathLength;
748 
749  // t = z for very small step
750  zPathLength = geomStepLength;
751  tPathLength = geomStepLength;
752  if(geomStepLength < tlimitminfix) return tPathLength;
753 
754  // recalculation
755  if((geomStepLength > lambda1*tausmall) && !insideskin)
756  {
757  if(par1 < 0.)
758  tPathLength = -lambda1*G4Log(1.-geomStepLength/lambda1) ;
759  else
760  {
761  if(par1*par3*geomStepLength < 1.)
762  tPathLength = (1.-G4Exp(G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
763  else
765  }
766  }
767  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
768  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
769 
770  return tPathLength;
771 }
772 
773 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
774 //Total & first transport x sections for e-/e+ generated from ELSEPA code
775 
777 {
778  G4String filename = "XSECTIONS.dat";
779 
780  char* path = getenv("G4LEDATA");
781  if (!path)
782  {
783  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()","em0006",
785  "Environment variable G4LEDATA not defined");
786  return;
787  }
788 
789  G4String pathString(path);
790  G4String dirFile = pathString + "/msc_GS/" + filename;
791  std::ifstream infile(dirFile, std::ios::in);
792  if( !infile.is_open()) {
794  ed << "Data file <" + dirFile + "> is not opened!";
795  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
796  "em0003",FatalException,ed);
797  return;
798  }
799 
800  // Read parameters from tables and take logarithms
801  G4double aRead;
802  for(G4int i=0 ; i<106 ;i++){
803  infile >> aRead;
804  if(infile.fail()) {
806  ed << "Error reading <" + dirFile + "> loop #1 i= " << i;
807  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
808  "em0003",FatalException,ed);
809  return;
810  } else {
811  if(aRead > 0.0) { aRead = G4Log(aRead); }
812  else { aRead = 0.0; }
813  }
814  ener[i]=aRead;
815  }
816  for(G4int j=0;j<103;j++){
817  for(G4int i=0;i<106;i++){
818  infile >> aRead;
819  if(infile.fail()) {
821  ed << "Error reading <" + dirFile + "> loop #2 j= " << j
822  << "; i= " << i;
823  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
824  "em0003",FatalException,ed);
825  return;
826  } else {
827  if(aRead > 0.0) { aRead = G4Log(aRead); }
828  else { aRead = 0.0; }
829  }
830  TCSE[j][i]=aRead;
831  }
832  }
833  for(G4int j=0;j<103;j++){
834  for(G4int i=0;i<106;i++){
835  infile >> aRead;
836  if(infile.fail()) {
838  ed << "Error reading <" + dirFile + "> loop #3 j= " << j
839  << "; i= " << i;
840  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
841  "em0003",FatalException,ed);
842  return;
843  } else {
844  if(aRead > 0.0) { aRead = G4Log(aRead); }
845  else { aRead = 0.0; }
846  }
847  FTCSE[j][i]=aRead;
848  }
849  }
850  for(G4int j=0;j<103;j++){
851  for(G4int i=0;i<106;i++){
852  infile >> aRead;
853  if(infile.fail()) {
855  ed << "Error reading <" + dirFile + "> loop #4 j= " << j
856  << "; i= " << i;
857  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
858  "em0003",FatalException,ed);
859  return;
860  } else {
861  if(aRead > 0.0) { aRead = G4Log(aRead); }
862  else { aRead = 0.0; }
863  }
864  TCSP[j][i]=aRead;
865  }
866  }
867  for(G4int j=0;j<103;j++){
868  for(G4int i=0;i<106;i++){
869  infile >> aRead;
870  if(infile.fail()) {
872  ed << "Error reading <" + dirFile + "> loop #5 j= " << j
873  << "; i= " << i;
874  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
875  "em0003",FatalException,ed);
876  return;
877  } else {
878  if(aRead > 0.0) { aRead = G4Log(aRead); }
879  else { aRead = 0.0; }
880  }
881  FTCSP[j][i]=aRead;
882  }
883  }
884 
885  infile.close();
886 }
887 
888 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4GoudsmitSaundersonMscModel(const G4String &nam="GoudsmitSaunderson")
void CalculateIntegrals(const G4ParticleDefinition *, G4double, G4double, G4double &, G4double &)
G4double facgeom
Definition: G4VMscModel.hh:178
virtual G4double ComputeTrueStepLength(G4double geomStepLength)
static const double MeV
Definition: G4SIunits.hh:193
virtual void StartTracking(G4Track *)
Definition: G4VEmModel.cc:284
static G4LossTableManager * Instance()
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double dtrl
Definition: G4VMscModel.hh:181
const G4DynamicParticle * GetDynamicParticle() const
G4double facrange
Definition: G4VMscModel.hh:177
G4StepStatus GetStepStatus() const
G4bool samplez
Definition: G4VMscModel.hh:189
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:247
G4double skin
Definition: G4VMscModel.hh:180
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
const G4ParticleDefinition * particle
G4double SampleTheta(G4double, G4double, G4double)
G4ParticleDefinition * GetDefinition() const
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
G4StepStatus
Definition: G4StepStatus.hh:49
G4StepPoint * GetPreStepPoint() const
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:309
static const G4double scrA[76 *11]
virtual G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep)
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetPosition() const
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:289
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:90
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
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:346
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
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 const double eV
Definition: G4SIunits.hh:194
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4double facsafety
Definition: G4VMscModel.hh:179
void SetParticle(const G4ParticleDefinition *p)
G4double GetDEDX(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:274
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:426
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double, G4double, G4double)
G4double GetSafety() const
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:187
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static const double TeV
Definition: G4SIunits.hh:197
void SampleCosineTheta(G4double, G4double, G4double &, G4double &)
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
static const double barn
Definition: G4SIunits.hh:95
static const double keV
Definition: G4SIunits.hh:195
virtual G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety)
double G4double
Definition: G4Types.hh:76
static const double mm
Definition: G4SIunits.hh:102
const G4MaterialCutsCouple * currentCouple
const G4Material * GetMaterial() const
virtual G4double ComputeGeomPathLength(G4double truePathLength)