Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 79188 2014-02-20 09:22:48Z 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 
93 G4double G4GoudsmitSaundersonMscModel::ener[] = {-1.};
94 G4double G4GoudsmitSaundersonMscModel::TCSE[103][106] ;
95 G4double G4GoudsmitSaundersonMscModel::FTCSE[103][106] ;
96 G4double G4GoudsmitSaundersonMscModel::TCSP[103][106] ;
97 G4double G4GoudsmitSaundersonMscModel::FTCSP[103][106] ;
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101  : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV)
102 {
103  currentKinEnergy=currentRange=skindepth=par1=par2=par3
104  =zPathLength=truePathLength
105  =tausmall=taulim=tlimit=charge=lambdalimit=tPathLength=lambda0=lambda1
106  =lambda11=mass=0.0;
107  currentMaterialIndex = -1;
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;
112  tlimitminfix=1.e-6*mm;stepmin=tlimitminfix;lambdalimit=1.*mm;smallstep=1.e10;
113  theManager=G4LossTableManager::Instance();
114  inside=false;insideskin=false;
115  samplez=false;
116  firstStep = true;
117 
118  GSTable = new G4GoudsmitSaundersonTable();
119 
120  if(ener[0] < 0.0){
121  G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
122  LoadELSEPAXSections();
123  }
124 }
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
127 {
128  delete GSTable;
129 }
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132  const G4DataVector&)
133 {
134  skindepth=skin*stepmin;
135  SetParticle(p);
136  fParticleChange = GetParticleChangeForMSC(p);
137 }
138 
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 
141 G4double
143  G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
144 {
145  G4double kinEnergy = kineticEnergy;
146  if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
147  if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
148 
149  G4double cs(0.0), cs0(0.0);
150  CalculateIntegrals(p,Z,kinEnergy,cs0,cs);
151 
152  return cs;
153 }
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
155 
158 {
159  fDisplacement.set(0.0,0.0,0.0);
160  G4double kineticEnergy = currentKinEnergy;
161  //dynParticle->GetKineticEnergy();
162  if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)||
163  (tPathLength/tausmall < lambda1)) { return fDisplacement; }
164 
165  ///////////////////////////////////////////
166  // Effective energy
167  G4double eloss = 0.0;
168  if (tPathLength > currentRange*dtrl) {
169  eloss = kineticEnergy -
170  GetEnergy(particle,currentRange-tPathLength,currentCouple);
171  } else {
172  eloss = tPathLength*GetDEDX(particle,kineticEnergy,currentCouple);
173  }
174  /*
175  G4double ttau = kineticEnergy/electron_mass_c2;
176  G4double ttau2 = ttau*ttau;
177  G4double epsilonpp = eloss/kineticEnergy;
178  G4double cst1 = epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
179  kineticEnergy *= (1 - cst1);
180  */
181  kineticEnergy -= 0.5*eloss;
182 
183  ///////////////////////////////////////////
184  // additivity rule for mixture and compound xsection's
185  const G4Material* mat = currentCouple->GetMaterial();
186  const G4ElementVector* theElementVector = mat->GetElementVector();
187  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
188  G4int nelm = mat->GetNumberOfElements();
189  G4double s0(0.0), s1(0.0);
190  lambda0 = 0.0;
191  for(G4int i=0;i<nelm;i++)
192  {
193  CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,s0,s1);
194  lambda0 += (theAtomNumDensityVector[i]*s0);
195  }
196  if(lambda0>0.0) { lambda0 =1./lambda0; }
197 
198  // Newton-Raphson root's finding method of scrA from:
199  // Sig1(PWA)/Sig0(PWA)=g1=2*scrA*((1+scrA)*log(1+1/scrA)-1)
200  G4double g1=0.0;
201  if(lambda1>0.0) { g1 = lambda0/lambda1; }
202 
203  G4double logx0,x1,delta;
204  G4double x0=g1*0.5;
205  // V.Ivanchenko added limit of the loop
206  for(G4int i=0;i<1000;++i)
207  {
208  logx0 = G4Log(1.+1./x0);
209  x1 = x0-(x0*((1.+x0)*logx0-1.0)-g1*0.5)/( (1.+2.*x0)*logx0-2.0);
210 
211  // V.Ivanchenko cut step size of iterative procedure
212  if(x1 < 0.0) { x1 = 0.5*x0; }
213  else if(x1 > 2*x0) { x1 = 2*x0; }
214  else if(x1 < 0.5*x0) { x1 = 0.5*x0; }
215  delta = std::fabs( x1 - x0 );
216  x0 = x1;
217  if(delta < 1.0e-3*x1) { break;}
218  }
219  G4double scrA = x1;
220 
221  G4double lambdan=0.;
222 
223  if(lambda0>0.0) { lambdan=tPathLength/lambda0; }
224  if(lambdan<=1.0e-12) { return fDisplacement; }
225 
226  //G4cout << "E(eV)= " << kineticEnergy/eV << " L0= " << lambda0
227  // << " L1= " << lambda1 << G4endl;
228 
229  G4double Qn1 = lambdan *g1;//2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
230  G4double Qn12 = 0.5*Qn1;
231 
232  G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
233  G4double cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
234  G4double us=0.0,vs=0.0,ws=1.0,wss=0.,x_coord=0.0,y_coord=0.0,z_coord=1.0;
235 
236  G4double epsilon1=G4UniformRand();
237  G4double expn = G4Exp(-lambdan);
238 
239  if(epsilon1<expn)// no scattering
240  { return fDisplacement; }
241  else if((epsilon1<((1.+lambdan)*expn))||(lambdan<1.))//single or plural scattering (Rutherford DCS's)
242  {
243  G4double xi=G4UniformRand();
244  xi= 2.*scrA*xi/(1.-xi + scrA);
245  if(xi<0.)xi=0.;
246  else if(xi>2.)xi=2.;
247  ws=(1. - xi);
248  wss=std::sqrt(xi*(2.-xi));
249  G4double phi0=CLHEP::twopi*G4UniformRand();
250  us=wss*cos(phi0);
251  vs=wss*sin(phi0);
252  }
253  else // multiple scattering
254  {
255  // Ref.2 subsection 4.4 "The best solution found"
256  // Sample first substep scattering angle
257  SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
258  G4double phi1 = CLHEP::twopi*G4UniformRand();
259  cosPhi1 = cos(phi1);
260  sinPhi1 = sin(phi1);
261 
262  // Sample second substep scattering angle
263  SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
264  G4double phi2 = CLHEP::twopi*G4UniformRand();
265  cosPhi2 = cos(phi2);
266  sinPhi2 = sin(phi2);
267 
268  // Overall scattering direction
269  us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
270  vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
271  ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2;
272  G4double sqrtA=sqrt(scrA);
273  if(acos(ws)<sqrtA)//small angle approximation for theta less than screening angle
274  {
275  G4int i=0;
276  do{i++;
277  ws=1.+Qn12*G4Log(G4UniformRand());
278  }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
279  if(i>=19)ws=cos(sqrtA);
280  wss=std::sqrt((1.-ws*ws));
281  us=wss*std::cos(phi1);
282  vs=wss*std::sin(phi1);
283  }
284  }
285 
286  //G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
287  G4ThreeVector newDirection(us,vs,ws);
288  newDirection.rotateUz(oldDirection);
289  fParticleChange->ProposeMomentumDirection(newDirection);
290 
291  // corresponding to error less than 1% in the exact formula of <z>
292  if(Qn1<0.02) { z_coord = 1.0 - Qn1*(0.5 - Qn1/6.); }
293  else { z_coord = (1.-G4Exp(-Qn1))/Qn1; }
294  G4double rr = zPathLength*std::sqrt((1.- z_coord*z_coord)/(1.-ws*ws));
295  x_coord = rr*us;
296  y_coord = rr*vs;
297 
298  // displacement is computed relatively to the end point
299  z_coord -= 1.0;
300  z_coord *= zPathLength;
301  /*
302  G4cout << "G4GS::SampleSecondaries: e(MeV)= " << kineticEnergy
303  << " sinTheta= " << sqrt(1.0 - ws*ws)
304  << " trueStep(mm)= " << tPathLength
305  << " geomStep(mm)= " << zPathLength
306  << G4endl;
307  */
308 
309  fDisplacement.set(x_coord,y_coord,z_coord);
310  fDisplacement.rotateUz(oldDirection);
311 
312  return fDisplacement;
313 }
314 
315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
316 
317 void
318 G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
319  G4double &cost, G4double &sint)
320 {
321  G4double r1,tet,xi=0.;
322  G4double Qn1 = 2.* lambdan;
323  if(scrA < 10.) { Qn1 *= scrA*((1.+scrA)*G4Log(1.+1./scrA)-1.); }
324  else { Qn1*= (1.0 - 0.5/scrA - 0.5/(scrA*scrA)) ; }
325  if (Qn1<0.001)
326  {
327  do{
328  r1=G4UniformRand();
329  xi=-0.5*Qn1*G4Log(G4UniformRand());
330  tet=acos(1.-xi);
331  }while(tet*r1*r1>sin(tet));
332  }
333  else if(Qn1>0.5) { xi=2.*G4UniformRand(); }//isotropic distribution
334  else{ xi=2.*(GSTable->SampleTheta(lambdan,scrA,G4UniformRand()));}
335 
336 
337  if(xi<0.)xi=0.;
338  else if(xi>2.)xi=2.;
339 
340  cost=(1. - xi);
341  sint=sqrt(xi*(2.-xi));
342 
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
346 // Polynomial log-log interpolation of Lambda0 and Lambda1 between 100 eV - 1 GeV
347 // linear log-log extrapolation between 1 GeV - 100 TeV
348 
349 void
350 G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z,
351  G4double kinEnergy,G4double &Sig0,
352  G4double &Sig1)
353 {
354  G4double x1,x2,y1,y2,acoeff,bcoeff;
355  G4double kineticE = kinEnergy;
356  if(kineticE<lowKEnergy)kineticE=lowKEnergy;
357  if(kineticE>highKEnergy)kineticE=highKEnergy;
358  kineticE /= eV;
359  G4double logE=G4Log(kineticE);
360 
361  G4int iZ = G4int(Z);
362  if(iZ > 103) iZ = 103;
363 
364  G4int enerInd=0;
365  for(G4int i=0;i<105;i++)
366  {
367  if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;}
368  }
369 
370  if(p==G4Electron::Electron())
371  {
372  if(kineticE<=1.0e+9)//Interpolation of the form y=ax²+b
373  {
374  x1=ener[enerInd];
375  x2=ener[enerInd+1];
376  y1=TCSE[iZ-1][enerInd];
377  y2=TCSE[iZ-1][enerInd+1];
378  acoeff=(y2-y1)/(x2*x2-x1*x1);
379  bcoeff=y2-acoeff*x2*x2;
380  Sig0=acoeff*logE*logE+bcoeff;
381  Sig0 =G4Exp(Sig0);
382  y1=FTCSE[iZ-1][enerInd];
383  y2=FTCSE[iZ-1][enerInd+1];
384  acoeff=(y2-y1)/(x2*x2-x1*x1);
385  bcoeff=y2-acoeff*x2*x2;
386  Sig1=acoeff*logE*logE+bcoeff;
387  Sig1=G4Exp(Sig1);
388  }
389  else //Interpolation of the form y=ax+b
390  {
391  x1=ener[104];
392  x2=ener[105];
393  y1=TCSE[iZ-1][104];
394  y2=TCSE[iZ-1][105];
395  Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
396  Sig0=G4Exp(Sig0);
397  y1=FTCSE[iZ-1][104];
398  y2=FTCSE[iZ-1][105];
399  Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
400  Sig1=G4Exp(Sig1);
401  }
402  }
403  if(p==G4Positron::Positron())
404  {
405  if(kinEnergy<=1.0e+9)
406  {
407  x1=ener[enerInd];
408  x2=ener[enerInd+1];
409  y1=TCSP[iZ-1][enerInd];
410  y2=TCSP[iZ-1][enerInd+1];
411  acoeff=(y2-y1)/(x2*x2-x1*x1);
412  bcoeff=y2-acoeff*x2*x2;
413  Sig0=acoeff*logE*logE+bcoeff;
414  Sig0 =G4Exp(Sig0);
415  y1=FTCSP[iZ-1][enerInd];
416  y2=FTCSP[iZ-1][enerInd+1];
417  acoeff=(y2-y1)/(x2*x2-x1*x1);
418  bcoeff=y2-acoeff*x2*x2;
419  Sig1=acoeff*logE*logE+bcoeff;
420  Sig1=G4Exp(Sig1);
421  }
422  else
423  {
424  x1=ener[104];
425  x2=ener[105];
426  y1=TCSP[iZ-1][104];
427  y2=TCSP[iZ-1][105];
428  Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
429  Sig0 =G4Exp(Sig0);
430  y1=FTCSP[iZ-1][104];
431  y2=FTCSP[iZ-1][105];
432  Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
433  Sig1=G4Exp(Sig1);
434  }
435  }
436 
437  Sig0 *= barn;
438  Sig1 *= barn;
439 
440 }
441 
442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
443 
445 {
446  SetParticle(track->GetDynamicParticle()->GetDefinition());
447  firstStep = true;
448  inside = false;
449  insideskin = false;
450  tlimit = geombig;
451 }
452 
453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
454 //t->g->t step transformations taken from Ref.6
455 
456 G4double
458  G4double& currentMinimalStep)
459 {
460  tPathLength = currentMinimalStep;
461  const G4DynamicParticle* dp = track.GetDynamicParticle();
462  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
463  G4StepStatus stepStatus = sp->GetStepStatus();
464  currentCouple = track.GetMaterialCutsCouple();
465  SetCurrentCouple(currentCouple);
466  currentMaterialIndex = currentCouple->GetIndex();
467  currentKinEnergy = dp->GetKineticEnergy();
468  currentRange = GetRange(particle,currentKinEnergy,currentCouple);
469 
470  lambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
471 
472  // stop here if small range particle
473  if(inside || tPathLength < tlimitminfix) {
474  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
475  }
476  if(tPathLength > currentRange) tPathLength = currentRange;
477 
478  G4double presafety = sp->GetSafety();
479 
480  //G4cout << "G4GS::StepLimit tPathLength= "
481  // <<tPathLength<<" safety= " << presafety
482  // << " range= " <<currentRange<< " lambda= "<<lambda1
483  // << " Alg: " << steppingAlgorithm <<G4endl;
484 
485  // far from geometry boundary
486  if(currentRange < presafety)
487  {
488  inside = true;
489  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
490  }
491 
492  // standard version
493  //
495  {
496  //compute geomlimit and presafety
497  G4double geomlimit = ComputeGeomLimit(track, presafety, tPathLength);
498 
499  // is far from boundary
500  if(currentRange <= presafety)
501  {
502  inside = true;
503  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
504  }
505 
506  smallstep += 1.;
507  insideskin = false;
508 
509  if(firstStep || stepStatus == fGeomBoundary)
510  {
511  rangeinit = currentRange;
512  if(firstStep) smallstep = 1.e10;
513  else smallstep = 1.;
514 
515  //define stepmin here (it depends on lambda!)
516  //rough estimation of lambda_elastic/lambda_transport
517  G4double rat = currentKinEnergy/MeV ;
518  rat = 1.e-3/(rat*(10.+rat)) ;
519  //stepmin ~ lambda_elastic
520  stepmin = rat*lambda1;
521  skindepth = skin*stepmin;
522  //define tlimitmin
523  tlimitmin = 10.*stepmin;
524  if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
525 
526  //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
527  // << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
528  // constraint from the geometry
529  if((geomlimit < geombig) && (geomlimit > geommin))
530  {
531  if(stepStatus == fGeomBoundary)
532  tgeom = geomlimit/facgeom;
533  else
534  tgeom = 2.*geomlimit/facgeom;
535  }
536  else
537  tgeom = geombig;
538 
539  }
540 
541  //step limit
542  tlimit = facrange*rangeinit;
543  if(tlimit < facsafety*presafety)
544  tlimit = facsafety*presafety;
545 
546  //lower limit for tlimit
547  if(tlimit < tlimitmin) tlimit = tlimitmin;
548 
549  if(tlimit > tgeom) tlimit = tgeom;
550 
551  //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
552  // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
553 
554  // shortcut
555  if((tPathLength < tlimit) && (tPathLength < presafety) &&
556  (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
557  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
558 
559  // step reduction near to boundary
560  if(smallstep < skin)
561  {
562  tlimit = stepmin;
563  insideskin = true;
564  }
565  else if(geomlimit < geombig)
566  {
567  if(geomlimit > skindepth)
568  {
569  if(tlimit > geomlimit-0.999*skindepth)
570  tlimit = geomlimit-0.999*skindepth;
571  }
572  else
573  {
574  insideskin = true;
575  if(tlimit > stepmin) tlimit = stepmin;
576  }
577  }
578 
579  if(tlimit < stepmin) tlimit = stepmin;
580 
581  if(tPathLength > tlimit) tPathLength = tlimit;
582 
583  }
584  // for 'normal' simulation with or without magnetic field
585  // there no small step/single scattering at boundaries
586  else if(steppingAlgorithm == fUseSafety)
587  {
588  // compute presafety again if presafety <= 0 and no boundary
589  // i.e. when it is needed for optimization purposes
590  if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
591  presafety = ComputeSafety(sp->GetPosition(),tPathLength);
592 
593  // is far from boundary
594  if(currentRange < presafety)
595  {
596  inside = true;
597  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
598  }
599 
600  if(firstStep || stepStatus == fGeomBoundary)
601  {
602  rangeinit = currentRange;
603  fr = facrange;
604  // 9.1 like stepping for e+/e- only (not for muons,hadrons)
605  if(mass < masslimite)
606  {
607  if(lambda1 > currentRange)
608  rangeinit = lambda1;
609  if(lambda1 > lambdalimit)
610  fr *= 0.75+0.25*lambda1/lambdalimit;
611  }
612 
613  //lower limit for tlimit
614  G4double rat = currentKinEnergy/MeV ;
615  rat = 1.e-3/(rat*(10.+rat)) ;
616  tlimitmin = 10.*lambda1*rat;
617  if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
618  }
619  //step limit
620  tlimit = fr*rangeinit;
621 
622  if(tlimit < facsafety*presafety)
623  tlimit = facsafety*presafety;
624 
625  //lower limit for tlimit
626  if(tlimit < tlimitmin) tlimit = tlimitmin;
627 
628  if(tPathLength > tlimit) tPathLength = tlimit;
629  }
630 
631  // version similar to 7.1 (needed for some experiments)
632  else
633  {
634  if (stepStatus == fGeomBoundary)
635  {
636  if (currentRange > lambda1) tlimit = facrange*currentRange;
637  else tlimit = facrange*lambda1;
638 
639  if(tlimit < tlimitmin) tlimit = tlimitmin;
640  if(tPathLength > tlimit) tPathLength = tlimit;
641  }
642  }
643  //G4cout << "tPathLength= " << tPathLength
644  // << " currentMinimalStep= " << currentMinimalStep << G4endl;
645  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
646 }
647 
648 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
649 // taken from Ref.6
651 {
652  firstStep = false;
653  par1 = -1. ;
654  par2 = par3 = 0. ;
655 
656  // do the true -> geom transformation
657  zPathLength = tPathLength;
658 
659  // z = t for very small tPathLength
660  if(tPathLength < tlimitminfix) { return zPathLength; }
661 
662  // this correction needed to run MSC with eIoni and eBrem inactivated
663  // and makes no harm for a normal run
664  if(tPathLength > currentRange)
665  { tPathLength = currentRange; }
666 
667  G4double tau = tPathLength/lambda1 ;
668 
669  if ((tau <= tausmall) || insideskin) {
670  zPathLength = tPathLength;
671  if(zPathLength > lambda1) { zPathLength = lambda1; }
672  return zPathLength;
673  }
674 
675  G4double zmean = tPathLength;
676  if (tPathLength < currentRange*dtrl) {
677  if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
678  else zmean = lambda1*(1.-G4Exp(-tau));
679  } else if(currentKinEnergy < mass || tPathLength == currentRange) {
680  par1 = 1./currentRange ;
681  par2 = 1./(par1*lambda1) ;
682  par3 = 1.+par2 ;
683  if(tPathLength < currentRange)
684  zmean = (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3) ;
685  else
686  zmean = 1./(par1*par3) ;
687  } else {
688  G4double T1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
689 
690  lambda11 = GetTransportMeanFreePath(particle,T1);
691 
692  par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
693  par2 = 1./(par1*lambda1) ;
694  par3 = 1.+par2 ;
695  zmean = (1.-G4Exp(par3*G4Log(lambda11/lambda1)))/(par1*par3) ;
696  }
697 
698  zPathLength = zmean ;
699  // sample z
700  if(samplez) {
701 
702  const G4double ztmax = 0.99;
703  G4double zt = zmean/tPathLength ;
704 
705  if (tPathLength > stepmin && zt < ztmax) {
706 
707  G4double u,cz1;
708  if(zt >= 0.333333333) {
709 
710  G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
711  cz1 = 1.+cz ;
712  G4double u0 = cz/cz1 ;
713  G4double grej ;
714  do {
715  u = G4Exp(G4Log(G4UniformRand())/cz1) ;
716  grej = exp(cz*G4Log(u/u0))*(1.-u)/(1.-u0) ;
717  } while (grej < G4UniformRand()) ;
718 
719  } else {
720  cz1 = 1./zt-1.;
721  u = 1.-G4Exp(G4Log(G4UniformRand())/cz1) ;
722  }
723  zPathLength = tPathLength*u ;
724  }
725  }
726  if(zPathLength > lambda1) zPathLength = lambda1;
727  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda1 << G4endl;
728 
729  return zPathLength;
730 }
731 
732 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
733 // taken from Ref.6
734 G4double
736 {
737  // step defined other than transportation
738  if(geomStepLength == zPathLength && tPathLength <= currentRange)
739  return tPathLength;
740 
741  // t = z for very small step
742  zPathLength = geomStepLength;
743  tPathLength = geomStepLength;
744  if(geomStepLength < tlimitminfix) return tPathLength;
745 
746  // recalculation
747  if((geomStepLength > lambda1*tausmall) && !insideskin)
748  {
749  if(par1 < 0.)
750  tPathLength = -lambda1*G4Log(1.-geomStepLength/lambda1) ;
751  else
752  {
753  if(par1*par3*geomStepLength < 1.)
754  tPathLength = (1.-G4Exp(G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
755  else
756  tPathLength = currentRange;
757  }
758  }
759  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
760  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
761 
762  return tPathLength;
763 }
764 
765 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
766 //Total & first transport x sections for e-/e+ generated from ELSEPA code
767 
768 void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
769 {
770  G4String filename = "XSECTIONS.dat";
771 
772  char* path = getenv("G4LEDATA");
773  if (!path)
774  {
775  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()","em0006",
777  "Environment variable G4LEDATA not defined");
778  return;
779  }
780 
781  G4String pathString(path);
782  G4String dirFile = pathString + "/msc_GS/" + filename;
783  std::ifstream infile(dirFile, std::ios::in);
784  if( !infile.is_open()) {
786  ed << "Data file <" + dirFile + "> is not opened!";
787  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
788  "em0003",FatalException,ed);
789  return;
790  }
791 
792  // Read parameters from tables and take logarithms
793  G4double aRead;
794  for(G4int i=0 ; i<106 ;i++){
795  infile >> aRead;
796  if(infile.fail()) {
798  ed << "Error reading <" + dirFile + "> loop #1 i= " << i;
799  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
800  "em0003",FatalException,ed);
801  return;
802  } else {
803  if(aRead > 0.0) { aRead = G4Log(aRead); }
804  else { aRead = 0.0; }
805  }
806  ener[i]=aRead;
807  }
808  for(G4int j=0;j<103;j++){
809  for(G4int i=0;i<106;i++){
810  infile >> aRead;
811  if(infile.fail()) {
813  ed << "Error reading <" + dirFile + "> loop #2 j= " << j
814  << "; i= " << i;
815  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
816  "em0003",FatalException,ed);
817  return;
818  } else {
819  if(aRead > 0.0) { aRead = G4Log(aRead); }
820  else { aRead = 0.0; }
821  }
822  TCSE[j][i]=aRead;
823  }
824  }
825  for(G4int j=0;j<103;j++){
826  for(G4int i=0;i<106;i++){
827  infile >> aRead;
828  if(infile.fail()) {
830  ed << "Error reading <" + dirFile + "> loop #3 j= " << j
831  << "; i= " << i;
832  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
833  "em0003",FatalException,ed);
834  return;
835  } else {
836  if(aRead > 0.0) { aRead = G4Log(aRead); }
837  else { aRead = 0.0; }
838  }
839  FTCSE[j][i]=aRead;
840  }
841  }
842  for(G4int j=0;j<103;j++){
843  for(G4int i=0;i<106;i++){
844  infile >> aRead;
845  if(infile.fail()) {
847  ed << "Error reading <" + dirFile + "> loop #4 j= " << j
848  << "; i= " << i;
849  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
850  "em0003",FatalException,ed);
851  return;
852  } else {
853  if(aRead > 0.0) { aRead = G4Log(aRead); }
854  else { aRead = 0.0; }
855  }
856  TCSP[j][i]=aRead;
857  }
858  }
859  for(G4int j=0;j<103;j++){
860  for(G4int i=0;i<106;i++){
861  infile >> aRead;
862  if(infile.fail()) {
864  ed << "Error reading <" + dirFile + "> loop #5 j= " << j
865  << "; i= " << i;
866  G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
867  "em0003",FatalException,ed);
868  return;
869  } else {
870  if(aRead > 0.0) { aRead = G4Log(aRead); }
871  else { aRead = 0.0; }
872  }
873  FTCSP[j][i]=aRead;
874  }
875  }
876 
877  infile.close();
878 }
879 
880 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void set(double x, double y, double z)
G4GoudsmitSaundersonMscModel(const G4String &nam="GoudsmitSaunderson")
G4double facgeom
Definition: G4VMscModel.hh:177
virtual G4double ComputeTrueStepLength(G4double geomStepLength)
static G4LossTableManager * Instance()
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
G4double dtrl
Definition: G4VMscModel.hh:180
const G4DynamicParticle * GetDynamicParticle() const
G4double facrange
Definition: G4VMscModel.hh:176
const char * p
Definition: xmltok.h:285
G4StepStatus GetStepStatus() const
G4bool samplez
Definition: G4VMscModel.hh:188
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:246
G4double skin
Definition: G4VMscModel.hh:179
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4double SampleTheta(G4double, G4double, G4double)
G4ParticleDefinition * GetDefinition() const
const G4Step * GetStep() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
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:308
virtual G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep)
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:204
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetPosition() const
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:288
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:89
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
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
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)
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4double facsafety
Definition: G4VMscModel.hh:178
G4double GetDEDX(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:273
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:419
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double, G4double, G4double)
G4double GetSafety() const
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:186
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
virtual G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety)
double G4double
Definition: G4Types.hh:76
const G4Material * GetMaterial() const
virtual G4double ComputeGeomPathLength(G4double truePathLength)