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